Merge pull request #97 from rubisco-sfa/alt_calm

ADD: CALM active layer thickness

Author: nocollier

ilamb3/configure/ilamb.yaml

@@ -323,6 +323,21 @@ Cryosphere:
       analyses:
       - Area Comparison
 
+  Active Layer Thickness:
+
+    CALM:
+      sources:
+        active_layer_thickness: active_layer_thickness/CALM/CALM_1990-2024.nc
+      transforms:
+      - select_time:
+          vmin: 1990-01
+          vmax: 2023-01
+      - select_lat:
+          vmin: 30.
+          vmax: 90.
+      - active_layer_thickness:
+          method: slater2013
+
 Radiation and Energy Cycle:
 
   Surface Upward SW Radiation:

ilamb3/dataset.py

@@ -3,6 +3,7 @@
 import warnings
 from typing import Any, Literal
 
+import cftime as cf
 import numpy as np
 import xarray as xr
 from scipy.interpolate import NearestNDInterpolator
@@ -1050,3 +1051,30 @@ def shift_time_by_years(ds: xr.Dataset, years: int) -> xr.Dataset:
         for t in ds[time_dim]
     ]
     return ds
+
+
+def convert_year_to_datetime(ds: xr.Dataset) -> xr.Dataset:
+    """
+    Return the dataset with a integer 'years' dimension converted to a
+    CF-compliant time.
+    """
+    if "year" not in ds.coords:
+        raise ValueError("Function expects a 'year' coordinate with integer years.")
+    # We are about to reassign the 'time' dimension, let's remove everything
+    # with 'time' in it. We are assuming that it is no longer relevant.
+    if "time" in ds.coords:
+        ds = ds.drop_vars([v for v in ds if "time" in ds[v].dims]).drop_dims("time")
+    years = ds["year"]
+    ds = ds.rename_dims(year="time").rename_vars(year="time")
+    ds["time"] = [cf.DatetimeNoLeap(y, 7, 1) for y in years]
+    ds["time_bnds"] = (
+        ("time", "nb"),
+        np.array(
+            [
+                [cf.DatetimeNoLeap(y, 1, 1) for y in years],
+                [cf.DatetimeNoLeap(y + 1, 1, 1) for y in years],
+            ]
+        ).T,
+    )
+    ds["time"].attrs["bounds"] = "time_bnds"
+    return ds

ilamb3/registry/ilamb.txt

@@ -1,4 +1,4 @@
-active_layer_thickness/CALM/CALM.nc sha1:fe084ab58a81911e3739c8bb3fad210b3863af27
+active_layer_thickness/CALM/CALM_1990-2024.nc sha1:d874d6efb07ceeb135c30607729653cf0d0c6651
 biomass/ESACCI/biomass.nc sha1:2eeac984738fc17e8b6cdf0564fff192a4a08718
 biomass/GEOCARBON/biomass.nc sha1:369235319774f273e376a9e7e3ffd665e5d1c83d
 biomass/NBCD2000/biomass_0.5x0.5.nc sha1:6354c20dfb61f97454ba0fa70023ee3e959c0c77

ilamb3/tests/test_transform.py

@@ -41,6 +41,38 @@ def gen_msftmz_dset(seed: int = 1):
     )
 
 
+def gen_permafrost_dset(seed: int = 1):
+    rs = np.random.RandomState(seed)
+    coords = {}
+    # time, depth, lat, lon
+    coords["time"] = [
+        cf.DatetimeNoLeap(2000 + int(m / 12), (m % 12) + 1, 15) for m in range(3 * 12)
+    ]
+    coords["depth"] = 0.5 * (
+        np.linspace(0, 5, 10 + 1)[1:] + np.linspace(0, 5, 10 + 1)[:-1]
+    )
+    coords["lat"] = 0.5 * (
+        np.linspace(60, 90, 10 + 1)[1:] + np.linspace(60, 90, 10 + 1)[:-1]
+    )
+    coords["lon"] = 0.5 * (
+        np.linspace(-180, 180, 20 + 1)[1:] + np.linspace(-180, 180, 20 + 1)[:-1]
+    )
+    dims = [c for c in coords]
+    ds = xr.Dataset(
+        data_vars={
+            "tsl": xr.DataArray(
+                rs.rand(*[len(coords[d]) for d in dims]) * 20 - 18,
+                coords=coords,
+                dims=dims,
+                name="tsl",
+                attrs={"units": "degC"},
+            ),
+        }
+    )
+    ds = ds.cf.add_bounds("depth")
+    return ds
+
+
 def gen_expression_dset(
     expr: str,
     var_meta: dict[str, dict[str, float | str]],
@@ -100,6 +132,7 @@ def gen_expression_dset(
     "depth_gradient": generate_test_dset(
         "thetao", "K", nyear=1, nlat=2, nlon=4, ndepth=10
     ),
+    "active_layer_thickness": gen_permafrost_dset(),
 }
 
 
@@ -112,12 +145,17 @@ def gen_expression_dset(
         ("select_depth", {"value": 0}, "thetao", 9.43861150676275),
         ("select_depth", {"vmin": 1, "vmax": 40}, "thetao", 9.983843647875275),
         ("depth_gradient", {}, "thetao", -0.005550578833866464),
+        (
+            "active_layer_thickness",
+            {"method": "slater2013"},
+            "active_layer_thickness",
+            2.224852071005917,
+        ),
     ],
 )
 def test_transform(name, kwargs, out, value):
     transform = ALL_TRANSFORMS[name](**kwargs)
     ds = transform(DATA[name])
-    print(ds[out].mean().values)  # this value should be the correct one
     assert np.allclose(value, ds[out].mean().values)
 
 

ilamb3/transform/__init__.py

@@ -2,13 +2,14 @@
 from ilamb3.transform.expression import expression
 from ilamb3.transform.gradient import depth_gradient
 from ilamb3.transform.ohc import ocean_heat_content
-from ilamb3.transform.permafrost import permafrost_extent
+from ilamb3.transform.permafrost import active_layer_thickness, permafrost_extent
 from ilamb3.transform.runoff_sensitivity import runoff_sensitivity
 from ilamb3.transform.select import select_depth, select_lat, select_lon, select_time
 from ilamb3.transform.soilmoisture import soil_moisture_to_vol_fraction
 from ilamb3.transform.stratification_index import stratification_index
 
 ALL_TRANSFORMS = {
+    "active_layer_thickness": active_layer_thickness,
     "depth_gradient": depth_gradient,
     "expression": expression,
     "msftmz_to_rapid": msftmz_to_rapid,

ilamb3/transform/permafrost.py

@@ -33,6 +33,10 @@ def _permafrost_extent_slater2013(
     ds["permafrost_extent"] = (frozen & frozen.shift(year=-1, fill_value=False)).any(
         dim=depth_dim
     )
+    ds["permafrost_extent"].attrs = {
+        "long_name": "Permafrost extent based on slater2013",
+        "units": "1",
+    }
 
     # the active layer is then the non-frozen thickness over permafrosted areas
     depth_bnds = ds[ds.cf.bounds[depth_dim][0]]
@@ -42,19 +46,29 @@ def _permafrost_extent_slater2013(
     ds["active_layer_thickness"] = xr.where(
         ds["permafrost_extent"], (~frozen * layer_thickness).sum(dim=depth_dim), np.nan
     )
+    ds["active_layer_thickness"].attrs = {
+        "long_name": "Active layer thickness based on slater2013",
+        "units": (
+            ds[depth_dim].attrs["units"] if "units" in ds[depth_dim].attrs else "m"
+        ),
+    }
 
     return ds
 
 
 EXTENT_METHODS = {"slater2013": _permafrost_extent_slater2013}
+ALT_METHODS = {"slater2013": _permafrost_extent_slater2013}
 
 
 class permafrost_extent(ILAMBTransform):
     """
-    An ILAMB transform for .
+    An ILAMB transform for estimating permafrost from layered soil temperature
+    (`tsl`).
 
     Parameters
     ----------
+    method: str
+        The method to use to estimate extent from `tsl`.
     """
 
     def __init__(self, method: Literal["slater2013"] = "slater2013", **kwargs: Any):
@@ -75,12 +89,43 @@ def __call__(self, ds: xr.Dataset) -> xr.Dataset:
         if not set(self.required_variables()).issubset(ds):
             return ds
         ds = EXTENT_METHODS[self.method](ds)
-        ds["permafrost_extent"].attrs = {
-            "long_name": f"Permafrost extent based on {self.method}",
-            "units": "1",
-        }
         ds = ds["permafrost_extent"].any(dim="year").squeeze().to_dataset()
         ds["permafrost_extent"] = xr.where(
             ds["permafrost_extent"] > 0, ds["permafrost_extent"], np.nan
         )
+        ds = ds.drop_vars("tsl")
+        return ds
+
+
+class active_layer_thickness(ILAMBTransform):
+    """
+    An ILAMB transform for estimating active layer thickness from layered soil
+    temperature (`tsl`).
+
+    Parameters
+    ----------
+    method: str
+        The method to use to estimate active layer thickness from `tsl`.
+    """
+
+    def __init__(self, method: Literal["slater2013"] = "slater2013", **kwargs: Any):
+        self.method = method
+
+    def required_variables(self) -> list[str]:
+        """
+        Return the required variables for the transform.
+        """
+        return ["tsl"]
+
+    def __call__(self, ds: xr.Dataset) -> xr.Dataset:
+        """
+        Compute
+        """
+        if "active_layer_thickness" in ds:
+            return ds
+        if not set(self.required_variables()).issubset(ds):
+            return ds
+        ds = ALT_METHODS[self.method](ds)
+        ds = ds.drop_vars("tsl")
+        ds = dset.convert_year_to_datetime(ds)
         return ds