From d8605e1ea6a26812399074a8c7abef4117654c86 Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Wed, 19 Jul 2023 14:06:18 -0700 Subject: [PATCH 01/55] Update calendar --- .../lib/lib_variability_across_timescales.py | 4 ++-- 1 file changed, 2 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 5f54506c7..2eaaf0673 100644 --- a/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py +++ b/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py @@ -130,7 +130,7 @@ def ClimAnom(d, ntd, syr, eyr, cal): # Year segment nyr = eyr - syr + 1 - if "gregorian" in cal: + if "gregorian" in cal or "standard" in cal: ndy = 366 ldy = 31 dseg = np.zeros((nyr, ndy, ntd, d.shape[1], d.shape[2]), dtype=float) @@ -171,7 +171,7 @@ def ClimAnom(d, ntd, syr, eyr, cal): # Anomaly anom = np.array([]) - if "gregorian" in cal: + if "gregorian" in cal or "standard" in cal: for iyr, year in enumerate(range(syr, eyr + 1)): yrtmp = d.sel(time=slice(str(year) + "-01-01 00:00:00",str(year) + "-12-" + str(ldy) + " 23:59:59")) From 5c0f295be001d0b11fe8cafb341b0fbf98b7c260 Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Tue, 22 Aug 2023 13:38:13 -0700 Subject: [PATCH 02/55] Add resolution --- .../lib/lib_variability_across_timescales.py | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 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 5f54506c7..c1fcd5e0c 100644 --- a/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py +++ b/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py @@ -17,7 +17,7 @@ # ================================================================================== def precip_variability_across_timescale( - file, syr, eyr, dfrq, mip, dat, var, fac, nperseg, noverlap, outdir, cmec + file, syr, eyr, dfrq, mip, dat, var, fac, nperseg, noverlap, res, outdir, cmec ): """ Regridding -> Anomaly -> Power spectra -> Domain&Frequency average -> Write @@ -38,7 +38,7 @@ 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 = Regrid2deg(do, var)*float(fac) + rgtmp = RegridHoriz(do, var, res)*float(fac) if iyr == syr: drg = copy.deepcopy(rgtmp) else: @@ -97,7 +97,7 @@ def precip_variability_across_timescale( # ================================================================================== -def Regrid2deg(d, var): +def RegridHoriz(d, var, res): """ Regrid to 2deg (180lon*90lat) horizontal resolution Input @@ -106,7 +106,14 @@ def Regrid2deg(d, var): Output - drg: xCDAT variable with 2deg horizontal resolution """ - tgrid = grid.create_uniform_grid(-89, 89, 2.0, 0.0, 358., 2.0) + start_lat=-90.+res/2. + start_lon=0. + end_lat = 90.-res/2. + end_lon = 360.-res + nlat = ((end_lat - start_lat) * 1./res) + 1 + nlon = ((end_lon - start_lon) * 1./res) + 1 + + tgrid = grid.create_uniform_grid(start_lat,end_lat,res,start_lon,end_lon,res) drg = d.regridder.horizontal(var, tgrid, tool="xesmf", method="conservative_normed",periodic=True)[var] print("Complete regridding from", d[var].shape, "to", drg.shape) From 1bbeef4165d8794ecd7e14b2656ed2a761023a19 Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Tue, 22 Aug 2023 13:39:24 -0700 Subject: [PATCH 03/55] Add resolution --- .../variability_across_timescales_PS_driver.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) 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 c606be499..5e981a2fb 100644 --- a/pcmdi_metrics/precip_variability/variability_across_timescales_PS_driver.py +++ b/pcmdi_metrics/precip_variability/variability_across_timescales_PS_driver.py @@ -24,10 +24,12 @@ fac = param.fac nperseg = param.nperseg noverlap = param.noverlap +res = param.res print(modpath) print(mod) print(prd) print(nperseg, noverlap) +print(res) # Get flag for CMEC output cmec = param.cmec @@ -57,5 +59,5 @@ syr = prd[0] eyr = prd[1] precip_variability_across_timescale( - file_list, syr, eyr, dfrq, mip, dat, var, fac, nperseg, noverlap, outdir, cmec + file_list, syr, eyr, dfrq, mip, dat, var, fac, nperseg, noverlap, res, outdir, cmec ) From e4e59c473b3f0a4d75136e4edfec749e8f9ff498 Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Tue, 22 Aug 2023 13:42:01 -0700 Subject: [PATCH 04/55] Add res to arguments --- pcmdi_metrics/precip_variability/lib/argparse_functions.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/pcmdi_metrics/precip_variability/lib/argparse_functions.py b/pcmdi_metrics/precip_variability/lib/argparse_functions.py index 98ce33225..8643dc105 100644 --- a/pcmdi_metrics/precip_variability/lib/argparse_functions.py +++ b/pcmdi_metrics/precip_variability/lib/argparse_functions.py @@ -59,6 +59,13 @@ def AddParserArgument(P): P.add_argument( "--ref", type=str, dest="ref", default=None, help="reference data path" ) + P.add_argument( + "--res", + type=int, + dest="res", + default=2, + help="Horizontal resolution [degree] for interpolation (lon, lat)", + ) P.add_argument( "--cmec", dest="cmec", From 434a202bce1e0ae5dd8acc9945d8bb54abd0a16e Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Tue, 22 Aug 2023 13:48:08 -0700 Subject: [PATCH 05/55] update file names --- .../lib/lib_variability_across_timescales.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 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 c1fcd5e0c..7d0b68da2 100644 --- a/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py +++ b/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py @@ -44,7 +44,9 @@ def precip_variability_across_timescale( else: drg = xr.concat([drg, rgtmp], dim='time') print(iyr, drg.shape) - + nlon = str(len(drg.lon)) + nlat = str(len(drg.lat)) + f.close() # Anomaly @@ -61,7 +63,7 @@ def precip_variability_across_timescale( # Domain & Frequency average psdmfm_forced = Avg_PS_DomFrq(ps, freqs, ntd, dat, mip, "forced") # Write data (nc file) - outfilename = "PS_pr." + str(dfrq) + "_regrid.180x90_" + dat + ".nc" + outfilename = "PS_pr." + str(dfrq) + "_regrid." + nlon + "x" + nlat + "_" + dat + ".nc" custom_dataset = xr.merge([freqs, ps, rn, sig95]) custom_dataset.to_netcdf(path=os.path.join(outdir(output_type="diagnostic_results"), outfilename)) @@ -70,7 +72,7 @@ def precip_variability_across_timescale( # Domain & Frequency average psdmfm_unforced = Avg_PS_DomFrq(ps, freqs, ntd, dat, mip, "unforced") # Write data (nc file) - outfilename = "PS_pr." + str(dfrq) + "_regrid.180x90_" + dat + "_unforced.nc" + outfilename = "PS_pr." + str(dfrq) + "_regrid." + nlon + "x" + nlat + "_" + dat + "_unforced.nc" custom_dataset = xr.merge([freqs, ps, rn, sig95]) custom_dataset.to_netcdf(path=os.path.join(outdir(output_type="diagnostic_results"), outfilename)) @@ -80,7 +82,7 @@ def precip_variability_across_timescale( psdmfm["RESULTS"][dat]["unforced"] = psdmfm_unforced outfilename = ( - "PS_pr." + str(dfrq) + "_regrid.180x90_area.freq.mean_" + dat + ".json" + "PS_pr." + str(dfrq) + "_regrid." + nlon + "x" + nlat + "_area.freq.mean_" + dat + ".json" ) JSON = pcmdi_metrics.io.base.Base( outdir(output_type="metrics_results"), outfilename From aedfbdc8a0b9e76ca1d5116dc16bb4f02afea4f4 Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Tue, 22 Aug 2023 14:00:45 -0700 Subject: [PATCH 06/55] add bounding box --- .../lib/lib_variability_across_timescales.py | 29 ++++++++++++++++++- 1 file changed, 28 insertions(+), 1 deletion(-) 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 7d0b68da2..8376ad08f 100644 --- a/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py +++ b/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py @@ -17,7 +17,7 @@ # ================================================================================== def precip_variability_across_timescale( - file, syr, eyr, dfrq, mip, dat, var, fac, nperseg, noverlap, res, outdir, cmec + file, syr, eyr, dfrq, mip, dat, var, fac, nperseg, noverlap, res, lat_range, lon_range, outdir, cmec ): """ Regridding -> Anomaly -> Power spectra -> Domain&Frequency average -> Write @@ -39,6 +39,8 @@ def precip_variability_across_timescale( # Regridding rgtmp = RegridHoriz(do, var, res)*float(fac) + if len(lat_range) > 0 and len(lon_range) > 0: + rgtmp = CropLatLon(rgtmp, lat_range, lon_range) if iyr == syr: drg = copy.deepcopy(rgtmp) else: @@ -121,6 +123,31 @@ def RegridHoriz(d, var, res): print("Complete regridding from", d[var].shape, "to", drg.shape) return drg +# ================================================================================== +def CropLatLon(d,lat_range,lon_range): + """ + Select a subgrid of the dataset defined by lat1, lat2, lon1, and lon2. + Input + - d: xCDAT variable + - lat_range: list of floats + - lon_range: list of floats + Output + - dnew: xCDAT variable selected over region of interest + """ + lat1 = lat_range[0] + lat2 = lat_range[1] + lon1 = lon_range[0] + lon2 = lon_range[1] + + try: + dnew = d.sel(lat=slice(lat1,lat2), lon=slice(lon1,lon2)) + + except Exception as e: + print("Error:",e) + print("Could not select lat/lon box",lat1,lat2,lon1,lon2) + dnew = d + + return dnew # ================================================================================== def ClimAnom(d, ntd, syr, eyr, cal): From 4af3b02b42ae712f527f9628614be93e74ab8313 Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Tue, 22 Aug 2023 14:17:06 -0700 Subject: [PATCH 07/55] Add regions specs --- .../variability_across_timescales_PS_driver.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) 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 5e981a2fb..809bf3eab 100644 --- a/pcmdi_metrics/precip_variability/variability_across_timescales_PS_driver.py +++ b/pcmdi_metrics/precip_variability/variability_across_timescales_PS_driver.py @@ -25,6 +25,7 @@ nperseg = param.nperseg noverlap = param.noverlap res = param.res +regions_specs = param.regions_specs print(modpath) print(mod) print(prd) @@ -59,5 +60,5 @@ syr = prd[0] eyr = prd[1] precip_variability_across_timescale( - file_list, syr, eyr, dfrq, mip, dat, var, fac, nperseg, noverlap, res, outdir, cmec + file_list, syr, eyr, dfrq, mip, dat, var, fac, nperseg, noverlap, res, regions_specs, outdir, cmec ) From c20e975937afa6ccb12778dfb8cadca5ba257348 Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Tue, 22 Aug 2023 14:25:40 -0700 Subject: [PATCH 08/55] Use regions specs for region --- .../lib/lib_variability_across_timescales.py | 26 ++++++------------- 1 file changed, 8 insertions(+), 18 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 8376ad08f..2262df7f5 100644 --- a/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py +++ b/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py @@ -17,7 +17,7 @@ # ================================================================================== def precip_variability_across_timescale( - file, syr, eyr, dfrq, mip, dat, var, fac, nperseg, noverlap, res, lat_range, lon_range, outdir, cmec + file, syr, eyr, dfrq, mip, dat, var, fac, nperseg, noverlap, res, regions_specs, outdir, cmec ): """ Regridding -> Anomaly -> Power spectra -> Domain&Frequency average -> Write @@ -39,8 +39,8 @@ def precip_variability_across_timescale( # Regridding rgtmp = RegridHoriz(do, var, res)*float(fac) - if len(lat_range) > 0 and len(lon_range) > 0: - rgtmp = CropLatLon(rgtmp, lat_range, lon_range) + if regions_specs is not None or bool(regions_specs): + rgtmp = CropLatLon(rgtmp, regions_specs) if iyr == syr: drg = copy.deepcopy(rgtmp) else: @@ -124,28 +124,18 @@ def RegridHoriz(d, var, res): return drg # ================================================================================== -def CropLatLon(d,lat_range,lon_range): +def CropLatLon(d, regions_specs): """ - Select a subgrid of the dataset defined by lat1, lat2, lon1, and lon2. + Select a subgrid of the dataset defined by the regions_specs dictionary. Input - d: xCDAT variable - - lat_range: list of floats - - lon_range: list of floats + - regions_specs: a dictionary Output - dnew: xCDAT variable selected over region of interest """ - lat1 = lat_range[0] - lat2 = lat_range[1] - lon1 = lon_range[0] - lon2 = lon_range[1] + region_name = list(regions_specs.keys())[0] - try: - dnew = d.sel(lat=slice(lat1,lat2), lon=slice(lon1,lon2)) - - except Exception as e: - print("Error:",e) - print("Could not select lat/lon box",lat1,lat2,lon1,lon2) - dnew = d + dnew = region_subset(d, regions_specs, region=region_name) return dnew From 93ecfaf2bff3a00dc99bc58bb4e30a671aafb966 Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Tue, 22 Aug 2023 14:26:10 -0700 Subject: [PATCH 09/55] add import --- .../precip_variability/lib/lib_variability_across_timescales.py | 1 + 1 file changed, 1 insertion(+) 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 2262df7f5..8f8f535d7 100644 --- a/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py +++ b/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py @@ -10,6 +10,7 @@ from scipy.stats import chi2 import pcmdi_metrics +from pcmdi_metrics.io import region_subset import xarray as xr import xcdat From b66ff65b8c871a316cd7b97ee24235d5fd4e98d2 Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Tue, 22 Aug 2023 14:39:12 -0700 Subject: [PATCH 10/55] add regions_specs --- .../precip_variability/lib/argparse_functions.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/pcmdi_metrics/precip_variability/lib/argparse_functions.py b/pcmdi_metrics/precip_variability/lib/argparse_functions.py index 8643dc105..285ba1127 100644 --- a/pcmdi_metrics/precip_variability/lib/argparse_functions.py +++ b/pcmdi_metrics/precip_variability/lib/argparse_functions.py @@ -66,6 +66,14 @@ def AddParserArgument(P): default=2, help="Horizontal resolution [degree] for interpolation (lon, lat)", ) + P.add_argument( + "--regions_specs", + type=ast.literal_eval, + dest="regions_specs", + help="Provide a single custom region", + default=None, + required=False, + ) P.add_argument( "--cmec", dest="cmec", From 90177edd8831a379916c7335056a2bc413c0793c Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Tue, 22 Aug 2023 15:26:32 -0700 Subject: [PATCH 11/55] Add new functions --- pcmdi_metrics/precip_variability/lib/__init__.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pcmdi_metrics/precip_variability/lib/__init__.py b/pcmdi_metrics/precip_variability/lib/__init__.py index b1979aa4d..2ce38a093 100644 --- a/pcmdi_metrics/precip_variability/lib/__init__.py +++ b/pcmdi_metrics/precip_variability/lib/__init__.py @@ -4,7 +4,8 @@ ClimAnom, Powerspectrum, RedNoiseSignificanceLevel, - Regrid2deg, + RegridHoriz, + CropLatLon, lag1_autocorrelation, prdday_to_frqidx, precip_variability_across_timescale, From 30f08c5828dadea13b4be370734ab221efde529b Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Tue, 22 Aug 2023 15:27:32 -0700 Subject: [PATCH 12/55] add import for regions_specs --- pcmdi_metrics/precip_variability/lib/argparse_functions.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/pcmdi_metrics/precip_variability/lib/argparse_functions.py b/pcmdi_metrics/precip_variability/lib/argparse_functions.py index 285ba1127..9ab94b017 100644 --- a/pcmdi_metrics/precip_variability/lib/argparse_functions.py +++ b/pcmdi_metrics/precip_variability/lib/argparse_functions.py @@ -1,3 +1,5 @@ +import ast + def AddParserArgument(P): P.add_argument( "--mip", type=str, dest="mip", default=None, help="cmip5, cmip6 or other mip" From da35fe9a810fddaae78a6a44fba512bdf50605c4 Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Wed, 23 Aug 2023 13:42:51 -0700 Subject: [PATCH 13/55] address xesmf zeros issue --- .../lib/lib_variability_across_timescales.py | 9 +++++++++ 1 file changed, 9 insertions(+) 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 8f8f535d7..f197c923d 100644 --- a/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py +++ b/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py @@ -120,6 +120,15 @@ def RegridHoriz(d, var, res): tgrid = grid.create_uniform_grid(start_lat,end_lat,res,start_lon,end_lon,res) drg = d.regridder.horizontal(var, tgrid, tool="xesmf", method="conservative_normed",periodic=True)[var] + + # Workaround for regional grid cases where areas outside + # original grid might be populated with 0's due to xesmf. + drg = drg.where( + (drg.lat >= d.lat.min()) & + (drg.lat <= d.lat.max()) & + (drg.lon >= d.lon.min()) & + (drg.lon <= d.lon.max()) + ) print("Complete regridding from", d[var].shape, "to", drg.shape) return drg From 74cd62773fd787fae3b01518d8c2edeb7d76416d Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Wed, 23 Aug 2023 16:10:44 -0700 Subject: [PATCH 14/55] handle change in lon axis --- .../precip_variability/lib/lib_variability_across_timescales.py | 2 ++ 1 file changed, 2 insertions(+) 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 1f310f440..c60de7014 100644 --- a/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py +++ b/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py @@ -123,6 +123,8 @@ def RegridHoriz(d, var, res): # Workaround for regional grid cases where areas outside # original grid might be populated with 0's due to xesmf. + if d.lon.min() < 0: + d = xc.swap_lon_axis(d, to=(0, 360)) drg = drg.where( (drg.lat >= d.lat.min()) & (drg.lat <= d.lat.max()) & From cc58283a3f6a016bfa08958af5311941dd349bf1 Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Wed, 23 Aug 2023 16:16:59 -0700 Subject: [PATCH 15/55] update call --- .../precip_variability/lib/lib_variability_across_timescales.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 c60de7014..1707a7603 100644 --- a/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py +++ b/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py @@ -124,7 +124,7 @@ def RegridHoriz(d, var, res): # Workaround for regional grid cases where areas outside # original grid might be populated with 0's due to xesmf. if d.lon.min() < 0: - d = xc.swap_lon_axis(d, to=(0, 360)) + d = xcdat.swap_lon_axis(d, to=(0, 360)) drg = drg.where( (drg.lat >= d.lat.min()) & (drg.lat <= d.lat.max()) & From 38ee1d87869125c6a0cfa9d72b0054197db13b27 Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Mon, 11 Sep 2023 11:38:36 -0700 Subject: [PATCH 16/55] Use unmapped_to_nan in regridding --- .../lib/lib_variability_across_timescales.py | 15 +-------------- 1 file changed, 1 insertion(+), 14 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 1707a7603..d64c64251 100644 --- a/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py +++ b/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py @@ -115,22 +115,9 @@ def RegridHoriz(d, var, res): start_lon=0. end_lat = 90.-res/2. end_lon = 360.-res - nlat = ((end_lat - start_lat) * 1./res) + 1 - nlon = ((end_lon - start_lon) * 1./res) + 1 tgrid = grid.create_uniform_grid(start_lat,end_lat,res,start_lon,end_lon,res) - drg = d.regridder.horizontal(var, tgrid, tool="xesmf", method="conservative_normed",periodic=True)[var] - - # Workaround for regional grid cases where areas outside - # original grid might be populated with 0's due to xesmf. - if d.lon.min() < 0: - d = xcdat.swap_lon_axis(d, to=(0, 360)) - drg = drg.where( - (drg.lat >= d.lat.min()) & - (drg.lat <= d.lat.max()) & - (drg.lon >= d.lon.min()) & - (drg.lon <= d.lon.max()) - ) + drg = d.regridder.horizontal(var, tgrid, tool="xesmf", method="conservative_normed",periodic=True,unmapped_to_nan=True)[var] print("Complete regridding from", d[var].shape, "to", drg.shape) return drg From c79d7ee82953febf9d21e76e32f3ff5383b85f1f Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Mon, 11 Sep 2023 11:50:10 -0700 Subject: [PATCH 17/55] Keep dataset for cropping --- .../lib/lib_variability_across_timescales.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 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 d64c64251..9a3d41fa7 100644 --- a/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py +++ b/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py @@ -39,9 +39,10 @@ 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 = RegridHoriz(do, var, res)*float(fac) + rgtmp_ds = RegridHoriz(do, res) if regions_specs is not None or bool(regions_specs): rgtmp = CropLatLon(rgtmp, regions_specs) + rgtmp = rgtmp_ds[var]*float(fac) if iyr == syr: drg = copy.deepcopy(rgtmp) else: @@ -102,7 +103,7 @@ def precip_variability_across_timescale( # ================================================================================== -def RegridHoriz(d, var, res): +def RegridHoriz(d, res): """ Regrid to 2deg (180lon*90lat) horizontal resolution Input @@ -117,7 +118,7 @@ def RegridHoriz(d, var, res): end_lon = 360.-res tgrid = grid.create_uniform_grid(start_lat,end_lat,res,start_lon,end_lon,res) - drg = d.regridder.horizontal(var, tgrid, tool="xesmf", method="conservative_normed",periodic=True,unmapped_to_nan=True)[var] + drg = d.regridder.horizontal(var, tgrid, tool="xesmf", method="conservative_normed",periodic=True,unmapped_to_nan=True) print("Complete regridding from", d[var].shape, "to", drg.shape) return drg From 9bd3f0c64201e3aca169f87cb3553f8faa5981aa Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Mon, 11 Sep 2023 11:52:17 -0700 Subject: [PATCH 18/55] Regridder needs var --- .../lib/lib_variability_across_timescales.py | 4 ++-- 1 file changed, 2 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 9a3d41fa7..613800724 100644 --- a/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py +++ b/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py @@ -39,7 +39,7 @@ 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, res) + rgtmp_ds = RegridHoriz(do, var, res) if regions_specs is not None or bool(regions_specs): rgtmp = CropLatLon(rgtmp, regions_specs) rgtmp = rgtmp_ds[var]*float(fac) @@ -103,7 +103,7 @@ def precip_variability_across_timescale( # ================================================================================== -def RegridHoriz(d, res): +def RegridHoriz(d, var, res): """ Regrid to 2deg (180lon*90lat) horizontal resolution Input From 686e515fdbd4253f27cdc22a3aef23b77cb19e42 Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Mon, 11 Sep 2023 11:54:35 -0700 Subject: [PATCH 19/55] add var call --- .../precip_variability/lib/lib_variability_across_timescales.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 613800724..634503638 100644 --- a/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py +++ b/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py @@ -120,7 +120,7 @@ def RegridHoriz(d, var, res): tgrid = grid.create_uniform_grid(start_lat,end_lat,res,start_lon,end_lon,res) drg = d.regridder.horizontal(var, tgrid, tool="xesmf", method="conservative_normed",periodic=True,unmapped_to_nan=True) - print("Complete regridding from", d[var].shape, "to", drg.shape) + print("Complete regridding from", d[var].shape, "to", drg[var].shape) return drg # ================================================================================== From d007ffa5d45adaa5ae0278b70ac139ad91583d9d Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Mon, 11 Sep 2023 11:56:17 -0700 Subject: [PATCH 20/55] Use rgtmp_ds variable --- .../precip_variability/lib/lib_variability_across_timescales.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 634503638..5b5abb50f 100644 --- a/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py +++ b/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py @@ -41,7 +41,7 @@ def precip_variability_across_timescale( # Regridding rgtmp_ds = RegridHoriz(do, var, res) if regions_specs is not None or bool(regions_specs): - rgtmp = CropLatLon(rgtmp, regions_specs) + rgtmp_ds = CropLatLon(rgtmp_ds, regions_specs) rgtmp = rgtmp_ds[var]*float(fac) if iyr == syr: drg = copy.deepcopy(rgtmp) From 453869435ebb66bbb8a690deaaaf7c514ba3dea4 Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Mon, 11 Sep 2023 13:51:15 -0700 Subject: [PATCH 21/55] Regional metrics in JSON --- .../lib/lib_variability_across_timescales.py | 20 +++++++++++++------ 1 file changed, 14 insertions(+), 6 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 5b5abb50f..81dbe5487 100644 --- a/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py +++ b/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py @@ -331,7 +331,7 @@ def RedNoiseSignificanceLevel(nu, rn, p=0.050): # ================================================================================== -def Avg_PS_DomFrq(d, frequency, ntd, dat, mip, frc): +def Avg_PS_DomFrq(d, frequency, ntd, dat, mip, frc, regions_specs): """ Domain and Frequency average of spectral power Input @@ -341,6 +341,7 @@ def Avg_PS_DomFrq(d, frequency, ntd, dat, mip, frc): - dat: model name - mip: mip name - frc: forced or unforced + - regions_specs: dictionary defining custom region Output - psdmfm: Domain and Frequency averaged of spectral power (json) """ @@ -358,6 +359,8 @@ def Avg_PS_DomFrq(d, frequency, ntd, dat, mip, frc): "Ocean_50S30S", "Land_50S30S", ] + if regions_specs: + domains = list(regions_specs.keys()) if ntd == 8: # 3-hourly data frqs_forced = ["semi-diurnal", "diurnal", "semi-annual", "annual"] @@ -389,9 +392,9 @@ def Avg_PS_DomFrq(d, frequency, ntd, dat, mip, frc): for dom in domains: psdmfm[dom] = {} - if "Ocean" in dom: + if "Ocean" in dom or regions_specs[dom].get("value",-1)==0: dmask = d.where(mask==0) - elif "Land" in dom: + elif "Land" in dom or regions_specs[dom].get("value",-1)==1: dmask = d.where(mask==1) else: dmask = d @@ -402,12 +405,17 @@ def Avg_PS_DomFrq(d, frequency, ntd, dat, mip, frc): if "50S50N" in dom: am = dmask.sel(lat=slice(-50, 50)).spatial.average("ps", axis=["X", "Y"], weights='generate')["ps"] - if "30N50N" in dom: + elif "30N50N" in dom: am = dmask.sel(lat=slice(30, 50)).spatial.average("ps", axis=["X", "Y"], weights='generate')["ps"] - if "30S30N" in dom: + elif "30S30N" in dom: am = dmask.sel(lat=slice(-30, 30)).spatial.average("ps", axis=["X", "Y"], weights='generate')["ps"] - if "50S30S" in dom: + elif "50S30S" in dom: am = dmask.sel(lat=slice(-50, -30)).spatial.average("ps", axis=["X", "Y"], weights='generate')["ps"] + else: + # Custom region + lats=[regions_specs[dom]["domain"]["latitude"]] + lons=[regions_specs[dom]["domain"]["longitude"]] + am = dmask.sel(lat=slice(lats[0],lats[1]),lon=slice(lons[0],lons[1])).spatial.average("ps", axis=["X", "Y"], weights='generate')["ps"] am = np.array(am) for frq in frqs: From de031f3a0ef2311d4c4d4ff51b643c858867d961 Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Mon, 11 Sep 2023 14:20:40 -0700 Subject: [PATCH 22/55] Add regions_specs to function call --- .../lib/lib_variability_across_timescales.py | 4 ++-- 1 file changed, 2 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 81dbe5487..e3e039eb6 100644 --- a/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py +++ b/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py @@ -65,7 +65,7 @@ def precip_variability_across_timescale( # Power spectum of total freqs, ps, rn, sig95 = Powerspectrum(drg, nperseg, noverlap) # Domain & Frequency average - psdmfm_forced = Avg_PS_DomFrq(ps, freqs, ntd, dat, mip, "forced") + psdmfm_forced = Avg_PS_DomFrq(ps, freqs, ntd, dat, mip, "forced", regions_specs) # Write data (nc file) outfilename = "PS_pr." + str(dfrq) + "_regrid." + nlon + "x" + nlat + "_" + dat + ".nc" custom_dataset = xr.merge([freqs, ps, rn, sig95]) @@ -74,7 +74,7 @@ def precip_variability_across_timescale( # Power spectum of anomaly freqs, ps, rn, sig95 = Powerspectrum(anom, nperseg, noverlap) # Domain & Frequency average - psdmfm_unforced = Avg_PS_DomFrq(ps, freqs, ntd, dat, mip, "unforced") + psdmfm_unforced = Avg_PS_DomFrq(ps, freqs, ntd, dat, mip, "unforced", regions_specs) # Write data (nc file) outfilename = "PS_pr." + str(dfrq) + "_regrid." + nlon + "x" + nlat + "_" + dat + "_unforced.nc" custom_dataset = xr.merge([freqs, ps, rn, sig95]) From dd37085c9964d709eea231e8395b8873438d8c22 Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Mon, 11 Sep 2023 15:10:37 -0700 Subject: [PATCH 23/55] not a list --- .../lib/lib_variability_across_timescales.py | 4 ++-- 1 file changed, 2 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 e3e039eb6..28f5e1ec4 100644 --- a/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py +++ b/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py @@ -413,8 +413,8 @@ def Avg_PS_DomFrq(d, frequency, ntd, dat, mip, frc, regions_specs): am = dmask.sel(lat=slice(-50, -30)).spatial.average("ps", axis=["X", "Y"], weights='generate')["ps"] else: # Custom region - lats=[regions_specs[dom]["domain"]["latitude"]] - lons=[regions_specs[dom]["domain"]["longitude"]] + lats=regions_specs[dom]["domain"]["latitude"] + lons=regions_specs[dom]["domain"]["longitude"] am = dmask.sel(lat=slice(lats[0],lats[1]),lon=slice(lons[0],lons[1])).spatial.average("ps", axis=["X", "Y"], weights='generate')["ps"] am = np.array(am) From d6ea6d69110a64163367be966baea06e81b52038 Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Mon, 11 Sep 2023 16:29:26 -0700 Subject: [PATCH 24/55] Fix region metrics selection --- .../lib/lib_variability_across_timescales.py | 7 +++---- 1 file changed, 3 insertions(+), 4 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 28f5e1ec4..a1d4ae863 100644 --- a/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py +++ b/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py @@ -412,10 +412,9 @@ def Avg_PS_DomFrq(d, frequency, ntd, dat, mip, frc, regions_specs): elif "50S30S" in dom: am = dmask.sel(lat=slice(-50, -30)).spatial.average("ps", axis=["X", "Y"], weights='generate')["ps"] else: - # Custom region - lats=regions_specs[dom]["domain"]["latitude"] - lons=regions_specs[dom]["domain"]["longitude"] - am = dmask.sel(lat=slice(lats[0],lats[1]),lon=slice(lons[0],lons[1])).spatial.average("ps", axis=["X", "Y"], weights='generate')["ps"] + # Custom region. Don't need to slice again because d only covered this region + am = dmask.spatial.average("ps", axis=["X", "Y"], weights='generate')["ps"] + am = np.array(am) for frq in frqs: From dcc9fa2f94d3ec135b32e9af8ab1d7ad92291482 Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Mon, 18 Sep 2023 13:21:32 -0700 Subject: [PATCH 25/55] Add region cropping --- .../lib/lib_variability_across_timescales.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) 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 a1d4ae863..f79267be1 100644 --- a/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py +++ b/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py @@ -18,7 +18,8 @@ # ================================================================================== def precip_variability_across_timescale( - file, syr, eyr, dfrq, mip, dat, var, fac, nperseg, noverlap, res, regions_specs, outdir, cmec + file, syr, eyr, dfrq, mip, dat, var, fac, nperseg, noverlap, res, regions_specs, \ + outdir, cmec, fshp, feature, attr ): """ Regridding -> Anomaly -> Power spectra -> Domain&Frequency average -> Write @@ -40,6 +41,8 @@ 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): rgtmp_ds = CropLatLon(rgtmp_ds, regions_specs) rgtmp = rgtmp_ds[var]*float(fac) From 617c4d07f74796e6fdfc973176d39b0d0f911ccd Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Mon, 18 Sep 2023 13:22:56 -0700 Subject: [PATCH 26/55] Add params for shpfile crop --- .../variability_across_timescales_PS_driver.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) 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 809bf3eab..87894180c 100644 --- a/pcmdi_metrics/precip_variability/variability_across_timescales_PS_driver.py +++ b/pcmdi_metrics/precip_variability/variability_across_timescales_PS_driver.py @@ -26,6 +26,9 @@ noverlap = param.noverlap res = param.res regions_specs = param.regions_specs +fshp = param.region_file +feature = param.feature +attr = param.attr print(modpath) print(mod) print(prd) @@ -60,5 +63,6 @@ syr = prd[0] eyr = prd[1] precip_variability_across_timescale( - file_list, syr, eyr, dfrq, mip, dat, var, fac, nperseg, noverlap, res, regions_specs, outdir, cmec + file_list, syr, eyr, dfrq, mip, dat, var, fac, nperseg, noverlap, res, regions_specs, \ + outdir, cmec, fshp, feature, attr ) From d22b6cbc6b4161c33d1ef89764a76c8da239bf85 Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Mon, 18 Sep 2023 14:04:10 -0700 Subject: [PATCH 27/55] Add args for region from shpfile --- .../lib/argparse_functions.py | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/pcmdi_metrics/precip_variability/lib/argparse_functions.py b/pcmdi_metrics/precip_variability/lib/argparse_functions.py index 9ab94b017..88ba7fea0 100644 --- a/pcmdi_metrics/precip_variability/lib/argparse_functions.py +++ b/pcmdi_metrics/precip_variability/lib/argparse_functions.py @@ -90,6 +90,24 @@ def AddParserArgument(P): action="store_false", help="Do not save CMEC format metrics JSON", ) + P.add_argument( + "--region_file", + dest="region_file", + help="File containing vector region of interest." + default=None + ) + P.add_argument( + "--feature", + dest="feature", + help="Feature in vectorized region." + default=None + ) + P.add_argument( + "--attr", + dest="attr", + help="Attribute containing feature in vectorized region." + default=None + ) P.set_defaults(cmec=False) return P From 370d6b223b21cd0444a9f96470e23930676e8db6 Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Mon, 18 Sep 2023 14:22:41 -0700 Subject: [PATCH 28/55] add import statement --- .../precip_variability/lib/lib_variability_across_timescales.py | 1 + 1 file changed, 1 insertion(+) 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 f79267be1..1540f85e8 100644 --- a/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py +++ b/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py @@ -11,6 +11,7 @@ import pcmdi_metrics from pcmdi_metrics.io import region_subset +from pcmdi_metrics.io.region_from_file import region_from_file import xarray as xr import xcdat From c320b303e935a8ede96e47b7964cf28198e193e9 Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Mon, 18 Sep 2023 14:23:37 -0700 Subject: [PATCH 29/55] Fix commas --- pcmdi_metrics/precip_variability/lib/argparse_functions.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/pcmdi_metrics/precip_variability/lib/argparse_functions.py b/pcmdi_metrics/precip_variability/lib/argparse_functions.py index 88ba7fea0..59376454e 100644 --- a/pcmdi_metrics/precip_variability/lib/argparse_functions.py +++ b/pcmdi_metrics/precip_variability/lib/argparse_functions.py @@ -93,19 +93,19 @@ def AddParserArgument(P): P.add_argument( "--region_file", dest="region_file", - help="File containing vector region of interest." + help="File containing vector region of interest.", default=None ) P.add_argument( "--feature", dest="feature", - help="Feature in vectorized region." + help="Feature in vectorized region.", default=None ) P.add_argument( "--attr", dest="attr", - help="Attribute containing feature in vectorized region." + help="Attribute containing feature in vectorized region.", default=None ) P.set_defaults(cmec=False) From 14a61f27757cf877a369cdcc397910fc1e14b0f8 Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Mon, 14 Aug 2023 16:16:55 -0700 Subject: [PATCH 30/55] 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 31/55] 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 32/55] 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 f11756e10cbf774c74889ceef86d3fa4ce9a689f Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Fri, 6 Oct 2023 10:37:57 -0700 Subject: [PATCH 33/55] Remove region cropping into regridding func --- .../lib/lib_variability_across_timescales.py | 6 +++--- 1 file changed, 3 insertions(+), 3 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 aa829dd23..c54f0dfeb 100644 --- a/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py +++ b/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py @@ -41,11 +41,9 @@ 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) + rgtmp_ds = RegridHoriz(do, var, res, regions_specs) 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): - rgtmp_ds = CropLatLon(rgtmp_ds, regions_specs) rgtmp = rgtmp_ds[var]*float(fac) if iyr == syr: drg = copy.deepcopy(rgtmp) @@ -122,6 +120,8 @@ def RegridHoriz(d, var, res): end_lon = 360.-res tgrid = grid.create_uniform_grid(start_lat,end_lat,res,start_lon,end_lon,res) + if regions_specs is not None: + tgrid = CropLatLon(tgrid, regions_specs) drg = d.regridder.horizontal(var, tgrid, tool="xesmf", method="conservative_normed",periodic=True,unmapped_to_nan=True) print("Complete regridding from", d[var].shape, "to", drg[var].shape) From fbed73256caa698c2d3a6f0ff27f853d9ed00857 Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Tue, 20 Feb 2024 10:22:06 -0800 Subject: [PATCH 34/55] 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) # ================================================================================== From 214a10e76647ae309d5da6af225173d49d23cf28 Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Fri, 1 Mar 2024 11:26:31 -0800 Subject: [PATCH 35/55] clean up and precommit --- pcmdi_metrics/io/__init__.py | 2 +- pcmdi_metrics/io/region_from_file.py | 4 +- .../lib/argparse_functions.py | 10 +- .../lib/lib_variability_across_timescales.py | 120 +++++++++++++----- ...variability_across_timescales_PS_driver.py | 24 +++- 5 files changed, 110 insertions(+), 50 deletions(-) diff --git a/pcmdi_metrics/io/__init__.py b/pcmdi_metrics/io/__init__.py index fa3fbf9d6..6f8e77148 100644 --- a/pcmdi_metrics/io/__init__.py +++ b/pcmdi_metrics/io/__init__.py @@ -5,6 +5,7 @@ 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 + from .xcdat_dataset_io import ( # noqa get_axis_list, get_data_list, @@ -22,4 +23,3 @@ get_time_key, select_subset, ) - diff --git a/pcmdi_metrics/io/region_from_file.py b/pcmdi_metrics/io/region_from_file.py index ff49aa277..bcb5144fe 100644 --- a/pcmdi_metrics/io/region_from_file.py +++ b/pcmdi_metrics/io/region_from_file.py @@ -1,5 +1,3 @@ -import shapely - import geopandas as gpd import regionmask @@ -14,7 +12,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) diff --git a/pcmdi_metrics/precip_variability/lib/argparse_functions.py b/pcmdi_metrics/precip_variability/lib/argparse_functions.py index f64d753c2..e71a34292 100644 --- a/pcmdi_metrics/precip_variability/lib/argparse_functions.py +++ b/pcmdi_metrics/precip_variability/lib/argparse_functions.py @@ -1,5 +1,6 @@ import ast + def AddParserArgument(P): P.add_argument( "--mip", type=str, dest="mip", default=None, help="cmip5, cmip6 or other mip" @@ -92,19 +93,16 @@ def AddParserArgument(P): "--region_file", dest="region_file", help="File containing vector region of interest.", - default=None + default=None, ) P.add_argument( - "--feature", - dest="feature", - help="Feature in vectorized region.", - default=None + "--feature", dest="feature", help="Feature in vectorized region.", default=None ) P.add_argument( "--attr", dest="attr", help="Attribute containing feature in vectorized region.", - default=None + default=None, ) P.set_defaults(cmec=False) 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 1656cc1cd..985d0a492 100644 --- a/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py +++ b/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py @@ -11,23 +11,38 @@ from scipy.stats import chi2 from xcdat.regridder import grid -import pcmdi_metrics from pcmdi_metrics.io import region_subset -from pcmdi_metrics.io.region_from_file import region_from_file from pcmdi_metrics.io.base import Base +from pcmdi_metrics.io.region_from_file import region_from_file from pcmdi_metrics.utils import create_land_sea_mask # ================================================================================== def precip_variability_across_timescale( - file, syr, eyr, dfrq, mip, dat, var, fac, nperseg, noverlap, res, regions_specs, \ - outdir, cmec, fshp, feature, attr + file, + syr, + eyr, + dfrq, + mip, + dat, + var, + fac, + nperseg, + noverlap, + res, + regions_specs, + outdir, + cmec, + fshp, + feature, + attr, ): """ Regridding -> Anomaly -> Power spectra -> Domain&Frequency average -> Write """ from time import perf_counter - t1_start = perf_counter() + + t1_start = perf_counter() psdmfm = {"RESULTS": {}} try: @@ -36,10 +51,10 @@ def precip_variability_across_timescale( f = xcdat.open_mfdataset(file, decode_times=False) cal = f.time.calendar # Add any calendar fixes here - cal = cal.replace("-","_") + cal = cal.replace("-", "_") f.time.attrs["calendar"] = cal f = xcdat.decode_time(f) - cal = f.time.encoding["calendar"] + cal = f.time.encoding["calendar"] if "360" in cal: ldy = 30 else: @@ -58,19 +73,19 @@ 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() + t2_start = perf_counter() rgtmp_ds = CropLatLon(rgtmp_ds, regions_specs) t2_stop = perf_counter() - print("Elapsed time cropping",t2_stop-t2_start) + 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) - rgtmp = rgtmp_ds[var]*float(fac) + rgtmp_ds = region_from_file(rgtmp_ds, fshp, attr, feature) + rgtmp = rgtmp_ds[var] * float(fac) if iyr == syr: drg = copy.deepcopy(rgtmp) else: - drg = xr.concat([drg, rgtmp], dim='time') - print(iyr, drg.shape) + drg = xr.concat([drg, rgtmp], dim="time") + print(iyr, drg.shape) nlon = str(len(drg.lon)) nlat = str(len(drg.lat)) f.close() @@ -89,7 +104,9 @@ def precip_variability_across_timescale( # Domain & Frequency average psdmfm_forced = Avg_PS_DomFrq(ps, freqs, ntd, dat, mip, "forced", regions_specs) # Write data (nc file) - outfilename = "PS_pr." + str(dfrq) + "_regrid." + nlon + "x" + nlat + "_" + dat + ".nc" + outfilename = ( + "PS_pr." + str(dfrq) + "_regrid." + nlon + "x" + nlat + "_" + dat + ".nc" + ) custom_dataset = xr.merge([freqs, ps, rn, sig95]) custom_dataset.to_netcdf( path=os.path.join( @@ -102,7 +119,17 @@ def precip_variability_across_timescale( # Domain & Frequency average psdmfm_unforced = Avg_PS_DomFrq(ps, freqs, ntd, dat, mip, "unforced", regions_specs) # Write data (nc file) - outfilename = "PS_pr." + str(dfrq) + "_regrid." + nlon + "x" + nlat + "_" + dat + "_unforced.nc" + outfilename = ( + "PS_pr." + + str(dfrq) + + "_regrid." + + nlon + + "x" + + nlat + + "_" + + dat + + "_unforced.nc" + ) custom_dataset = xr.merge([freqs, ps, rn, sig95]) custom_dataset.to_netcdf( path=os.path.join( @@ -116,7 +143,15 @@ def precip_variability_across_timescale( psdmfm["RESULTS"][dat]["unforced"] = psdmfm_unforced outfilename = ( - "PS_pr." + str(dfrq) + "_regrid." + nlon + "x" + nlat + "_area.freq.mean_" + dat + ".json" + "PS_pr." + + str(dfrq) + + "_regrid." + + nlon + + "x" + + nlat + + "_area.freq.mean_" + + dat + + ".json" ) JSON = Base(outdir.replace("%(output_type)", "metrics_results"), outfilename) JSON.write( @@ -129,11 +164,11 @@ 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) + print("Elapsed time cropping", t1_stop - t1_start) # ================================================================================== -def RegridHoriz(d, var, res): +def RegridHoriz(d, var, res, regions_specs=None): """ Regrid to 2deg (180lon*90lat) horizontal resolution Input @@ -143,20 +178,28 @@ def RegridHoriz(d, var, res): - drg: xCDAT variable with 2deg horizontal resolution """ - start_lat=-90.+res/2. - start_lon=0. - end_lat = 90.-res/2. - end_lon = 360.-res + start_lat = -90.0 + res / 2.0 + start_lon = 0.0 + end_lat = 90.0 - res / 2.0 + end_lon = 360.0 - res - tgrid = grid.create_uniform_grid(start_lat,end_lat,res,start_lon,end_lon,res) + tgrid = grid.create_uniform_grid(start_lat, end_lat, res, start_lon, end_lon, res) if regions_specs is not None: tgrid = CropLatLon(tgrid, regions_specs) - drg = d.regridder.horizontal(var, tgrid, tool="xesmf", method="conservative_normed",periodic=True,unmapped_to_nan=True) - + drg = d.regridder.horizontal( + var, + tgrid, + tool="xesmf", + method="conservative_normed", + periodic=True, + unmapped_to_nan=True, + ) + print("Complete regridding from", d[var].shape, "to", drg[var].shape) return drg + # ================================================================================== def CropLatLon(d, regions_specs): """ @@ -173,6 +216,7 @@ def CropLatLon(d, regions_specs): return dnew + # ================================================================================== def ClimAnom(d, ntd, syr, eyr, cal): """ @@ -445,10 +489,10 @@ def Avg_PS_DomFrq(d, frequency, ntd, dat, mip, frc, regions_specs): for dom in domains: psdmfm[dom] = {} - if "Ocean" in dom or regions_specs[dom].get("value",-1)==0: - dmask = d.where(mask==0) - elif "Land" in dom or regions_specs[dom].get("value",-1)==1: - dmask = d.where(mask==1) + if "Ocean" in dom or regions_specs[dom].get("value", -1) == 0: + dmask = d.where(mask == 0) + elif "Land" in dom or regions_specs[dom].get("value", -1) == 1: + dmask = d.where(mask == 1) else: dmask = d @@ -457,16 +501,24 @@ def Avg_PS_DomFrq(d, frequency, ntd, dat, mip, frc, regions_specs): dmask = dmask.bounds.add_bounds(axis="Y") if "50S50N" in dom: - am = dmask.sel(lat=slice(-50, 50)).spatial.average("ps", axis=["X", "Y"], weights='generate')["ps"] + am = dmask.sel(lat=slice(-50, 50)).spatial.average( + "ps", axis=["X", "Y"], weights="generate" + )["ps"] elif "30N50N" in dom: - am = dmask.sel(lat=slice(30, 50)).spatial.average("ps", axis=["X", "Y"], weights='generate')["ps"] + am = dmask.sel(lat=slice(30, 50)).spatial.average( + "ps", axis=["X", "Y"], weights="generate" + )["ps"] elif "30S30N" in dom: - am = dmask.sel(lat=slice(-30, 30)).spatial.average("ps", axis=["X", "Y"], weights='generate')["ps"] + am = dmask.sel(lat=slice(-30, 30)).spatial.average( + "ps", axis=["X", "Y"], weights="generate" + )["ps"] elif "50S30S" in dom: - am = dmask.sel(lat=slice(-50, -30)).spatial.average("ps", axis=["X", "Y"], weights='generate')["ps"] + am = dmask.sel(lat=slice(-50, -30)).spatial.average( + "ps", axis=["X", "Y"], weights="generate" + )["ps"] else: # Custom region. Don't need to slice again because d only covered this region - am = dmask.spatial.average("ps", axis=["X", "Y"], weights='generate')["ps"] + am = dmask.spatial.average("ps", axis=["X", "Y"], weights="generate")["ps"] am = np.array(am) 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 68ed416be..5a9b4cdae 100644 --- a/pcmdi_metrics/precip_variability/variability_across_timescales_PS_driver.py +++ b/pcmdi_metrics/precip_variability/variability_across_timescales_PS_driver.py @@ -1,7 +1,4 @@ #!/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 @@ -65,6 +62,21 @@ syr = prd[0] eyr = prd[1] precip_variability_across_timescale( - file_list, syr, eyr, dfrq, mip, dat, var, fac, nperseg, noverlap, res, regions_specs, \ - outdir_template, cmec, fshp, feature, attr - ) + file_list, + syr, + eyr, + dfrq, + mip, + dat, + var, + fac, + nperseg, + noverlap, + res, + regions_specs, + outdir_template, + cmec, + fshp, + feature, + attr, +) From e37bafd975f9820f9fb37861043c0ff38375fc38 Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Fri, 1 Mar 2024 11:28:30 -0800 Subject: [PATCH 36/55] remove timing --- .../lib/lib_variability_across_timescales.py | 8 -------- 1 file changed, 8 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 985d0a492..c900150ed 100644 --- a/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py +++ b/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py @@ -40,9 +40,6 @@ def precip_variability_across_timescale( """ Regridding -> Anomaly -> Power spectra -> Domain&Frequency average -> Write """ - from time import perf_counter - - t1_start = perf_counter() psdmfm = {"RESULTS": {}} try: @@ -73,10 +70,7 @@ 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) @@ -163,8 +157,6 @@ 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) # ================================================================================== From 157adeac018a924a6e07905c50bc110bf67c98a2 Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Fri, 1 Mar 2024 11:32:26 -0800 Subject: [PATCH 37/55] fix white space --- share/DefArgsCIA.json | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/share/DefArgsCIA.json b/share/DefArgsCIA.json index cd33a055d..8507f33ba 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 4df78facf89723f03b2cd26d91b5678e75187858 Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Tue, 12 Mar 2024 14:48:43 -0700 Subject: [PATCH 38/55] add figure scripts --- .../calc_ps_area.freq.mean.regional.py | 140 +++++++++++++++ .../calc_ps_area.mean.regional.py | 98 +++++++++++ .../calc_ps_freq.mean.regional.py | 165 ++++++++++++++++++ 3 files changed, 403 insertions(+) create mode 100644 pcmdi_metrics/precip_variability/scripts_pcmdi/calc_ps_area.freq.mean.regional.py create mode 100644 pcmdi_metrics/precip_variability/scripts_pcmdi/calc_ps_area.mean.regional.py create mode 100644 pcmdi_metrics/precip_variability/scripts_pcmdi/calc_ps_freq.mean.regional.py diff --git a/pcmdi_metrics/precip_variability/scripts_pcmdi/calc_ps_area.freq.mean.regional.py b/pcmdi_metrics/precip_variability/scripts_pcmdi/calc_ps_area.freq.mean.regional.py new file mode 100644 index 000000000..bc71f4a45 --- /dev/null +++ b/pcmdi_metrics/precip_variability/scripts_pcmdi/calc_ps_area.freq.mean.regional.py @@ -0,0 +1,140 @@ +import numpy as np +import datetime +import pickle +import json +import copy +import os + +def prdday_to_frq3hridx(prdday,frequency): + frq3hr=1./(float(prdday)*8.) + idx = (np.abs(frequency-frq3hr)).argmin() + return int(idx) + +def prdday_to_frq1didx(prdday,frequency): + frq24hr=1./(float(prdday)) + idx = (np.abs(frequency-frq24hr)).argmin() + return int(idx) + +#-------------------- +# User settings here +#-------------------- +# "3hr" or "day" +hr='' +# Input file name (pickle .pkl output from calc_ps_area.mean.py) +fname = '' +# Output directory (not including version) +outdir = '' + +#----------------------- +# The rest of the script +#----------------------- +ver = datetime.datetime.now().strftime("v%Y%m%d") + +if hr == 'day': + frqs_forced=['semi-annual','annual'] + frqs_unforced=['synoptic','sub-seasonal','seasonal-annual','interannual'] +elif hr == '3hr': + frqs_forced=['semi-diurnal','diurnal','semi-annual','annual'] + frqs_unforced=['sub-daily','synoptic','sub-seasonal','seasonal-annual','interannual'] + +infile = open(fname, 'rb') +psdm = pickle.load(infile) + +psdmfm=copy.deepcopy(psdm) +for frc in psdm.keys(): + if (frc=='forced'): + frqs=frqs_forced + elif (frc=='unforced'): + frqs=frqs_unforced + print(frc) + for mip in psdm[frc].keys(): + print(mip) + for dat in psdm[frc][mip].keys(): + frequency=np.array(psdm[frc][mip][dat]['freqs']) + del psdm[frc][mip][dat]['freqs'] + del psdmfm[frc][mip][dat]['freqs'] + + for var in psdm[frc][mip][dat].keys(): + print(dat,var) + for idm, dom in enumerate(psdm[frc][mip][dat][var].keys()): + am=np.array(psdm[frc][mip][dat][var][dom]) + del psdmfm[frc][mip][dat][var][dom] + psdmfm[frc][mip][dat][var][dom] = {} + print(dom) + for frq in frqs: + print(frq) + if (frq=='semi-diurnal'): # pr=0.5day + idx=prdday_to_frq1didx(0.5,frequency) + amfm = am[idx] + elif (frq=='diurnal'): # pr=1day + idx=prdday_to_frq1didx(1,frequency) + amfm = am[idx] + if (frq=='semi-annual'): # 180day= Date: Tue, 12 Mar 2024 14:52:42 -0700 Subject: [PATCH 39/55] run precommit --- .../calc_ps_area.freq.mean.regional.py | 200 ++++++------ .../calc_ps_area.mean.regional.py | 115 +++---- .../calc_ps_freq.mean.regional.py | 284 ++++++++++-------- 3 files changed, 329 insertions(+), 270 deletions(-) diff --git a/pcmdi_metrics/precip_variability/scripts_pcmdi/calc_ps_area.freq.mean.regional.py b/pcmdi_metrics/precip_variability/scripts_pcmdi/calc_ps_area.freq.mean.regional.py index bc71f4a45..2d1aabbe7 100644 --- a/pcmdi_metrics/precip_variability/scripts_pcmdi/calc_ps_area.freq.mean.regional.py +++ b/pcmdi_metrics/precip_variability/scripts_pcmdi/calc_ps_area.freq.mean.regional.py @@ -1,140 +1,152 @@ -import numpy as np +import copy import datetime -import pickle import json -import copy import os +import pickle + +import numpy as np -def prdday_to_frq3hridx(prdday,frequency): - frq3hr=1./(float(prdday)*8.) - idx = (np.abs(frequency-frq3hr)).argmin() + +def prdday_to_frq3hridx(prdday, frequency): + frq3hr = 1.0 / (float(prdday) * 8.0) + idx = (np.abs(frequency - frq3hr)).argmin() return int(idx) -def prdday_to_frq1didx(prdday,frequency): - frq24hr=1./(float(prdday)) - idx = (np.abs(frequency-frq24hr)).argmin() + +def prdday_to_frq1didx(prdday, frequency): + frq24hr = 1.0 / (float(prdday)) + idx = (np.abs(frequency - frq24hr)).argmin() return int(idx) -#-------------------- + +# -------------------- # User settings here -#-------------------- +# -------------------- # "3hr" or "day" -hr='' +hr = "" # Input file name (pickle .pkl output from calc_ps_area.mean.py) -fname = '' +fname = "" # Output directory (not including version) -outdir = '' +outdir = "" -#----------------------- +# ----------------------- # The rest of the script -#----------------------- +# ----------------------- ver = datetime.datetime.now().strftime("v%Y%m%d") -if hr == 'day': - frqs_forced=['semi-annual','annual'] - frqs_unforced=['synoptic','sub-seasonal','seasonal-annual','interannual'] -elif hr == '3hr': - frqs_forced=['semi-diurnal','diurnal','semi-annual','annual'] - frqs_unforced=['sub-daily','synoptic','sub-seasonal','seasonal-annual','interannual'] +if hr == "day": + frqs_forced = ["semi-annual", "annual"] + frqs_unforced = ["synoptic", "sub-seasonal", "seasonal-annual", "interannual"] +elif hr == "3hr": + frqs_forced = ["semi-diurnal", "diurnal", "semi-annual", "annual"] + frqs_unforced = [ + "sub-daily", + "synoptic", + "sub-seasonal", + "seasonal-annual", + "interannual", + ] -infile = open(fname, 'rb') +infile = open(fname, "rb") psdm = pickle.load(infile) -psdmfm=copy.deepcopy(psdm) +psdmfm = copy.deepcopy(psdm) for frc in psdm.keys(): - if (frc=='forced'): - frqs=frqs_forced - elif (frc=='unforced'): - frqs=frqs_unforced + if frc == "forced": + frqs = frqs_forced + elif frc == "unforced": + frqs = frqs_unforced print(frc) for mip in psdm[frc].keys(): print(mip) for dat in psdm[frc][mip].keys(): - frequency=np.array(psdm[frc][mip][dat]['freqs']) - del psdm[frc][mip][dat]['freqs'] - del psdmfm[frc][mip][dat]['freqs'] + frequency = np.array(psdm[frc][mip][dat]["freqs"]) + del psdm[frc][mip][dat]["freqs"] + del psdmfm[frc][mip][dat]["freqs"] for var in psdm[frc][mip][dat].keys(): - print(dat,var) + print(dat, var) for idm, dom in enumerate(psdm[frc][mip][dat][var].keys()): - am=np.array(psdm[frc][mip][dat][var][dom]) + am = np.array(psdm[frc][mip][dat][var][dom]) del psdmfm[frc][mip][dat][var][dom] psdmfm[frc][mip][dat][var][dom] = {} print(dom) for frq in frqs: print(frq) - if (frq=='semi-diurnal'): # pr=0.5day - idx=prdday_to_frq1didx(0.5,frequency) + if frq == "semi-diurnal": # pr=0.5day + idx = prdday_to_frq1didx(0.5, frequency) amfm = am[idx] - elif (frq=='diurnal'): # pr=1day - idx=prdday_to_frq1didx(1,frequency) + elif frq == "diurnal": # pr=1day + idx = prdday_to_frq1didx(1, frequency) amfm = am[idx] - if (frq=='semi-annual'): # 180day= Date: Tue, 4 Jun 2024 11:47:28 -0700 Subject: [PATCH 40/55] Update __init__.py --- pcmdi_metrics/io/__init__.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/pcmdi_metrics/io/__init__.py b/pcmdi_metrics/io/__init__.py index 9addbd94b..3e2b7d7c9 100644 --- a/pcmdi_metrics/io/__init__.py +++ b/pcmdi_metrics/io/__init__.py @@ -3,9 +3,6 @@ from .string_constructor import StringConstructor, fill_template # noqa # isort:skip from . import base # noqa 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 from .xcdat_dataset_io import ( # noqa # isort:skip da_to_ds, From 807dab90369ce4b00364fae8311fa37e1e9f3f82 Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Mon, 17 Jun 2024 13:41:18 -0700 Subject: [PATCH 41/55] update lat lon vars --- .../lib/lib_variability_across_timescales.py | 35 +++++++++++++++---- 1 file changed, 29 insertions(+), 6 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 c900150ed..bc6dbc37b 100644 --- a/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py +++ b/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py @@ -80,8 +80,11 @@ def precip_variability_across_timescale( else: drg = xr.concat([drg, rgtmp], dim="time") print(iyr, drg.shape) - nlon = str(len(drg.lon)) - nlat = str(len(drg.lat)) + + xvar = lib.find_lon(drg) + yvar = lib.find_lat(drg) + nlon = str(len(xvar)) + nlat = str(len(yvar)) f.close() # Anomaly @@ -158,6 +161,25 @@ def precip_variability_across_timescale( if cmec: JSON.write_cmec(indent=4, separators=(",", ": ")) +# ================================================================================== +def find_lon(ds): + for key in ds.coords: + if key in ["lon", "longitude"]: + return key + for key in ds.keys(): + if key in ["lon", "longitude"]: + return key + return None + +# ================================================================================== +def find_lat(ds): + for key in ds.coords: + if key in ["lat", "latitude"]: + return key + for key in ds.keys(): + if key in ["lat", "latitude"]: + return key + return None # ================================================================================== def RegridHoriz(d, var, res, regions_specs=None): @@ -476,6 +498,7 @@ def Avg_PS_DomFrq(d, frequency, ntd, dat, mip, frc, regions_specs): # generate land sea mask mask = create_land_sea_mask(d[0]) + yvar = lib.find_lat(drg) psdmfm = {} for dom in domains: @@ -493,19 +516,19 @@ def Avg_PS_DomFrq(d, frequency, ntd, dat, mip, frc, regions_specs): dmask = dmask.bounds.add_bounds(axis="Y") if "50S50N" in dom: - am = dmask.sel(lat=slice(-50, 50)).spatial.average( + am = dmask.sel({yvar: slice(-50, 50)}).spatial.average( "ps", axis=["X", "Y"], weights="generate" )["ps"] elif "30N50N" in dom: - am = dmask.sel(lat=slice(30, 50)).spatial.average( + am = dmask.sel({yvar: slice(30, 50)}).spatial.average( "ps", axis=["X", "Y"], weights="generate" )["ps"] elif "30S30N" in dom: - am = dmask.sel(lat=slice(-30, 30)).spatial.average( + am = dmask.sel({yvar: slice(-30, 30)}).spatial.average( "ps", axis=["X", "Y"], weights="generate" )["ps"] elif "50S30S" in dom: - am = dmask.sel(lat=slice(-50, -30)).spatial.average( + am = dmask.sel({yvar: slice(-50, -30)}).spatial.average( "ps", axis=["X", "Y"], weights="generate" )["ps"] else: From 5974ccaf1098642d7b84b6da8163d7313f97d1d7 Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Mon, 17 Jun 2024 13:46:34 -0700 Subject: [PATCH 42/55] move check --- .../lib/lib_variability_across_timescales.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 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 bc6dbc37b..8b4c33928 100644 --- a/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py +++ b/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py @@ -41,6 +41,13 @@ def precip_variability_across_timescale( Regridding -> Anomaly -> Power spectra -> Domain&Frequency average -> Write """ psdmfm = {"RESULTS": {}} + + if dfrq == "day": + ntd = 1 + elif dfrq == "3hr": + ntd = 8 + else: + sys.exit("ERROR: dfrq " + dfrq + " is not defined!") try: f = xcdat.open_mfdataset(file) @@ -88,12 +95,6 @@ def precip_variability_across_timescale( f.close() # Anomaly - if dfrq == "day": - ntd = 1 - elif dfrq == "3hr": - ntd = 8 - else: - sys.exit("ERROR: dfrq " + dfrq + " is not defined!") clim, anom = ClimAnom(drg, ntd, syr, eyr, cal) # Power spectum of total From 6510a2656a1fddbafd1fca2878864e3f79b27dc0 Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Mon, 17 Jun 2024 13:56:46 -0700 Subject: [PATCH 43/55] run precommit --- .../lib/lib_variability_across_timescales.py | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 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 8b4c33928..4f51cbd6b 100644 --- a/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py +++ b/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py @@ -41,7 +41,7 @@ def precip_variability_across_timescale( Regridding -> Anomaly -> Power spectra -> Domain&Frequency average -> Write """ psdmfm = {"RESULTS": {}} - + if dfrq == "day": ntd = 1 elif dfrq == "3hr": @@ -87,9 +87,9 @@ def precip_variability_across_timescale( else: drg = xr.concat([drg, rgtmp], dim="time") print(iyr, drg.shape) - - xvar = lib.find_lon(drg) - yvar = lib.find_lat(drg) + + xvar = find_lon(drg) + yvar = find_lat(drg) nlon = str(len(xvar)) nlat = str(len(yvar)) f.close() @@ -162,6 +162,7 @@ def precip_variability_across_timescale( if cmec: JSON.write_cmec(indent=4, separators=(",", ": ")) + # ================================================================================== def find_lon(ds): for key in ds.coords: @@ -172,6 +173,7 @@ def find_lon(ds): return key return None + # ================================================================================== def find_lat(ds): for key in ds.coords: @@ -182,6 +184,7 @@ def find_lat(ds): return key return None + # ================================================================================== def RegridHoriz(d, var, res, regions_specs=None): """ @@ -499,7 +502,7 @@ def Avg_PS_DomFrq(d, frequency, ntd, dat, mip, frc, regions_specs): # generate land sea mask mask = create_land_sea_mask(d[0]) - yvar = lib.find_lat(drg) + yvar = find_lat(d) psdmfm = {} for dom in domains: From 6630229bd632c24b64fa13afeb6252da33e55dd3 Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Mon, 17 Jun 2024 14:59:25 -0700 Subject: [PATCH 44/55] fix regino --- .../lib/lib_variability_across_timescales.py | 19 +++++++++++++++---- 1 file changed, 15 insertions(+), 4 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 4f51cbd6b..8367a6a3f 100644 --- a/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py +++ b/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py @@ -16,6 +16,8 @@ from pcmdi_metrics.io.region_from_file import region_from_file from pcmdi_metrics.utils import create_land_sea_mask +import resource + # ================================================================================== def precip_variability_across_timescale( @@ -40,6 +42,7 @@ def precip_variability_across_timescale( """ Regridding -> Anomaly -> Power spectra -> Domain&Frequency average -> Write """ + print("Start:",resource.getrusage(resource.RUSAGE_SELF).ru_maxrss) psdmfm = {"RESULTS": {}} if dfrq == "day": @@ -87,6 +90,7 @@ def precip_variability_across_timescale( else: drg = xr.concat([drg, rgtmp], dim="time") print(iyr, drg.shape) + print("Mem:",resource.getrusage(resource.RUSAGE_SELF).ru_maxrss) xvar = find_lon(drg) yvar = find_lat(drg) @@ -230,7 +234,7 @@ def CropLatLon(d, regions_specs): """ region_name = list(regions_specs.keys())[0] - dnew = region_subset(d, regions_specs, region=region_name) + dnew = region_subset(d, region_name, None, regions_specs, False) return dnew @@ -249,7 +253,8 @@ def ClimAnom(d, ntd, syr, eyr, cal): - clim: climatology (climatological diurnal and annual cycles) - anom: anomaly departure from the climatological diurnal and annual cycles """ - + print("ClimAnom Start:",resource.getrusage(resource.RUSAGE_SELF).ru_maxrss) + # Year segment nyr = eyr - syr + 1 if "gregorian" in cal or "standard" in cal: @@ -301,7 +306,7 @@ def ClimAnom(d, ntd, syr, eyr, cal): # Climatology clim = np.nanmean(dseg, axis=0) - + print("ClimAnom post-mean:",resource.getrusage(resource.RUSAGE_SELF).ru_maxrss) # Anomaly anom = np.array([]) if "gregorian" in cal or "standard" in cal: @@ -325,8 +330,11 @@ def ClimAnom(d, ntd, syr, eyr, cal): anom = np.append(anom, (dseg[iyr] - clim)) # Reahape and Dimension information + print("Starting reshape") clim = np.reshape(clim, (ndy * ntd, d.shape[1], d.shape[2])) + print("ClimAnom post clim reshape:",resource.getrusage(resource.RUSAGE_SELF).ru_maxrss) anom = np.reshape(anom, (d.shape[0], d.shape[1], d.shape[2])) + print("ClimAnom post Anom reshape:",resource.getrusage(resource.RUSAGE_SELF).ru_maxrss) climtime = pd.period_range( start="0001-01-01", periods=clim.shape[0], freq=str(24 / ntd) + "H" @@ -354,6 +362,7 @@ def Powerspectrum(d, nperseg, noverlap): - rn: Rednoise - sig95: 95% rednoise confidence level """ + print("Powerspectrum Start:",resource.getrusage(resource.RUSAGE_SELF).ru_maxrss) # Fill missing date using interpolation dnp = np.array(d) dfm = np.zeros((d.shape[0], d.shape[1], d.shape[2]), dtype=float) @@ -364,9 +373,11 @@ def Powerspectrum(d, nperseg, noverlap): dfm[:, iy, ix] = np.array(ypfm) # Calculate power spectrum + print("Pre-welch:",resource.getrusage(resource.RUSAGE_SELF).ru_maxrss) freqs, psd = signal.welch( dfm, scaling="spectrum", nperseg=nperseg, noverlap=noverlap, axis=0 ) + print("Post-welch:",resource.getrusage(resource.RUSAGE_SELF).ru_maxrss) # Signigicance test of power spectra (from J. Lee's MOV code) nps = max( @@ -380,7 +391,7 @@ def Powerspectrum(d, nperseg, noverlap): rn[:, iy, ix] = rednoise(psd[:, iy, ix], len(freqs), r1) nu = 2 * nps sig95[:, iy, ix] = RedNoiseSignificanceLevel(nu, rn[:, iy, ix]) - + print("Post-sig:",resource.getrusage(resource.RUSAGE_SELF).ru_maxrss) # Decorate arrays with dimensions # axisfrq = np.arange(len(freqs)) axisfrq = range(len(freqs)) From 5b197c2bf63b03a1c25ad1a7c9fcdcb2e4bf7687 Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Tue, 25 Jun 2024 15:00:20 -0700 Subject: [PATCH 45/55] fix lon issues --- .../lib/lib_variability_across_timescales.py | 59 ++++++------------- 1 file changed, 19 insertions(+), 40 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 8367a6a3f..da8f90b2e 100644 --- a/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py +++ b/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py @@ -16,8 +16,6 @@ from pcmdi_metrics.io.region_from_file import region_from_file from pcmdi_metrics.utils import create_land_sea_mask -import resource - # ================================================================================== def precip_variability_across_timescale( @@ -42,7 +40,6 @@ def precip_variability_across_timescale( """ Regridding -> Anomaly -> Power spectra -> Domain&Frequency average -> Write """ - print("Start:",resource.getrusage(resource.RUSAGE_SELF).ru_maxrss) psdmfm = {"RESULTS": {}} if dfrq == "day": @@ -68,34 +65,22 @@ def precip_variability_across_timescale( ldy = 31 print(dat, cal) print("syr, eyr:", syr, eyr) - for iyr in range(syr, eyr + 1): - print(iyr) - 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) - 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) - else: - drg = xr.concat([drg, rgtmp], dim="time") - print(iyr, drg.shape) - print("Mem:",resource.getrusage(resource.RUSAGE_SELF).ru_maxrss) - xvar = find_lon(drg) - yvar = find_lat(drg) + # Regridding + xvar = find_lon(f) + yvar = find_lat(f) nlon = str(len(xvar)) nlat = str(len(yvar)) + if np.min(f[xvar]) < 0: + lon_range = [-180.,180.] + else: + lon_range = [0.,360.] + rgtmp = RegridHoriz(f, var, res, regions_specs, lon_range) + if fshp is not None: + print("Cropping from shapefile") + rgtmp = region_from_file(rgtmp, fshp, attr, feature) + drg = rgtmp[var] * float(fac) + f.close() # Anomaly @@ -190,7 +175,7 @@ def find_lat(ds): # ================================================================================== -def RegridHoriz(d, var, res, regions_specs=None): +def RegridHoriz(d, var, res, regions_specs=None, lon_range=[0.,360.]): """ Regrid to 2deg (180lon*90lat) horizontal resolution Input @@ -201,9 +186,9 @@ def RegridHoriz(d, var, res, regions_specs=None): """ start_lat = -90.0 + res / 2.0 - start_lon = 0.0 + start_lon = lon_range[0] end_lat = 90.0 - res / 2.0 - end_lon = 360.0 - res + end_lon = lon_range[1] - res tgrid = grid.create_uniform_grid(start_lat, end_lat, res, start_lon, end_lon, res) if regions_specs is not None: @@ -253,7 +238,6 @@ def ClimAnom(d, ntd, syr, eyr, cal): - clim: climatology (climatological diurnal and annual cycles) - anom: anomaly departure from the climatological diurnal and annual cycles """ - print("ClimAnom Start:",resource.getrusage(resource.RUSAGE_SELF).ru_maxrss) # Year segment nyr = eyr - syr + 1 @@ -306,7 +290,7 @@ def ClimAnom(d, ntd, syr, eyr, cal): # Climatology clim = np.nanmean(dseg, axis=0) - print("ClimAnom post-mean:",resource.getrusage(resource.RUSAGE_SELF).ru_maxrss) + # Anomaly anom = np.array([]) if "gregorian" in cal or "standard" in cal: @@ -330,11 +314,8 @@ def ClimAnom(d, ntd, syr, eyr, cal): anom = np.append(anom, (dseg[iyr] - clim)) # Reahape and Dimension information - print("Starting reshape") clim = np.reshape(clim, (ndy * ntd, d.shape[1], d.shape[2])) - print("ClimAnom post clim reshape:",resource.getrusage(resource.RUSAGE_SELF).ru_maxrss) anom = np.reshape(anom, (d.shape[0], d.shape[1], d.shape[2])) - print("ClimAnom post Anom reshape:",resource.getrusage(resource.RUSAGE_SELF).ru_maxrss) climtime = pd.period_range( start="0001-01-01", periods=clim.shape[0], freq=str(24 / ntd) + "H" @@ -362,7 +343,7 @@ def Powerspectrum(d, nperseg, noverlap): - rn: Rednoise - sig95: 95% rednoise confidence level """ - print("Powerspectrum Start:",resource.getrusage(resource.RUSAGE_SELF).ru_maxrss) + # Fill missing date using interpolation dnp = np.array(d) dfm = np.zeros((d.shape[0], d.shape[1], d.shape[2]), dtype=float) @@ -373,11 +354,9 @@ def Powerspectrum(d, nperseg, noverlap): dfm[:, iy, ix] = np.array(ypfm) # Calculate power spectrum - print("Pre-welch:",resource.getrusage(resource.RUSAGE_SELF).ru_maxrss) freqs, psd = signal.welch( dfm, scaling="spectrum", nperseg=nperseg, noverlap=noverlap, axis=0 ) - print("Post-welch:",resource.getrusage(resource.RUSAGE_SELF).ru_maxrss) # Signigicance test of power spectra (from J. Lee's MOV code) nps = max( @@ -391,7 +370,7 @@ def Powerspectrum(d, nperseg, noverlap): rn[:, iy, ix] = rednoise(psd[:, iy, ix], len(freqs), r1) nu = 2 * nps sig95[:, iy, ix] = RedNoiseSignificanceLevel(nu, rn[:, iy, ix]) - print("Post-sig:",resource.getrusage(resource.RUSAGE_SELF).ru_maxrss) + # Decorate arrays with dimensions # axisfrq = np.arange(len(freqs)) axisfrq = range(len(freqs)) From 43f8de5004241cee253fe3f6796105be6666d0da Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Tue, 25 Jun 2024 15:02:36 -0700 Subject: [PATCH 46/55] fix lon --- .../lib/lib_variability_across_timescales.py | 13 ++++--------- 1 file changed, 4 insertions(+), 9 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 da8f90b2e..75689f24a 100644 --- a/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py +++ b/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py @@ -1,4 +1,3 @@ -import copy import math import os import sys @@ -59,10 +58,6 @@ def precip_variability_across_timescale( f.time.attrs["calendar"] = cal f = xcdat.decode_time(f) cal = f.time.encoding["calendar"] - if "360" in cal: - ldy = 30 - else: - ldy = 31 print(dat, cal) print("syr, eyr:", syr, eyr) @@ -72,9 +67,9 @@ def precip_variability_across_timescale( nlon = str(len(xvar)) nlat = str(len(yvar)) if np.min(f[xvar]) < 0: - lon_range = [-180.,180.] + lon_range = [-180.0, 180.0] else: - lon_range = [0.,360.] + lon_range = [0.0, 360.0] rgtmp = RegridHoriz(f, var, res, regions_specs, lon_range) if fshp is not None: print("Cropping from shapefile") @@ -175,7 +170,7 @@ def find_lat(ds): # ================================================================================== -def RegridHoriz(d, var, res, regions_specs=None, lon_range=[0.,360.]): +def RegridHoriz(d, var, res, regions_specs=None, lon_range=[0.0, 360.0]): """ Regrid to 2deg (180lon*90lat) horizontal resolution Input @@ -238,7 +233,7 @@ def ClimAnom(d, ntd, syr, eyr, cal): - clim: climatology (climatological diurnal and annual cycles) - anom: anomaly departure from the climatological diurnal and annual cycles """ - + # Year segment nyr = eyr - syr + 1 if "gregorian" in cal or "standard" in cal: From cd1327fa9546a1b3b46ee098bf12477a6bff4369 Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Thu, 1 Aug 2024 13:50:04 -0700 Subject: [PATCH 47/55] update syntax --- .../precip_variability/lib/lib_variability_across_timescales.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 75689f24a..23a2097cf 100644 --- a/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py +++ b/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py @@ -214,7 +214,7 @@ def CropLatLon(d, regions_specs): """ region_name = list(regions_specs.keys())[0] - dnew = region_subset(d, region_name, None, regions_specs, False) + dnew = region_subset(d, region_name, regions_specs=regions_specs) return dnew From a81c78533598c4dcaf87ae4ee4d72e8606f254c4 Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Thu, 1 Aug 2024 14:05:29 -0700 Subject: [PATCH 48/55] Update params --- docs/metrics_precip-variability.rst | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/docs/metrics_precip-variability.rst b/docs/metrics_precip-variability.rst index 6655ef14a..10374c892 100644 --- a/docs/metrics_precip-variability.rst +++ b/docs/metrics_precip-variability.rst @@ -62,6 +62,10 @@ Options available to set in the parameter file include: * **nperseg**: Length of segment in power spectra. * **noverlap**: Length of overlap between segments in power spectra. * **ref**: Reference data path. +* **regions_specs**: Dictionary containing region bounding box. Uses format {"your region name": {"domain": {"latitude": (min, max), "longitude": (min, max)}}}. Min and max should be replaced with the values that define the region. +* **region_file**: Path to a shapefile containing vector region outline. Must use with **attr** and **feature** parameters. +* **attr**: Attribute used to identify region in shapefile (eg, column of attribute table). For example, "COUNTRY" in a shapefile of countries. +* **feature**: Unique feature value of the region that occurs in the attribute given by "--attr". Must match only one geometry in the shapefile. An example is a feature called "USA" under the attribute "COUNTRY". * **cmec**: Set to True to output CMEC formatted JSON. Metric From 172d8f21fa11b5cece953410583c90161ac82898 Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Thu, 1 Aug 2024 14:13:55 -0700 Subject: [PATCH 49/55] remove files --- .../calc_ps_area.freq.mean.regional.py | 152 ------------- .../calc_ps_area.mean.regional.py | 103 --------- .../calc_ps_freq.mean.regional.py | 207 ------------------ 3 files changed, 462 deletions(-) delete mode 100644 pcmdi_metrics/precip_variability/scripts_pcmdi/calc_ps_area.freq.mean.regional.py delete mode 100644 pcmdi_metrics/precip_variability/scripts_pcmdi/calc_ps_area.mean.regional.py delete mode 100644 pcmdi_metrics/precip_variability/scripts_pcmdi/calc_ps_freq.mean.regional.py diff --git a/pcmdi_metrics/precip_variability/scripts_pcmdi/calc_ps_area.freq.mean.regional.py b/pcmdi_metrics/precip_variability/scripts_pcmdi/calc_ps_area.freq.mean.regional.py deleted file mode 100644 index 2d1aabbe7..000000000 --- a/pcmdi_metrics/precip_variability/scripts_pcmdi/calc_ps_area.freq.mean.regional.py +++ /dev/null @@ -1,152 +0,0 @@ -import copy -import datetime -import json -import os -import pickle - -import numpy as np - - -def prdday_to_frq3hridx(prdday, frequency): - frq3hr = 1.0 / (float(prdday) * 8.0) - idx = (np.abs(frequency - frq3hr)).argmin() - return int(idx) - - -def prdday_to_frq1didx(prdday, frequency): - frq24hr = 1.0 / (float(prdday)) - idx = (np.abs(frequency - frq24hr)).argmin() - return int(idx) - - -# -------------------- -# User settings here -# -------------------- -# "3hr" or "day" -hr = "" -# Input file name (pickle .pkl output from calc_ps_area.mean.py) -fname = "" -# Output directory (not including version) -outdir = "" - -# ----------------------- -# The rest of the script -# ----------------------- -ver = datetime.datetime.now().strftime("v%Y%m%d") - -if hr == "day": - frqs_forced = ["semi-annual", "annual"] - frqs_unforced = ["synoptic", "sub-seasonal", "seasonal-annual", "interannual"] -elif hr == "3hr": - frqs_forced = ["semi-diurnal", "diurnal", "semi-annual", "annual"] - frqs_unforced = [ - "sub-daily", - "synoptic", - "sub-seasonal", - "seasonal-annual", - "interannual", - ] - -infile = open(fname, "rb") -psdm = pickle.load(infile) - -psdmfm = copy.deepcopy(psdm) -for frc in psdm.keys(): - if frc == "forced": - frqs = frqs_forced - elif frc == "unforced": - frqs = frqs_unforced - print(frc) - for mip in psdm[frc].keys(): - print(mip) - for dat in psdm[frc][mip].keys(): - frequency = np.array(psdm[frc][mip][dat]["freqs"]) - del psdm[frc][mip][dat]["freqs"] - del psdmfm[frc][mip][dat]["freqs"] - - for var in psdm[frc][mip][dat].keys(): - print(dat, var) - for idm, dom in enumerate(psdm[frc][mip][dat][var].keys()): - am = np.array(psdm[frc][mip][dat][var][dom]) - del psdmfm[frc][mip][dat][var][dom] - psdmfm[frc][mip][dat][var][dom] = {} - print(dom) - for frq in frqs: - print(frq) - if frq == "semi-diurnal": # pr=0.5day - idx = prdday_to_frq1didx(0.5, frequency) - amfm = am[idx] - elif frq == "diurnal": # pr=1day - idx = prdday_to_frq1didx(1, frequency) - amfm = am[idx] - if frq == "semi-annual": # 180day= Date: Thu, 1 Aug 2024 14:15:58 -0700 Subject: [PATCH 50/55] Add res --- docs/metrics_precip-variability.rst | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/metrics_precip-variability.rst b/docs/metrics_precip-variability.rst index 10374c892..7b5b08028 100644 --- a/docs/metrics_precip-variability.rst +++ b/docs/metrics_precip-variability.rst @@ -62,6 +62,7 @@ Options available to set in the parameter file include: * **nperseg**: Length of segment in power spectra. * **noverlap**: Length of overlap between segments in power spectra. * **ref**: Reference data path. +* **res**: Target resolution in degrees. * **regions_specs**: Dictionary containing region bounding box. Uses format {"your region name": {"domain": {"latitude": (min, max), "longitude": (min, max)}}}. Min and max should be replaced with the values that define the region. * **region_file**: Path to a shapefile containing vector region outline. Must use with **attr** and **feature** parameters. * **attr**: Attribute used to identify region in shapefile (eg, column of attribute table). For example, "COUNTRY" in a shapefile of countries. From bacb89145c01405cefda2321c84a4ca9e9e39270 Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Mon, 26 Aug 2024 11:48:01 -0700 Subject: [PATCH 51/55] fix regionspsecs --- .../lib/lib_variability_across_timescales.py | 11 +++++++++-- 1 file changed, 9 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 23a2097cf..76c8619d0 100644 --- a/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py +++ b/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py @@ -446,6 +446,13 @@ def Avg_PS_DomFrq(d, frequency, ntd, dat, mip, frc, regions_specs): Output - psdmfm: Domain and Frequency averaged of spectral power (json) """ + + def val_from_rs(regions_specs, dom): + if regions_specs is not None: + if dom in regions_specs: + return regions_specs[dom].get("value", -1) + return -1 + domains = [ "Total_50S50N", "Ocean_50S50N", @@ -493,9 +500,9 @@ def Avg_PS_DomFrq(d, frequency, ntd, dat, mip, frc, regions_specs): for dom in domains: psdmfm[dom] = {} - if "Ocean" in dom or regions_specs[dom].get("value", -1) == 0: + if "Ocean" in dom or val_from_rs(regions_specs, dom) == 0: dmask = d.where(mask == 0) - elif "Land" in dom or regions_specs[dom].get("value", -1) == 1: + elif "Land" in dom or val_from_rs(regions_specs, dom) == 1: dmask = d.where(mask == 1) else: dmask = d From e094b867c422eb318ca807d31b5c8a70c6e60f78 Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Mon, 26 Aug 2024 14:41:34 -0700 Subject: [PATCH 52/55] bug fixes --- .../lib/lib_variability_across_timescales.py | 22 +++++++++++++++++-- 1 file changed, 20 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 76c8619d0..cf359835a 100644 --- a/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py +++ b/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py @@ -61,11 +61,24 @@ def precip_variability_across_timescale( print(dat, cal) print("syr, eyr:", syr, eyr) + if "360" in cal: + f = f.sel( + time=slice( + str(syr) + "-01-01 00:00:00", + str(eyr) + "-12-" + str(30) + " 23:59:59", + ) + ) + else: + f = f.sel( + time=slice( + str(syr) + "-01-01 00:00:00", + str(eyr) + "-12-" + str(31) + " 23:59:59", + ) + ) + # Regridding xvar = find_lon(f) yvar = find_lat(f) - nlon = str(len(xvar)) - nlat = str(len(yvar)) if np.min(f[xvar]) < 0: lon_range = [-180.0, 180.0] else: @@ -75,6 +88,8 @@ def precip_variability_across_timescale( print("Cropping from shapefile") rgtmp = region_from_file(rgtmp, fshp, attr, feature) drg = rgtmp[var] * float(fac) + nlon = str(len(drg[xvar])) + nlat = str(len(drg[yvar])) f.close() @@ -296,6 +311,8 @@ def ClimAnom(d, ntd, syr, eyr, cal): str(year) + "-12-" + str(ldy) + " 23:59:59", ) ) + print(str(year) + "-01-01 00:00:00") + print(str(year) + "-12-" + str(ldy) + " 23:59:59") if yrtmp.shape[0] == 365 * ntd: anom = np.append( @@ -303,6 +320,7 @@ def ClimAnom(d, ntd, syr, eyr, cal): (np.delete(dseg[iyr], 59, axis=0) - np.delete(clim, 59, axis=0)), ) else: + print("other") anom = np.append(anom, (dseg[iyr] - clim)) else: for iyr, year in enumerate(range(syr, eyr + 1)): From 7b60bb7e8ae0ee8cf00a5958bf041e2251f183d8 Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Wed, 28 Aug 2024 15:13:29 -0700 Subject: [PATCH 53/55] bug fixes --- .../lib/lib_variability_across_timescales.py | 43 ++++++++++++++----- 1 file changed, 32 insertions(+), 11 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 cf359835a..0fd53adab 100644 --- a/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py +++ b/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py @@ -59,7 +59,7 @@ def precip_variability_across_timescale( f = xcdat.decode_time(f) cal = f.time.encoding["calendar"] print(dat, cal) - print("syr, eyr:", syr, eyr) + print(syr, eyr) if "360" in cal: f = f.sel( @@ -99,6 +99,17 @@ def precip_variability_across_timescale( # Power spectum of total freqs, ps, rn, sig95 = Powerspectrum(drg, nperseg, noverlap) # Domain & Frequency average + if fshp and regions_specs is None: + # Set up the regions_specs to cover the whole earth; areas outside + # the shapefile region are NaN, and will not be included. + regions_specs = { + feature: { + "domain": { + "latitude": (-90, 90), + "longitude": (lon_range[0], lon_range[1]), + } + } + } psdmfm_forced = Avg_PS_DomFrq(ps, freqs, ntd, dat, mip, "forced", regions_specs) # Write data (nc file) outfilename = ( @@ -311,8 +322,6 @@ def ClimAnom(d, ntd, syr, eyr, cal): str(year) + "-12-" + str(ldy) + " 23:59:59", ) ) - print(str(year) + "-01-01 00:00:00") - print(str(year) + "-12-" + str(ldy) + " 23:59:59") if yrtmp.shape[0] == 365 * ntd: anom = np.append( @@ -320,7 +329,6 @@ def ClimAnom(d, ntd, syr, eyr, cal): (np.delete(dseg[iyr], 59, axis=0) - np.delete(clim, 59, axis=0)), ) else: - print("other") anom = np.append(anom, (dseg[iyr] - clim)) else: for iyr, year in enumerate(range(syr, eyr + 1)): @@ -465,10 +473,10 @@ def Avg_PS_DomFrq(d, frequency, ntd, dat, mip, frc, regions_specs): - psdmfm: Domain and Frequency averaged of spectral power (json) """ - def val_from_rs(regions_specs, dom): - if regions_specs is not None: - if dom in regions_specs: - return regions_specs[dom].get("value", -1) + def val_from_rs(regions_specs, reg_dom): + if regions_specs: + if reg_dom in regions_specs: + return regions_specs[reg_dom].get("value", -1) return -1 domains = [ @@ -518,9 +526,9 @@ def val_from_rs(regions_specs, dom): for dom in domains: psdmfm[dom] = {} - if "Ocean" in dom or val_from_rs(regions_specs, dom) == 0: + if "Ocean_" in dom or val_from_rs(regions_specs, dom) == 0: dmask = d.where(mask == 0) - elif "Land" in dom or val_from_rs(regions_specs, dom) == 1: + elif "Land_" in dom or val_from_rs(regions_specs, dom) == 1: dmask = d.where(mask == 1) else: dmask = d @@ -551,6 +559,19 @@ def val_from_rs(regions_specs, dom): am = np.array(am) + def check_nan(data): + if isinstance(data, list): + data2 = [] + for item in data: + if np.isnan(item): + data2.append(str(data)) + else: + data2.append(data) + return data2 + if np.isnan(data): + return str(data) + return data + for frq in frqs: if frq == "semi-diurnal": # pr=0.5day idx = prdday_to_frqidx(0.5, frequency, ntd) @@ -584,7 +605,7 @@ def val_from_rs(regions_specs, dom): elif frq == "interannual": # 365day= Date: Wed, 18 Sep 2024 14:49:30 -0700 Subject: [PATCH 54/55] add example --- .../Demo/Demo_7_precip_variability.ipynb | 97 ++++++++++++++++++- 1 file changed, 92 insertions(+), 5 deletions(-) diff --git a/doc/jupyter/Demo/Demo_7_precip_variability.ipynb b/doc/jupyter/Demo/Demo_7_precip_variability.ipynb index 4e5652d20..ff57033b8 100644 --- a/doc/jupyter/Demo/Demo_7_precip_variability.ipynb +++ b/doc/jupyter/Demo/Demo_7_precip_variability.ipynb @@ -586,7 +586,9 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "tags": [] + }, "source": [ "## Precipitation Variability Metric" ] @@ -785,19 +787,104 @@ "print(json.dumps(metric, indent=2))" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Regional metrics" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The precipitation variability metrics have a set of default regions, but users can instead define a single spatial region to compute metrics over. There are two ways to do this.\n", + "\n", + "1. Use the `regions_specs` parameter to define a latitude/longitude box. \n", + "Parameter file example:\n", + "```\n", + "regions_specs={\"CONUS\": {\"domain\": {\"latitude\": (24.7, 49.4), \"longitude\": (235.22, 293.08)}}}\n", + "```\n", + "\n", + "2. Use a shapefile to define a region. Users must provide the path to the shapefile along with the attribute/feature pair that defines the region. \n", + "Parameter file example:\n", + "```\n", + "region_file=\"CONUS.shp\" # Shapefile path\n", + "attr=\"NAME\" # An attribute in the shapefile\n", + "feature=\"CONUS\" # A unique feature name that can be \n", + " # found under the \"attr\" attribute\n", + "```" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "First, we generate a simple shapefile for use in this example. The shapefile contains one feature, a box that defines the CONUS region." + ] + }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], - "source": [] + "source": [ + "from shapely import Polygon\n", + "import geo pandas as gpd\n", + "import pandas as pd\n", + "\n", + "coords = ((233.,22.),(233.,50.),(294.,50.),(294.,22))\n", + "df = pd.DataFrame({\"Region\": [\"CONUS\"], \"Coords\": [Polygon(coords)]})\n", + "gdf = gpd.GeoDataFrame(df, geometry=\"Coords\", crs=\"EPSG:4326\")\n", + "gdf.to_file(demo_output_directory+'/shp/CONUS.shp')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Add the information for this shapefile to the variability_across_timescales_PS_driver.py run command." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%%bash -s \"$demo_output_directory\"\n", + "variability_across_timescales_PS_driver.py -p basic_precip_variability_param.py \\\n", + "--region_file $1/shp/CONUS.shp \\\n", + "--attr 'Region' \\\n", + "--feature 'CONUS' \\\n", + "--results_dir $1/precip_variability/region_ex" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The metrics output will look different than the default example. Metrics will only be produced for a single region that we defined in this shapefile." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "output_path = os.path.join(demo_output_directory,\"precip_variability/region_ex/PS_pr.day_regrid.180x90_area.freq.mean_GISS-E2-H.r6i1p1.json\")\n", + "with open(output_path) as f:\n", + " metric = json.load(f)[\"RESULTS\"]\n", + "print(json.dumps(metric, indent=2))" + ] } ], "metadata": { "kernelspec": { - "display_name": "Python [conda env:pmp_v20220110] *", + "display_name": "pmp_pv_dev", "language": "python", - "name": "conda-env-pmp_v20220110-py" + "name": "pmp_pv_dev" }, "language_info": { "codemirror_mode": { @@ -809,7 +896,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.9.10" + "version": "3.10.12" } }, "nbformat": 4, From d56d3734a7957f55176a87b3ad1577e976ba9bc7 Mon Sep 17 00:00:00 2001 From: Ana Ordonez Date: Wed, 18 Sep 2024 16:12:16 -0700 Subject: [PATCH 55/55] update region ex --- .../Demo/Demo_7_precip_variability.ipynb | 518 +++++++++--------- 1 file changed, 268 insertions(+), 250 deletions(-) diff --git a/doc/jupyter/Demo/Demo_7_precip_variability.ipynb b/doc/jupyter/Demo/Demo_7_precip_variability.ipynb index ff57033b8..2996603ef 100644 --- a/doc/jupyter/Demo/Demo_7_precip_variability.ipynb +++ b/doc/jupyter/Demo/Demo_7_precip_variability.ipynb @@ -66,10 +66,14 @@ " [--fac FAC]\n", " [--nperseg NPERSEG]\n", " [--noverlap NOVERLAP]\n", - " [--ref REF] [--cmec]\n", - " [--no_cmec]\n", + " [--ref REF] [--res RES]\n", + " [--regions_specs REGIONS_SPECS]\n", + " [--cmec] [--no_cmec]\n", + " [--region_file REGION_FILE]\n", + " [--feature FEATURE]\n", + " [--attr ATTR]\n", "\n", - "optional arguments:\n", + "options:\n", " -h, --help show this help message and exit\n", " --parameters PARAMETERS, -p PARAMETERS\n", " --diags OTHER_PARAMETERS [OTHER_PARAMETERS ...], -d OTHER_PARAMETERS [OTHER_PARAMETERS ...]\n", @@ -91,26 +95,25 @@ " --noverlap NOVERLAP length of overlap between segments in power spectra\n", " (default: None)\n", " --ref REF reference data path (default: None)\n", + " --res RES Horizontal resolution [degree] for interpolation (lon,\n", + " lat) (default: 2)\n", + " --regions_specs REGIONS_SPECS\n", + " Provide a single custom region (default: None)\n", " --cmec Use to save CMEC format metrics JSON (default: False)\n", - " --no_cmec Do not save CMEC format metrics JSON (default: False)\n" + " --no_cmec Do not save CMEC format metrics JSON (default: False)\n", + " --region_file REGION_FILE\n", + " File containing vector region of interest. (default:\n", + " None)\n", + " --feature FEATURE Feature in vectorized region. (default: None)\n", + " --attr ATTR Attribute containing feature in vectorized region.\n", + " (default: None)\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ - "/home/ahn6/anaconda3/envs/pmp_v20220110/lib/python3.9/site-packages/pkg_resources/__init__.py:116: PkgResourcesDeprecationWarning: v3.0-46-ge6f562d is an invalid version and will not be supported in a future release\n", - " warnings.warn(\n", - "/home/ahn6/anaconda3/envs/pmp_v20220110/lib/python3.9/site-packages/pkg_resources/__init__.py:116: PkgResourcesDeprecationWarning: v1.2.1-1163-g66b2186 is an invalid version and will not be supported in a future release\n", - " warnings.warn(\n", - "/home/ahn6/anaconda3/envs/pmp_v20220110/lib/python3.9/site-packages/pkg_resources/__init__.py:116: PkgResourcesDeprecationWarning: v1.2.1-1165-g43e256e is an invalid version and will not be supported in a future release\n", - " warnings.warn(\n", - "/home/ahn6/anaconda3/envs/pmp_v20220110/lib/python3.9/site-packages/pkg_resources/__init__.py:116: PkgResourcesDeprecationWarning: v1.2.1-640-gc4401ed is an invalid version and will not be supported in a future release\n", - " warnings.warn(\n", - "/home/ahn6/anaconda3/envs/pmp_v20220110/lib/python3.9/site-packages/pkg_resources/__init__.py:116: PkgResourcesDeprecationWarning: v2.4.0-23-g4fa5682 is an invalid version and will not be supported in a future release\n", - " warnings.warn(\n", - "/home/ahn6/anaconda3/envs/pmp_v20220110/lib/python3.9/site-packages/pkg_resources/__init__.py:116: PkgResourcesDeprecationWarning: v2.5.1-10-g93c30ce is an invalid version and will not be supported in a future release\n", - " warnings.warn(\n" + "[WARNING] yaksa: 10 leaked handle pool objects\n" ] } ], @@ -145,8 +148,8 @@ "mod = \"pr_day_GISS-E2-H_historical_r6i1p1_*.nc\"\n", "var = \"pr\"\n", "frq = \"day\"\n", - "modpath = 'demo_data/CMIP5_demo_timeseries/historical/atmos/day/pr/'\n", - "results_dir = 'demo_output/precip_variability/GISS-E2-H/'\n", + "modpath = 'demo_data_tmp/CMIP5_demo_timeseries/historical/atmos/day/pr/'\n", + "results_dir = 'demo_output_tmp/precip_variability/GISS-E2-H/'\n", "prd = [2000,2005] # analysis period\n", "fac = 86400 # factor to make unit of [mm/day]\n", "\n", @@ -187,39 +190,33 @@ "scrolled": true }, "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "INFO::2024-09-18 15:56::pcmdi_metrics:: Results saved to a json file: /home/ordonez4/git/pcmdi_metrics/doc/jupyter/Demo/demo_output_tmp/precip_variability/GISS-E2-H/PS_pr.day_regrid.180x90_area.freq.mean_GISS-E2-H.r6i1p1.json\n", + "2024-09-18 15:56:49,164 [INFO]: base.py(write:422) >> Results saved to a json file: /home/ordonez4/git/pcmdi_metrics/doc/jupyter/Demo/demo_output_tmp/precip_variability/GISS-E2-H/PS_pr.day_regrid.180x90_area.freq.mean_GISS-E2-H.r6i1p1.json\n", + "2024-09-18 15:56:49,164 [INFO]: base.py(write:422) >> Results saved to a json file: /home/ordonez4/git/pcmdi_metrics/doc/jupyter/Demo/demo_output_tmp/precip_variability/GISS-E2-H/PS_pr.day_regrid.180x90_area.freq.mean_GISS-E2-H.r6i1p1.json\n" + ] + }, { "name": "stdout", "output_type": "stream", "text": [ - "demo_data/CMIP5_demo_timeseries/historical/atmos/day/pr/\n", + "demo_data_tmp/CMIP5_demo_timeseries/historical/atmos/day/pr/\n", "pr_day_GISS-E2-H_historical_r6i1p1_*.nc\n", "[2000, 2005]\n", "730 365\n", - "demo_output/precip_variability/GISS-E2-H/\n", - "demo_output/precip_variability/GISS-E2-H/\n", - "demo_output/precip_variability/GISS-E2-H/\n", + "2\n", + "demo_output_tmp/precip_variability/GISS-E2-H/\n", + "demo_output_tmp/precip_variability/GISS-E2-H/\n", + "demo_output_tmp/precip_variability/GISS-E2-H/\n", + "['demo_data_tmp/CMIP5_demo_timeseries/historical/atmos/day/pr/pr_day_GISS-E2-H_historical_r6i1p1_20000101-20051231.nc']\n", "GISS-E2-H.r6i1p1\n", - "['demo_data/CMIP5_demo_timeseries/historical/atmos/day/pr/pr_day_GISS-E2-H_historical_r6i1p1_20000101-20051231.nc']\n", + "['demo_data_tmp/CMIP5_demo_timeseries/historical/atmos/day/pr/pr_day_GISS-E2-H_historical_r6i1p1_20000101-20051231.nc']\n", "GISS-E2-H.r6i1p1 365_day\n", - "syr, eyr: 2000 2005\n", - "2000\n", - "Complete regridding from (365, 90, 144) to (365, 90, 180)\n", - "2000 (365, 90, 180)\n", - "2001\n", - "Complete regridding from (365, 90, 144) to (365, 90, 180)\n", - "2001 (730, 90, 180)\n", - "2002\n", - "Complete regridding from (365, 90, 144) to (365, 90, 180)\n", - "2002 (1095, 90, 180)\n", - "2003\n", - "Complete regridding from (365, 90, 144) to (365, 90, 180)\n", - "2003 (1460, 90, 180)\n", - "2004\n", - "Complete regridding from (365, 90, 144) to (365, 90, 180)\n", - "2004 (1825, 90, 180)\n", - "2005\n", - "Complete regridding from (365, 90, 144) to (365, 90, 180)\n", - "2005 (2190, 90, 180)\n", + "2000 2005\n", + "Complete regridding from (2190, 90, 144) to (2190, 90, 180)\n", "Complete calculating climatology and anomaly for calendar of 365_day\n", "Complete power spectra (segment: 730 nps: 5.0 )\n", "Complete domain and frequency average of spectral power\n", @@ -231,20 +228,7 @@ "name": "stderr", "output_type": "stream", "text": [ - "/home/ahn6/anaconda3/envs/pmp_v20220110/lib/python3.9/site-packages/pkg_resources/__init__.py:116: PkgResourcesDeprecationWarning: v3.0-46-ge6f562d is an invalid version and will not be supported in a future release\n", - " warnings.warn(\n", - "/home/ahn6/anaconda3/envs/pmp_v20220110/lib/python3.9/site-packages/pkg_resources/__init__.py:116: PkgResourcesDeprecationWarning: v1.2.1-1163-g66b2186 is an invalid version and will not be supported in a future release\n", - " warnings.warn(\n", - "/home/ahn6/anaconda3/envs/pmp_v20220110/lib/python3.9/site-packages/pkg_resources/__init__.py:116: PkgResourcesDeprecationWarning: v1.2.1-1165-g43e256e is an invalid version and will not be supported in a future release\n", - " warnings.warn(\n", - "/home/ahn6/anaconda3/envs/pmp_v20220110/lib/python3.9/site-packages/pkg_resources/__init__.py:116: PkgResourcesDeprecationWarning: v1.2.1-640-gc4401ed is an invalid version and will not be supported in a future release\n", - " warnings.warn(\n", - "/home/ahn6/anaconda3/envs/pmp_v20220110/lib/python3.9/site-packages/pkg_resources/__init__.py:116: PkgResourcesDeprecationWarning: v2.4.0-23-g4fa5682 is an invalid version and will not be supported in a future release\n", - " warnings.warn(\n", - "/home/ahn6/anaconda3/envs/pmp_v20220110/lib/python3.9/site-packages/pkg_resources/__init__.py:116: PkgResourcesDeprecationWarning: v2.5.1-10-g93c30ce is an invalid version and will not be supported in a future release\n", - " warnings.warn(\n", - "INFO::2023-02-23 11:51::pcmdi_metrics:: Results saved to a json file: /home/ahn6/PCMDI/pcmdi_metrics/branch/900_msa_precip_variability_demo/pcmdi_metrics/doc/jupyter/Demo/demo_output/precip_variability/GISS-E2-H/PS_pr.day_regrid.180x90_area.freq.mean_GISS-E2-H.r6i1p1.json\n", - "2023-02-23 11:51:11,427 [INFO]: base.py(write:237) >> Results saved to a json file: /home/ahn6/PCMDI/pcmdi_metrics/branch/900_msa_precip_variability_demo/pcmdi_metrics/doc/jupyter/Demo/demo_output/precip_variability/GISS-E2-H/PS_pr.day_regrid.180x90_area.freq.mean_GISS-E2-H.r6i1p1.json\n" + "[WARNING] yaksa: 10 leaked handle pool objects\n" ] } ], @@ -274,9 +258,9 @@ "name": "stdout", "output_type": "stream", "text": [ - "PS_pr.day_regrid.180x90_area.freq.mean_GISS-E2-H.r6i1p1.json\n", - "PS_pr.day_regrid.180x90_GISS-E2-H.r6i1p1.nc\n", - "PS_pr.day_regrid.180x90_GISS-E2-H.r6i1p1_unforced.nc\n" + "PS_pr.day_regrid.180x90_area.freq.mean_GISS-E2-H.r6i1p1.json\r\n", + "PS_pr.day_regrid.180x90_GISS-E2-H.r6i1p1.nc\r\n", + "PS_pr.day_regrid.180x90_GISS-E2-H.r6i1p1_unforced.nc\r\n" ] } ], @@ -304,36 +288,36 @@ " \"GISS-E2-H.r6i1p1\": {\n", " \"forced\": {\n", " \"Land_30N50N\": {\n", - " \"annual\": 1.154266721510137,\n", - " \"semi-annual\": 0.3692551744241903\n", + " \"annual\": 1.153948602189096,\n", + " \"semi-annual\": 0.3675381067314767\n", " },\n", " \"Land_30S30N\": {\n", - " \"annual\": 6.8655960795131294,\n", - " \"semi-annual\": 1.1969126049181855\n", + " \"annual\": 6.8509958100746555,\n", + " \"semi-annual\": 1.1945015595812805\n", " },\n", " \"Land_50S30S\": {\n", - " \"annual\": 0.7829928891110198,\n", - " \"semi-annual\": 0.33398811326967975\n", + " \"annual\": 0.8090939740005696,\n", + " \"semi-annual\": 0.34297346148163804\n", " },\n", " \"Land_50S50N\": {\n", - " \"annual\": 4.803117924524398,\n", - " \"semi-annual\": 0.8989181591887316\n", + " \"annual\": 4.793570167683052,\n", + " \"semi-annual\": 0.8971106124805638\n", " },\n", " \"Ocean_30N50N\": {\n", - " \"annual\": 1.4467988289024762,\n", - " \"semi-annual\": 0.37232162338162866\n", + " \"annual\": 1.450126151318265,\n", + " \"semi-annual\": 0.3738726067518909\n", " },\n", " \"Ocean_30S30N\": {\n", - " \"annual\": 4.568654517465613,\n", - " \"semi-annual\": 1.5044899979603008\n", + " \"annual\": 4.561426422605001,\n", + " \"semi-annual\": 1.5069884231014545\n", " },\n", " \"Ocean_50S30S\": {\n", - " \"annual\": 0.5918242629787758,\n", - " \"semi-annual\": 0.1927211439124904\n", + " \"annual\": 0.5890515819402276,\n", + " \"semi-annual\": 0.19150748548003316\n", " },\n", " \"Ocean_50S50N\": {\n", - " \"annual\": 3.3099973296409195,\n", - " \"semi-annual\": 1.0764366904440072\n", + " \"annual\": 3.3050864193776026,\n", + " \"semi-annual\": 1.0780758057454556\n", " },\n", " \"Total_30N50N\": {\n", " \"annual\": 1.3110986682307972,\n", @@ -354,52 +338,52 @@ " },\n", " \"unforced\": {\n", " \"Land_30N50N\": {\n", - " \"interannual\": 0.11027519561502575,\n", - " \"seasonal-annual\": 0.15034540027412407,\n", - " \"sub-seasonal\": 0.13605316533174094,\n", - " \"synoptic\": 0.06277267344693233\n", + " \"interannual\": 0.11025112312631524,\n", + " \"seasonal-annual\": 0.1502570664470804,\n", + " \"sub-seasonal\": 0.13618888930844514,\n", + " \"synoptic\": 0.06327297649960462\n", " },\n", " \"Land_30S30N\": {\n", - " \"interannual\": 0.31262787817502063,\n", - " \"seasonal-annual\": 0.30924523792899317,\n", - " \"sub-seasonal\": 0.243897137677461,\n", - " \"synoptic\": 0.07552274806731542\n", + " \"interannual\": 0.31535297942347246,\n", + " \"seasonal-annual\": 0.31179854291318526,\n", + " \"sub-seasonal\": 0.2477967897126997,\n", + " \"synoptic\": 0.076484979080103\n", " },\n", " \"Land_50S30S\": {\n", - " \"interannual\": 0.1527755501819138,\n", - " \"seasonal-annual\": 0.2041066189679213,\n", - " \"sub-seasonal\": 0.17203311677473293,\n", - " \"synoptic\": 0.07133658627473073\n", + " \"interannual\": 0.16178541870984994,\n", + " \"seasonal-annual\": 0.21589364787265297,\n", + " \"sub-seasonal\": 0.1847557860658534,\n", + " \"synoptic\": 0.07524240453524904\n", " },\n", " \"Land_50S50N\": {\n", - " \"interannual\": 0.24222527413177844,\n", - " \"seasonal-annual\": 0.25493923345176356,\n", - " \"sub-seasonal\": 0.20701576396034696,\n", - " \"synoptic\": 0.07136923241812429\n", + " \"interannual\": 0.24443780233759468,\n", + " \"seasonal-annual\": 0.25718039033896883,\n", + " \"sub-seasonal\": 0.2102202999468355,\n", + " \"synoptic\": 0.07234360585017287\n", " },\n", " \"Ocean_30N50N\": {\n", - " \"interannual\": 0.13240454583620145,\n", - " \"seasonal-annual\": 0.17549307553001822,\n", - " \"sub-seasonal\": 0.15428702961801613,\n", - " \"synoptic\": 0.09824822890957895\n", + " \"interannual\": 0.13265625643216272,\n", + " \"seasonal-annual\": 0.1758330640897642,\n", + " \"sub-seasonal\": 0.15435681112427357,\n", + " \"synoptic\": 0.0981749977902816\n", " },\n", " \"Ocean_30S30N\": {\n", - " \"interannual\": 0.6531055721473543,\n", - " \"seasonal-annual\": 0.6376662281993245,\n", - " \"sub-seasonal\": 0.43458496427824367,\n", - " \"synoptic\": 0.11441802205292609\n", + " \"interannual\": 0.6539803811119562,\n", + " \"seasonal-annual\": 0.6385364543707663,\n", + " \"sub-seasonal\": 0.43424291163472306,\n", + " \"synoptic\": 0.11428977945404156\n", " },\n", " \"Ocean_50S30S\": {\n", - " \"interannual\": 0.09839333071364982,\n", - " \"seasonal-annual\": 0.13364245158376373,\n", - " \"sub-seasonal\": 0.12036603745194574,\n", - " \"synoptic\": 0.06906461136148313\n", + " \"interannual\": 0.09747609150424397,\n", + " \"seasonal-annual\": 0.13244482423836793,\n", + " \"sub-seasonal\": 0.11915711328928573,\n", + " \"synoptic\": 0.06874014945078849\n", " },\n", " \"Ocean_50S50N\": {\n", - " \"interannual\": 0.46681841351291453,\n", - " \"seasonal-annual\": 0.4697785382523058,\n", - " \"sub-seasonal\": 0.3309051790489744,\n", - " \"synoptic\": 0.10250804393167526\n", + " \"interannual\": 0.467278699215871,\n", + " \"seasonal-annual\": 0.4701741107777076,\n", + " \"sub-seasonal\": 0.33044759093021675,\n", + " \"synoptic\": 0.10233245216785025\n", " },\n", " \"Total_30N50N\": {\n", " \"interannual\": 0.12213915511604374,\n", @@ -471,81 +455,33 @@ "execution_count": 8, "metadata": {}, "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "INFO::2024-09-18 16:08::pcmdi_metrics:: Results saved to a json file: /home/ordonez4/git/pcmdi_metrics/doc/jupyter/Demo/demo_output_tmp/precip_variability/GPCP-1-3/PS_pr.day_regrid.180x90_area.freq.mean_GPCP-1-3.json\n", + "2024-09-18 16:08:28,962 [INFO]: base.py(write:422) >> Results saved to a json file: /home/ordonez4/git/pcmdi_metrics/doc/jupyter/Demo/demo_output_tmp/precip_variability/GPCP-1-3/PS_pr.day_regrid.180x90_area.freq.mean_GPCP-1-3.json\n", + "2024-09-18 16:08:28,962 [INFO]: base.py(write:422) >> Results saved to a json file: /home/ordonez4/git/pcmdi_metrics/doc/jupyter/Demo/demo_output_tmp/precip_variability/GPCP-1-3/PS_pr.day_regrid.180x90_area.freq.mean_GPCP-1-3.json\n" + ] + }, { "name": "stdout", "output_type": "stream", "text": [ - "demo_data/obs4MIPs_PCMDI_daily/NASA-JPL/GPCP-1-3/day/pr/gn/latest/\n", + "demo_data_tmp/obs4MIPs_PCMDI_daily/NASA-JPL/GPCP-1-3/day/pr/gn/latest/\n", "pr_day_GPCP-1-3_PCMDI_gn_19961002-20170101.nc\n", "[1997, 2016]\n", "730 365\n", - "demo_output/precip_variability/GPCP-1-3/\n", - "demo_output/precip_variability/GPCP-1-3/\n", - "demo_output/precip_variability/GPCP-1-3/\n", + "2\n", + "demo_output_tmp/precip_variability/GPCP-1-3/\n", + "demo_output_tmp/precip_variability/GPCP-1-3/\n", + "demo_output_tmp/precip_variability/GPCP-1-3/\n", + "['demo_data_tmp/obs4MIPs_PCMDI_daily/NASA-JPL/GPCP-1-3/day/pr/gn/latest/pr_day_GPCP-1-3_PCMDI_gn_19961002-20170101.nc']\n", "GPCP-1-3\n", - "['demo_data/obs4MIPs_PCMDI_daily/NASA-JPL/GPCP-1-3/day/pr/gn/latest/pr_day_GPCP-1-3_PCMDI_gn_19961002-20170101.nc']\n", + "['demo_data_tmp/obs4MIPs_PCMDI_daily/NASA-JPL/GPCP-1-3/day/pr/gn/latest/pr_day_GPCP-1-3_PCMDI_gn_19961002-20170101.nc']\n", "GPCP-1-3 gregorian\n", - "syr, eyr: 1997 2016\n", - "1997\n", - "Complete regridding from (365, 180, 360) to (365, 90, 180)\n", - "1997 (365, 90, 180)\n", - "1998\n", - "Complete regridding from (365, 180, 360) to (365, 90, 180)\n", - "1998 (730, 90, 180)\n", - "1999\n", - "Complete regridding from (365, 180, 360) to (365, 90, 180)\n", - "1999 (1095, 90, 180)\n", - "2000\n", - "Complete regridding from (366, 180, 360) to (366, 90, 180)\n", - "2000 (1461, 90, 180)\n", - "2001\n", - "Complete regridding from (365, 180, 360) to (365, 90, 180)\n", - "2001 (1826, 90, 180)\n", - "2002\n", - "Complete regridding from (365, 180, 360) to (365, 90, 180)\n", - "2002 (2191, 90, 180)\n", - "2003\n", - "Complete regridding from (365, 180, 360) to (365, 90, 180)\n", - "2003 (2556, 90, 180)\n", - "2004\n", - "Complete regridding from (366, 180, 360) to (366, 90, 180)\n", - "2004 (2922, 90, 180)\n", - "2005\n", - "Complete regridding from (365, 180, 360) to (365, 90, 180)\n", - "2005 (3287, 90, 180)\n", - "2006\n", - "Complete regridding from (365, 180, 360) to (365, 90, 180)\n", - "2006 (3652, 90, 180)\n", - "2007\n", - "Complete regridding from (365, 180, 360) to (365, 90, 180)\n", - "2007 (4017, 90, 180)\n", - "2008\n", - "Complete regridding from (366, 180, 360) to (366, 90, 180)\n", - "2008 (4383, 90, 180)\n", - "2009\n", - "Complete regridding from (365, 180, 360) to (365, 90, 180)\n", - "2009 (4748, 90, 180)\n", - "2010\n", - "Complete regridding from (365, 180, 360) to (365, 90, 180)\n", - "2010 (5113, 90, 180)\n", - "2011\n", - "Complete regridding from (365, 180, 360) to (365, 90, 180)\n", - "2011 (5478, 90, 180)\n", - "2012\n", - "Complete regridding from (366, 180, 360) to (366, 90, 180)\n", - "2012 (5844, 90, 180)\n", - "2013\n", - "Complete regridding from (365, 180, 360) to (365, 90, 180)\n", - "2013 (6209, 90, 180)\n", - "2014\n", - "Complete regridding from (365, 180, 360) to (365, 90, 180)\n", - "2014 (6574, 90, 180)\n", - "2015\n", - "Complete regridding from (365, 180, 360) to (365, 90, 180)\n", - "2015 (6939, 90, 180)\n", - "2016\n", - "Complete regridding from (366, 180, 360) to (366, 90, 180)\n", - "2016 (7305, 90, 180)\n", + "1997 2016\n", + "Complete regridding from (7305, 180, 360) to (7305, 90, 180)\n", "Complete calculating climatology and anomaly for calendar of gregorian\n", "Complete power spectra (segment: 730 nps: 19.0 )\n", "Complete domain and frequency average of spectral power\n", @@ -557,20 +493,7 @@ "name": "stderr", "output_type": "stream", "text": [ - "/home/ahn6/anaconda3/envs/pmp_v20220110/lib/python3.9/site-packages/pkg_resources/__init__.py:116: PkgResourcesDeprecationWarning: v3.0-46-ge6f562d is an invalid version and will not be supported in a future release\n", - " warnings.warn(\n", - "/home/ahn6/anaconda3/envs/pmp_v20220110/lib/python3.9/site-packages/pkg_resources/__init__.py:116: PkgResourcesDeprecationWarning: v1.2.1-1163-g66b2186 is an invalid version and will not be supported in a future release\n", - " warnings.warn(\n", - "/home/ahn6/anaconda3/envs/pmp_v20220110/lib/python3.9/site-packages/pkg_resources/__init__.py:116: PkgResourcesDeprecationWarning: v1.2.1-1165-g43e256e is an invalid version and will not be supported in a future release\n", - " warnings.warn(\n", - "/home/ahn6/anaconda3/envs/pmp_v20220110/lib/python3.9/site-packages/pkg_resources/__init__.py:116: PkgResourcesDeprecationWarning: v1.2.1-640-gc4401ed is an invalid version and will not be supported in a future release\n", - " warnings.warn(\n", - "/home/ahn6/anaconda3/envs/pmp_v20220110/lib/python3.9/site-packages/pkg_resources/__init__.py:116: PkgResourcesDeprecationWarning: v2.4.0-23-g4fa5682 is an invalid version and will not be supported in a future release\n", - " warnings.warn(\n", - "/home/ahn6/anaconda3/envs/pmp_v20220110/lib/python3.9/site-packages/pkg_resources/__init__.py:116: PkgResourcesDeprecationWarning: v2.5.1-10-g93c30ce is an invalid version and will not be supported in a future release\n", - " warnings.warn(\n", - "INFO::2023-02-23 11:58::pcmdi_metrics:: Results saved to a json file: /home/ahn6/PCMDI/pcmdi_metrics/branch/900_msa_precip_variability_demo/pcmdi_metrics/doc/jupyter/Demo/demo_output/precip_variability/GPCP-1-3/PS_pr.day_regrid.180x90_area.freq.mean_GPCP-1-3.json\n", - "2023-02-23 11:58:48,155 [INFO]: base.py(write:237) >> Results saved to a json file: /home/ahn6/PCMDI/pcmdi_metrics/branch/900_msa_precip_variability_demo/pcmdi_metrics/doc/jupyter/Demo/demo_output/precip_variability/GPCP-1-3/PS_pr.day_regrid.180x90_area.freq.mean_GPCP-1-3.json\n" + "[WARNING] yaksa: 10 leaked handle pool objects\n" ] } ], @@ -616,13 +539,20 @@ "name": "stdout", "output_type": "stream", "text": [ - "reference: demo_output/precip_variability/GPCP-1-3/PS_pr.day_regrid.180x90_area.freq.mean_GPCP-1-3.json\n", - "modpath: demo_output/precip_variability/GISS-E2-H/\n", - "outdir: demo_output/precip_variability/ratio/\n", - "['demo_output/precip_variability/GISS-E2-H/PS_pr.day_regrid.180x90_area.freq.mean_GISS-E2-H.r6i1p1.json']\n", + "reference: demo_output_tmp/precip_variability/GPCP-1-3/PS_pr.day_regrid.180x90_area.freq.mean_GPCP-1-3.json\n", + "modpath: demo_output_tmp/precip_variability/GISS-E2-H/\n", + "outdir: demo_output_tmp/precip_variability/ratio/\n", + "['demo_output_tmp/precip_variability/GISS-E2-H/PS_pr.day_regrid.180x90_area.freq.mean_GISS-E2-H.r6i1p1.json']\n", "Complete GISS-E2-H.r6i1p1\n", "Complete all\n" ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "[WARNING] yaksa: 10 leaked handle pool objects\n" + ] } ], "source": [ @@ -653,36 +583,36 @@ " \"GISS-E2-H.r6i1p1\": {\n", " \"forced\": {\n", " \"Land_30N50N\": {\n", - " \"annual\": 1.6223030227997506,\n", - " \"semi-annual\": 1.873331686130227\n", + " \"annual\": 1.6279984575673894,\n", + " \"semi-annual\": 1.867057373601494\n", " },\n", " \"Land_30S30N\": {\n", - " \"annual\": 1.3409912459551663,\n", - " \"semi-annual\": 1.3385919476532728\n", + " \"annual\": 1.3338114720532706,\n", + " \"semi-annual\": 1.333350181560781\n", " },\n", " \"Land_50S30S\": {\n", - " \"annual\": 1.1582631259388922,\n", - " \"semi-annual\": 1.903328778893799\n", + " \"annual\": 1.164227264547647,\n", + " \"semi-annual\": 1.9246852085563568\n", " },\n", " \"Land_50S50N\": {\n", - " \"annual\": 1.3568447816315299,\n", - " \"semi-annual\": 1.3967541356262723\n", + " \"annual\": 1.3503132388688357,\n", + " \"semi-annual\": 1.391749566706111\n", " },\n", " \"Ocean_30N50N\": {\n", - " \"annual\": 1.0571429112202069,\n", - " \"semi-annual\": 0.8535214354376027\n", + " \"annual\": 1.0524861972773814,\n", + " \"semi-annual\": 0.8517712548298377\n", " },\n", " \"Ocean_30S30N\": {\n", - " \"annual\": 1.4932022320513534,\n", - " \"semi-annual\": 1.817141396507603\n", + " \"annual\": 1.499118828822202,\n", + " \"semi-annual\": 1.8222593026548162\n", " },\n", " \"Ocean_50S30S\": {\n", - " \"annual\": 1.4346150163209932,\n", - " \"semi-annual\": 1.053929465464535\n", + " \"annual\": 1.4363958284724372,\n", + " \"semi-annual\": 1.0484119422307991\n", " },\n", " \"Ocean_50S50N\": {\n", - " \"annual\": 1.4578241823817903,\n", - " \"semi-annual\": 1.6866782169880241\n", + " \"annual\": 1.4625476582104198,\n", + " \"semi-annual\": 1.6902905191733497\n", " },\n", " \"Total_30N50N\": {\n", " \"annual\": 1.2324909366302752,\n", @@ -703,52 +633,52 @@ " },\n", " \"unforced\": {\n", " \"Land_30N50N\": {\n", - " \"interannual\": 1.3925193177807649,\n", - " \"seasonal-annual\": 1.4578033501404122,\n", - " \"sub-seasonal\": 1.2833073189255852,\n", - " \"synoptic\": 0.9599889206339726\n", + " \"interannual\": 1.3879961062215058,\n", + " \"seasonal-annual\": 1.4543733420466998,\n", + " \"sub-seasonal\": 1.2722446114532955,\n", + " \"synoptic\": 0.9550314725762122\n", " },\n", " \"Land_30S30N\": {\n", - " \"interannual\": 1.5881700737299942,\n", - " \"seasonal-annual\": 1.3863862394635518,\n", - " \"sub-seasonal\": 1.0245268449368474,\n", - " \"synoptic\": 0.6287695564515763\n", + " \"interannual\": 1.568479703495435,\n", + " \"seasonal-annual\": 1.3855140760562272,\n", + " \"sub-seasonal\": 1.0320215218679585,\n", + " \"synoptic\": 0.6344408069821\n", " },\n", " \"Land_50S30S\": {\n", - " \"interannual\": 1.2117987802714973,\n", - " \"seasonal-annual\": 1.439080108544409,\n", - " \"sub-seasonal\": 1.0621359559625092,\n", - " \"synoptic\": 0.6461197014112282\n", + " \"interannual\": 1.273480429665751,\n", + " \"seasonal-annual\": 1.4835739537712782,\n", + " \"sub-seasonal\": 1.1166071488025653,\n", + " \"synoptic\": 0.6682326701057775\n", " },\n", " \"Land_50S50N\": {\n", - " \"interannual\": 1.543064523906278,\n", - " \"seasonal-annual\": 1.400904149842401,\n", - " \"sub-seasonal\": 1.0699981481996261,\n", - " \"synoptic\": 0.6950528347881917\n", + " \"interannual\": 1.5292151952175095,\n", + " \"seasonal-annual\": 1.4013209418868053,\n", + " \"sub-seasonal\": 1.076210914944252,\n", + " \"synoptic\": 0.6996958985943696\n", " },\n", " \"Ocean_30N50N\": {\n", - " \"interannual\": 0.7064700849633351,\n", - " \"seasonal-annual\": 0.6481146462601631,\n", - " \"sub-seasonal\": 0.6149584603478464,\n", - " \"synoptic\": 0.6871119360750948\n", + " \"interannual\": 0.7043783826080335,\n", + " \"seasonal-annual\": 0.6455934192259553,\n", + " \"sub-seasonal\": 0.6137724411737419,\n", + " \"synoptic\": 0.6863874501625437\n", " },\n", " \"Ocean_30S30N\": {\n", - " \"interannual\": 1.249329338868364,\n", - " \"seasonal-annual\": 1.5507541726899665,\n", - " \"sub-seasonal\": 1.1807618126643655,\n", - " \"synoptic\": 1.0941063854826194\n", + " \"interannual\": 1.250341415643576,\n", + " \"seasonal-annual\": 1.5516779450827425,\n", + " \"sub-seasonal\": 1.1798960241814673,\n", + " \"synoptic\": 1.0953812575228667\n", " },\n", " \"Ocean_50S30S\": {\n", - " \"interannual\": 0.8611006199438214,\n", - " \"seasonal-annual\": 0.8486520728964435,\n", - " \"sub-seasonal\": 0.7681080333270882,\n", - " \"synoptic\": 0.6800227284451652\n", + " \"interannual\": 0.8539632674914027,\n", + " \"seasonal-annual\": 0.8423603608480983,\n", + " \"sub-seasonal\": 0.7618579649944118,\n", + " \"synoptic\": 0.6782173179198747\n", " },\n", " \"Ocean_50S50N\": {\n", - " \"interannual\": 1.1919874028769966,\n", - " \"seasonal-annual\": 1.3887066381499753,\n", - " \"sub-seasonal\": 1.0768347202958415,\n", - " \"synoptic\": 0.9425899855273115\n", + " \"interannual\": 1.192306567530281,\n", + " \"seasonal-annual\": 1.388569073500126,\n", + " \"sub-seasonal\": 1.075457264771947,\n", + " \"synoptic\": 0.9428927883322024\n", " },\n", " \"Total_30N50N\": {\n", " \"interannual\": 0.8901421815356266,\n", @@ -798,7 +728,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "The precipitation variability metrics have a set of default regions, but users can instead define a single spatial region to compute metrics over. There are two ways to do this.\n", + "The precipitation variability metrics have a set of default regions. However, users can instead define a single spatial region to compute metrics over. There are two ways to do this.\n", "\n", "1. Use the `regions_specs` parameter to define a latitude/longitude box. \n", "Parameter file example:\n", @@ -813,29 +743,39 @@ "attr=\"NAME\" # An attribute in the shapefile\n", "feature=\"CONUS\" # A unique feature name that can be \n", " # found under the \"attr\" attribute\n", - "```" + "```\n", + "\n", + "Both options can be used at the same time. In that case, the area defined by regions_specs is applied first and can be used to trim down very large, high resolution datasets. Then the metrics are computed for the area defined by the shapefile region." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "First, we generate a simple shapefile for use in this example. The shapefile contains one feature, a box that defines the CONUS region." + "### Region example\n", + "First, we generate a simple shapefile for use in this demo. The shapefile contains one feature, a box that defines the CONUS region." ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "from shapely import Polygon\n", - "import geo pandas as gpd\n", + "import geopandas as gpd\n", "import pandas as pd\n", "\n", + "# Define region box\n", "coords = ((233.,22.),(233.,50.),(294.,50.),(294.,22))\n", + "\n", + "# Add to pandas dataframe, then convert to geopandas dataframe\n", "df = pd.DataFrame({\"Region\": [\"CONUS\"], \"Coords\": [Polygon(coords)]})\n", "gdf = gpd.GeoDataFrame(df, geometry=\"Coords\", crs=\"EPSG:4326\")\n", + "\n", + "# Create the output location\n", + "if not os.path.exists(demo_output_directory+\"/shp\"):\n", + " os.mkdir(demo_output_directory+\"/shp\")\n", "gdf.to_file(demo_output_directory+'/shp/CONUS.shp')" ] }, @@ -848,9 +788,55 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 14, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/home/ordonez4/miniconda3/envs/pmp_dev/lib/python3.10/site-packages/pcmdi_metrics/precip_variability/lib/lib_variability_across_timescales.py:313: RuntimeWarning: Mean of empty slice\n", + " clim = np.nanmean(dseg, axis=0)\n", + "INFO::2024-09-18 16:11::pcmdi_metrics:: Results saved to a json file: /home/ordonez4/git/pcmdi_metrics/doc/jupyter/Demo/demo_output_tmp/precip_variability/region_ex/PS_pr.day_regrid.180x90_area.freq.mean_GISS-E2-H.r6i1p1.json\n", + "2024-09-18 16:11:40,113 [INFO]: base.py(write:422) >> Results saved to a json file: /home/ordonez4/git/pcmdi_metrics/doc/jupyter/Demo/demo_output_tmp/precip_variability/region_ex/PS_pr.day_regrid.180x90_area.freq.mean_GISS-E2-H.r6i1p1.json\n", + "2024-09-18 16:11:40,113 [INFO]: base.py(write:422) >> Results saved to a json file: /home/ordonez4/git/pcmdi_metrics/doc/jupyter/Demo/demo_output_tmp/precip_variability/region_ex/PS_pr.day_regrid.180x90_area.freq.mean_GISS-E2-H.r6i1p1.json\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "demo_data_tmp/CMIP5_demo_timeseries/historical/atmos/day/pr/\n", + "pr_day_GISS-E2-H_historical_r6i1p1_*.nc\n", + "[2000, 2005]\n", + "730 365\n", + "2\n", + "demo_output_tmp/precip_variability/region_ex\n", + "demo_output_tmp/precip_variability/region_ex\n", + "demo_output_tmp/precip_variability/region_ex\n", + "['demo_data_tmp/CMIP5_demo_timeseries/historical/atmos/day/pr/pr_day_GISS-E2-H_historical_r6i1p1_20000101-20051231.nc']\n", + "GISS-E2-H.r6i1p1\n", + "['demo_data_tmp/CMIP5_demo_timeseries/historical/atmos/day/pr/pr_day_GISS-E2-H_historical_r6i1p1_20000101-20051231.nc']\n", + "GISS-E2-H.r6i1p1 365_day\n", + "2000 2005\n", + "Complete regridding from (2190, 90, 144) to (2190, 90, 180)\n", + "Cropping from shapefile\n", + "Reading region from file.\n", + "Complete calculating climatology and anomaly for calendar of 365_day\n", + "Complete power spectra (segment: 730 nps: 5.0 )\n", + "Complete domain and frequency average of spectral power\n", + "Complete power spectra (segment: 730 nps: 5.0 )\n", + "Complete domain and frequency average of spectral power\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "[WARNING] yaksa: 10 leaked handle pool objects\n" + ] + } + ], "source": [ "%%bash -s \"$demo_output_directory\"\n", "variability_across_timescales_PS_driver.py -p basic_precip_variability_param.py \\\n", @@ -869,22 +855,54 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 15, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "{\n", + " \"GISS-E2-H.r6i1p1\": {\n", + " \"forced\": {\n", + " \"CONUS\": {\n", + " \"annual\": 1.2011870574080201,\n", + " \"semi-annual\": 0.380975826207154\n", + " }\n", + " },\n", + " \"unforced\": {\n", + " \"CONUS\": {\n", + " \"interannual\": 0.1521909521737256,\n", + " \"seasonal-annual\": 0.20428410514869913,\n", + " \"sub-seasonal\": 0.20652699240276465,\n", + " \"synoptic\": 0.10360220715481439\n", + " }\n", + " }\n", + " }\n", + "}\n" + ] + } + ], "source": [ "output_path = os.path.join(demo_output_directory,\"precip_variability/region_ex/PS_pr.day_regrid.180x90_area.freq.mean_GISS-E2-H.r6i1p1.json\")\n", "with open(output_path) as f:\n", " metric = json.load(f)[\"RESULTS\"]\n", "print(json.dumps(metric, indent=2))" ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { "kernelspec": { - "display_name": "pmp_pv_dev", + "display_name": "Python [conda env:pmp_dev] *", "language": "python", - "name": "pmp_pv_dev" + "name": "conda-env-pmp_dev-py" }, "language_info": { "codemirror_mode": {