staged

pcmdi_metrics/monsoon_wang/__init__.py

@@ -0,0 +1,2 @@
+from .monsoon_precip_index_fncs import mpd, mpi_skill_scores  # noqa
+from .monsoon_wang_driver import create_monsoon_wang_parser, monsoon_wang_runner  # noqa

pcmdi_metrics/monsoon_wang/basic_monsoon_wang_param.py

@@ -0,0 +1,57 @@
+#
+#  OPTIONS ARE SET BY USER IN THIS FILE AS INDICATED BELOW BY:
+#
+#
+
+# LIST OF MODEL VERSIONS TO BE TESTED
+
+#  trouble model ICON-ESM-LR'
+
+modnames = [
+    "ACCESS-CM2",
+    "ACCESS-ESM1-5",
+    "AWI-ESM-1-1-LR",
+    "BCC-CSM2-MR",
+    "BCC-ESM1",
+    "CAMS-CSM1-0",
+    "CanESM5",
+    "CAS-ESM2-0",
+    "CESM2",
+    "CESM2-WACCM",
+    "CIESM",
+    "CMCC-ESM2",
+    "E3SM-1-0",
+    "E3SM-1-1",
+    "E3SM-2-0",
+    "EC-Earth3",
+    "FGOALS-f3-L",
+    "GFDL-CM4",
+    "GFDL-ESM4",
+    "GISS-E2-1-G",
+    "GISS-E2-1-H",
+    "IPSL-CM5A2-INCA",
+    "IPSL-CM6A-LR",
+    "KIOST-ESM",
+    "MIROC6",
+    "MPI-ESM1-2-LR",
+    "MRI-ESM2-0",
+    "NESM3",
+    "NorESM2-LM",
+    "TaiESM1",
+]
+
+
+# ROOT PATH FOR MODELS CLIMATOLOGIES
+test_data_path = "/p/user_pub/pmp/pmp_results/pmp_v1.1.2/diagnostic_results/CMIP_CLIMS/cmip6/historical/v20230823/pr/cmip6.historical.%(model).r1i1p1f1.mon.pr.198101-200512.AC.v20230823.nc"
+
+# ROOT PATH FOR OBSERVATIONS
+reference_data_path = "/p/user_pub/PCMDIobs/obs4MIPs_legacy/PCMDIobs2_clims/atmos/pr/TRMM-3B43v-7/pr_mon_TRMM-3B43v-7_BE_gn_199801-201712.v20200421.AC.nc"
+
+# DIRECTORY WHERE TO PUT RESULTS
+results_dir = "$OUTPUT_DIR$/monsoon_wang"
+
+# Threshold
+threshold = 2.5 / 86400
+
+# monsoon domain mask based on observations
+obs_mask = True

pcmdi_metrics/monsoon_wang/monsoon_precip_index_fncs.py

@@ -0,0 +1,153 @@
+from typing import Union
+
+import numpy as np
+import xarray as xr
+
+#  SEASONAL RANGE - USING ANNUAL CYCLE CLIMATOLGIES 0=Jan, 11=Dec
+
+
+def da_to_ds(d: Union[xr.Dataset, xr.DataArray], var: str = "variable") -> xr.Dataset:
+    """Convert xarray DataArray to Dataset
+
+    Parameters
+    ----------
+    d : Union[xr.Dataset, xr.DataArray]
+        Input dataArray. If dataset is given, no process will be done
+    var : str, optional
+        Name of dataArray, by default "variable"
+
+    Returns
+    -------
+    xr.Dataset
+        xarray Dataset
+
+    Raises
+    ------
+    TypeError
+        Raised when given input is not xarray based variables
+    """
+    if isinstance(d, xr.Dataset):
+        return d.copy()
+    elif isinstance(d, xr.DataArray):
+        return d.to_dataset(name=var).bounds.add_missing_bounds().copy()
+    else:
+        raise TypeError(
+            "Input must be an instance of either xarrary.DataArray or xarrary.Dataset"
+        )
+
+
+def regrid(da_in, da_grid, data_var="pr"):
+    ds_in = da_to_ds(da_in, data_var)
+    ds_grid = da_to_ds(da_grid, data_var)
+
+    ds_out = ds_in.regridder.horizontal(data_var, ds_grid, tool="regrid2")
+    da_out = ds_out[data_var]
+
+    return da_out
+
+
+def compute_season(data, season_indices, weights):
+    out = np.ma.zeros(data.shape[1:], dtype=data.dtype)
+    N = 0
+    for i in season_indices:
+        out += data[i] * weights[i]
+        N += weights[i]
+    return out / N
+
+
+def mpd(data):
+    """Monsoon precipitation intensity and annual range calculation
+    .. describe:: Input
+        *  data
+            * Assumes climatology array with 12 times step first one January
+    """
+    months_length = [
+        31.0,
+        28.0,
+        31.0,
+        30.0,
+        31.0,
+        30.0,
+        31.0,
+        31.0,
+        30.0,
+        31.0,
+        30.0,
+        31.0,
+    ]
+    mjjas = compute_season(data, [4, 5, 6, 7, 8], months_length)
+    ndjfm = compute_season(data, [10, 11, 0, 1, 2], months_length)
+    ann = compute_season(data, list(range(12)), months_length)
+
+    data_map = data.isel(time=0)
+
+    annrange = mjjas - ndjfm
+
+    da_annrange = xr.DataArray(annrange, coords=data_map.coords, dims=data_map.dims)
+    da_annrange = da_annrange.where(da_annrange.lat >= 0, da_annrange * -1)
+
+    mpi = np.divide(da_annrange.values, ann, where=ann.astype(bool))
+
+    da_mpi = xr.DataArray(mpi, coords=data_map.coords, dims=data_map.dims)
+
+    return da_annrange, da_mpi
+
+
+def mpi_skill_scores(annrange_mod_dom, annrange_obs_dom, threshold=2.5 / 86400.0):
+    """Monsoon precipitation index skill score calculation
+    see Wang et al., doi:10.1007/s00382-010-0877-0
+
+      .. describe:: Input
+
+          *  annrange_mod_dom
+
+              * Model Values Range (summer - winter)
+
+          *  annrange_obs_dom
+
+              * Observations Values Range (summer - winter)
+
+          *  threshold [default is 2.5/86400.]
+
+              * threshold in same units as inputs
+    """
+    mt = np.ma.greater(annrange_mod_dom, threshold)
+    ot = np.ma.greater(annrange_obs_dom, threshold)
+
+    hitmap = mt * ot  # only where both  mt and ot are True
+    hit = float(hitmap.sum())
+
+    xor = np.ma.logical_xor(mt, ot)
+    missmap = xor * ot
+    missed = float(missmap.sum())
+
+    falarmmap = xor * mt
+    falarm = float(falarmmap.sum())
+
+    if (hit + missed + falarm) > 0.0:
+        score = hit / (hit + missed + falarm)
+    else:
+        score = 1.0e20
+
+    xr_hitmap = xr.DataArray(
+        hitmap,
+        name="hitmap",
+        coords={"lat": annrange_mod_dom.lat, "lon": annrange_mod_dom.lon},
+        dims=["lat", "lon"],
+    )
+
+    xr_missmap = xr.DataArray(
+        missmap,
+        name="missmap",
+        coords={"lat": annrange_mod_dom.lat, "lon": annrange_mod_dom.lon},
+        dims=["lat", "lon"],
+    )
+
+    xr_falarmmap = xr.DataArray(
+        falarmmap,
+        name="falarmmap",
+        coords={"lat": annrange_mod_dom.lat, "lon": annrange_mod_dom.lon},
+        dims=["lat", "lon"],
+    )
+
+    return hit, missed, falarm, score, xr_hitmap, xr_missmap, xr_falarmmap