From 14a61f27757cf877a369cdcc397910fc1e14b0f8 Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Mon, 14 Aug 2023 16:16:55 -0700 Subject: [PATCH 1/4] add cal --- .../variability_across_timescales_PS_driver.py | 1 + 1 file changed, 1 insertion(+) diff --git a/pcmdi_metrics/precip_variability/variability_across_timescales_PS_driver.py b/pcmdi_metrics/precip_variability/variability_across_timescales_PS_driver.py index 87894180c..83385ed21 100644 --- a/pcmdi_metrics/precip_variability/variability_across_timescales_PS_driver.py +++ b/pcmdi_metrics/precip_variability/variability_across_timescales_PS_driver.py @@ -50,6 +50,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: From 454e05eac9ffcb72e4d16a1fe1b0ef07b2a3fa9b Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Mon, 25 Sep 2023 13:41:47 -0700 Subject: [PATCH 2/4] tests --- pcmdi_metrics/io/__init__.py | 1 + .../lib/lib_variability_across_timescales.py | 6 ++++-- .../variability_across_timescales_PS_driver.py | 4 +++- share/DefArgsCIA.json | 2 +- 4 files changed, 9 insertions(+), 4 deletions(-) diff --git a/pcmdi_metrics/io/__init__.py b/pcmdi_metrics/io/__init__.py index 9822a2637..784876fbc 100644 --- a/pcmdi_metrics/io/__init__.py +++ b/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 diff --git a/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py b/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py index 1540f85e8..eba058ca0 100644 --- a/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py +++ b/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py @@ -42,10 +42,12 @@ def precip_variability_across_timescale( # Regridding rgtmp_ds = RegridHoriz(do, var, res) - if fshp is not None: - rgtmp_ds = region_from_file(rgtmp_ds,fshp,feature,attr) if regions_specs is not None or bool(regions_specs): + print("Cropping region from bounds") rgtmp_ds = CropLatLon(rgtmp_ds, regions_specs) + if fshp is not None: + 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) diff --git a/pcmdi_metrics/precip_variability/variability_across_timescales_PS_driver.py b/pcmdi_metrics/precip_variability/variability_across_timescales_PS_driver.py index 83385ed21..dcdee102a 100644 --- a/pcmdi_metrics/precip_variability/variability_across_timescales_PS_driver.py +++ b/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 diff --git a/share/DefArgsCIA.json b/share/DefArgsCIA.json index 8507f33ba..cd33a055d 100644 --- a/share/DefArgsCIA.json +++ b/share/DefArgsCIA.json @@ -163,4 +163,4 @@ ], "help":"A list of variables to be processed" } -} +} \ No newline at end of file From 9d721ac290443623556561d923ceb7d5839aeea0 Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Mon, 25 Sep 2023 13:42:42 -0700 Subject: [PATCH 3/4] tests --- pcmdi_metrics/io/region_from_file.py | 28 ++++++++++++++++++++++++++++ 1 file changed, 28 insertions(+) create mode 100644 pcmdi_metrics/io/region_from_file.py diff --git a/pcmdi_metrics/io/region_from_file.py b/pcmdi_metrics/io/region_from_file.py new file mode 100644 index 000000000..f50262bb7 --- /dev/null +++ b/pcmdi_metrics/io/region_from_file.py @@ -0,0 +1,28 @@ +import shapely +import geopandas as gpd +import regionmask +import xarray as xr +import xcdat +import fiona + +def region_from_file(data,rgn_path,attr,feature): + # Return data masked from a feature in input file. + # Arguments: + # data: xcdat dataset + # feature: str, name of region + # rgn_path: str, path to file + # attr: str, attribute name + + lon = data["lon"].data + lat = data["lat"].data + + print("Reading region from file:",rgn_path) + regions_df = gpd.read_file(rgn_path) + regions = regionmask.from_geopandas(regions_df) + mask = regions.mask(lon, lat) + # Can't match mask by name, rather index of name + val = list(regions_df[attr]).index(feature) + masked_data = data.where(mask == val) + + return masked_data + From fbed73256caa698c2d3a6f0ff27f853d9ed00857 Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Tue, 20 Feb 2024 10:22:06 -0800 Subject: [PATCH 4/4] fix calendar --- .../lib/lib_variability_across_timescales.py | 18 ++++++++++++++++-- 1 file changed, 16 insertions(+), 2 deletions(-) diff --git a/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py b/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py index eba058ca0..5294802ca 100644 --- a/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py +++ b/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 @@ -44,7 +53,10 @@ def precip_variability_across_timescale( 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: print("Cropping from shapefile") rgtmp_ds = region_from_file(rgtmp_ds,fshp,attr,feature) @@ -106,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) # ==================================================================================