Merge pull request #1022 from PCMDI/new_monsoon_lee1043_clean_up

 New monsoon sperber code in xCDAT -- PR cleaned up

Author: lee1043

doc/jupyter/Demo/Demo_2b_monsoon_sperber.ipynb

@@ -2,3 +2,4 @@
 from .calc_metrics import sperber_metrics  # noqa
 from .divide_chunks import divide_chunks, divide_chunks_advanced, interp1d  # noqa
 from .model_land_only import model_land_only  # noqa
+from .lib_monsoon_sperber import pick_year_last_day, tree

pcmdi_metrics/monsoon_sperber/driver_monsoon_sperber.py

@@ -53,7 +53,10 @@ def AddParserArgument(P):
     P.add_argument(
         "--meyear", dest="meyear", type=int, help="End year for model data set"
     )
-    P.add_argument("--modnames", type=list, default=None, help="List of models")
+    P.add_argument("--modnames", type=str, default=None, help="List of models")
+    P.add_argument(
+        "--list_monsoon_regions", type=str, default=None, help="List of regions"
+    )
     P.add_argument(
         "-r",
         "--realization",
@@ -95,6 +98,23 @@ def AddParserArgument(P):
         default=True,
         help="Option for update existing JSON file: True (i.e., update) (default) / False (i.e., overwrite)",
     )
+    # CMEC
+    P.add_argument(
+        "--cmec",
+        dest="cmec",
+        default=False,
+        action="store_true",
+        help="Use to save CMEC format metrics JSON",
+    )
+    P.add_argument(
+        "--no_cmec",
+        dest="cmec",
+        default=False,
+        action="store_false",
+        help="Do not save CMEC format metrics JSON",
+    )
+    P.set_defaults(cmec=False)
+
     return P
 
 

pcmdi_metrics/monsoon_sperber/lib/__init__.py

@@ -8,34 +8,45 @@
 
 Drafted: Jiwoo Lee, 2018-07
 Revised: Jiwoo Lee, 2019-05
+Revised: Bo Dong, 2023-12
 
 Note: Code for picking onset/decay index inspired by
 https://stackoverflow.com/questions/2236906/first-python-list-index-greater-than-x
 """
 
-import MV2
-
 
 def sperber_metrics(d, region, debug=False):
     """d: input, 1d array of cumulative pentad time series"""
     # Convert accumulation to fractional accumulation; normalize by sum
     d_sum = d[-1]
+
     # Normalize
-    frac_accum = MV2.divide(d, d_sum)
+    frac_accum = d / d_sum
+
     # Stat 1: Onset
     onset_index = next(i for i, v in enumerate(frac_accum) if v >= 0.2)
+    i = onset_index
+    v = frac_accum[i]
+
+    if debug:
+        print("i = , ", i, "  v =  ", v)
+
     # Stat 2: Decay
     if region == "GoG":
         decay_threshold = 0.6
     else:
         decay_threshold = 0.8
+
     decay_index = next(i for i, v in enumerate(frac_accum) if v >= decay_threshold)
+
     # Stat 3: Slope
     slope = (frac_accum[decay_index] - frac_accum[onset_index]) / float(
         decay_index - onset_index
     )
+
     # Stat 4: Duration
     duration = decay_index - onset_index + 1
+
     # Calc done, return result as dic
     return {
         "frac_accum": frac_accum,

pcmdi_metrics/monsoon_sperber/lib/argparse_functions.py

@@ -1,5 +1,3 @@
-from __future__ import print_function
-
 import sys
 
 import numpy as np
@@ -14,7 +12,7 @@
 
 def divide_chunks(data, n):
     # looping till length data
-    for i in range(0, len(data), n):
+    for i in range(0, data.time.shape[0], n):
         yield data[i : i + n]
 
 
@@ -24,43 +22,51 @@ def divide_chunks(data, n):
 
 def divide_chunks_advanced(data, n, debug=False):
     # Double check first date should be Jan 1 (except for SH monsoon)
-    tim = data.getTime()
-    calendar = tim.calendar
-    month = tim.asComponentTime()[0].month
-    day = tim.asComponentTime()[0].day
+
+    tim = data.time.dt
+    month = tim.month[0]
+    day = tim.day[0]
+    month = month.values
+    day = day.values
+    calendar = "gregorian"
+
     if debug:
         print("debug: first day of year is " + str(month) + "/" + str(day))
+
     if month not in [1, 7] or day != 1:
         sys.exit(
             "error: first day of year time series is " + str(month) + "/" + str(day)
         )
 
     # Check number of days in given year
-    nday = len(data)
+    nday = data.time.shape[0]
 
     if nday in [365, 360]:
         # looping till length data
         for i in range(0, nday, n):
             yield data[i : i + n]
+
     elif nday == 366:
         # until leap year day detected
         for i in range(0, nday, n):
             # Check if leap year date included
             leap_detect = False
             for ii in range(i, i + n):
-                date = data.getTime().asComponentTime()[ii]
-                month = date.month
-                day = date.day
+                date = data.time.dt
+                month = date.month[ii]
+                day = date.day[ii]
                 if month == 2 and day > 28:
                     if debug:
                         print("debug: leap year detected:", month, "/", day)
                     leap_detect = True
+
             if leap_detect:
                 yield data[i : i + n + 1]
                 tmp = i + n + 1
                 break
             else:
                 yield data[i : i + n]
+
         # after leap year day passed
         if leap_detect:
             for i in range(tmp, nday, n):
@@ -76,6 +82,7 @@ def divide_chunks_advanced(data, n, debug=False):
         # looping till length data
         for i in range(0, nday, n):
             yield data[i : i + n]
+
     else:
         sys.exit("error: number of days in year is " + str(nday))
 

pcmdi_metrics/monsoon_sperber/lib/calc_metrics.py

@@ -0,0 +1,22 @@
+from collections import defaultdict
+
+import xcdat as xc
+
+
+def tree():
+    return defaultdict(tree)
+
+
+def pick_year_last_day(ds):
+    eday = 31
+    try:
+        time_key = xc.axis.get_dim_keys(ds, axis="T")
+        if "calendar" in ds[time_key].attrs.keys():
+            if "360" in ds[time_key]["calendar"]:
+                eday = 30
+        else:
+            if "360" in ds[time_key][0].values.item().calendar:
+                eday = 30
+    except Exception:
+        pass
+    return eday

pcmdi_metrics/monsoon_sperber/lib/divide_chunks.py

@@ -1,80 +1,47 @@
 import cartopy.crs as ccrs
-import genutil
 import matplotlib.pyplot as plt
-import MV2
+import numpy as np
 
 
 def model_land_only(model, model_timeseries, lf, debug=False):
     # -------------------------------------------------
     # Mask out over ocean grid
     # - - - - - - - - - - - - - - - - - - - - - - - - -
+
     if debug:
-        plot_map(model_timeseries[0], "_".join(["test", model, "beforeMask.png"]))
+        # plot_map(model_timeseries[0], "_".join(["test", model, "beforeMask.png"]))
         print("debug: plot for beforeMask done")
 
     # Check land fraction variable to see if it meet criteria
     # (0 for ocean, 100 for land, no missing value)
-    lat_c = lf.getAxis(0)
-    lon_c = lf.getAxis(1)
-    lf_id = lf.id
-
-    lf = MV2.array(lf.filled(0.0))
-
-    lf.setAxis(0, lat_c)
-    lf.setAxis(1, lon_c)
-    lf.id = lf_id
 
-    if float(MV2.max(lf)) == 1.0:
-        lf = MV2.multiply(lf, 100.0)
-
-    # Matching dimension
-    if debug:
-        print("debug: match dimension in model_land_only")
-    model_timeseries, lf_timeConst = genutil.grower(model_timeseries, lf)
-
-    # Conserve axes
-    time_c = model_timeseries.getAxis(0)
-    lat_c2 = model_timeseries.getAxis(1)
-    lon_c2 = model_timeseries.getAxis(2)
+    if np.max(lf) == 1.0:
+        lf = lf * 100.0
 
     opt1 = False
 
     if opt1:  # Masking out partial ocean grids as well
         # Mask out ocean even fractional (leave only pure ocean grid)
-        model_timeseries_masked = MV2.masked_where(lf_timeConst < 100, model_timeseries)
+        model_timeseries_masked = model_timeseries.where(lf > 0 & lf < 100)
+
     else:  # Mask out only full ocean grid & use weighting for partial ocean grid
-        model_timeseries_masked = MV2.masked_where(
-            lf_timeConst == 0, model_timeseries
-        )  # mask out pure ocean grids
+        model_timeseries_masked = model_timeseries.where(lf > 0)
+
         if model == "EC-EARTH":
             # Mask out over 90% land grids for models those consider river as
             # part of land-sea fraction. So far only 'EC-EARTH' does..
-            model_timeseries_masked = MV2.masked_where(
-                lf_timeConst < 90, model_timeseries
-            )
-        lf2 = MV2.divide(lf, 100.0)
-        model_timeseries, lf2_timeConst = genutil.grower(
-            model_timeseries, lf2
-        )  # Matching dimension
-        model_timeseries_masked = MV2.multiply(
-            model_timeseries_masked, lf2_timeConst
-        )  # consider land fraction like as weighting
-
-    # Make sure to have consistent axes
-    model_timeseries_masked.setAxis(0, time_c)
-    model_timeseries_masked.setAxis(1, lat_c2)
-    model_timeseries_masked.setAxis(2, lon_c2)
+            model_timeseries_masked = model_timeseries.where(lf > 90)
 
     if debug:
-        plot_map(model_timeseries_masked[0], "_".join(["test", model, "afterMask.png"]))
+        # plot_map(model_timeseries_masked[0], "_".join(["test", model, "afterMask.png"]))
         print("debug: plot for afterMask done")
 
     return model_timeseries_masked
 
 
 def plot_map(data, filename):
-    lons = data.getLongitude()
-    lats = data.getLatitude()
+    lons = data["lon"]
+    lats = data["lat"]
     ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=180))
     ax.contourf(lons, lats, data, transform=ccrs.PlateCarree(), cmap="viridis")
     ax.coastlines()

pcmdi_metrics/monsoon_sperber/lib/lib_monsoon_sperber.py

@@ -0,0 +1,76 @@
+import datetime
+import os
+
+# =================================================
+# Background Information
+# -------------------------------------------------
+mip = "cmip5"
+exp = "historical"
+frequency = "da"
+realm = "atm"
+
+# =================================================
+# Miscellaneous
+# -------------------------------------------------
+update_json = False
+debug = False
+# debug = True
+
+# list_monsoon_regions = ["AIR", "AUS", "Sahel", "GoG", "NAmo", "SAmo"]
+list_monsoon_regions = ["AUS"]
+# =================================================
+# Observation
+# -------------------------------------------------
+reference_data_name = "GPCP-1-3"
+reference_data_path = "/p/user_pub/PCMDIobs/obs4MIPs/NASA-GSFC/GPCP-1DD-CDR-v1-3/day/pr/1x1/latest/pr_day_GPCP-1DD-CDR-v1-3_PCMDIFROGS_1x1_19961001-20201231.nc"
+reference_data_lf_path = (
+    "/work/lee1043/DATA/LandSeaMask_1x1_NCL/NCL_LandSeaMask_rewritten.nc"  # noqa
+)
+
+varOBS = "pr"
+ObsUnitsAdjust = (True, "multiply", 86400.0)  # kg m-2 s-1 to mm day-1
+
+osyear = 1998
+oeyear = 1999
+
+includeOBS = True
+
+# =================================================
+# Models
+# -------------------------------------------------
+modpath = "/work/lee1043/ESGF/xmls/cmip5/historical/day/pr/cmip5.%(model).%(exp).%(realization).day.pr.xml"
+modpath_lf = "/work/lee1043/ESGF/xmls/cmip5/historical/fx/sftlf/cmip5.%(model).historical.r0i0p0.fx.sftlf.xml"
+
+# /p/css03/scratch/published-older/cmip5/output1/CSIRO-BOM/ACCESS1-0/historical/day/atmos/day/r1i1p1/v4/pr/pr_day_ACCESS1-0_historical_r1i1p1_19750101-19991231.nc
+
+# modnames = ['ACCESS1-0', 'ACCESS1-3', 'BCC-CSM1-1', 'BCC-CSM1-1-M', 'BNU-ESM', 'CanCM4', 'CanESM2', 'CCSM4', 'CESM1-BGC', 'CESM1-CAM5', 'CESM1-FASTCHEM', 'CMCC-CESM', 'CMCC-CM', 'CMCC-CMS', 'CNRM-CM5', 'CSIRO-Mk3-6-0', 'EC-EARTH', 'FGOALS-g2', 'GFDL-CM3', 'GFDL-ESM2G', 'GFDL-ESM2M', 'GISS-E2-H', 'GISS-E2-R', 'HadGEM2-AO', 'HadGEM2-CC', 'HadGEM2-ES', 'INMCM4', 'IPSL-CM5A-LR', 'IPSL-CM5A-MR', 'IPSL-CM5B-LR', 'MIROC-ESM', 'MIROC-ESM-CHEM', 'MIROC4h', 'MIROC5', 'MPI-ESM-MR', 'MPI-ESM-P', 'MRI-CGCM3', 'MRI-ESM1', 'NorESM1-M']  # noqa
+
+modnames = ["ACCESS1-0"]
+
+realization = "r1i1p1"
+# realization = '*'
+
+varModel = "pr"
+ModUnitsAdjust = (True, "multiply", 86400.0)  # kg m-2 s-1 to mm day-1
+units = "mm/d"
+
+msyear = 1998
+meyear = 1999
+
+# =================================================
+# Output
+# -------------------------------------------------
+# pmprdir = "/p/user_pub/pmp/pmp_results/pmp_v1.1.2"
+pmprdir = "/p/user_pub/climate_work/dong12/PMP_result/"
+case_id = "{:v%Y%m%d}".format(datetime.datetime.now())
+
+if debug:
+    pmprdir = "/p/user_pub/climate_work/dong12/PMP_result/"
+    case_id = "{:v%Y%m%d-%H%M}".format(datetime.datetime.now())
+
+results_dir = os.path.join(
+    pmprdir, "%(output_type)", "monsoon", "monsoon_sperber", mip, exp, case_id
+)
+
+nc_out = True  # Write output in NetCDF
+plot = True  # Create map graphics

pcmdi_metrics/monsoon_sperber/lib/model_land_only.py

@@ -27,8 +27,8 @@
 varOBS = "pr"
 ObsUnitsAdjust = (True, "multiply", 86400.0)  # kg m-2 s-1 to mm day-1
 
-osyear = 1996
-oeyear = 2016
+osyear = 1998
+oeyear = 1999
 
 includeOBS = True
 
@@ -49,7 +49,7 @@
 ModUnitsAdjust = (True, "multiply", 86400.0)  # kg m-2 s-1 to mm day-1
 units = "mm/d"
 
-msyear = 1961
+msyear = 1998
 meyear = 1999
 
 # =================================================