diff --git a/m2plots/time-series/README.txt b/m2plots/time-series/README.txt new file mode 100644 index 0000000..ac0cb6d --- /dev/null +++ b/m2plots/time-series/README.txt @@ -0,0 +1 @@ +This directory contains codes that are used to generate the monthly time series figures that are posted on the MERRA-2 Annual Cycle instance of FLUID. Additionally, latestyear_spaghetti_ncaregions_statD.py will produce a daily version of the figures for temperature, however, due to the additional computational burden of daily data points, statistics for the mean and percentiles must be pre-computed ahead of time. diff --git a/m2plots/time-series/latestyear_daily.py b/m2plots/time-series/latestyear_daily.py deleted file mode 100755 index 3ed9933..0000000 --- a/m2plots/time-series/latestyear_daily.py +++ /dev/null @@ -1,81 +0,0 @@ -#!/usr/bin/env python - -import numpy as np -import os.path -import xarray as xr -import metpy.calc as mpcalc -from metpy.units import units -import matplotlib.pyplot as plt - - -####INPUT parameters#### -variablename='T2MMAX' -varlongname='Daily Max 2 m Temperature' -units='K' -collection='statD_2d_slv_Nx' -lat1=29.5 -lat2=32 -lon1=-98 -lon2=-95 -#lat1=25 -#lat2=50 -#lon1=-125 -#lon2=-65 -region='SE Texas' -endyear=2024 -endmonth=5 -endday= -landonly=1 -ncaregion=0 -unitconversion=1 - -####LOAD DATA#### -DS = xr.open_mfdataset('/discover/nobackup/projects/gmao/merra2/data/products/MERRA2_all/Y*/M05/MERRA2.' + collection + '.*.nc4') -ncaregions=xr.open_dataset('/discover/nobackup/acollow/MERRA2/NCA_regs_MERRA-2.nc') -m2constants=xr.open_dataset('/discover/nobackup/projects/gmao/merra2/data/products/MERRA2_all/MERRA2.const_2d_asm_Nx.00000000.nc4') -subset=DS[variablename].sel(lon=slice(lon1,lon2),lat=slice(lat1,lat2)) -DS2024 = xr.open_mfdataset('/discover/nobackup/projects/gmao/merra2/data/products/MERRA2_400/Y2024/M05/MERRA2.' + collection + '.*.nc4') -subset2024=DS2024[variablename].sel(lon=slice(lon1,lon2),lat=slice(lat1,lat2)) - - -if landonly==1: - land=m2constants.FRLAND+m2constants.FRLANDICE - land_subset=land.sel(lon=slice(lon1,lon2),lat=slice(lat1,lat2)).squeeze(['time'],drop=True) - subset=subset.where(land_subset>0.3) - -if ncaregion>0: - nca_subset=ncaregions['regs05'].sel(lon=slice(lon1,lon2),lat=slice(lat1,lat2)) - subset=subset.where(nca_subset==ncaregion) - -####Get area average#### -weights=np.cos(np.deg2rad(subset.lat)) -#normalizedweights=(weights - np.min(weights)) / (np.max(weights) - np.min(weights)) -subset_weighted=subset.weighted(weights) -weighted_mean = unitconversion*subset_weighted.mean(("lon", "lat")) - -####Compute Stats#### -climo=weighted_mean.groupby("time.dayofyear").mean() -minimum=weighted_mean.groupby("time.dayofyear").min() -maximum=weighted_mean.groupby("time.dayofyear").max() -pctl15=weighted_mean.groupby('time.dayofyear').reduce(np.nanpercentile, dim='time', q=15) -pctl85=weighted_mean.groupby('time.dayofyear').reduce(np.nanpercentile, dim='time', q=85) -#print(weighted_mean.sel(time=slice(str(endyear) + "-01-01", str(endyear) + "-" + str(endmonth) + "-01"))) - -####Generate Figure#### -xaxis=np.arange(1,32,1) -fig, ax = plt.subplots() -line1=ax.plot(np.arange(1,endmonth+1,1),weighted_mean.sel(time=slice(str(endyear) + "-01-01", str(endyear) + "-" + str(endmonth) + "-01")),'r',label=str(endyear)) -line2=ax.plot(xaxis,climo,'k',label="Climo Mean") -line3=ax.fill_between(xaxis,pctl15,pctl85,color='lightgray',label="15th-85th Percentile") -ax.plot(xaxis,minimum,'k',linewidth=0.5) -line4=ax.plot(xaxis,maximum,'k',linewidth=0.5,label="Min/Max") -plt.ylabel(varlongname + ' (' + units + ')') -plt.xlabel('May') -ax.legend([str(endyear),"Climo Mean","15th-85th Percentile","Min/Max"]) -#plt.xticks(ticks=np.arange(1,13,1), labels=['Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec']) -plt.xlim([1,31]) -#plt.title(region + ' (Lon = ' + str(lon1) + ' to ' + str(lon2) + ', Lat = ' + str(lat1) + ' to ' + str(lat2) + ')') -plt.title(region) -plt.show() - -fig.savefig('%s_%s_%4d.png'%(variablename,region,endmonth)) diff --git a/m2plots/time-series/latestyear_spaghetti.py b/m2plots/time-series/latestyear_spaghetti.py deleted file mode 100755 index 839d404..0000000 --- a/m2plots/time-series/latestyear_spaghetti.py +++ /dev/null @@ -1,76 +0,0 @@ -#!/usr/bin/env python - -import numpy as np -import os.path -import xarray as xr -import metpy.calc as mpcalc -from metpy.units import units -import matplotlib.pyplot as plt - - -####INPUT parameters#### -variablename='T2M' -varlongname='2 m Temperature' -units='K' -collection='tavgM_2d_slv_Nx' -lat1=-90 -lat2=90 -lon1=-180 -lon2=179.625 -#lat1=25 -#lat2=50 -#lon1=-125 -#lon2=-65 -region='Sudan' -endyear=2023 -endmonth=2 -landonly=0 -ncaregion=0 -unitconversion=1 - -####LOAD DATA#### -DS = xr.open_mfdataset('/discover/nobackup/projects/gmao/merra2/data/products/MERRA2_all/Y*/M*/MERRA2.' + collection + '.*.nc4') -ncaregions=xr.open_dataset('/discover/nobackup/acollow/MERRA2/NCA_regs_MERRA-2.nc') -m2constants=xr.open_dataset('/discover/nobackup/projects/gmao/merra2/data/products/MERRA2_all/MERRA2.const_2d_asm_Nx.00000000.nc4') -subset=DS[variablename].sel(lon=slice(lon1,lon2),lat=slice(lat1,lat2)) - -if landonly==1: - land=m2constants.FRLAND+m2constants.FRLANDICE - land_subset=land.sel(lon=slice(lon1,lon2),lat=slice(lat1,lat2)).squeeze(['time'],drop=True) - subset=subset.where(land_subset>0.3) - -if ncaregion>0: - nca_subset=ncaregions['regs05'].sel(lon=slice(lon1,lon2),lat=slice(lat1,lat2)) - subset=subset.where(nca_subset==ncaregion) - -####Get area average#### -weights=np.cos(np.deg2rad(subset.lat)) -#normalizedweights=(weights - np.min(weights)) / (np.max(weights) - np.min(weights)) -subset_weighted=subset.weighted(weights) -weighted_mean = unitconversion*subset_weighted.mean(("lon", "lat")) - -####Compute Stats#### -climo=weighted_mean.groupby("time.month").mean() -minimum=weighted_mean.groupby("time.month").min() -maximum=weighted_mean.groupby("time.month").max() -pctl15=weighted_mean.groupby('time.month').reduce(np.nanpercentile, dim='time', q=15) -pctl85=weighted_mean.groupby('time.month').reduce(np.nanpercentile, dim='time', q=85) -#print(weighted_mean.sel(time=slice(str(endyear) + "-01-01", str(endyear) + "-" + str(endmonth) + "-01"))) - -####Generate Figure#### -xaxis=np.arange(1,13,1) -fig, ax = plt.subplots() -line1=ax.plot(np.arange(1,endmonth+1,1),weighted_mean.sel(time=slice(str(endyear) + "-01-01", str(endyear) + "-" + str(endmonth) + "-01")),'r',label=str(endyear)) -line2=ax.plot(xaxis,climo,'k',label="Climo Mean") -line3=ax.fill_between(xaxis,pctl15,pctl85,color='lightgray',label="15th-85th Percentile") -ax.plot(xaxis,minimum,'k',linewidth=0.5) -line4=ax.plot(xaxis,maximum,'k',linewidth=0.5,label="Min/Max") -plt.ylabel(varlongname + ' (' + units + ')') -ax.legend([str(endyear),"Climo Mean","15th-85th Percentile","Min/Max"]) -plt.xticks(ticks=np.arange(1,13,1), labels=['Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec']) -plt.xlim([1,12]) -#plt.title(region + ' (Lon = ' + str(lon1) + ' to ' + str(lon2) + ', Lat = ' + str(lat1) + ' to ' + str(lat2) + ')') -plt.title(region) -plt.show() - -fig.savefig('%s_%s_%4d.png'%(variablename,region,endyear)) diff --git a/m2plots/time-series/latestyear_spaghetti_ncaregions_precip.py b/m2plots/time-series/latestyear_spaghetti_ncaregions_precip.py deleted file mode 100755 index 7da687c..0000000 --- a/m2plots/time-series/latestyear_spaghetti_ncaregions_precip.py +++ /dev/null @@ -1,72 +0,0 @@ -#!/usr/bin/env python - -import numpy as np -import os.path -import xarray as xr -import metpy.calc as mpcalc -from metpy.units import units -import matplotlib.pyplot as plt - - -####INPUT parameters#### -variablename='PRECTOT' -varlongname='Precipitation' -units='mm/day' -collection='tavgM_2d_flx_Nx' -lat1=25 -lat2=50 -lon1=-125 -lon2=-65 -region='Southeast United States' -endyear=2005 -endmonth=12 -landonly=0 -ncaregion=2 -unitconversion=86400 - -####LOAD DATA#### -DS = xr.open_mfdataset('/discover/nobackup/projects/gmao/merra2/data/products/MERRA2_all/Y*/M*/MERRA2.' + collection + '.*.nc4') -ncaregions=xr.open_dataset('/discover/nobackup/acollow/MERRA2/NCA_regs_MERRA-2.nc') -m2constants=xr.open_dataset('/discover/nobackup/projects/gmao/merra2/data/products/MERRA2_all/MERRA2.const_2d_asm_Nx.00000000.nc4') -subset=DS[variablename].sel(lon=slice(lon1,lon2),lat=slice(lat1,lat2)) - -if landonly==1: - land=m2constants.FRLAND+m2constants.FRLANDICE - land_subset=land.sel(lon=slice(lon1,lon2),lat=slice(lat1,lat2)).squeeze(['time'],drop=True) - subset=subset.where(land_subset>0.3) - -if ncaregion>0: - nca_subset=ncaregions['regs05'].sel(lon=slice(lon1,lon2),lat=slice(lat1,lat2)) - subset=subset.where(nca_subset==ncaregion) - -####Get area average#### -weights=np.cos(np.deg2rad(subset.lat)) -#normalizedweights=(weights - np.min(weights)) / (np.max(weights) - np.min(weights)) -subset_weighted=subset.weighted(weights) -weighted_mean = unitconversion*subset_weighted.mean(("lon", "lat")) - -####Compute Stats#### -climo=weighted_mean.groupby("time.month").mean() -minimum=weighted_mean.groupby("time.month").min() -maximum=weighted_mean.groupby("time.month").max() -pctl15=weighted_mean.groupby('time.month').reduce(np.nanpercentile, dim='time', q=15) -pctl85=weighted_mean.groupby('time.month').reduce(np.nanpercentile, dim='time', q=85) -#print(weighted_mean.sel(time=slice(str(endyear) + "-01-01", str(endyear) + "-" + str(endmonth) + "-01"))) - -####Generate Figure#### -xaxis=np.arange(1,13,1) -fig, ax = plt.subplots() -line1=ax.plot(np.arange(1,endmonth+1,1),weighted_mean.sel(time=slice(str(endyear) + "-01-01", str(endyear) + "-" + str(endmonth) + "-01")),'r',label=str(endyear)) -line2=ax.plot(xaxis,climo,'k',label="Climo Mean") -line3=ax.fill_between(xaxis,pctl15,pctl85,color='lightgray',label="15th-85th Percentile") -ax.plot(xaxis,minimum,'k',linewidth=0.5) -line4=ax.plot(xaxis,maximum,'k',linewidth=0.5,label="Min/Max") -plt.ylabel(varlongname + ' (' + units + ')') -ax.legend([str(endyear),"Climo Mean","15th-85th Percentile","Min/Max"]) -plt.xticks(ticks=np.arange(1,13,1), labels=['Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec']) -plt.xlim([1,12]) -#plt.title(region + ' (Lon = ' + str(lon1) + ' to ' + str(lon2) + ', Lat = ' + str(lat1) + ' to ' + str(lat2) + ')') -plt.title(region) -plt.show() - -fig.savefig('%s_%s_%4d.png'%(variablename,region,endyear)) diff --git a/m2plots/time-series/latestyear_spaghetti_ncaregions_statD.py b/m2plots/time-series/latestyear_spaghetti_ncaregions_statD.py new file mode 100755 index 0000000..da650d4 --- /dev/null +++ b/m2plots/time-series/latestyear_spaghetti_ncaregions_statD.py @@ -0,0 +1,229 @@ +#!/usr/bin/env python + +import numpy as np +import os.path +import xarray as xr +import metpy.calc as mpcalc +from metpy.units import units +import matplotlib.pyplot as plt +import yaml +import warnings +warnings.filterwarnings('ignore') +import sys +import datetime +import netCDF4 + +#Required command line inputs in order are 1) the variable to be plotted, 2) the region, 3) the year to be highlighted, and 4) the ending month of the year to be highlighted. +#Example usage: ./latestyear_spaghetti_ncaregions.py t2m ne 2023 12 +yamlkey_var=sys.argv[1] +yamlkey_reg=sys.argv[2] +endyear=int(sys.argv[3]) +endmonth=int(sys.argv[4]) + +dt=datetime.datetime(endyear,endmonth,11) +tt=dt.timetuple().tm_yday + +####INPUT parameters#### +#yamlkey_var='t2m' +#yamlkey_reg='ne' +#endyear=2023 +#endmonth=12 +######################## + +NCA_map = """ +ne: + region: 'Northeast United States' + regionshortname: 'neus' + regionnumber: 1 + landonly: 1 + lat1: 24 + lat2: 50 + lon1: -125 + lon2: -65 +se: + region: 'Southeast United States' + regionshortname: 'seus' + regionnumber: 2 + landonly: 1 + lat1: 24 + lat2: 50 + lon1: -125 + lon2: -65 +mw: + region: 'Midwest United States' + regionshortname: 'mwus' + regionnumber: 3 + landonly: 1 + lat1: 24 + lat2: 50 + lon1: -125 + lon2: -65 +ngp: + region: 'Northern Great Plains' + regionshortname: 'ngpus' + regionnumber: 4 + landonly: 1 + lat1: 24 + lat2: 50 + lon1: -125 + lon2: -65 +sgp: + region: 'Southern Great Plains' + regionshortname: 'sgpus' + regionnumber: 5 + landonly: 1 + lat1: 24 + lat2: 50 + lon1: -125 + lon2: -65 +nw: + region: 'Northwest United States' + regionshortname: 'nwus' + regionnumber: 6 + landonly: 1 + lat1: 24 + lat2: 50 + lon1: -125 + lon2: -65 +sw: + region: 'Southwest United States' + regionshortname: 'swus' + regionnumber: 7 + landonly: 1 + lat1: 24 + lat2: 50 + lon1: -125 + lon2: -65 +conus: + region: 'CONUS' + regionshortname: 'conus' + regionnumber: 0 + landonly: 1 + lat1: 24 + lat2: 50 + lon1: -125 + lon2: -65 +global: + region: 'Global' + regionshortname: 'global' + regionnumber: 0 + landonly: 0 + lat1: -90 + lat2: 90 + lon1: -180 + lon2: 179.375 +globalland: + region: 'Global Land' + regionshortname: 'globalland' + regionnumber: 0 + landonly: 1 + lat1: -90 + lat2: 90 + lon1: -180 + lon2: 179.375 +globalland6060: + region: 'Land 60N to 60S' + regionshortname: 'globalland6060' + regionnumber: 0 + landonly: 1 + lat1: -60 + lat2: 60 + lon1: -180 + lon2: 179.375 + +""" + +variable_map = """ +t2mmean: + variablename: 'T2MMEAN' + varlongname: '2 m Mean Temperature' + units: 'K' + unitconversion: 1 + collection: 'statD_2d_slv_Nx' + +t2mmax: + variablename: 'T2MMAX' + varlongname: '2 m Maximum Temperature' + units: 'K' + unitconversion: 1 + collection: 'statD_2d_slv_Nx' + +t2mmin: + variablename: 'T2MMIN' + varlongname: '2 m Minimum Temperature' + units: 'K' + unitconversion: 1 + collection: 'statD_2d_slv_Nx' + + +""" + +####Import yaml info#### +var = yaml.safe_load(variable_map) +region = yaml.safe_load(NCA_map) + + +####LOAD DATA#### +DS = xr.open_mfdataset('/discover/nobackup/projects/gmao/merra2/data/products/MERRA2_400/Y'+str(endyear)+'/M0?/MERRA2_400.' + var[yamlkey_var]['collection'] + '.2*.nc4') +lon1=region[yamlkey_reg]['lon1'] +lon2=region[yamlkey_reg]['lon2'] +lat1=region[yamlkey_reg]['lat1'] +lat2=region[yamlkey_reg]['lat2'] + + +####Subset for Selected Region#### +subset=DS[var[yamlkey_var]['variablename']].sel(lon=slice(lon1,lon2),lat=slice(lat1,lat2)) +if region[yamlkey_reg]['landonly']==1: + m2constants=xr.open_dataset('/discover/nobackup/projects/gmao/merra2/data/products/MERRA2_all/MERRA2.const_2d_asm_Nx.00000000.nc4') + land=m2constants.FRLAND+m2constants.FRLANDICE + land_subset=land.sel(lon=slice(lon1,lon2),lat=slice(lat1,lat2)).squeeze(['time'],drop=True) + subset=subset.where(land_subset>0.3) + +if region[yamlkey_reg]['regionnumber']>0 & region[yamlkey_reg]['regionnumber']<8: + ncaregions=xr.open_dataset('/discover/nobackup/acollow/MERRA2/NCA_regs_MERRA-2.nc') + nca_subset=ncaregions['regs05'].sel(lon=slice(lon1,lon2),lat=slice(lat1,lat2)) + subset=subset.where(nca_subset==region[yamlkey_reg]['regionnumber']) +elif region[yamlkey_reg]['regionnumber']==8: + nca_subset=ncaregions['regs05'].sel(lon=slice(lon1,lon2),lat=slice(lat1,lat2)) + print(max(nca_subset)) + subset=subset.where(nca_subset>0 & nca_subset<8) + +####Get area average#### +weights=np.cos(np.deg2rad(subset.lat)) +subset_weighted=subset.weighted(weights) +weighted_mean = var[yamlkey_var]['unitconversion']*subset_weighted.mean(("lon", "lat")) + +####Compute Stats#### +#stats_subset=weighted_mean.sel(time=slice("1980-01-01","2023-12-01")) +#climo=stats_subset.groupby("time.month").mean() +#minimum=stats_subset.groupby("time.month").min() +#maximum=stats_subset.groupby("time.month").max() +#pctl15=stats_subset.groupby('time.month').reduce(np.nanpercentile, dim='time', q=15) +#pctl85=stats_subset.groupby('time.month').reduce(np.nanpercentile, dim='time', q=85) +#print(weighted_mean.sel(time=slice(str(endyear) + "-01-01", str(endyear) + "-" + str(endmonth) + "-01"))) + +#### Load stats +f1 = netCDF4.Dataset(yamlkey_var + '_'+ yamlkey_reg +'_stats.nc','r'); +climo=f1.variables['tmean'][0:365] +minimum=f1.variables['tmin'][0:365] +maximum=f1.variables['tmax'][0:365] +pctl15=f1.variables['tmeanpc15'][0:365] +pctl85=f1.variables['tmeanpc85'][0:365] + +####Generate Figure#### +xaxis=np.arange(1,366,1) +fig, ax = plt.subplots() +ax.tick_params(axis='both', which='major', labelsize=14) +line1=ax.plot(np.arange(1,tt+1,1),weighted_mean,'r',label=str(endyear)) +line2=ax.plot(xaxis,climo,'k',label="Climo Mean") +line3=ax.fill_between(xaxis,pctl15,pctl85,color='lightgray',label="15th-85th Percentile") +ax.plot(xaxis,minimum,'k',linewidth=0.5) +line4=ax.plot(xaxis,maximum,'k',linewidth=0.5,label="Min/Max") +plt.ylabel(var[yamlkey_var]['varlongname'] + ' (' + var[yamlkey_var]['units'] + ')', fontsize=14) +ax.legend([str(endyear),"Climo Mean (1980-2023)","15th-85th Percentile","Min/Max"]) +plt.xticks(ticks=[15,47,76,107,137,167,197,227,260,290,320,350], labels=['Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec']) +plt.xlim([1,366]) +plt.title(region[yamlkey_reg]['region'], fontsize=14, fontweight='bold') +plt.show() + +#fig.savefig('%s_%s_%4d.png'%(var[yamlkey_var]['variablename'],region[yamlkey_reg]['regionshortname'],endyear)) diff --git a/m2plots/time-series/latestyear_spaghetti_singlepoint.py b/m2plots/time-series/latestyear_spaghetti_singlepoint.py deleted file mode 100755 index c9274d5..0000000 --- a/m2plots/time-series/latestyear_spaghetti_singlepoint.py +++ /dev/null @@ -1,38 +0,0 @@ -#!/usr/bin/env python - -import numpy as np -import os.path -import xarray as xr -import metpy.calc as mpcalc -from metpy.units import units -import matplotlib.pyplot as plt - - -####INPUT parameters#### -variablename='T2M' -collection='tavgM_2d_slv_Nx' -latinput=39.3 -loninput=-76.6 -location='Baltimore, MD' -endyear=2023 -endmonth=2 - -####LOAD DATA#### -DS = xr.open_mfdataset('/discover/nobackup/projects/gmao/merra2/data/products/MERRA2_all/Y*/M*/MERRA2.' + collection + '.*.nc4') -spatialsubset=DS.sel(lon=loninput,lat=latinput,method="nearest") -climosubset=spatialsubset.sel(time=slice("1991-01-01", "2020-12-01")) -data=climosubset.T2M - -####Compute Stats#### -varmean=data.groupby("time.month").mean() -pctl85=data.groupby("time.month").quantile(0.85) -pctl15=data.groupby("time.month").quantile(0.15) - - -####Create Figure#### -xaxis=np.arange(1,13,1) -fig, ax = plt.subplots() -ax.fill_between(xaxis,pctl15,pctl85,color='lightgray') -ax.plot(xaxis,varmean) -ax.plot(np.arange(1,endmonth+1,1),spatialsubset.T2M.sel(time=slice(str(endyear) + "-01-01", str(endyear) + "-" + str(endmonth) + "-01"))) -plt.show() diff --git a/m2plots/time-series/latestyear_spaghetti_sudan.py b/m2plots/time-series/latestyear_spaghetti_sudan.py deleted file mode 100755 index 502854f..0000000 --- a/m2plots/time-series/latestyear_spaghetti_sudan.py +++ /dev/null @@ -1,76 +0,0 @@ -#!/usr/bin/env python - -import numpy as np -import os.path -import xarray as xr -import metpy.calc as mpcalc -from metpy.units import units -import matplotlib.pyplot as plt - - -####INPUT parameters#### -variablename='T2MMAX' -varlongname='Maximum 2 m Temperature' -units='K' -collection='statM_2d_slv_Nx' -lat1=3 -lat2=22 -lon1=22 -lon2=38 -#lat1=25 -#lat2=50 -#lon1=-125 -#lon2=-65 -region='Sudan' -endyear=2024 -endmonth=3 -landonly=0 -ncaregion=0 -unitconversion=1 - -####LOAD DATA#### -DS = xr.open_mfdataset('/discover/nobackup/projects/gmao/merra2/data/products/MERRA2_400/Y*/M*/MERRA2_400.' + collection + '.*.nc4') -ncaregions=xr.open_dataset('/discover/nobackup/acollow/MERRA2/NCA_regs_MERRA-2.nc') -m2constants=xr.open_dataset('/discover/nobackup/projects/gmao/merra2/data/products/MERRA2_all/MERRA2.const_2d_asm_Nx.00000000.nc4') -subset=DS[variablename].sel(lon=slice(lon1,lon2),lat=slice(lat1,lat2)) - -if landonly==1: - land=m2constants.FRLAND+m2constants.FRLANDICE - land_subset=land.sel(lon=slice(lon1,lon2),lat=slice(lat1,lat2)).squeeze(['time'],drop=True) - subset=subset.where(land_subset>0.3) - -if ncaregion>0: - nca_subset=ncaregions['regs05'].sel(lon=slice(lon1,lon2),lat=slice(lat1,lat2)) - subset=subset.where(nca_subset==ncaregion) - -####Get area average#### -weights=np.cos(np.deg2rad(subset.lat)) -#normalizedweights=(weights - np.min(weights)) / (np.max(weights) - np.min(weights)) -subset_weighted=subset.weighted(weights) -weighted_mean = unitconversion*subset_weighted.mean(("lon", "lat")) - -####Compute Stats#### -climo=weighted_mean.groupby("time.month").mean() -minimum=weighted_mean.groupby("time.month").min() -maximum=weighted_mean.groupby("time.month").max() -pctl15=weighted_mean.groupby('time.month').reduce(np.nanpercentile, dim='time', q=15) -pctl85=weighted_mean.groupby('time.month').reduce(np.nanpercentile, dim='time', q=85) -#print(weighted_mean.sel(time=slice(str(endyear) + "-01-01", str(endyear) + "-" + str(endmonth) + "-01"))) - -####Generate Figure#### -xaxis=np.arange(1,13,1) -fig, ax = plt.subplots() -line1=ax.plot(np.arange(1,endmonth+1,1),weighted_mean.sel(time=slice(str(endyear) + "-01-01", str(endyear) + "-" + str(endmonth+1) + "-01")),'r',label=str(endyear)) -line2=ax.plot(xaxis,climo,'k',label="Climo Mean") -line3=ax.fill_between(xaxis,pctl15,pctl85,color='lightgray',label="15th-85th Percentile") -ax.plot(xaxis,minimum,'k',linewidth=0.5) -line4=ax.plot(xaxis,maximum,'k',linewidth=0.5,label="Min/Max") -plt.ylabel(varlongname + ' (' + units + ')') -ax.legend([str(endyear),"Climo Mean","15th-85th Percentile","Min/Max"]) -plt.xticks(ticks=np.arange(1,13,1), labels=['Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec']) -plt.xlim([1,12]) -#plt.title(region + ' (Lon = ' + str(lon1) + ' to ' + str(lon2) + ', Lat = ' + str(lat1) + ' to ' + str(lat2) + ')') -plt.title(region) -plt.show() - -fig.savefig('%s_%s_%4d.png'%(variablename,region,endyear)) diff --git a/m2stats/NCAindices/monthly/README.txt b/m2stats/NCAindices/monthly/README.txt new file mode 100644 index 0000000..8882a8b --- /dev/null +++ b/m2stats/NCAindices/monthly/README.txt @@ -0,0 +1,2 @@ +The codes in this directory are used to generate the extremes detection indices files that are distributed through GES DISC and posted as imagery on FLUID, in addtion to the code that is used to compute monthly wind speeds for FLUID. +Each code uses an input of "yyyy mm climoperiod", where climoperiod is currently 19902020. The codes are too hardwired at the moment for MERRA-2 and precomputed percentiles, which will hopefully be modified in future versions. diff --git a/m2stats/NCAindices/seasonal/get5dayprecip_seasonal.bash b/m2stats/NCAindices/seasonal/get5dayprecip_seasonal.bash new file mode 100755 index 0000000..3a41e3d --- /dev/null +++ b/m2stats/NCAindices/seasonal/get5dayprecip_seasonal.bash @@ -0,0 +1,62 @@ +#!/bin/bash + + +yr2=$1 +ssnnum=$2 + +if [ ${ssnnum} -eq 1 ]; then + m1=oct + m2=feb + yr1=`expr ${yr2} - 1` +elif [ ${ssnnum} -eq 2 ]; then + m1=nov + m2=mar + yr1=`expr ${yr2} - 1` +elif [ ${ssnnum} -eq 3 ]; then + m1=dec + m2=apr + yr1=`expr ${yr2} - 1` +elif [ ${ssnnum} -eq 4 ]; then + m1=jan + m2=may + yr1=$yr2 +elif [ ${ssnnum} -eq 5 ]; then + m1=feb + m2=jun + yr1=$yr2 +elif [ ${ssnnum} -eq 6 ]; then + m1=mar + m2=jul + yr1=$yr2 +elif [ ${ssnnum} -eq 7 ]; then + m1=apr + m2=aug + yr1=$yr2 +elif [ ${ssnnum} -eq 8 ]; then + m1=may + m2=sep + yr1=$yr2 +elif [ ${ssnnum} -eq 9 ]; then + m1=jun + m2=oct + yr1=$yr2 +elif [ ${ssnnum} -eq 10 ]; then + m1=jul + m2=nov + yr1=$yr2 +elif [ ${ssnnum} -eq 11 ]; then + m1=aug + m2=dec + yr1=$yr2 +elif [ ${ssnnum} -eq 12 ]; then + m1=sep + m2=jan + yr1=${yr2} + yr2=`expr ${yr2} + 1` +else + echo "Incorrect Season Selection" +fi + +/discover/nobackup/projects/gmao/share/dasilva/opengrads/Contents/lats4d.sh -i /discover/nobackup/acollow/MERRA2/opendap/tavg1_2d_flx_Nx -o MERRA2.prectot.daily.uptodate -time 00:30Z01${m1}${yr1} 00:30Z03${m2}${yr2} 24 -func "mean(@,t-12,t+11)*86400" -vars prectot -mxtimes 400000 +wait +/discover/nobackup/projects/gmao/share/dasilva/opengrads/Contents/lats4d.sh -i MERRA2.prectot.daily.uptodate.nc -o MERRA2.prectot.5daytotal.uptodate -time 0:30Z01${m1}${yr1} 0:30Z01${m2}${yr2} 1 -func "sum(@,t-2,t+2)" -vars prectot -mxtimes 366