###################################################################### CHJ ##### ## Name : plot_fv3lam_his2d_bndr.py ## Language : Python 3.7 ## Usage : Plot boundary of domain from fv3_history2d.nc (output) ## Input files : fv3_history2d.nc ## NOAA/NWS/NCEP/EMC ## History =============================== ## V000: 2020/12/24: Chan-Hoo Jeon : Preliminary version ###################################################################### CHJ ##### import os, sys import matplotlib matplotlib.use('Agg') import matplotlib.pyplot as plt import matplotlib.colors as colors import matplotlib.lines as mlines import numpy as np from netCDF4 import Dataset #import xarray as xr from mpl_toolkits.axes_grid1 import make_axes_locatable # HPC machine ('hera','orion') machine='hera' print(' You are on', machine) plt.switch_backend('agg') # Global variables ======================================== CHJ ===== # ..... Case-dependent input :: should be changed case-by-case ..... # ****** # INPUT # ****** # Path to the directory where the input NetCDF file is located. dnm_in="/scratch2/NCEPDEV/fv3-cam/Benjamin.Blake/SAR_paper/sar_2019080800/" # grid file name fnm_input='phyf000.nc' # Variables #vars_his=["ugrd","vgrd"] vars_his=["prateb_ave"] # number of the target vertical level (only for 3-D fields) #ilvl=4 # Time level of plotting ( 0 < prt_tlvl < maximum-1) prt_tlvl=0 # Length of boundary (unit: number of grid points) len_bndr=41 # tick interval for short sides dtick_s=4 # tick interval for long sides dtick_l=200 # Domain name: domain='HRRR' # Grid resolution ('C96'/'C768'): res='C768' # Grid type ('ESG'/'GFDL') gtype='GFDL' # GFDL grid-refinement ratio (for ESG grid, refine=0) if gtype=='ESG': refine=0 elif gtype=='GFDL': refine=3 # ******* # OUTPUT # ******* # Path to directory if machine=='hera': out_fig_dir="/scratch2/NCEPDEV/fv3-cam/Benjamin.Blake/SAR_paper/" elif machine=='orion': out_fig_dir="/work/noaa/fv3-cam/chjeon/tools/Fig/" else: sys.exit('ERROR: path to output directory is not set !!!') # basic forms of title and file name if gtype=='ESG': out_title_base='FV3::HIS2D::'+domain+'(ESG)::'+res+'::' out_fname_base='fv3lam_out_his2d_'+domain+'_esg_'+res+'_' elif gtype=='GFDL': out_title_base='PHY::'+domain+'(GFDL)::'+res+'(x'+str(refine)+')'+'::' out_fname_base='fv3lam_out_his2d_'+domain+'_gfdl_'+res+'_' # Machine-specific mfdataset arguments if machine=='hera': mfdt_kwargs={'parallel':False} elif machine=='orion': mfdt_kwargs={'parallel':False,'combine':'by_coords'} else: mfdt_kwargs={'parallel':False} # Main part (will be called at the end) ======================= CHJ ===== def main(): # ============================================================= CHJ ===== global his2d print(' ===== OUTPUT: history2d =======================================') # open the data file fname=os.path.join(dnm_in,fnm_input) try: his2d=Dataset(fname,'r') except: raise Exception('Could NOT find the file',fname) print(his2d) # Variables for svar in vars_his: his_plot(svar) # ===== plot ================================================== CHJ ===== def his_plot(svar): # ============================================================= CHJ ===== global out_title,out_fname,cs_min,cs_max,sfld2d global cs_cmap,lb_ext,tick_ln,tick_wd,tlb_sz,n_rnd global nts,nys,nxs print(' ===== '+svar+' ===== history2d ===============================') # Extract data array # prate=np.ma.masked_invalid(his2d[svar].data) # (nts,nys,nxs)=prate.shape # print(' time+2D: nts=',nts,' nys=',nys,' nxs=',nxs) # if prt_tlvl>=nts: # sys.exit('ERROR: prt_tlvl >= max. time level !!!') list = ['00','03','06','09','12'] for t in list: fnm_input='phyf0'+t+'.nc' fname=os.path.join(dnm_in,fnm_input) try: his2d=Dataset(fname,'r') except: raise Exception('Could NOT find the file',fname) prate=np.ma.masked_invalid(his2d.variables[svar][:]) (nts,nys,nxs)=prate.shape prate2d=prate[0,:,:] # Rotate the array 180 degrees prate2d=np.fliplr(prate2d) prate2d=np.flipud(prate2d) if t == '00': sfld2d = prate2d * 3 * 3600 * 0.0393701 else: sfld2d = sfld2d + (prate2d * 3 * 3600 * 0.0393701) out_title=out_title_base+svar+'::time level='+str(prt_tlvl) out_fname=out_fname_base+svar cs_cmap='viridis' lb_ext='max' tick_ln=1.5 tick_wd=0.45 tlb_sz=5 n_rnd=5 cmap_range='fixed' print(' Plotting field=',svar) # Max and Min of the field fmax=np.max(sfld2d) fmin=np.min(sfld2d) print(' fld_max=',fmax) print(' flx_min=',fmin) # Make the colormap range symmetry print(' cmap range=',cmap_range) if cmap_range=='symmetry': tmp_cmp=max(abs(fmax),abs(fmin)) cs_min=round(-tmp_cmp,n_rnd) cs_max=round(tmp_cmp,n_rnd) elif cmap_range=='round': cs_min=round(fmin,n_rnd) cs_max=round(fmax,n_rnd) elif cmap_range=='real': cs_min=fmin cs_max=fmax elif cmap_range=='fixed': cs_min=0.0 cs_max=10.0 else: sys.exit('ERROR: wrong colormap-range flag !!!') print(' cs_max=',cs_max) print(' cs_min=',cs_min) print(' Plotting the entire domain ... ') plot_domain(svar) print(' Plotting boundaries ... ') plot_bndry(svar) # Plot: Boundary ========================================== CHJ ===== def plot_bndry(svar): # ========================================================= CHJ ===== # Plot field fig=plt.figure(figsize=(4,4)) # fig=plt.figure(figsize=(5,4)) grid=plt.GridSpec(4,4,wspace=0.1,hspace=0.1) # grid=plt.GridSpec(4,5,wspace=0.1,hspace=0.1) # center # ax_c=fig.add_subplot(grid[1:-1,1:-2]) ax_c=fig.add_subplot(grid[1:-1,1:-1]) # left ax_l=fig.add_subplot(grid[1:-1,0]) # right ax_r=fig.add_subplot(grid[1:-1,3]) # colorbar # ax_cb=fig.add_axes([0.8,0.2,0.02,0.6]) # top # ax_t=fig.add_subplot(grid[0,0:4]) ax_t=fig.add_subplot(grid[0,:]) # bottom # ax_b=fig.add_subplot(grid[3,0:4]) ax_b=fig.add_subplot(grid[3,:]) tick_ln=1.5 tick_wd=0.45 tlb_sz=4.5 len_bndr_m1=len_bndr-1 # tick: bottom tick_b_x=np.arange(0,nxs,dtick_l) tick_lbl_b_x=np.arange(0,nxs+1,dtick_l) tick_b_y=np.arange(0,len_bndr,dtick_s) tick_lbl_b_y=np.arange(0,len_bndr+1,dtick_s) # tick: top tick_t_x=np.arange(0,nxs,dtick_l) tick_lbl_t_x=np.arange(0,nxs+1,dtick_l) tick_t_y=np.arange(0,len_bndr,dtick_s) tick_lbl_t_y=np.arange(nys-len_bndr_m1,nys+1,dtick_s) # tick: left tick_l_x=np.arange(0,len_bndr,dtick_s) tick_lbl_l_x=np.arange(0,len_bndr+1,dtick_s) tick_l_y=np.arange(0,nys,dtick_l) tick_lbl_l_y=np.arange(0,nys+1,dtick_l) # tick: right tick_r_x=np.arange(0,len_bndr,dtick_s) tick_lbl_r_x=np.arange(nxs-len_bndr_m1,nxs+1,dtick_s) tick_r_y=np.arange(0,nys,dtick_l) tick_lbl_r_y=np.arange(0,nys+1,dtick_l) # Use standard colormap for precip clevs = [0.01,0.1,0.25,0.5,0.75,1,1.25,1.5,1.75,2,2.5,3,4,5,6] colorlist = ['chartreuse','limegreen','green','blue','dodgerblue','deepskyblue','cyan','mediumpurple','mediumorchid','darkmagenta','darkred','crimson','orangered','goldenrod'] cs_cmap = matplotlib.colors.ListedColormap(colorlist) cs_cmap.set_under('white',alpha=0.) cs_cmap.set_over('gold') norm = matplotlib.colors.BoundaryNorm(clevs, cs_cmap.N) # bottom (South) ax_b.pcolormesh(sfld2d[:len_bndr-1,:],cmap=cs_cmap,rasterized=True,vmin=0.01,norm=norm) # ax_b.pcolormesh(sfld2d[:len_bndr-1,:],cmap=cs_cmap,rasterized=True,vmin=cs_min,vmax=cs_max) ax_b.tick_params(direction='out',length=tick_ln,width=tick_wd,labelsize=tlb_sz) ax_b.set_xticks(tick_b_x) ax_b.set_xticklabels(tick_lbl_b_x) ax_b.set_yticks(tick_b_y) ax_b.set_yticklabels(tick_lbl_b_y) # top (North) ax_t.pcolormesh(sfld2d[-len_bndr_m1:,:],cmap=cs_cmap,rasterized=True,vmin=0.01,norm=norm) # ax_t.pcolormesh(sfld2d[-len_bndr_m1:,:],cmap=cs_cmap,rasterized=True,vmin=cs_min,vmax=cs_max) ax_t.tick_params(direction='out',length=tick_ln,width=tick_wd,labelsize=tlb_sz, bottom=False,labelbottom=False,top=True,labeltop=False) ax_t.set_xticks(tick_t_x) ax_t.set_xticklabels(tick_lbl_t_x) ax_t.set_yticks(tick_t_y) ax_t.set_yticklabels(tick_lbl_t_y) # left (West) ax_l.pcolormesh(sfld2d[len_bndr:-len_bndr,:len_bndr-1],cmap=cs_cmap,rasterized=True,vmin=0.01,norm=norm) # ax_l.pcolormesh(sfld2d[len_bndr:-len_bndr,:len_bndr-1],cmap=cs_cmap,rasterized=True,vmin=cs_min,vmax=cs_max) ax_l.tick_params(direction='out',length=tick_ln,width=tick_wd,labelsize=tlb_sz, labelbottom=False) ax_l.set_xticks(tick_l_x) ax_l.set_xticklabels(tick_lbl_l_x) ax_l.set_yticks(tick_l_y) ax_l.set_yticklabels(tick_lbl_l_y) # right (East) cs_r = ax_r.pcolormesh(sfld2d[len_bndr:-len_bndr,-len_bndr_m1:],cmap=cs_cmap,rasterized=True,vmin=0.01,norm=norm) # ax_r.pcolormesh(sfld2d[len_bndr:-len_bndr,-len_bndr_m1:],cmap=cs_cmap,rasterized=True,vmin=cs_min,vmax=cs_max) ax_r.tick_params(direction='out',length=tick_ln,width=tick_wd, left=False,labelleft=False,labelbottom=False,right=True) ax_r.set_xticks(tick_r_x) ax_r.set_xticklabels(tick_lbl_r_x) ax_r.set_yticks(tick_r_y) ax_r.set_yticklabels(tick_lbl_r_y) # center ax_c.axis([0,10,0,10]) ax_c.axis('off') ax_c.text(5,0.05,'South',fontsize=tlb_sz,horizontalalignment='center', verticalalignment='center') ax_c.text(5,9.95,'North',fontsize=tlb_sz,horizontalalignment='center', verticalalignment='center') ax_c.text(0.05,5,'West',fontsize=tlb_sz,horizontalalignment='center', verticalalignment='center',rotation='vertical') ax_c.text(9.95,5,'East',fontsize=tlb_sz,horizontalalignment='center', verticalalignment='center',rotation='vertical') ax_c.text(5,7,'d) FV3 LAM',fontsize=tlb_sz+1.5,color='black', horizontalalignment='center',verticalalignment='center') ax_c.text(5,6,'accumulated precipitation (in)',fontsize=tlb_sz+1.5,color='black', horizontalalignment='center',verticalalignment='center') ax_c.text(5,5,'Forecast hour 12',fontsize=tlb_sz+1.5,color='black', horizontalalignment='center',verticalalignment='center') # ax_c.text(5,4,'Boundary cells='+str(len_bndr),fontsize=tlb_sz+1.5,color='black', # horizontalalignment='center',verticalalignment='center') # cbar=plt.colorbar(cs_r,cax=ax_cb,extend=lb_ext,pad=0.1,ticks=clevs) # cs_r.cmap.set_under('white',alpha=0.) # cs_r.cmap.set_over('gold') # cbar.ax.tick_params(labelsize=tlb_sz) # cbar.ax.set_yticklabels([0.01,0.1,0.25,0.5,0.75,1,1.25,1.5,1.75,2,2.5,3,4,5,6]) # cbar.set_label('in',fontsize=tlb_sz) out_fname_bndr=out_fname+'_bndr' # Output figure ndpi=300 out_file(out_fname_bndr,ndpi) # Plot: Entire domain ===================================== CHJ ===== def plot_domain(svar): # ========================================================= CHJ ===== # tick_x=np.arange(0,nxs,dtick_l) # tick_lbl_x=np.arange(1,nxs+1,dtick_l) # tick_y=np.arange(0,nys,dtick_l) # tick_lbl_y=np.arange(1,nys+1,dtick_l) tick_x=np.arange(140,221,20) tick_lbl_x=np.arange(140,221,20) tick_y=np.arange(0,41,5) tick_lbl_y=np.arange(0,41,5) fig,ax=plt.subplots(1,1,figsize=(3.5,2.85)) ax.set_title('a) FV3 LAM accumulated precipitation (in) - Forecast hour 12', fontsize=tlb_sz+1) # Use standard colormap for precip clevs = [0.01,0.1,0.25,0.5,0.75,1,1.25,1.5,1.75,2,2.5,3,4,5,6] colorlist = ['chartreuse','limegreen','green','blue','dodgerblue','deepskyblue','cyan','mediumpurple','mediumorchid','darkmagenta','darkred','crimson','orangered','goldenrod'] cs_cmap = matplotlib.colors.ListedColormap(colorlist) norm = matplotlib.colors.BoundaryNorm(clevs, cs_cmap.N) cs=ax.pcolormesh(sfld2d,cmap=cs_cmap,rasterized=True,vmin=0.01,norm=norm) # cs=ax.pcolormesh(sfld2d,cmap=cs_cmap,rasterized=True,vmin=cs_min,vmax=cs_max) ax.tick_params(direction='out',length=tick_ln,width=tick_wd,labelsize=tlb_sz) ax.set_xlim([130,230]) ax.set_ylim([0,40]) ax.set_xticks(tick_x) ax.set_xticklabels(tick_lbl_x) ax.set_yticks(tick_y) ax.set_yticklabels(tick_lbl_y) # divider=make_axes_locatable(ax) # ax_cb=divider.new_horizontal(size="3%",pad=0.1,axes_class=plt.Axes) # fig.add_axes(ax_cb) # cbar=plt.colorbar(cs,cax=ax_cb,extend=lb_ext,ticks=clevs) # cs.cmap.set_under('white',alpha=0.) # cs.cmap.set_over('gold') # cbar.ax.tick_params(labelsize=tlb_sz) # cbar.ax.set_yticklabels([0.01,0.1,0.25,0.5,0.75,1,1.25,1.5,1.75,2,2.5,3,4,5,6]) # cbar.set_label('in',fontsize=tlb_sz) out_fname_domain=out_fname+'_domain' # Output figure ndpi=300 out_file(out_fname_domain,ndpi) # Output file ============================================= CHJ ===== def out_file(out_fname,ndpi): # ========================================================= CHJ ===== # Output figure plt.savefig(out_fig_dir+out_fname+'.png',dpi=ndpi,bbox_inches='tight') plt.close('all') # Main call ================================================ CHJ ===== if __name__=='__main__': main()