Merge pull request #94 from dgergel/add_precip_unit_conversion

Add precip unit conversion

Author: brews

HISTORY.rst

@@ -6,6 +6,7 @@ History
 0.X.X (XXXX-XX-XX)
 ------------------
 * Add ``include-quantiles`` flag to ``apply_qdm`` to allow for including quantile information in bias corrected output. (PR #95, @dgergel)
+* Add precipitation unit conversion to ``standardize_gcm``. (PR #94, @dgergel)
 * Add ``astype`` argument to ``regrid``. (PR #92, @brews)
 * Make ``dodola`` container's default CMD. (PR #90, @brews)
 * Add ``dask-kubernetes``, ``distributed`` to container environment. (PR #90, @brews)

dodola/core.py

@@ -307,6 +307,20 @@ def standardize_gcm(ds, leapday_removal=True):
     )
 
     # Cleanup time.
+
+    # if variable is precip, need to update units to mm day-1
+    if "pr" in ds_cleaned.variables:
+        # units should be kg/m2/s in CMIP6 output
+        if ds_cleaned["pr"].units == "kg m-2 s-1":
+            # convert to mm/day
+            mmday_conversion = 24 * 60 * 60
+            ds_cleaned["pr"] = ds_cleaned["pr"] * mmday_conversion
+            # update units attribute
+            ds_cleaned["pr"].attrs["units"] = "mm day-1"
+        else:
+            # we want this to fail, as pr units are something we don't expect
+            raise ValueError("check units: pr units attribute is not kg m-2 s-1")
+
     if leapday_removal:
         # if calendar is just integers, xclim cannot understand it
         if ds.time.dtype == "int64":

dodola/tests/test_services.py

@@ -44,7 +44,7 @@ def _datafactory(x, start_time="1950-01-01", variable_name="fakevariable"):
     return out
 
 
-def _gcmfactory(x, start_time="1950-01-01"):
+def _gcmfactory(x, gcm_variable="fakevariable", start_time="1950-01-01"):
     """Populate xr.Dataset with synthetic GCM data for testing
     that includes extra dimensions and leap days to be removed.
     """
@@ -58,7 +58,7 @@ def _gcmfactory(x, start_time="1950-01-01"):
 
     out = xr.Dataset(
         {
-            "fakevariable": (
+            gcm_variable: (
                 ["time", "lon", "lat", "member_id"],
                 x[:, np.newaxis, np.newaxis, np.newaxis],
             )
@@ -74,6 +74,7 @@ def _gcmfactory(x, start_time="1950-01-01"):
             "time_bnds": (["time", "bnds"], np.ones((len(x), 2))),
         },
     )
+
     return out
 
 
@@ -462,6 +463,35 @@ def test_clean_cmip6():
     assert "time_bnds" not in ds_cleaned.coords
 
 
+@pytest.mark.parametrize(
+    "gcm_variable", [pytest.param("tasmax"), pytest.param("tasmin"), pytest.param("pr")]
+)
+def test_cmip6_precip_unitconversion(gcm_variable):
+    """Tests that precip units are converted in CMIP6 cleanup if variable is precip"""
+    # Setup input data
+    n = 1500  # need over four years of daily data
+    ts = np.sin(np.linspace(-10 * np.pi, 10 * np.pi, n)) * 0.5
+    ds_gcm = _gcmfactory(ts, gcm_variable=gcm_variable, start_time="1950-01-01")
+
+    if gcm_variable == "pr":
+        # assign units to typical GCM pr units so they can be cleaned
+        ds_gcm["pr"].attrs["units"] = "kg m-2 s-1"
+
+    in_url = "memory://test_clean_cmip6/an/input/path.zarr"
+    out_url = "memory://test_clean_cmip6/an/output/path.zarr"
+    repository.write(in_url, ds_gcm)
+
+    clean_cmip6(in_url, out_url, leapday_removal=True)
+    ds_cleaned = repository.read(out_url)
+
+    assert "height" not in ds_cleaned.coords
+    assert "member_id" not in ds_cleaned.coords
+    assert "time_bnds" not in ds_cleaned.coords
+
+    if "pr" in ds_cleaned.variables:
+        assert ds_cleaned["pr"].units == "mm day-1"
+
+
 def test_remove_leapdays():
     """Test that leapday removal service removes leap days"""
     # Setup input data