Merge pull request #1065 from PCMDI/1064_precip_var_ao

Join branches for precip variability

pcmdi_metrics/io/__init__.py

@@ -4,3 +4,4 @@
 from .base import MV2Json  # noqa
 from .default_regions_define import load_regions_specs  # noqa
 from .default_regions_define import region_subset  # noqa
+from .region_from_file import region_from_file

pcmdi_metrics/io/region_from_file.py

@@ -1,8 +1,11 @@
+import shapely
+
 import geopandas as gpd
 import regionmask
 import xarray as xr
 import xcdat
 
+
 def region_from_file(data,rgn_path,attr,feature):
     # Return data masked from a feature in input file.
     # Arguments:
@@ -13,7 +16,7 @@ def region_from_file(data,rgn_path,attr,feature):
 
     lon = data["lon"].data
     lat = data["lat"].data
-
+    
     print("Reading region from file.")
     try:
         regions_df = gpd.read_file(rgn_path)

pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py

@@ -25,10 +25,19 @@ def precip_variability_across_timescale(
     """
     Regridding -> Anomaly -> Power spectra -> Domain&Frequency average -> Write
     """
-
+    from time import perf_counter
+    t1_start = perf_counter() 
     psdmfm = {"RESULTS": {}}
 
-    f = xcdat.open_mfdataset(file)
+    try:
+        f = xcdat.open_mfdataset(file)
+    except ValueError:
+        f = xcdat.open_mfdataset(file, decode_times=False)
+        cal = f.time.calendar
+        # Add any calendar fixes here
+        cal = cal.replace("-","_")
+        f.time.attrs["calendar"] = cal
+        f = xcdat.decode_time(f)
     cal = f.time.encoding["calendar"]   
     if "360" in cal:
         ldy = 30
@@ -41,9 +50,16 @@ def precip_variability_across_timescale(
         do = f.sel(time=slice(str(iyr) + "-01-01 00:00:00",str(iyr) + "-12-" + str(ldy) + " 23:59:59"))
 
         # Regridding
-        rgtmp_ds = RegridHoriz(do, var, res, regions_specs)
+        rgtmp_ds = RegridHoriz(do, var, res)
+        if regions_specs is not None or bool(regions_specs):
+            print("Cropping region from bounds")
+            t2_start = perf_counter() 
+            rgtmp_ds = CropLatLon(rgtmp_ds, regions_specs)
+            t2_stop = perf_counter()
+            print("Elapsed time cropping",t2_stop-t2_start)
         if fshp is not None:
-            rgtmp_ds = region_from_file(rgtmp_ds,fshp,feature,attr)
+            print("Cropping from shapefile")
+            rgtmp_ds = region_from_file(rgtmp_ds,fshp,attr,feature)
         rgtmp = rgtmp_ds[var]*float(fac)
         if iyr == syr:
             drg = copy.deepcopy(rgtmp)
@@ -102,6 +118,8 @@ def precip_variability_across_timescale(
     )
     if cmec:
         JSON.write_cmec(indent=4, separators=(",", ": "))
+    t1_stop = perf_counter()
+    print("Elapsed time cropping",t1_stop-t1_start)
 
 
 # ==================================================================================

pcmdi_metrics/precip_variability/variability_across_timescales_PS_driver.py

@@ -1,5 +1,7 @@
 #!/usr/bin/env python
-
+#import faulthandler; faulthandler.enable()
+import geopandas
+from pcmdi_metrics.io.region_from_file import region_from_file
 import glob
 import os
 
@@ -47,6 +49,7 @@
 
 # Check data in advance
 file_list = sorted(glob.glob(os.path.join(modpath, mod)))
+print(file_list)
 if mip == "obs":
     dat = file_list[0].split("/")[-1].split("_")[2]
 else:

share/DefArgsCIA.json

@@ -163,4 +163,4 @@
         ],
         "help":"A list of variables to be processed"
     }
-}
+}