From 8703cd92faba0ddf5643e72ee46e4e6bc203f8f3 Mon Sep 17 00:00:00 2001 From: Nathan Collier Date: Thu, 25 May 2023 11:48:22 -0400 Subject: [PATCH 1/4] added initial confrontatoin --- src/ILAMB/ConfTerrestrialCoupling.py | 69 ++++++++++++++++++++++++++++ src/ILAMB/Scoreboard.py | 4 +- 2 files changed, 72 insertions(+), 1 deletion(-) create mode 100644 src/ILAMB/ConfTerrestrialCoupling.py diff --git a/src/ILAMB/ConfTerrestrialCoupling.py b/src/ILAMB/ConfTerrestrialCoupling.py new file mode 100644 index 00000000..702d387f --- /dev/null +++ b/src/ILAMB/ConfTerrestrialCoupling.py @@ -0,0 +1,69 @@ +""".""" +from copy import deepcopy + +import numpy as np +import xarray as xr + +from .Confrontation import Confrontation +from .Variable import Variable + + +def compute_coupling_index(dset: xr.Dataset, control: str, respond: str) -> xr.Dataset: + """Return the dataset with the coupling index added. + + Parameters + ---------- + dset + The dataset containing DataArray's labeled `control` and `respond` + control + The name of the controlling variable + respond + The name of the responding variable, the index will be in the units of + this variable. + """ + assert control in dset + assert respond in dset + ssn = dset.sel(time=dset["time.season"] == "JJA") + dset["ci"] = xr.cov(ssn[control], ssn[respond], dim="time") / ssn[control].std( + dim="time" + ) + dset["ci"].attrs["units"] = dset[respond].attrs["units"] + return dset + + +class ConfTerrestrialCoupling(Confrontation): + def stageData(self, m): + # eventually we need mrsos_source and hfls_source + obs_ds = xr.open_dataset(self.source) + obs_ds = compute_coupling_index(obs_ds, "mrsos", "hfls") + + # for backwards compatibility, now convert to ILAMB object + tmin = obs_ds["time"].min() + tmax = obs_ds["time"].max() + tbnds = np.asarray( + [ + [ + ((tmin.dt.year - 1850) * 365) + tmin.dt.dayofyear, + ((tmax.dt.year - 1850) * 365) + tmax.dt.dayofyear, + ] + ], + dtype=float, + ) + obs = Variable( + name="ci", + data=np.ma.masked_array( + obs_ds["ci"].to_numpy(), mask=np.zeros_like(obs_ds["ci"]) + ).reshape((1, -1)), + unit=obs_ds["ci"].attrs["units"], + ndata=len(obs_ds["SITE"]), + lat=obs_ds["lat"].to_numpy(), + lon=obs_ds["lon"].to_numpy(), + time=tbnds.mean(axis=1), + time_bnds=tbnds, + ) + + # just for testing + mod = deepcopy(obs) + mod.data += obs.data.max() * 0.1 * (np.random.rand(*mod.data.shape) - 0.5) + + return obs, mod diff --git a/src/ILAMB/Scoreboard.py b/src/ILAMB/Scoreboard.py index 5cdb3f3c..6d4290a0 100644 --- a/src/ILAMB/Scoreboard.py +++ b/src/ILAMB/Scoreboard.py @@ -13,6 +13,7 @@ from .ConfSoilCarbon import ConfSoilCarbon from .ConfUncertainty import ConfUncertainty from .ConfBurntArea import ConfBurntArea +from .ConfTerrestrialCoupling import ConfTerrestrialCoupling try: from .ConfUSGS import ConfUSGS except: @@ -328,7 +329,8 @@ def _loadScores(node): "ConfSoilCarbon" : ConfSoilCarbon, "ConfUncertainty" : ConfUncertainty, "ConfBurntArea" : ConfBurntArea, - "ConfUSGS" : ConfUSGS } + "ConfUSGS" : ConfUSGS , + "ConfTerrestrialCoupling" : ConfTerrestrialCoupling } class Scoreboard(): """ From e1ec59fbf9ec2db89ae8338f6c534700a839e350 Mon Sep 17 00:00:00 2001 From: Nathan Collier Date: Thu, 25 May 2023 11:48:37 -0400 Subject: [PATCH 2/4] small ilamblib fixes --- src/ILAMB/ilamblib.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/ILAMB/ilamblib.py b/src/ILAMB/ilamblib.py index 8b6c7a8c..7f37e6eb 100644 --- a/src/ILAMB/ilamblib.py +++ b/src/ILAMB/ilamblib.py @@ -912,7 +912,7 @@ def Score(var,normalizer): name = name.replace("iav" ,"iav_score") np.seterr(over='ignore',under='ignore') data = np.exp(-np.abs(var.data/normalizer.data)) - data[data<1e-16] = 0. + data = (data<1e-16)*0 + (data>=1e-16)*data np.seterr(over='raise',under='raise') return Variable(name = name, data = data, @@ -1100,7 +1100,7 @@ def AnalysisMeanStateSites(ref,com,**keywords): else: msg = f"[{name}] Bias scored using Collier2018" logger.info(msg) - bias_score_map = Score(bias,crms) + bias_score_map = Score(bias,crms if REF.time.size > 1 else REF_timeint) if not skip_rmse: cCOM = Variable(name = "centralized %s" % COM.name, unit = COM.unit, From a7c94c29dc75868787adf0d728fd3aa3d4db9821 Mon Sep 17 00:00:00 2001 From: Nathan Collier Date: Fri, 26 May 2023 11:41:19 -0400 Subject: [PATCH 3/4] readin model results --- src/ILAMB/ConfTerrestrialCoupling.py | 47 +++++++++++++++++++++++++--- 1 file changed, 43 insertions(+), 4 deletions(-) diff --git a/src/ILAMB/ConfTerrestrialCoupling.py b/src/ILAMB/ConfTerrestrialCoupling.py index 702d387f..c4e60fe2 100644 --- a/src/ILAMB/ConfTerrestrialCoupling.py +++ b/src/ILAMB/ConfTerrestrialCoupling.py @@ -62,8 +62,47 @@ def stageData(self, m): time_bnds=tbnds, ) - # just for testing - mod = deepcopy(obs) - mod.data += obs.data.max() * 0.1 * (np.random.rand(*mod.data.shape) - 0.5) - + # load the model result + mod_ds = {} + units = {} + lat = xr.DataArray(obs.lat, dims="SITE") + lon = xr.DataArray(obs.lon, dims="SITE") + for vname in ["mrsos", "hfls"]: + var = xr.open_mfdataset(m.variables[vname]) + cal = var["time"].values[0].__class__ + t0 = tmin.dt.day + tf = tmax.dt.day + if cal.__name__ == "Datetime360Day": + if t0 == 31: + t0 = 30 + if tf == 31: + tf = 30 + units[vname] = var[vname].attrs["units"] + var = var.sel( + time=slice( + cal(tmin.dt.year, tmin.dt.month, t0), + cal(tmax.dt.year, tmax.dt.month, tf), + ), + ) + var = var.sel(lat=lat, lon=lon, method="nearest") + var.load() + mod_ds[vname] = (["time", "SITE"], var[vname].data) + mod_ds = xr.Dataset( + data_vars=mod_ds, + coords={"time": var["time"].data, "SITE": var["SITE"].data}, + ) + for var, unit in units.items(): + mod_ds[var].attrs["units"] = unit + mod_ds = compute_coupling_index(mod_ds, "mrsos", "hfls") + data = np.ma.masked_invalid(mod_ds["ci"].to_numpy()) + mod = Variable( + name="ci", + data=data.reshape((1,) + data.shape), + unit=mod_ds["ci"].attrs["units"], + lat=lat.to_numpy(), + lon=lon.to_numpy(), + ndata=len(lat), + time=tbnds.mean(axis=1), + time_bnds=tbnds, + ) return obs, mod From 93a0e3fc086a74cb5fc693bc4317df7e961babc9 Mon Sep 17 00:00:00 2001 From: Nathan Collier Date: Wed, 31 May 2023 22:48:03 -0400 Subject: [PATCH 4/4] working checkin --- src/ILAMB/ConfTerrestrialCoupling.py | 103 ++++++++++++++++++--------- 1 file changed, 71 insertions(+), 32 deletions(-) diff --git a/src/ILAMB/ConfTerrestrialCoupling.py b/src/ILAMB/ConfTerrestrialCoupling.py index c4e60fe2..71edd5e5 100644 --- a/src/ILAMB/ConfTerrestrialCoupling.py +++ b/src/ILAMB/ConfTerrestrialCoupling.py @@ -1,45 +1,81 @@ """.""" -from copy import deepcopy +import os import numpy as np import xarray as xr +from typing import Tuple +import datetime from .Confrontation import Confrontation from .Variable import Variable -def compute_coupling_index(dset: xr.Dataset, control: str, respond: str) -> xr.Dataset: - """Return the dataset with the coupling index added. +def coupling_index(control: xr.DataArray, respond: xr.DataArray) -> xr.DataArray: + """Return the coupling index in units of the responding variable. Parameters ---------- - dset - The dataset containing DataArray's labeled `control` and `respond` control - The name of the controlling variable + The controlling variable respond - The name of the responding variable, the index will be in the units of + The responding variable, the index will be in the units of this variable. + """ - assert control in dset - assert respond in dset - ssn = dset.sel(time=dset["time.season"] == "JJA") - dset["ci"] = xr.cov(ssn[control], ssn[respond], dim="time") / ssn[control].std( - dim="time" - ) - dset["ci"].attrs["units"] = dset[respond].attrs["units"] - return dset + control, response = xr.align(control, respond, join="override") + control = control.sel(time=control["time.season"] == "JJA") + respond = respond.sel(time=respond["time.season"] == "JJA") + cov = xr.cov(control, respond, dim="time") + std = control.std(dim="time") + cov.load() + std.load() + coupling = cov/std + coupling.attrs = {"long_name": "Coupling index", "units": respond.attrs["units"]} + coupling.load() + return coupling + + +def max_time_bounds(*dss): + """.""" + t0 = [] + tf = [] + for ds in dss: + if "time" not in ds: + continue + time = ds["time"].attrs["bounds"] if "bounds" in ds["time"].attrs else "time" + time = ds[time] if time in ds else ds["time"] + t0.append(time.min()) + tf.append(time.max()) + t0 = max(t0) + tf = min(tf) + return t0, tf class ConfTerrestrialCoupling(Confrontation): + def __init__(self, **keywords): + super(ConfTerrestrialCoupling, self).__init__(**keywords) + for srcname in ["mrsos_source", "hfss_source"]: + src = self.keywords.get(srcname, None) + if src is None: + continue + self.keywords[srcname] = os.path.join(os.environ["ILAMB_ROOT"], src) + def stageData(self, m): - # eventually we need mrsos_source and hfls_source - obs_ds = xr.open_dataset(self.source) - obs_ds = compute_coupling_index(obs_ds, "mrsos", "hfls") + + # read in reference data and find maximal overlap + mrsos_obs = xr.open_dataset( + self.keywords.get("mrsos_source", self.source), chunks=dict(time=1800) + ) + hfss_obs = xr.open_dataset( + self.keywords.get("hfss_source", self.source), chunks=dict(time=1800) + ) + tmin, tmax = max_time_bounds(mrsos_obs, hfss_obs) + if len(mrsos_obs['time']) != len(hfss_obs['time']): + mrsos_obs = mrsos_obs.sel(time=slice(tmin, tmax)) + hfss_obs = hfss_obs.sel(time=slice(tmin, tmax)) + ci_obs = coupling_index(mrsos_obs["mrsos"], hfss_obs["hfss"]) # for backwards compatibility, now convert to ILAMB object - tmin = obs_ds["time"].min() - tmax = obs_ds["time"].max() tbnds = np.asarray( [ [ @@ -49,25 +85,28 @@ def stageData(self, m): ], dtype=float, ) + ndata = len(ci_obs['SITE']) if 'SITE' in ci_obs.dims else None + data = np.ma.masked_invalid(ci_obs.to_numpy()) + if len(data.mask)==1: data.mask = np.zeros_like(data) + data = data.reshape((1,)+data.shape) obs = Variable( name="ci", - data=np.ma.masked_array( - obs_ds["ci"].to_numpy(), mask=np.zeros_like(obs_ds["ci"]) - ).reshape((1, -1)), - unit=obs_ds["ci"].attrs["units"], - ndata=len(obs_ds["SITE"]), - lat=obs_ds["lat"].to_numpy(), - lon=obs_ds["lon"].to_numpy(), + data=data, + unit=ci_obs.attrs["units"], + ndata=ndata, + lat=mrsos_obs["lat"].to_numpy(), + lon=mrsos_obs["lon"].to_numpy(), time=tbnds.mean(axis=1), time_bnds=tbnds, ) - + print(obs) + # load the model result mod_ds = {} units = {} lat = xr.DataArray(obs.lat, dims="SITE") lon = xr.DataArray(obs.lon, dims="SITE") - for vname in ["mrsos", "hfls"]: + for vname in ["mrsos", "hfss"]: var = xr.open_mfdataset(m.variables[vname]) cal = var["time"].values[0].__class__ t0 = tmin.dt.day @@ -93,12 +132,12 @@ def stageData(self, m): ) for var, unit in units.items(): mod_ds[var].attrs["units"] = unit - mod_ds = compute_coupling_index(mod_ds, "mrsos", "hfls") - data = np.ma.masked_invalid(mod_ds["ci"].to_numpy()) + ci_mod = coupling_index(mod_ds["mrsos"], mod_ds["hfss"]) + data = np.ma.masked_invalid(ci_mod.to_numpy()) mod = Variable( name="ci", data=data.reshape((1,) + data.shape), - unit=mod_ds["ci"].attrs["units"], + unit=ci_mod.attrs["units"], lat=lat.to_numpy(), lon=lon.to_numpy(), ndata=len(lat),