converted to xcdat

Author: bosup

.gitignore

@@ -90,4 +90,9 @@ ENV/
 # Rope project settings
 .ropeproject
 
+# debug_data
+pcmdi_metrics/monsoon_wang/*.nc
+pcmdi_metrics/monsoon_wang/debug_regions_specs.py
+pcmdi_metrics/monsoon_wang/test_param.py
+
 test/

pcmdi_metrics/io/default_regions_define.py

@@ -1,9 +1,9 @@
-import xcdat as xc
 from typing import Union
 
 import xarray as xr
+import xcdat as xc
 
-#from pcmdi_metrics.io import da_to_ds, get_longitude, select_subset
+# from pcmdi_metrics.io import da_to_ds, get_longitude, select_subset
 from .xcdat_dataset_io import da_to_ds, get_longitude, select_subset
 
 
@@ -51,7 +51,7 @@ def load_regions_specs():
         # South American Monsoon
         "SAMM": {"domain": {"latitude": (-45.0, 0.0), "longitude": (240.0, 330.0)}},
         # North African Monsoon
-        #"NAFM": {"domain": {"latitude": (0.0, 45.0), "longitude": (310.0, 60.0)}},
+        # "NAFM": {"domain": {"latitude": (0.0, 45.0), "longitude": (310.0, 60.0)}},
         "NAFM": {"domain": {"latitude": (0.0, 45.0), "longitude": (-50, 60.0)}},
         # South African Monsoon
         "SAFM": {"domain": {"latitude": (-45.0, 0.0), "longitude": (0.0, 90.0)}},
@@ -77,7 +77,7 @@ def load_regions_specs():
     return regions_specs
 
 
-#def region_subset(ds, regions_specs, region=None):
+# def region_subset(ds, regions_specs, region=None):
 #    """
 #    d: xarray.Dataset
 #    regions_specs: dict
@@ -132,6 +132,7 @@ def load_regions_specs():
 #
 #    return ds
 
+
 def region_subset(
     ds: Union[xr.Dataset, xr.DataArray],
     region: str,

pcmdi_metrics/io/xcdat_dataset_io.py

@@ -3,7 +3,6 @@
 import xarray as xr
 import xcdat as xc
 
-
 # Internal function
 
 

pcmdi_metrics/monsoon_sperber/driver_monsoon_sperber.py

@@ -46,17 +46,18 @@
 from shutil import copyfile
 
 import matplotlib
-matplotlib.use('Agg')
-#import matplotlib.pyplot as plt
-from matplotlib import pyplot as plt
+
+matplotlib.use("Agg")
 import numpy as np
 import pandas as pd
 import xarray as xr
 import xcdat as xc
+# import matplotlib.pyplot as plt
+from matplotlib import pyplot as plt
 
 import pcmdi_metrics
 from pcmdi_metrics import resources
-from pcmdi_metrics.io import load_regions_specs, region_subset
+from pcmdi_metrics.io import load_regions_specs, region_subset, xcdat_open
 from pcmdi_metrics.mean_climate.lib import pmp_parser
 from pcmdi_metrics.monsoon_sperber.lib import (
     AddParserArgument,
@@ -67,7 +68,6 @@
     sperber_metrics,
 )
 from pcmdi_metrics.utils import create_land_sea_mask, fill_template
-from pcmdi_metrics.io import xcdat_open
 
 
 def tree():
@@ -92,7 +92,7 @@ def pick_year_last_day(ds):
 # =================================================
 # Hard coded options... will be moved out later
 # -------------------------------------------------
-#list_monsoon_regions = ["AIR", "AUS", "Sahel", "GoG", "NAmo", "SAmo"]
+# list_monsoon_regions = ["AIR", "AUS", "Sahel", "GoG", "NAmo", "SAmo"]
 list_monsoon_regions = ["AIR", "Sahel"]
 # list_monsoon_regions = ["all"]
 
@@ -293,7 +293,16 @@ def pick_year_last_day(ds):
                 modpath(model=model, exp=exp, realization=realization, variable=var)
             )
             if debug:
-                print("model: ", model, "   exp: ", exp, "  realization: ", realization, "  variable: ", var)
+                print(
+                    "model: ",
+                    model,
+                    "   exp: ",
+                    exp,
+                    "  realization: ",
+                    realization,
+                    "  variable: ",
+                    var,
+                )
                 print("debug: model_path_list: ", model_path_list)
             # land fraction
             model_lf_path = modpath_lf(model=model)
@@ -314,18 +323,18 @@ def pick_year_last_day(ds):
                     ds_lf = xcdat_open(model_lf_path)
                 except Exception:
                     ds_lf = None
-        
-        if not ds_lf: 
+
+        if not ds_lf:
             lf_array = create_land_sea_mask(ds_lf, method="pcmdi")
             ds_lf = lf_array.to_dataset().compute()
             ds_lf = ds_lf.rename_vars({"lsmask": "sftlf"})
 
-	#  use pcmdi mask
-	#  lf_array = create_land_sea_mask(ds_lf, method="pcmdi")
-	#  ds_lf = lf_array.to_dataset().compute()
-	#  ds_lf = ds_lf.rename_vars({"lsmask": "sftlf"})
+        #  use pcmdi mask
+        #  lf_array = create_land_sea_mask(ds_lf, method="pcmdi")
+        #  ds_lf = lf_array.to_dataset().compute()
+        #  ds_lf = ds_lf.rename_vars({"lsmask": "sftlf"})
 
-        if model in [ "EC-EARTH" ]: #, "BNU-ESM" ]:
+        if model in ["EC-EARTH"]:  # , "BNU-ESM" ]:
             ds_lf = ds_lf.isel(lat=slice(None, None, -1))
         lf = ds_lf.sftlf.sel(lat=slice(-90, 90))  # land frac file must be global
 
@@ -350,10 +359,9 @@ def pick_year_last_day(ds):
                 # Get time coordinate information
                 print("model_path =   ", model_path)
 
-
                 dc = xcdat_open(model_path, decode_times=True)
-                dc['time'].attrs['axis'] = 'T'
-                dc['time'].attrs['standard_name'] = 'time'
+                dc["time"].attrs["axis"] = "T"
+                dc["time"].attrs["standard_name"] = "time"
                 dc = xr.decode_cf(dc, decode_times=True)
                 dc = dc.bounds.add_missing_bounds("X")
                 dc = dc.bounds.add_missing_bounds("Y")
@@ -363,7 +371,6 @@ def pick_year_last_day(ds):
                 c = xc.center_times(dc)
                 eday = pick_year_last_day(dc)
 
-
                 # Get starting and ending year and month
                 startYear = c.time.values[0].year
                 startMonth = c.time.values[0].month
@@ -493,7 +500,6 @@ def pick_year_last_day(ds):
                         d.values = d.values * 86400.0
                         d["units"] = units
 
-
                     # variable for over land only
                     d_land = model_land_only(model, d, lf, debug=debug)
 
@@ -523,7 +529,6 @@ def pick_year_last_day(ds):
                             d_sub_pr.values = d_sub_pr.values * 86400.0
                             d_sub_pr["units"] = units
 
-
                         else:
                             # land-only rainfall
 
@@ -544,11 +549,9 @@ def pick_year_last_day(ds):
                                 model, d_sub_pr, lf_sub, debug=debug
                             )
 
-
                             d_sub_pr.values = d_sub_pr.values * 86400.0
                             d_sub_pr["units"] = units
 
-
                         # Area average
 
                         ds_sub_pr = d_sub_pr.to_dataset().compute()
@@ -557,7 +560,6 @@ def pick_year_last_day(ds):
                         ds_sub_pr = ds_sub_pr.bounds.add_missing_bounds("Y")
                         ds_sub_pr = ds_sub_pr.bounds.add_missing_bounds("T")
 
-
                         if "lat_bnds" not in ds_sub_pr.variables:
                             lat_bnds = dc["lat_bnds"].sel(lat=ds_sub_pr["lat"])
                             ds_sub_pr["lat_bnds"] = lat_bnds
@@ -567,8 +569,6 @@ def pick_year_last_day(ds):
                         ).compute()
                         d_sub_aave = ds_sub_aave.pr
 
-
-
                         if debug:
                             print("debug: region:", region)
 
@@ -603,7 +603,6 @@ def pick_year_last_day(ds):
 
                                 d_sub_aave = xr.concat([part1, part2], dim="time")
 
-
                                 if debug:
                                     print(
                                         "debug: ",
@@ -638,7 +637,6 @@ def pick_year_last_day(ds):
                             coords={"time": time_coords},
                         )
 
-
                         if debug:
                             print(
                                 "debug: pentad_time_series length: ",
@@ -648,14 +646,14 @@ def pick_year_last_day(ds):
                         # Keep pentad time series length in consistent
                         ref_length = int(365 / n)
                         if len(pentad_time_series) < ref_length:
-
                             pentad_time_series = pentad_time_series.interp(
-                                time=pd.date_range(time_coords[0], time_coords[-1], periods=ref_length)
+                                time=pd.date_range(
+                                    time_coords[0], time_coords[-1], periods=ref_length
+                                )
                             )
 
                             time_coords = pentad_time_series.coords["time"]
 
-
                         pentad_time_series_cumsum = np.cumsum(pentad_time_series)
                         pentad_time_series = xr.DataArray(
                             pentad_time_series,
@@ -672,7 +670,6 @@ def pick_year_last_day(ds):
                         pentad_time_series_cumsum.attrs["units"] = str(d.units.values)
                         pentad_time_series_cumsum.coords["time"] = time_coords
 
-
                         if nc_out:
                             # Archive individual year time series in netCDF file
                             pentad_time_series.to_netcdf(file_path, mode="a")
@@ -819,9 +816,7 @@ def pick_year_last_day(ds):
 
                         # obs
                         ax[region].plot(
-                            np.array(
-                                dict_obs_composite[reference_data_name][region]
-                            ),
+                            np.array(dict_obs_composite[reference_data_name][region]),
                             c="blue",
                             label=reference_data_name,
                         )
@@ -836,9 +831,9 @@ def pick_year_last_day(ds):
                             ax[region].axvline(
                                 x=idx,
                                 ymin=0,
-                                ymax=dict_obs_composite[reference_data_name][
-                                    region
-                                ][idx].item(),
+                                ymax=dict_obs_composite[reference_data_name][region][
+                                    idx
+                                ].item(),
                                 c="blue",
                                 ls="--",
                             )

pcmdi_metrics/monsoon_sperber/lib/argparse_functions.py

@@ -54,7 +54,9 @@ def AddParserArgument(P):
         "--meyear", dest="meyear", type=int, help="End year for model data set"
     )
     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(
+        "--list_monsoon_regions", type=str, default=None, help="List of regions"
+    )
     P.add_argument(
         "-r",
         "--realization",

pcmdi_metrics/monsoon_sperber/lib/divide_chunks.py

@@ -30,7 +30,7 @@ def divide_chunks_advanced(data, n, debug=False):
     day = day.values
     calendar = "gregorian"
     if debug:
-        #print("month = ", month, "day = ", day)
+        # print("month = ", month, "day = ", day)
         print("debug: first day of year is " + str(month) + "/" + str(day))
     if month not in [1, 7] or day != 1:
         sys.exit(

pcmdi_metrics/monsoon_sperber/lib/model_land_only.py

@@ -9,7 +9,7 @@ def model_land_only(model, model_timeseries, lf, debug=False):
     # - - - - - - - - - - - - - - - - - - - - - - - - -
 
     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
@@ -33,7 +33,7 @@ def model_land_only(model, model_timeseries, lf, debug=False):
             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

pcmdi_metrics/monsoon_sperber/param/Bo_param.py

@@ -14,9 +14,9 @@
 # -------------------------------------------------
 update_json = False
 debug = False
-#debug = True
+# debug = True
 
-#list_monsoon_regions = ["AIR", "AUS", "Sahel", "GoG", "NAmo", "SAmo"]
+# list_monsoon_regions = ["AIR", "AUS", "Sahel", "GoG", "NAmo", "SAmo"]
 list_monsoon_regions = ["AUS"]
 # =================================================
 # Observation
@@ -41,7 +41,7 @@
 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
+# /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
 
@@ -60,7 +60,7 @@
 # =================================================
 # Output
 # -------------------------------------------------
-#pmprdir = "/p/user_pub/pmp/pmp_results/pmp_v1.1.2"
+# 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())