diff --git a/doc/jupyter/Demo/Demo_7_precip_variability.ipynb b/doc/jupyter/Demo/Demo_7_precip_variability.ipynb index 4e5652d20..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" ] } ], @@ -586,7 +509,9 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "tags": [] + }, "source": [ "## Precipitation Variability Metric" ] @@ -614,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": [ @@ -651,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", @@ -701,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", @@ -785,6 +717,179 @@ "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. 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", + "```\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", + "```\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": [ + "### 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": 13, + "metadata": {}, + "outputs": [], + "source": [ + "from shapely import Polygon\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')" + ] + }, + { + "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": 14, + "metadata": {}, + "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", + "--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": 15, + "metadata": {}, + "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, @@ -795,9 +900,9 @@ ], "metadata": { "kernelspec": { - "display_name": "Python [conda env:pmp_v20220110] *", + "display_name": "Python [conda env:pmp_dev] *", "language": "python", - "name": "conda-env-pmp_v20220110-py" + "name": "conda-env-pmp_dev-py" }, "language_info": { "codemirror_mode": { @@ -809,7 +914,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.9.10" + "version": "3.10.12" } }, "nbformat": 4, diff --git a/docs/metrics_precip-variability.rst b/docs/metrics_precip-variability.rst index 6655ef14a..7b5b08028 100644 --- a/docs/metrics_precip-variability.rst +++ b/docs/metrics_precip-variability.rst @@ -62,6 +62,11 @@ 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. +* **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 diff --git a/pcmdi_metrics/io/__init__.py b/pcmdi_metrics/io/__init__.py index 6414313fd..3e2b7d7c9 100644 --- a/pcmdi_metrics/io/__init__.py +++ b/pcmdi_metrics/io/__init__.py @@ -3,6 +3,7 @@ from .string_constructor import StringConstructor, fill_template # noqa # isort:skip from . import base # noqa from .base import MV2Json # noqa + from .xcdat_dataset_io import ( # noqa # isort:skip da_to_ds, get_axis_list, 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, diff --git a/pcmdi_metrics/precip_variability/lib/argparse_functions.py b/pcmdi_metrics/precip_variability/lib/argparse_functions.py index b447f02a0..e71a34292 100644 --- a/pcmdi_metrics/precip_variability/lib/argparse_functions.py +++ b/pcmdi_metrics/precip_variability/lib/argparse_functions.py @@ -1,3 +1,6 @@ +import ast + + def AddParserArgument(P): P.add_argument( "--mip", type=str, dest="mip", default=None, help="cmip5, cmip6 or other mip" @@ -57,6 +60,21 @@ 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( + "--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", @@ -71,6 +89,21 @@ 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 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 fd02ff1f2..0fd53adab 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 @@ -11,61 +10,111 @@ from scipy.stats import chi2 from xcdat.regridder import grid +from pcmdi_metrics.io import region_subset 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, 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 """ - psdmfm = {"RESULTS": {}} - f = xcdat.open_mfdataset(file) + 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) + 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"] + print(dat, cal) + print(syr, eyr) + if "360" in cal: - ldy = 30 + f = f.sel( + time=slice( + str(syr) + "-01-01 00:00:00", + str(eyr) + "-12-" + str(30) + " 23:59:59", + ) + ) else: - ldy = 31 - print(dat, cal) - print("syr, eyr:", syr, eyr) - for iyr in range(syr, eyr + 1): - print(iyr) - do = f.sel( + f = f.sel( time=slice( - str(iyr) + "-01-01 00:00:00", str(iyr) + "-12-" + str(ldy) + " 23:59:59" + str(syr) + "-01-01 00:00:00", + str(eyr) + "-12-" + str(31) + " 23:59:59", ) ) - # Regridding - rgtmp = Regrid2deg(do, var) * float(fac) - if iyr == syr: - drg = copy.deepcopy(rgtmp) - else: - drg = xr.concat([drg, rgtmp], dim="time") - print(iyr, drg.shape) + # Regridding + xvar = find_lon(f) + yvar = find_lat(f) + if np.min(f[xvar]) < 0: + lon_range = [-180.0, 180.0] + else: + lon_range = [0.0, 360.0] + 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) + nlon = str(len(drg[xvar])) + nlat = str(len(drg[yvar])) 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 freqs, ps, rn, sig95 = Powerspectrum(drg, nperseg, noverlap) # Domain & Frequency average - psdmfm_forced = Avg_PS_DomFrq(ps, freqs, ntd, dat, mip, "forced") + 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 = "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( @@ -76,9 +125,19 @@ 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.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( @@ -92,7 +151,15 @@ 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 = Base(outdir.replace("%(output_type)", "metrics_results"), outfilename) JSON.write( @@ -107,7 +174,29 @@ def precip_variability_across_timescale( # ================================================================================== -def Regrid2deg(d, var): +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, lon_range=[0.0, 360.0]): """ Regrid to 2deg (180lon*90lat) horizontal resolution Input @@ -116,15 +205,46 @@ def Regrid2deg(d, var): Output - drg: xCDAT variable with 2deg horizontal resolution """ - tgrid = grid.create_uniform_grid(-89, 89, 2.0, 0.0, 358.0, 2.0) + + start_lat = -90.0 + res / 2.0 + start_lon = lon_range[0] + end_lat = 90.0 - res / 2.0 + 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: + tgrid = CropLatLon(tgrid, regions_specs) drg = d.regridder.horizontal( - var, tgrid, tool="xesmf", method="conservative_normed", periodic=True - )[var] + var, + tgrid, + tool="xesmf", + method="conservative_normed", + periodic=True, + unmapped_to_nan=True, + ) + + print("Complete regridding from", d[var].shape, "to", drg[var].shape) - print("Complete regridding from", d[var].shape, "to", drg.shape) return drg +# ================================================================================== +def CropLatLon(d, regions_specs): + """ + Select a subgrid of the dataset defined by the regions_specs dictionary. + Input + - d: xCDAT variable + - regions_specs: a dictionary + Output + - dnew: xCDAT variable selected over region of interest + """ + region_name = list(regions_specs.keys())[0] + + dnew = region_subset(d, region_name, regions_specs=regions_specs) + + return dnew + + # ================================================================================== def ClimAnom(d, ntd, syr, eyr, cal): """ @@ -142,7 +262,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) @@ -194,7 +314,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( @@ -244,6 +364,7 @@ def Powerspectrum(d, nperseg, noverlap): - rn: Rednoise - sig95: 95% rednoise confidence level """ + # Fill missing date using interpolation dnp = np.array(d) dfm = np.zeros((d.shape[0], d.shape[1], d.shape[2]), dtype=float) @@ -337,7 +458,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 @@ -347,9 +468,17 @@ 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) """ + + 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 = [ "Total_50S50N", "Ocean_50S50N", @@ -364,6 +493,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,14 +520,15 @@ def Avg_PS_DomFrq(d, frequency, ntd, dat, mip, frc): # generate land sea mask mask = create_land_sea_mask(d[0]) + yvar = find_lat(d) psdmfm = {} for dom in domains: psdmfm[dom] = {} - if "Ocean" in dom: + if "Ocean_" in dom or val_from_rs(regions_specs, dom) == 0: dmask = d.where(mask == 0) - elif "Land" in dom: + elif "Land_" in dom or val_from_rs(regions_specs, dom) == 1: dmask = d.where(mask == 1) else: dmask = d @@ -406,23 +538,40 @@ def Avg_PS_DomFrq(d, frequency, ntd, dat, mip, frc): 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"] - if "30N50N" in dom: - am = dmask.sel(lat=slice(30, 50)).spatial.average( + elif "30N50N" in dom: + am = dmask.sel({yvar: slice(30, 50)}).spatial.average( "ps", axis=["X", "Y"], weights="generate" )["ps"] - if "30S30N" in dom: - am = dmask.sel(lat=slice(-30, 30)).spatial.average( + elif "30S30N" in dom: + am = dmask.sel({yvar: slice(-30, 30)}).spatial.average( "ps", axis=["X", "Y"], weights="generate" )["ps"] - if "50S30S" in dom: - am = dmask.sel(lat=slice(-50, -30)).spatial.average( + elif "50S30S" in dom: + am = dmask.sel({yvar: 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 = 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) @@ -456,7 +605,7 @@ def Avg_PS_DomFrq(d, frequency, ntd, dat, mip, frc): elif frq == "interannual": # 365day=