From ba37cf73d5586e3a9a59e0cb87bc01d25057a1b4 Mon Sep 17 00:00:00 2001 From: Meg Fowler Date: Wed, 23 Oct 2024 08:29:23 -0600 Subject: [PATCH] Change name of notebook --- ...strialCouplingIndex_VisualCompareObs.ipynb | 975 ++++++++++++++++++ 1 file changed, 975 insertions(+) create mode 100755 examples/nblibrary/lnd/Global_TerrestrialCouplingIndex_VisualCompareObs.ipynb diff --git a/examples/nblibrary/lnd/Global_TerrestrialCouplingIndex_VisualCompareObs.ipynb b/examples/nblibrary/lnd/Global_TerrestrialCouplingIndex_VisualCompareObs.ipynb new file mode 100755 index 0000000..826a2f0 --- /dev/null +++ b/examples/nblibrary/lnd/Global_TerrestrialCouplingIndex_VisualCompareObs.ipynb @@ -0,0 +1,975 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "99564fab-c321-4116-8229-b16eefa1536e", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "# Compute land-atmosphere coupling indices \n", + "This notebook takes in a series of CESM simulations, computes the land-atmosphere coupling index (CI; \n", + "terrestrial leg only currently), and plots those seasonal means.
\n", + "- Note: Built to use monthly output; ideally, CI should be based on daily data. \n", + "- Optional: Comparison against FLUXNET obs\n", + "

\n", + "Notebook created by mdfowler@ucar.edu; Last update: 2 Aug 2024 " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "750da831-1c5c-4b41-947e-a9e57a62a820", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "outputs": [], + "source": [ + "import os\n", + "import glob\n", + "import numpy as np\n", + "import xarray as xr\n", + "import datetime\n", + "from datetime import date, timedelta\n", + "import dask\n", + "import pandas as pd\n", + "import sys\n", + "\n", + "# Plotting utils\n", + "import matplotlib\n", + "import matplotlib.pyplot as plt\n", + "import cartopy\n", + "import cartopy.crs as ccrs\n", + "import uxarray as uxr" + ] + }, + { + "cell_type": "markdown", + "id": "774bd269-ce50-4449-b32f-83246b74b73c", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "## 1. Modify this section for each run" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7f8f2d17-c653-4ad1-9dc3-c49bf836ceb6", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [ + "parameters" + ] + }, + "outputs": [], + "source": [ + "## - - - - - - - - - - - - - - - - - - - - - -\n", + "## Settings for case locations + names\n", + "## - - - - - - - - - - - - - - - - - - - - - -\n", + "## Where observations are stored\n", + "# obsDir = '/glade/campaign/cgd/tss/people/mdfowler/FLUXNET2015/' ## Need to copy into CUPiD Data\n", + "\n", + "## Where CESM timeseries data is stored\n", + "CESM_output_dir = \"/glade/campaign/cesm/development/cross-wg/diagnostic_framework/CESM_output_for_testing/\"\n", + "\n", + "\n", + "## Full casenames that are present in CESM_output_dir and in individual filenames\n", + "# caseNames = [\n", + "# 'b.e23_alpha16b.BLT1850.ne30_t232.054',\n", + "# 'b.e30_beta02.BLT1850.ne30_t232.104',\n", + "# ]\n", + "case_name = \"b.e30_beta02.BLT1850.ne30_t232.104\"\n", + "\n", + "# clmFile_h = '.h0.'\n", + "\n", + "start_date = \"0001-01-01\"\n", + "end_date = \"0101-01-01\"\n", + "\n", + "## - - - - - - - - - - - - - - - - - - - - - -\n", + "## Optional settings for notebook\n", + "## - - - - - - - - - - - - - - - - - - - - - -\n", + "\n", + "## If comparison against FLUXNET desired\n", + "# fluxnet_comparison = True" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0014712f-d094-4dae-b583-740bf7a9789c", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "outputs": [], + "source": [ + "## - - - - - - - - - - - - - - - - - - - - - -\n", + "## Settings for computing coupling index\n", + "## - - - - - - - - - - - - - - - - - - - - - -\n", + "startYrs = [start_date.split(\"-\")[0]]\n", + "endYrs = [f\"{int(end_date.split('-')[0])-1:04d}\"]\n", + "\n", + "caseNames = [\n", + " case_name,\n", + " # base_case_name,\n", + "]\n", + "\n", + "shortNames = [case.split(\".\")[-1] for case in caseNames]" + ] + }, + { + "cell_type": "markdown", + "id": "d70024c7-0af2-48b9-9041-893f40e613ec", + "metadata": {}, + "source": [ + "## 2. Read in model data and compute terrestrial coupling index" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "304ce8d0-6aab-4fcb-9635-2a78e270f3c7", + "metadata": { + "editable": true, + "scrolled": true, + "slideshow": { + "slide_type": "" + }, + "tags": [ + "hide-input" + ] + }, + "outputs": [], + "source": [ + "\"\"\"\n", + "Inputs: xname -- controlling variable \n", + " yname -- responding variable\n", + " ds -- dataset containing xname and yname data \n", + " \n", + "This is pulled almost directly from Ahmed Tawfik's CI code here: \n", + " https://github.com/abtawfik/coupling-metrics/blob/new_version_1/src/comet/metrics/coupling_indices.py \n", + "\"\"\"\n", + "\n", + "\n", + "def compute_couplingIndex_cesm(xname, yname, xDS, yDS):\n", + " xday = xDS[xname].groupby(\"time.season\")\n", + " yday = yDS[yname].groupby(\"time.season\")\n", + "\n", + " # Get the covariance of the two (numerator in coupling index)\n", + " covarTerm = ((xday - xday.mean()) * (yday - yday.mean())).groupby(\n", + " \"time.season\"\n", + " ).sum() / xday.count()\n", + "\n", + " # Now compute the actual coupling index\n", + " couplingIndex = covarTerm / xday.std()\n", + "\n", + " return couplingIndex" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3c35ae8d-dfff-44b0-a854-dfc2b5b030c0", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [ + "hide-input" + ] + }, + "outputs": [], + "source": [ + "for iCase in range(len(caseNames)):\n", + " ## Check first if coupling index has already been created:\n", + " TCI_filePath = (\n", + " \"/glade/derecho/scratch/mdfowler/\"\n", + " + caseNames[0]\n", + " + \"_TerrestrialCouplingIndex_SHvsSM.nc\"\n", + " )\n", + "\n", + " if os.path.exists(TCI_filePath): # Use previously computed TCI\n", + " print(\"Using previously computed coupling index saved in file \", TCI_filePath)\n", + " else: # Compute TCI\n", + "\n", + " # Get list of necessary time series files\n", + " soilWater_file = np.sort(\n", + " glob.glob(\n", + " CESM_output_dir\n", + " + \"/\"\n", + " + caseNames[iCase]\n", + " + \"/lnd/proc/tseries/\"\n", + " + caseNames[iCase]\n", + " + clmFile_h\n", + " + \"SOILWATER_10CM.\"\n", + " + startYrs[iCase]\n", + " + \"??-\"\n", + " + endYrs[iCase]\n", + " + \"??.nc\"\n", + " )\n", + " )\n", + " if len(soilWater_file) == 0:\n", + " print(\"Soil moisture file not found!\")\n", + " elif len(soilWater_file) > 1:\n", + " print(\n", + " \"More than one file matches requested time period and case for soil moisture.\"\n", + " )\n", + " elif len(soilWater_file) == 1:\n", + " soilWater_DS = xr.open_dataset(soilWater_file[0], decode_times=True)\n", + "\n", + " sh_file = np.sort(\n", + " glob.glob(\n", + " CESM_output_dir\n", + " + \"/\"\n", + " + caseNames[iCase]\n", + " + \"/lnd/proc/tseries/\"\n", + " + caseNames[iCase]\n", + " + clmFile_h\n", + " + \"FSH_TO_COUPLER.\"\n", + " + startYrs[iCase]\n", + " + \"??-\"\n", + " + endYrs[iCase]\n", + " + \"??.nc\"\n", + " )\n", + " )\n", + " if len(sh_file) == 0:\n", + " print(\"Land-based SHFLX file not found!\")\n", + " elif len(sh_file) > 1:\n", + " print(\"More than one file matches requested time period and case for SH.\")\n", + " elif len(sh_file) == 1:\n", + " shflx_DS = xr.open_dataset(sh_file[0])\n", + "\n", + " # If years start at 0000, offset by 1700 years for analysis\n", + " yrOffset = 1850\n", + " if shflx_DS[\"time.year\"].values[0] < 1500:\n", + " shflx_DS[\"time\"] = shflx_DS.time + timedelta(days=yrOffset * 365)\n", + " if soilWater_DS[\"time.year\"].values[0] < 1500:\n", + " soilWater_DS[\"time\"] = soilWater_DS.time + timedelta(days=yrOffset * 365)\n", + " # Convert times to datetime for easier use\n", + " shflx_DS[\"time\"] = shflx_DS.indexes[\"time\"].to_datetimeindex()\n", + " soilWater_DS[\"time\"] = soilWater_DS.indexes[\"time\"].to_datetimeindex()\n", + "\n", + " # Add case ID (short name) to the DS\n", + " shflx_DS = shflx_DS.assign_coords({\"case\": shortNames[iCase]})\n", + " soilWater_DS = soilWater_DS.assign_coords({\"case\": shortNames[iCase]})\n", + "\n", + " ## Compute coupling index and save to netCDF file\n", + " ## - - - - - - - - - - - - - - - - - - - - - - - - -\n", + " xname = \"SOILWATER_10CM\" # Controlling variable\n", + " yname = \"FSH_TO_COUPLER\" # Responding variable\n", + "\n", + " xDS = soilWater_DS\n", + " yDS = shflx_DS\n", + "\n", + " couplingInd = compute_couplingIndex_cesm(xname, yname, xDS, yDS)\n", + "\n", + " filePath = (\n", + " \"/glade/derecho/scratch/mdfowler/\"\n", + " + caseNames[0]\n", + " + \"_TerrestrialCouplingIndex_SHvsSM.nc\"\n", + " )\n", + " couplingInd.to_netcdf(filePath)\n", + " print(\"File created: \", filePath)" + ] + }, + { + "cell_type": "markdown", + "id": "5f8fba2a-98d2-4e94-9d71-3b2625e16032", + "metadata": {}, + "source": [ + "### 1.1 Read in FLUXNET data if requested" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fdf20c82-5a01-4ab9-8881-d9388b1b2356", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [ + "hide-input" + ] + }, + "outputs": [], + "source": [ + "# --------------------------------------------------------\n", + "# Function to read requested variables from FLUXNET file.\n", + "# - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n", + "#\n", + "# Inputs: fileName --> Full path to FLUXNET data file\n", + "# varNames --> An array of variable names to be\n", + "# retrieved from said data file.\n", + "# NOTE: If you wish to retrieve *all*\n", + "# variables, pass the string 'ALL'.\n", + "#\n", + "# Outputs: fluxnetID --> ID string used to identify station\n", + "# fluxnetDS --> An x-array dataset containing the\n", + "# requested variables.\n", + "# Missing values will be set to NaN.\n", + "#\n", + "# --------------------------------------------------------\n", + "\n", + "\n", + "def readFLUXNET_var(fileName, varNames):\n", + " # Get ID of station\n", + " startID = fileName.find(\"FLX_\")\n", + " fluxnetID = fileName[startID + 4 : startID + 10]\n", + "\n", + " # If this is taking a long time or you just want to know where in the stations you are, uncomment print statement\n", + " # print('Reading in site - ', fluxnetID)\n", + "\n", + " # Read in CSV file containing data\n", + " dataDF = pd.read_csv(fileName)\n", + "\n", + " # Return ALL variables from dataDF if requested\n", + " if varNames == \"ALL\":\n", + " fluxnetDF = dataDF\n", + "\n", + " # Set any value that's missing to NaN (not -9999)\n", + " fluxnetDF = fluxnetDF.replace(-9999, np.nan)\n", + "\n", + " # If time has been requested, reformat to pandas date index\n", + " fluxnetDF[\"TIMESTAMP\"] = pd.to_datetime(\n", + " fluxnetDF[\"TIMESTAMP\"].values, format=\"%Y%m%d\"\n", + " )\n", + " fluxnetDF = fluxnetDF.set_index([\"TIMESTAMP\"])\n", + "\n", + " # Convert dataframe to Xarray Dataset (required to use coupling metrics toolbox)\n", + " # NOTE: since current implementation doesn't use the pre-formatted CoMeT, might not need this step now\n", + " fluxnetDS = fluxnetDF.to_xarray()\n", + "\n", + " # Reduce returned DF to contain only variables of interest\n", + " else:\n", + "\n", + " # Check that requested variables are available in specific file\n", + " errCount = 0 # Initialize flag for error\n", + " colNames = dataDF.columns.values # Available variables in file\n", + "\n", + " for iVar in range(len(varNames)): # Check each variable individually\n", + " if (varNames[iVar] in colNames) == False:\n", + " # Turn on print statement for more verbose output\n", + " # print('** ERROR: %13s not contained in file for %8s **' %(varNames[iVar], fluxnetID))\n", + "\n", + " # If any variable is not conatined in file, return a NaN\n", + " fluxnetDS = -999\n", + " errCount = errCount + 1\n", + "\n", + " # If all the variables *are* available, isolate those in DF and return that\n", + " if errCount == 0:\n", + " fluxnetDF = dataDF[varNames]\n", + "\n", + " # Set any value that's missing to NaN (not -999)\n", + " fluxnetDF = fluxnetDF.replace(-9999, np.nan)\n", + "\n", + " # If time has been requested, reformat to pandas make index\n", + " if (\"TIMESTAMP\" in varNames) == True:\n", + " fluxnetDF[\"TIMESTAMP\"] = pd.to_datetime(\n", + " fluxnetDF[\"TIMESTAMP\"].values, format=\"%Y%m%d\"\n", + " )\n", + " fluxnetDF = fluxnetDF.set_index([\"TIMESTAMP\"])\n", + "\n", + " # Convert dataframe to Xarray Dataset (required to use coupling metrics toolbox)\n", + " fluxnetDS = fluxnetDF.to_xarray()\n", + "\n", + " return (fluxnetID, fluxnetDS)\n", + "\n", + "\n", + "\"\"\"\n", + "Inputs: xname -- controlling variable \n", + " yname -- responding variable\n", + " ds -- dataset containing xname and yname data \n", + " \n", + "This is pulled almost directly from Ahmed Tawfik's CI code here: \n", + " https://github.com/abtawfik/coupling-metrics/blob/new_version_1/src/comet/metrics/coupling_indices.py \n", + "\"\"\"\n", + "\n", + "\n", + "def compute_couplingIndex_FLUXNET(xname, yname, ds):\n", + " xday = ds[xname].groupby(\"TIMESTAMP.season\")\n", + " yday = ds[yname].groupby(\"TIMESTAMP.season\")\n", + "\n", + " # Get the covariance of the two (numerator in coupling index)\n", + " covarTerm = ((xday - xday.mean()) * (yday - yday.mean())).groupby(\n", + " \"TIMESTAMP.season\"\n", + " ).sum() / xday.count()\n", + "\n", + " # Now compute the actual coupling index\n", + " couplingIndex = covarTerm / xday.std()\n", + "\n", + " return couplingIndex" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9acb032f-f106-4783-8c57-8e0d3a44eb3e", + "metadata": { + "editable": true, + "scrolled": true, + "slideshow": { + "slide_type": "" + }, + "tags": [ + "hide-input", + "hide-output" + ] + }, + "outputs": [], + "source": [ + "if fluxnet_comparison == True:\n", + " # obsDir = '/glade/campaign/cgd/tss/people/mdfowler/FLUXNET2015/' ## Need to copy into CUPiD Data\n", + "\n", + " ## Metadata files\n", + " siteInfoFile = obsDir + \"SiteList.csv\"\n", + " siteInfoDF = pd.read_csv(siteInfoFile)\n", + "\n", + " metadataFile = obsDir + \"FLX_AA-Flx_BIF_ALL_20200501/FLX_AA-Flx_BIF_DD_20200501.csv\"\n", + " metadataDF = pd.read_csv(metadataFile)\n", + "\n", + " ## List of all station files\n", + " dataFiles = glob.glob(obsDir + \"FLX_*/*SUBSET_DD*\")\n", + "\n", + " # Set up a few empty arrays to save data into\n", + " terraCI_fluxnetConverted = np.full(\n", + " [len(dataFiles), 4], np.nan\n", + " ) # CI using kg/m2 soil water content [nStations, seasons]\n", + "\n", + " # Also save out some data on each station\n", + " startTime_fluxnet = np.zeros(len(dataFiles), dtype=\"datetime64[s]\")\n", + " endTime_fluxnet = np.zeros(len(dataFiles), dtype=\"datetime64[s]\")\n", + " lat_fluxnet = np.full([len(dataFiles)], np.nan)\n", + " lon_fluxnet = np.full([len(dataFiles)], np.nan)\n", + " SWCdepth = np.full([len(dataFiles)], np.nan)\n", + "\n", + " stationID = []\n", + " stationID_converted = []\n", + "\n", + " allStationID = []\n", + "\n", + " # Variables I want returned:\n", + " varNames = [\"TIMESTAMP\", \"H_F_MDS\", \"SWC_F_MDS_1\", \"SWC_F_MDS_1_QC\"]\n", + "\n", + " # Loop over each station (data file)\n", + " for iStation in range(len(dataFiles)):\n", + "\n", + " # Read in data\n", + " # ----------------------------------------------------------\n", + " fluxnetID, fluxnetDS = readFLUXNET_var(dataFiles[iStation], varNames)\n", + "\n", + " # Save lat and lon for this station\n", + " # ----------------------------------------------------------\n", + " indStation = int(np.where(fluxnetID == siteInfoDF[\"SITE_ID\"])[0][0])\n", + " lat_fluxnet[iStation] = siteInfoDF[\"LOCATION_LAT\"].values[indStation]\n", + " lon_fluxnet[iStation] = siteInfoDF[\"LOCATION_LONG\"].values[indStation]\n", + " allStationID.append(fluxnetID)\n", + "\n", + " # Check that there was data saved for this particular site:\n", + " # ----------------------------------------------------------\n", + " if type(fluxnetDS) == int:\n", + " print(\"No data for station: %8s\" % fluxnetID)\n", + "\n", + " elif (np.all(np.isnan(fluxnetDS[\"H_F_MDS\"])) == True) | (\n", + " np.all(np.isnan(fluxnetDS[\"SWC_F_MDS_1\"])) == True\n", + " ):\n", + " print(\"No data for station: %8s\" % fluxnetID)\n", + "\n", + " # If data is present:\n", + " # ----------------------------------------------------------\n", + " else:\n", + " # Only consider where data is actually present for selected vars\n", + " iReal = np.where(\n", + " (np.isfinite(fluxnetDS[\"SWC_F_MDS_1\"]) == True)\n", + " & (np.isfinite(fluxnetDS[\"H_F_MDS\"]) == True)\n", + " )[0]\n", + " fluxnetDS = fluxnetDS.isel(TIMESTAMP=iReal)\n", + "\n", + " stationID.append(fluxnetID)\n", + "\n", + " # Convert units from volumetric (%) to mass (kg/m2)\n", + " # -------------------------------------------------\n", + " # Step 1: Convert from % to fraction\n", + " fracSM = (fluxnetDS[\"SWC_F_MDS_1\"].values) / 100.0\n", + "\n", + " # Step 2: Need to use depth of obs in conversion\n", + " metaData_station = metadataDF[metadataDF.SITE_ID == fluxnetID]\n", + " iSWC = np.where(metaData_station.DATAVALUE == \"SWC_F_MDS_1\")[0]\n", + " # Some locations (5) have two depths\n", + " if len(iSWC) > 1:\n", + " for iDepth in range(len(iSWC)):\n", + " SWC_DF = metaData_station[iSWC[iDepth] : iSWC[iDepth] + 4]\n", + "\n", + " depth = np.asarray(\n", + " SWC_DF[SWC_DF.VARIABLE == \"VAR_INFO_HEIGHT\"].DATAVALUE.values[0]\n", + " ).astype(float)\n", + " depthDay = np.asarray(\n", + " SWC_DF[SWC_DF.VARIABLE == \"VAR_INFO_DATE\"].DATAVALUE.values[0]\n", + " ).astype(int)\n", + " depthDay = int(\n", + " str(depthDay)[:8]\n", + " ) # Some weird ones have time attached; don't want that\n", + " depthDay = pd.to_datetime(depthDay, format=\"%Y%m%d\")\n", + "\n", + " # Keep deepest level as the depth for station\n", + " if iDepth == 0:\n", + " SWCdepth[iStation] = depth\n", + " convertSM = fracSM * 1000.0 * np.abs(depth)\n", + " else:\n", + " # Use date as break point for getting kg/m2 SWC\n", + " # Eq: SWC_kgm2 = SWC_vol [m3/m3] * 1000 [kg/m3] * depth [m]\n", + " dateArr = pd.DatetimeIndex(fluxnetDS.TIMESTAMP.values)\n", + " iTime = int(np.where(dateArr == depthDay)[0][0])\n", + " convertSM[iTime::] = (fracSM[iTime::]) * 1000.0 * np.abs(depth)\n", + "\n", + " # Keep deepest level as the depth for station\n", + " if depth < SWCdepth[iStation]:\n", + " SWCdepth[iStation] = depth\n", + "\n", + " stationID_converted.append(fluxnetID)\n", + "\n", + " # If station only has one level recorded, things are a bit easier:\n", + " else:\n", + " SWC_DF = metaData_station[iSWC[0] : iSWC[0] + 4]\n", + " SWCdepth[iStation] = np.asarray(\n", + " SWC_DF[SWC_DF.VARIABLE == \"VAR_INFO_HEIGHT\"].DATAVALUE.values[0]\n", + " ).astype(float)\n", + " convertSM = fracSM * 1000.0 * np.abs(SWCdepth[iStation])\n", + "\n", + " stationID_converted.append(fluxnetID)\n", + "\n", + " # Save converted soil moisture to dataset\n", + " fluxnetDS[\"SWC_F_MDS_1_convert\"] = ((\"TIMESTAMP\"), convertSM)\n", + "\n", + " # Save first and last time used for computing CI\n", + " # ----------------------------------------------\n", + " startTime_fluxnet[iStation] = fluxnetDS[\"TIMESTAMP\"].values[0]\n", + " endTime_fluxnet[iStation] = fluxnetDS[\"TIMESTAMP\"].values[-1]\n", + "\n", + " # Compute terrestrial coupling metric\n", + " # -----------------------------------\n", + " terraLeg = compute_couplingIndex_FLUXNET(\n", + " \"SWC_F_MDS_1_convert\", \"H_F_MDS\", fluxnetDS\n", + " )\n", + "\n", + " # If there's less than one full year of data, don't use station\n", + " # (i.e., as long as all 4 seasons are defined, save values)\n", + " if np.shape(terraLeg)[0] == 4:\n", + " terraCI_fluxnetConverted[iStation, :] = terraLeg\n", + "\n", + " seasons = terraLeg.season\n", + "\n", + " ## Print some useful information\n", + " print(\n", + " \"Number of FLUXNET stations with CI calculated: %i\"\n", + " % len(np.where(np.isfinite(terraCI_fluxnetConverted[:, 1]) == True)[0])\n", + " )\n", + "\n", + " # How many months go into each calculation of CI for JJA?\n", + " nMonths = np.full([len(dataFiles)], np.nan)\n", + "\n", + " for iSt in range(len(dataFiles)):\n", + " if np.isfinite(terraCI_fluxnetConverted[iSt, 1]):\n", + " dateRange = pd.date_range(\n", + " start=startTime_fluxnet[iSt], end=endTime_fluxnet[iSt], freq=\"ME\"\n", + " )\n", + " nMonths[iSt] = len(\n", + " np.where((dateRange.month >= 6) & (dateRange.month <= 8))[0]\n", + " )\n", + "\n", + " print(\n", + " \"Minimum number of months used for JJA mean CI: %i \" % int(np.nanmin(nMonths))\n", + " )\n", + " print(\n", + " \"Maximum number of months used for JJA mean CI: %i \" % int(np.nanmax(nMonths))\n", + " )" + ] + }, + { + "cell_type": "markdown", + "id": "dd29ba9d-bbbb-4f59-9829-250018215b53", + "metadata": {}, + "source": [ + "*Make some choices on limiting which stations are used*\n", + "- Let's limit usage to depths less than 20 cm (arbitrary, but I don't want us using non-surface soil moisture for this application). This will eliminate 11 stations.\n", + "- It would also be good to put some time limits on this. So let's say the observations need to have at least 9 months of data for JJA means (3-years). Otherwise, set terraCI to np.nan again so we don't use it." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "13afbfb4-2040-403c-9d82-3e370f47ec5d", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [ + "hide-input" + ] + }, + "outputs": [], + "source": [ + "if fluxnet_comparison == True:\n", + "\n", + " # Get stations with SWC from below 20 cm (or equal to zero)\n", + " iLimit = np.where((SWCdepth == 0.0) | (SWCdepth < -0.2))[0]\n", + "\n", + " # Set the terrestrial leg of CI to missing so we don't consider those\n", + " terraCI_fluxnetConverted[iLimit, :] = np.nan\n", + "\n", + " print(\n", + " \"Number of FLUXNET stations to use with reasonable depths of SWC: %i\"\n", + " % len(np.where(np.isfinite(terraCI_fluxnetConverted[:, 1]) == True)[0])\n", + " )\n", + "\n", + " # Get stations with less than 9 months used for JJA terrestrial CI\n", + " iLimit = np.where(nMonths < 9)[0]\n", + "\n", + " # Set to missing so we don't consider stations with less than three years of data going into the average\n", + " terraCI_fluxnetConverted[iLimit, :] = np.nan\n", + "\n", + " print(\n", + " \"Number of FLUXNET stations to use with 3+ years of JJA data: %i\"\n", + " % len(np.where(np.isfinite(terraCI_fluxnetConverted[:, 1]) == True)[0])\n", + " )" + ] + }, + { + "cell_type": "markdown", + "id": "bc387253-cdd7-4a36-956b-8ce548e963bd", + "metadata": {}, + "source": [ + "## 2. Make plots" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "67253fd1-d2f7-45fe-a59f-215303b93c06", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [ + "hide-input" + ] + }, + "outputs": [], + "source": [ + "## Load coupling index with uxarray\n", + "gridFile = (\n", + " \"/glade/p/cesmdata/cseg/inputdata/share/meshes/ne30pg3_ESMFmesh_cdf5_c20211018.nc\"\n", + ")\n", + "uxgrid = uxr.open_grid(gridFile)\n", + "uxgrid" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "43206b67-1313-4b61-94ea-b50b85a3d50c", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [ + "hide-input" + ] + }, + "outputs": [], + "source": [ + "for iCase in range(len(caseNames)):\n", + " filePath = (\n", + " \"/glade/derecho/scratch/mdfowler/\"\n", + " + caseNames[iCase]\n", + " + \"_TerrestrialCouplingIndex_SHvsSM.nc\"\n", + " )\n", + " couplingIndex_case = uxr.open_dataset(gridFile, filePath)\n", + " # Rename the variable:\n", + " couplingIndex_case = couplingIndex_case.rename(\n", + " {\"__xarray_dataarray_variable__\": \"CouplingIndex\"}\n", + " )\n", + "\n", + " # Assign case coord\n", + " couplingIndex_case = couplingIndex_case.assign_coords(\n", + " {\"case\": couplingIndex_case.case.values}\n", + " )\n", + "\n", + " # Return all the cases in a single dataset\n", + " if iCase == 0:\n", + " couplingIndex_DS = couplingIndex_case\n", + " del couplingIndex_case\n", + " else:\n", + " couplingIndex_DS = uxr.concat([couplingIndex_DS, couplingIndex_case], \"case\")\n", + " del couplingIndex_case\n", + "\n", + "print(\"Coupling index is now ready to go\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e9e22cd0-870e-47cb-b5c7-2e57bbd9af16", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [ + "hide-input" + ] + }, + "outputs": [], + "source": [ + "def make_cmap(colors, position=None, bit=False):\n", + " \"\"\"\n", + " make_cmap takes a list of tuples which contain RGB values. The RGB\n", + " values may either be in 8-bit [0 to 255] (in which bit must be set to\n", + " True when called) or arithmetic [0 to 1] (default). make_cmap returns\n", + " a cmap with equally spaced colors.\n", + " Arrange your tuples so that the first color is the lowest value for the\n", + " colorbar and the last is the highest.\n", + " position contains values from 0 to 1 to dictate the location of each color.\n", + " \"\"\"\n", + "\n", + " import matplotlib as mpl\n", + " import numpy as np\n", + "\n", + " bit_rgb = np.linspace(0, 1, 256)\n", + " if position == None:\n", + " position = np.linspace(0, 1, len(colors))\n", + " else:\n", + " if len(position) != len(colors):\n", + " sys.exit(\"position length must be the same as colors\")\n", + " elif position[0] != 0 or position[-1] != 1:\n", + " sys.exit(\"position must start with 0 and end with 1\")\n", + "\n", + " if bit:\n", + " for i in range(len(colors)):\n", + " colors[i] = (\n", + " bit_rgb[colors[i][0]],\n", + " bit_rgb[colors[i][1]],\n", + " bit_rgb[colors[i][2]],\n", + " )\n", + "\n", + " cdict = {\"red\": [], \"green\": [], \"blue\": []}\n", + " for pos, color in zip(position, colors):\n", + " cdict[\"red\"].append((pos, color[0], color[0]))\n", + " cdict[\"green\"].append((pos, color[1], color[1]))\n", + " cdict[\"blue\"].append((pos, color[2], color[2]))\n", + "\n", + " cmap = mpl.colors.LinearSegmentedColormap(\"my_colormap\", cdict, 256)\n", + " return cmap" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6e3e7b6f-98a8-4e21-8378-adfd8b399823", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [ + "hide-input" + ] + }, + "outputs": [], + "source": [ + "### Create a list of RGB tuples for terrestrial leg (SM, SHFLX)\n", + "colorsList_SMvSHF = [\n", + " (124, 135, 181),\n", + " (107, 109, 161),\n", + " (51, 82, 120),\n", + " (49, 114, 127),\n", + " (97, 181, 89),\n", + " (200, 218, 102),\n", + " (255, 242, 116),\n", + " (238, 164, 58),\n", + "] # This example uses the 8-bit RGB\n", + "my_cmap_SMvSHF = make_cmap(colorsList_SMvSHF, bit=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "882b1c92-e5b9-4cd0-aa55-4e04aca3c50c", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [ + "hide-input" + ] + }, + "outputs": [], + "source": [ + "def plotTCI_case(seasonstr, caseSel=None):\n", + "\n", + " transform = ccrs.PlateCarree()\n", + " projection = ccrs.PlateCarree()\n", + "\n", + " if caseSel:\n", + " # create a Poly Array from a 1D slice of a face-centered data variable\n", + " collection = (\n", + " couplingIndex_DS[\"CouplingIndex\"]\n", + " .sel(season=seasonstr)\n", + " .isel(case=caseSel)\n", + " .to_polycollection()\n", + " )\n", + "\n", + " collection.set_transform(transform)\n", + " collection.set_antialiased(False)\n", + " collection.set_cmap(my_cmap_SMvSHF)\n", + " collection.set_clim(vmin=-20, vmax=5)\n", + "\n", + " fig, ax = plt.subplots(\n", + " 1,\n", + " 1,\n", + " figsize=(12, 5),\n", + " facecolor=\"w\",\n", + " constrained_layout=True,\n", + " subplot_kw=dict(projection=projection),\n", + " )\n", + "\n", + " ax.coastlines()\n", + " ax.add_collection(collection)\n", + " ax.set_global()\n", + " fig.colorbar(collection, label=\"Terrestrial Coupling Index ($W m^{-2}$)\")\n", + " ax.set_title(\n", + " seasonstr\n", + " + \" Coupling Index: \"\n", + " + str(couplingIndex_DS.case.isel(case=caseSel).values)\n", + " )\n", + "\n", + " plt.show()\n", + " plt.close()\n", + "\n", + " else:\n", + " # create a Poly Array from a 1D slice of a face-centered data variable\n", + " collection = (\n", + " couplingIndex_DS[\"CouplingIndex\"].sel(season=seasonstr).to_polycollection()\n", + " )\n", + "\n", + " collection.set_transform(transform)\n", + " collection.set_antialiased(False)\n", + " collection.set_cmap(my_cmap_SMvSHF)\n", + " collection.set_clim(vmin=-20, vmax=5)\n", + "\n", + " fig, ax = plt.subplots(\n", + " 1,\n", + " 1,\n", + " figsize=(12, 5),\n", + " facecolor=\"w\",\n", + " constrained_layout=True,\n", + " subplot_kw=dict(projection=projection),\n", + " )\n", + "\n", + " ax.coastlines()\n", + " ax.add_collection(collection)\n", + " ax.set_global()\n", + " fig.colorbar(collection, label=\"Terrestrial Coupling Index ($W m^{-2}$)\")\n", + " ax.set_title(\n", + " seasonstr + \" Coupling Index: \" + str(couplingIndex_DS.case.values)\n", + " )\n", + "\n", + " if fluxnet_comparison == True:\n", + " ## Add FLUXNET obs\n", + " iSeason = np.where(seasons == seasonstr)[0]\n", + " iStations = np.where(\n", + " np.isfinite(terraCI_fluxnetConverted[:, iSeason]) == True\n", + " )[0]\n", + " norm_CI = matplotlib.colors.Normalize(vmin=-20, vmax=5)\n", + "\n", + " ax.scatter(\n", + " lon_fluxnet[iStations],\n", + " lat_fluxnet[iStations],\n", + " c=terraCI_fluxnetConverted[iStations, iSeason],\n", + " cmap=my_cmap_SMvSHF,\n", + " norm=norm_CI,\n", + " edgecolor=\"k\",\n", + " s=30,\n", + " marker=\"o\",\n", + " transform=ccrs.PlateCarree(),\n", + " )\n", + "\n", + " plt.show()\n", + " plt.close()\n", + "\n", + " return" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d43b6525-aa30-47a4-8b66-ca70a805a13a", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [ + "hide-input" + ] + }, + "outputs": [], + "source": [ + "if len(caseNames) == 1:\n", + " plotTCI_case(\"JJA\", None)\n", + " plotTCI_case(\"DJF\", None)\n", + "else:\n", + " for iCase in range(len(caseNames)):\n", + " plotTCI_case(\"JJA\", iCase)\n", + " plotTCI_case(\"DJF\", iCase)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python [conda env:cupid-analysis]", + "language": "python", + "name": "conda-env-cupid-analysis-py" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.4" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}