diff --git a/parm/templates/template.plot_hofx.yaml b/parm/templates/template.plot_hofx.yaml index 1b23f1a6..fbf9b0ee 100644 --- a/parm/templates/template.plot_hofx.yaml +++ b/parm/templates/template.plot_hofx.yaml @@ -5,3 +5,4 @@ nbins: XXNBINS plottype: 'XXPLOTTYPE' title_fig: 'XXFIGTITLE' output_prefix: 'XXPREOUTFN' +machine: 'XXMACHINE' diff --git a/scripts/exlandda_plot_stats.sh b/scripts/exlandda_plot_stats.sh index a2ab8b61..208f32bb 100755 --- a/scripts/exlandda_plot_stats.sh +++ b/scripts/exlandda_plot_stats.sh @@ -32,6 +32,7 @@ sed -i -e "s/XXNBINS/${NBINS}/g" plot_hofx.yaml sed -i -e "s/XXPLOTTYPE/${PLOTTYPE}/g" plot_hofx.yaml sed -i -e "s/XXFIGTITLE/${FIGTITLE}/g" plot_hofx.yaml sed -i -e "s/XXPREOUTFN/${PREOUTFN}/g" plot_hofx.yaml +sed -i -e "s/XXMACHINE/${MACHINE}/g" plot_hofx.yaml ${USHlandda}/hofx_analysis_stats.py if [[ $? != 0 ]]; then diff --git a/ush/hofx_analysis_stats.py b/ush/hofx_analysis_stats.py index 0e9156d6..38c33ca3 100755 --- a/ush/hofx_analysis_stats.py +++ b/ush/hofx_analysis_stats.py @@ -8,6 +8,7 @@ from scipy.stats import norm import cartopy import cartopy.crs as ccrs +import cartopy.feature as cfeature import xarray as xr import matplotlib.ticker import matplotlib as mpl @@ -15,37 +16,37 @@ def get_obs_stats(fdir, plottype): global lat,lon - omb_ = [] - oma_ = [] - obs_ = [] - qc_ = [] - err_ = [] - lat_ = [] - lon_ = [] + omb_=[] + oma_=[] + obs_=[] + qc_ =[] + err_=[] + lat_=[] + lon_=[] for fname in os.listdir(fdir): print("=== File Name:",fname) - f = netCDF4.Dataset(fdir+'/'+fname) + f=netCDF4.Dataset(fdir+'/'+fname) # print("NETCDF:",f) - obs = f.groups['ObsValue'].variables['totalSnowDepth'] + obs=f.groups['ObsValue'].variables['totalSnowDepth'] # print("ObsValue:",obs) - ombg = f.groups['ombg'].variables['totalSnowDepth'] + ombg=f.groups['ombg'].variables['totalSnowDepth'] # print("OMBG:",ombg) - oman = f.groups['oman'].variables['totalSnowDepth'] + oman=f.groups['oman'].variables['totalSnowDepth'] # print("OMAN:",oman) - qc = f.groups['PreQC'].variables['totalSnowDepth'] + qc=f.groups['PreQC'].variables['totalSnowDepth'] # print("PreQC:",qc) - obstime = f.groups['MetaData'].variables['dateTime'] + obstime=f.groups['MetaData'].variables['dateTime'] # print("OBS_TIME:",obstime) - if plottype == 'histogram': - ombg_ = np.ma.masked_where(qc != 0, ombg) - ombg_ = np.ma.masked_where(ombg == 0, ombg_) - oman_ = np.ma.masked_where(qc != 0, oman) - oman_ = np.ma.masked_where(ombg == 0, oman_) - ombg = ombg_ - oman = oman_ - lat = f.groups['MetaData'].variables['latitude'] - lon = f.groups['MetaData'].variables['longitude'] + if plottype=='histogram': + ombg_=np.ma.masked_where(qc != 0, ombg) + ombg_=np.ma.masked_where(ombg == 0, ombg_) + oman_=np.ma.masked_where(qc != 0, oman) + oman_=np.ma.masked_where(ombg == 0, oman_) + ombg=ombg_ + oman=oman_ + lat=f.groups['MetaData'].variables['latitude'] + lon=f.groups['MetaData'].variables['longitude'] obs_.append(obs[:]) omb_.append(ombg[:]) @@ -53,88 +54,95 @@ def get_obs_stats(fdir, plottype): lat_.append(lat[:]) lon_.append(lon[:]) - total_omb = np.concatenate(omb_) - total_oma = np.concatenate(oma_) - total_obs = np.concatenate(obs_) - total_lat = np.concatenate(lat_) - total_lon = np.concatenate(lon_) + total_omb=np.concatenate(omb_) + total_oma=np.concatenate(oma_) + total_obs=np.concatenate(obs_) + total_lat=np.concatenate(lat_) + total_lon=np.concatenate(lon_) return total_oma,total_omb,total_lat,total_lon def plot_scatter(): print("===== PLOT: SCATTER =====") - field_mean = float("{:.2f}".format(np.mean(np.absolute(field)))) - field_std = float("{:.2f}".format(np.std(np.absolute(field)))) - field_max = float("{:.2f}".format(np.max(np.absolute(field)))) - field_min = float("{:.2f}".format(np.min(np.absolute(field)))) + if yaml_data['machine']=='hera': + cartopy.config['data_dir']='/scratch2/NCEPDEV/fv3-cam/Chan-hoo.Jeon/tools/NaturalEarth' + elif yaml_data['machine']=='orion': + cartopy.config['data_dir']='/home/chjeon/tools/NaturalEarth' + + field_mean=float("{:.2f}".format(np.mean(np.absolute(field)))) + field_std=float("{:.2f}".format(np.std(np.absolute(field)))) + field_max=float("{:.2f}".format(np.max(np.absolute(field)))) + field_min=float("{:.2f}".format(np.min(np.absolute(field)))) print("Mean |OMA|=",field_mean) print("STDV |OMA|=",field_std) print("Max |OMA|=",field_max) print("Min |OMA|=",field_min) - crs = ccrs.PlateCarree() - fig = plt.figure(figsize=(8,5)) - ax = plt.subplot(111, projection=crs) - ax.coastlines(resolution='110m') - norm = plt.Normalize(yaml_data['field_range'][0],yaml_data['field_range'][1]) - num_cmap = 25 - cmap_neg = mpl.colormaps['Blues_r'].resampled(num_cmap) - cmap_pos = mpl.colormaps['Reds'].resampled(num_cmap) - cmap_color = np.vstack((cmap_neg(np.linspace(0.1,0.7,num_cmap)),cmap_pos(np.linspace(0.2,0.8,num_cmap)))) - cmap_new = ListedColormap(cmap_color, name='BlueRed_rw') - sc = ax.scatter(lon, lat, c=field, s=1.5, cmap=cmap_new, transform=crs, norm=norm) - cbar = plt.colorbar(sc, orientation="horizontal", shrink=0.5, pad=0.05) - stitle = yaml_data['title_fig']+' \n '+'Mean |OMA| ='+str(field_mean)+', STDV |OMA| ='+str(field_std) + crs=ccrs.PlateCarree() + fig=plt.figure(figsize=(8,5)) + ax=plt.subplot(111, projection=crs) + coastline=cfeature.NaturalEarthFeature('physical','coastline','50m',edgecolor='black',facecolor='none', + linewidth=0.5,alpha=0.7) + ax.add_feature(coastline) + norm=plt.Normalize(yaml_data['field_range'][0],yaml_data['field_range'][1]) + num_cmap=25 + cmap_neg=mpl.colormaps['Blues_r'].resampled(num_cmap) + cmap_pos=mpl.colormaps['Reds'].resampled(num_cmap) + cmap_color=np.vstack((cmap_neg(np.linspace(0.1,0.7,num_cmap)),cmap_pos(np.linspace(0.2,0.8,num_cmap)))) + cmap_new=ListedColormap(cmap_color, name='BlueRed_rw') + sc=ax.scatter(lon, lat, c=field, s=1.5, cmap=cmap_new, transform=crs, norm=norm) + cbar=plt.colorbar(sc, orientation="horizontal", shrink=0.5, pad=0.05) + stitle=yaml_data['title_fig']+' \n '+'Mean |OMA| ='+str(field_mean)+', STDV |OMA| ='+str(field_std) plt.title(stitle) - output_fn = yaml_data['output_prefix']+"_scatter.png" + output_fn=yaml_data['output_prefix']+"_scatter.png" plt.savefig(output_fn,dpi=200,bbox_inches='tight') plt.close('all') def plot_histogram(): print("===== PLOT: HISTOGRAM =====") - field_mean = float("{:.2f}".format(np.mean(field))) - field_std = float("{:.2f}".format(np.std(field))) - field_max = float("{:.2f}".format(np.max(field))) - field_min = float("{:.2f}".format(np.min(field))) + field_mean=float("{:.2f}".format(np.mean(field))) + field_std=float("{:.2f}".format(np.std(field))) + field_max=float("{:.2f}".format(np.max(field))) + field_min=float("{:.2f}".format(np.min(field))) print("Mean OMA=",field_mean) print("STDV OMA=",field_std) print("Max OMA=",field_max) print("Min OMA=",field_min) - nbins = yaml_data['nbins'] + nbins=yaml_data['nbins'] - opt_xlimit = 'auto' - if opt_xlimit == 'auto': - fld_min = int(field_min) - fld_max = int(field_max) - xlimit = [fld_min,fld_max] + opt_xlimit='auto' + if opt_xlimit=='auto': + fld_min=int(field_min) + fld_max=int(field_max) + xlimit=[fld_min,fld_max] print("xlimit min=",fld_min) print("xlimit max=",fld_max) print("xlimit=",xlimit) else: - xlimit = yaml_data['field_range'] + xlimit=yaml_data['field_range'] plt.hist(field[:], bins=nbins, range=xlimit, density=True, color ="blue") - stitle = yaml_data['title_fig']+' \n '+'Mean(OMA) ='+str(field_mean)+', STDV(OMA) ='+str(field_std) + stitle=yaml_data['title_fig']+' \n '+'Mean(OMA) ='+str(field_mean)+', STDV(OMA) ='+str(field_std) plt.title(stitle) - output_fn = yaml_data['output_prefix']+"_histogram.png" + output_fn=yaml_data['output_prefix']+"_histogram.png" plt.savefig(output_fn,dpi=150,bbox_inches='tight') plt.close('all') if __name__ == '__main__': global field,yaml_data - yaml_file = "plot_hofx.yaml" + yaml_file="plot_hofx.yaml" with open(yaml_file, 'r') as f: - yaml_data = yaml.load(f, Loader=yaml.FullLoader) + yaml_data=yaml.load(f, Loader=yaml.FullLoader) f.close() print("YAML_DATA:",yaml_data) - oma,omb,lat,lon= get_obs_stats(yaml_data['hofx_files'],yaml_data['plottype']) - - if yaml_data['field_var'] == 'OMA': - field = oma - if yaml_data['field_var'] == 'OMB': - field = omb - if yaml_data['plottype'] == 'scatter' or yaml_data['plottype'] == 'both': + oma,omb,lat,lon=get_obs_stats(yaml_data['hofx_files'],yaml_data['plottype']) + if yaml_data['field_var']=='OMA': + field=oma + if yaml_data['field_var']=='OMB': + field=omb + + if yaml_data['plottype']=='scatter' or yaml_data['plottype']=='both': plot_scatter() - if yaml_data['plottype'] == 'histogram' or yaml_data['plottype'] == 'both': + if yaml_data['plottype']=='histogram' or yaml_data['plottype']=='both': plot_histogram()