An Exception was encountered at ‘In [7]’.
Compute land-atmosphere coupling indices#
This notebook takes in a series of CESM simulations, computes the land-atmosphere coupling index (CI;
terrestrial leg only currently), and plots those seasonal means.
Note: Built to use monthly output; ideally, CI should be based on daily data.
Optional: Comparison against FLUXNET obs
Notebook created by mdfowler@ucar.edu; Last update: 11 Dec 2024
import os
import glob
import numpy as np
import xarray as xr
from datetime import timedelta
import pandas as pd
import sys
# Plotting utils
import matplotlib
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import uxarray as uxr
1. Modify this section for each run#
## - - - - - - - - - - - - - - - - - - - - - -
## Settings for case locations + names
## - - - - - - - - - - - - - - - - - - - - - -
## Where observations are stored
# '/glade/campaign/cesm/development/cross-wg/diagnostic_framework/CUPiD_obs_data/lnd/analysis_datasets/ungridded/timeseries/FLUXNET2015/'
obs_data_dir = ""
obsDir = ""
CESM_output_dir = os.path.join(
os.path.sep,
"glade",
"campaign",
"cesm",
"development",
"cross-wg",
"diagnostic_framework",
"CESM_output_for_testing",
)
ts_dir = None
## Full casenames that are present in CESM_output_dir or ts_dir and in individual filenames
# caseNames = [
# 'b.e23_alpha16b.BLT1850.ne30_t232.054',
# 'b.e30_beta02.BLT1850.ne30_t232.104',
# ]
case_name = "b.e30_beta02.BLT1850.ne30_t232.104"
# clmFile_h = '.h0.'
start_date = "0001-01-01"
end_date = "0101-01-01"
## - - - - - - - - - - - - - - - - - - - - - -
## Optional settings for notebook
## - - - - - - - - - - - - - - - - - - - - - -
## If comparison against FLUXNET desired
# fluxnet_comparison = True
# Parameters
case_name = "b.e30_beta06.B1850C_LTso.ne30_t232_wgx3.192"
base_case_name = "b.e30_beta06.B1850C_LTso.ne30_t232_wgx3.188"
CESM_output_dir = "/glade/derecho/scratch/hannay/archive"
base_case_output_dir = "/glade/derecho/scratch/gmarques/archive"
start_date = "0002-01-01"
end_date = "0021-12-01"
base_start_date = "0002-01-01"
base_end_date = "0021-12-01"
obs_data_dir = (
"/glade/campaign/cesm/development/cross-wg/diagnostic_framework/CUPiD_obs_data"
)
ts_dir = None
lc_kwargs = {"threads_per_worker": 1}
serial = True
clmFile_h = "h0"
fluxnet_comparison = True
obsDir = "lnd/analysis_datasets/ungridded/timeseries/FLUXNET2015/"
subset_kwargs = {}
product = "/glade/work/hannay/CUPiD/examples/key_metrics/computed_notebooks//lnd/Global_TerrestrialCouplingIndex_VisualCompareObs.ipynb"
# Set some parameter defaults
if ts_dir is None:
ts_dir = CESM_output_dir
## - - - - - - - - - - - - - - - - - - - - - -
## Settings for computing coupling index
## - - - - - - - - - - - - - - - - - - - - - -
# Set up directory for aux output like coupling index file
if "SCRATCH" in os.environ:
cupid_temp = os.path.join(os.path.sep, os.environ["SCRATCH"], "CUPiD_scratch")
os.makedirs(cupid_temp, exist_ok=True)
else:
cupid_temp = "."
startYrs = [start_date.split("-")[0]]
endYrs = [f"{int(end_date.split('-')[0])-1:04d}"]
caseNames = [
case_name,
# base_case_name,
]
shortNames = [case.split(".")[-1] for case in caseNames]
2. Read in model data and compute terrestrial coupling index#
Execution using papermill encountered an exception here and stopped:
Soil moisture file not found!
Land-based SHFLX file not found!
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[7], line 60
58 # If years start at 0000, offset by 1700 years for analysis
59 yrOffset = 1850
---> 60 if shflx_DS["time.year"].values[0] < 1500:
61 shflx_DS["time"] = shflx_DS.time + timedelta(days=yrOffset * 365)
62 if soilWater_DS["time.year"].values[0] < 1500:
NameError: name 'shflx_DS' is not defined
2.1 Read in FLUXNET data if requested#
Make some choices on limiting which stations are used
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.
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.
3. Make plots#
def plotScatter(seasonstr, caseSel=None):
node_lats = uxgrid.face_lat.values
node_lons = uxgrid.face_lon.values
predictions = []
fig, axs = plt.subplots(1, 1, figsize=(8, 8))
iSeason = np.where(seasons == seasonstr)[0]
iStations = np.where(np.isfinite(terraCI_fluxnetConverted[:, iSeason]) == True)[0]
for iPoint in range(len(iStations)):
this_lon = lon_fluxnet[iStations[iPoint]] # lon1
this_lat = lat_fluxnet[iStations[iPoint]] # lat1
obs_point = np.array((this_lon, this_lat))
# Get subset of relevant points
i = np.where(
(node_lats >= (this_lat - 2))
& (node_lats <= (this_lat + 2))
& (node_lons >= (this_lon - 2))
& (node_lons <= (this_lon + 2))
)[0]
minDistance = 100
for iSelClose in range(len(i)):
# Find point in uxarray? Use euclidian distance
distance = np.linalg.norm(
obs_point - np.array((node_lons[i[iSelClose]], node_lats[i[iSelClose]]))
)
if (distance < minDistance) & (
np.isfinite(
couplingIndex_DS["CouplingIndex"]
.sel(season="JJA")
.values[i[iSelClose]]
)
== True
):
minDistance = distance
selLoc = i[iSelClose]
predictions = np.append(
predictions,
couplingIndex_DS["CouplingIndex"].sel(season=seasonstr).values[selLoc],
)
axs.plot(
terraCI_fluxnetConverted[iPoint, iSeason],
couplingIndex_DS["CouplingIndex"].sel(season=seasonstr).values[selLoc],
"bo",
alpha=0.5,
)
axs.set_xlabel("FLUXNET CI Value", fontsize=12)
axs.set_ylabel("CESM CI Value", fontsize=12)
axs.set_title(
"Individual station CI vs. nearest gridcell CI: " + seasonstr, fontsize=14
)
axs.set_xlim([-25, 25])
axs.set_ylim([-25, 25])
axs.plot(np.arange(-25, 26), np.arange(-25, 26), "k--")
rmse = np.sqrt(
((predictions - terraCI_fluxnetConverted[iStations, iSeason]) ** 2).mean()
)
axs.text(0.05, 0.95, "RMSE: " + str(rmse), transform=axs.transAxes)
return axs