pro plot_cfradial ; Path prefix path_prefix='/Volumes/' ;path_prefix='/uufs/chpc.utah.edu/common/home/' ; Hid - fine grid hdir='mace-group5/gpm/hid/update_20210728/HID_SS_20180209_12-00z/' hfile='OPOL-SS_rnkHID_F_20180209_140448.lat_lon.nc' ; Cfradial cdir='mace-group4/Capricorn2/CBand/cfradial/' cfile='9776HUB-PPIVol-20180209-140504-0000.cfradial.pyart.81.401.nc' ; raw rdir='mace-group4/Capricorn2/CBand/cfradial/' file='9776HUB-PPIVol-20180209-140504-0000.cfradial.nc' ; Get the date and time from the filename parts=strsplit(hfile,'_',/extract) hdate=parts[3] htime=parts[4] hhour=strmid(htime,0,2) hmin=strmid(htime,2,2) hsec=strmid(htime,4,2) ; Read the hid file fid=ncdf_open(path_prefix+hdir+hfile) vid=ncdf_varid(fid,'time') & ncdf_varget,fid,vid,time_hid xid=ncdf_varid(fid,'x') & ncdf_varget,fid,xid,x_m_hid xid=ncdf_varid(fid,'y') & ncdf_varget,fid,xid,y_m_hid xid=ncdf_varid(fid,'z') & ncdf_varget,fid,xid,z_m_hid xid=ncdf_varid(fid,'origin_latitude') & ncdf_varget,fid,xid,olat_hid xid=ncdf_varid(fid,'origin_longitude') & ncdf_varget,fid,xid,olon_hid xid=ncdf_varid(fid,'origin_altitude') & ncdf_varget,fid,xid,oalt_hid xid=ncdf_varid(fid,'point_latitude') & if xid ne -1 then ncdf_varget,fid,xid,lat_3d_hid xid=ncdf_varid(fid,'point_longitude') & if xid ne -1 then ncdf_varget,fid,xid,lon_3d_hid xid=ncdf_varid(fid,'point_altitude') & if xid ne -1 then ncdf_varget,fid,xid,alt_3d_hid xid=ncdf_varid(fid,'corrected_specific_differential_phase') & ncdf_varget,fid,xid,kdp_hid xid=ncdf_varid(fid,'cross_correlation_ratio') & ncdf_varget,fid,xid,rhohv_hid xid=ncdf_varid(fid,'log_reflectivity') & ncdf_varget,fid,xid,dbz_hid xid=ncdf_varid(fid,'log_zdr') & ncdf_varget,fid,xid,zdr_hid xid=ncdf_varid(fid,'HID_class_01') & ncdf_varget,fid,xid,hid xid=ncdf_varid(fid,'hid_score_DZ') & ncdf_varget,fid,xid,hid_dz ;drizzle xid=ncdf_varid(fid,'hid_score_RN') & ncdf_varget,fid,xid,hid_rn ;rain xid=ncdf_varid(fid,'hid_score_CR') & ncdf_varget,fid,xid,hid_cr ;ice crystals xid=ncdf_varid(fid,'hid_score_AG') & ncdf_varget,fid,xid,hid_ag ;aggregates xid=ncdf_varid(fid,'hid_score_WS') & ncdf_varget,fid,xid,hid_ws ;wet snow xid=ncdf_varid(fid,'hid_score_VI') & ncdf_varget,fid,xid,hid_vi ;vertical ice xid=ncdf_varid(fid,'hid_score_GL') & ncdf_varget,fid,xid,hid_gl ;graupel xid=ncdf_varid(fid,'hid_score_GH') & ncdf_varget,fid,xid,hid_gh ;graupel HD xid=ncdf_varid(fid,'hid_score_HA') & ncdf_varget,fid,xid,hid_ha ;hail xid=ncdf_varid(fid,'hid_score_HA') & ncdf_varget,fid,xid,hid_ha ;hail xid=ncdf_varid(fid,'hid_score_BD') & ncdf_varget,fid,xid,hid_bd ;big drops/melting hail ncdf_close,fid ; Read the pyart cfradial file fid=ncdf_open(path_prefix+cdir+cfile) vid=ncdf_varid(fid,'base_time') & ncdf_varget,fid,vid,base_time_cfr vid=ncdf_varid(fid,'time_offset') & ncdf_varget,fid,vid,time_offset_cfr xid=ncdf_varid(fid,'x') & ncdf_varget,fid,xid,x_m_cfr xid=ncdf_varid(fid,'y') & ncdf_varget,fid,xid,y_m_cfr xid=ncdf_varid(fid,'z') & ncdf_varget,fid,xid,z_m_cfr vid=ncdf_varid(fid,'point_latitude') & ncdf_varget,fid,vid,lat_3d_cfr vid=ncdf_varid(fid,'point_longitude') & ncdf_varget,fid,vid,lon_3d_cfr vid=ncdf_varid(fid,'point_altitude') & ncdf_varget,fid,vid,alt_3d_cfr vid=ncdf_varid(fid,'corrected_specific_differential_phase') & ncdf_varget,fid,vid,kdp_cfr vid=ncdf_varid(fid,'corrected_reflectivity') & ncdf_varget,fid,vid,dbz_cfr vid=ncdf_varid(fid,'cross_correlation_ratio') & ncdf_varget,fid,vid,rhohv_cfr vid=ncdf_varid(fid,'temperature') & ncdf_varget,fid,vid,temp_cfr vid=ncdf_varid(fid,'radar_echo_classification') & ncdf_varget,fid,vid,hid_cfr vid=ncdf_varid(fid,'corrected_differential_reflectivity') & ncdf_varget,fid,vid,zdr_cfr ncdf_close,fid ; Convert distances to km x_km_hid=x_m_hid/1000.0 y_km_hid=y_m_hid/1000.0 z_km_hid=z_m_hid/1000.0 x_km_cfr=x_m_cfr/1000.0 y_km_cfr=y_m_cfr/1000.0 z_km_cfr=z_m_cfr/1000.0 ; Choose the 2km height level r=where(z_km_hid le 2.0,c) idx_hid=r[-1] r=where(z_km_cfr le 2.0,c) idx_cfr=r[-1] ; Subset it to be the same region as the hid file (from Stephanie) r=where(x_km_hid lt 0 and x_km_hid ge min(x_km_hid),c) lat_3d_hid=lat_3d_hid[r,*,*] lon_3d_hid=lon_3d_hid[r,*,*] kdp_hid=kdp_hid[r,*,*] ;hid=hid[r,*,*] ;rhohv=rhohv[r,*,*] dbz_hid=dbz_hid[r,*,*] ;zdr=zdr[r,*,*] x_km_hid=x_km_hid[r] r=where(x_km_cfr lt 0 and x_km_cfr ge min(x_km_hid),c) lat_3d_cfr=lat_3d_cfr[r,*,*] lon_3d_cfr=lon_3d_cfr[r,*,*] kdp_cfr=kdp_cfr[r,*,*] hid_cfr=hid_cfr[r,*,*] rhohv_cfr=rhohv_cfr[r,*,*] dbz_cfr=dbz_cfr[r,*,*] zdr_cfr=zdr_cfr[r,*,*] temp_cfr=temp_cfr[r,*,*] x_km_cfr=x_km_cfr[r] r=where(y_km_cfr lt max(y_km_hid) and y_km_cfr ge min(y_km_hid),c) lat_3d_cfr=lat_3d_cfr[*,r,*] lon_3d_cfr=lon_3d_cfr[*,r,*] kdp_cfr=kdp_cfr[*,r,*] hid_cfr=hid_cfr[*,r,*] rhohv_cfr=rhohv_cfr[*,r,*] dbz_cfr=dbz_cfr[*,r,*] zdr_cfr=zdr_cfr[*,r,*] temp_cfr=temp_cfr[*,r,*] y_km_cfr=y_km_cfr[r] ; Find the lat,lon range of the plot area ulat_hid=max(lat_3d_hid) llat_hid=min(lat_3d_hid) llon_hid=min(lon_3d_hid) rlon_hid=max(lon_3d_hid) if llon_hid lt 0 then lcolon_hid=360.0+llon_hid else lcolon_hid=llon_hid if rlon_hid lt 0 then rcolon_hid=360.0+rlon_hid else rcolon_hid=rlon_hid ;************************* ; Temp image ;************************* ; Set up plot size pxdim=1000 & pydim=900 ; Position the plots xl=0.05 & xr=0.90 yb=0.05 & yt=0.85 sx=0.05 sy=0.13 numplots_x=4 numplots_y=2 position_plots,xl,xr,yb,yt,sx,sy,numplots_x,numplots_y,pos ; Colorbar position cbpos=pos ;cbpos[*,0]=pos[*,2]+0.02 ;vertical ;cbpos[*,2]=cbpos[*,0]+0.01 cbpos[*,1]=pos[*,3]+0.07 ;horizontal cbpos[*,3]=cbpos[*,1]+0.01 dx=0.01 & dy=0.01 ;Clear plot space pnum=0 p0=plot([0,1],[0,1],position=pos[pnum,*],/buffer,dimensions=[pxdim,pydim],axis_style=4,/nodata) d=p0.convertcoord([pos[pnum,0],pos[pnum,2]],[pos[pnum,1],pos[pnum,3]],/normal,/to_device) isx=(d[0,1,0]-d[0,0,0]) isy=(d[1,1,0]-d[1,0,0]) ; Top is the last color to scale 256 colors, 0-255 top_color=252 ; Colortable 0-252 253=white mytable=colortable(39,ncolors=254) ;254=hot pink ;gray=255 mytable=[mytable,transpose([238,18,137]),transpose([230,230,230])] mycbtable=mytable[0:top_color,*] ; drizzle rain ice crystals aggregates wet snow hidct=['white smoke','sky blue','medium blue','orange', 'pink', 'aqua',$ ;vertical ice LD Graupel HD Graupel Hail Big Drops/Melting Hail 'dark gray', 'lime green', 'yellow', 'red','magenta'] ;*** PLOT KDP-hid if 1 eq 0 then begin pnum=0 data_var=reform(kdp_hid[*,*,idx_hid]) lat_var=reform(lat_3d_hid[*,*,idx_hid]) lon_var=reform(lon_3d_hid[*,*,idx_hid]) height_var=z_km_hid[idx_hid] dmax=max(data_var) r=where(data_var ne -9999,c) dmin=min(data_var[r]) data_image=bytscl(data_var,top=top_color,max=dmax,min=dmin) r=where(data_var eq -9999,c) if c gt 0 then data_image[r]=253 r=where(data_var eq 0,c) if c gt 0 then data_image[r]=255 p0=image(data_image,x_km_hid,y_km_hid,/current,$ dimensions=[pxdim,pydim],position=pos[pnum,*],$ font_size=fsm,rgb_table=mytable) c0=contour(data_var,x_km_hid,y_km_hid,/current,position=pos[pnum,*],/nodata,$ /xstyle,/ystyle,xtickdir=1,ytickdir=1) c0=colorbar(rgb_table=mycbtable,range=[dmin,dmax],orientation=1,$ position=reform(cbpos[pnum,*]),font_size=10,tickdir=1,border_on=1,title='deg/km') t1=text(pos[pnum,0]+0*dx,pos[pnum,3]+1*dy,$ 'KDP '+string(height_var,format='(F5.2)')+'km',font_size=fs1) t1=text(pos[pnum,0]+0*dx,pos[pnum,3]+3*dy,$ hfile,font_size=10) endif ;*** PLOT KDP-cfr pnum=0 data_var=reform(kdp_cfr[*,*,idx_cfr]) lat_var=reform(lat_3d_cfr[*,*,idx_cfr]) lon_var=reform(lon_3d_cfr[*,*,idx_cfr]) height_var=z_km_cfr[idx_cfr] dmax=max(data_var) r=where(data_var ne -9999,c) dmin=min(data_var[r]) data_image=bytscl(data_var,top=top_color,max=dmax,min=dmin) r=where(data_var eq -9999,c) if c gt 0 then data_image[r]=253 r=where(data_var eq 0,c) if c gt 0 then data_image[r]=255 p0=image(data_image,x_km_cfr,y_km_cfr,/current,$ image_dimensions=[isx,isy],$ dimensions=[pxdim,pydim],position=pos[pnum,*],$ font_size=fsm,rgb_table=mytable) c0=contour(data_var,x_km_cfr,y_km_cfr,/current,position=pos[pnum,*],/nodata,$ /xstyle,/ystyle,xtickdir=1,ytickdir=1) c0=colorbar(rgb_table=mycbtable,range=[dmin,dmax],orientation=0,$ position=reform(cbpos[pnum,*]),font_size=10,tickdir=1,border_on=1,title='deg/km') t1=text(pos[pnum,0]+0*dx,pos[pnum,3]+1*dy,$ 'KDP '+string(height_var,format='(F5.2)')+'km',font_size=fs1) ;t1=text(pos[pnum,0]+0*dx,pos[pnum,3]+3*dy,cfile,font_size=10) ;*** PLOT dbz-cfr pnum=1 data_var=reform(dbz_cfr[*,*,idx_cfr]) lat_var=reform(lat_3d_cfr[*,*,idx_cfr]) lon_var=reform(lon_3d_cfr[*,*,idx_cfr]) height_var=z_km_cfr[idx_cfr] dmax=max(data_var) r=where(data_var ne -9999,c) dmin=min(data_var[r]) data_image=bytscl(data_var,top=top_color,max=dmax,min=dmin) r=where(data_var eq -9999,c) if c gt 0 then data_image[r]=253 r=where(data_var eq 0,c) if c gt 0 then data_image[r]=255 p0=image(data_image,x_km_cfr,y_km_cfr,/current,$ image_dimensions=[isx,isy],$ dimensions=[pxdim,pydim],position=pos[pnum,*],$ font_size=fsm,rgb_table=mytable) c0=contour(data_var,x_km_cfr,y_km_cfr,/current,position=pos[pnum,*],/nodata,$ /xstyle,/ystyle,xtickdir=1,ytickdir=1) c0=colorbar(rgb_table=mycbtable,range=[dmin,dmax],orientation=0,$ position=reform(cbpos[pnum,*]),font_size=10,tickdir=1,border_on=1,title='dBZ') t1=text(pos[pnum,0]+0*dx,pos[pnum,3]+1*dy,$ 'REFLECT '+string(height_var,format='(F5.2)')+'km',font_size=fs1) ;t1=text(pos[pnum,0]+0*dx,pos[pnum,3]+3*dy,cfile,font_size=10) ;*** PLOT temp-cfr pnum=2 data_var=reform(temp_cfr[*,*,idx_cfr]) lat_var=reform(lat_3d_cfr[*,*,idx_cfr]) lon_var=reform(lon_3d_cfr[*,*,idx_cfr]) height_var=z_km_cfr[idx_cfr] r=where(data_var ne -9999 and data_var lt 1e6,c) dmax=max(data_var[r]) dmin=min(data_var[r]) data_image=bytscl(data_var,top=top_color,max=dmax,min=dmin) r=where(data_var eq -9999,c) if c gt 0 then data_image[r]=253 r=where(data_var eq 0,c) if c gt 0 then data_image[r]=255 p0=image(data_image,x_km_cfr,y_km_cfr,/current,$ image_dimensions=[isx,isy],$ dimensions=[pxdim,pydim],position=pos[pnum,*],$ font_size=fsm,rgb_table=mytable) c0=contour(data_var,x_km_cfr,y_km_cfr,/current,position=pos[pnum,*],/nodata,$ /xstyle,/ystyle,xtickdir=1,ytickdir=1) c0=colorbar(rgb_table=mycbtable,range=[dmin,dmax],orientation=0,$ position=reform(cbpos[pnum,*]),font_size=10,tickdir=1,border_on=0,title='C') t1=text(pos[pnum,0]+0*dx,pos[pnum,3]+1*dy,$ 'TEMP '+string(height_var,format='(F5.2)')+'km',font_size=fs1) ;t1=text(pos[pnum,0]+0*dx,pos[pnum,3]+3*dy,cfile,font_size=10) ;*** PLOT hid-cfr pnum=3 data_var=reform(hid_cfr[*,*,idx_cfr]) lat_var=reform(lat_3d_cfr[*,*,idx_cfr]) lon_var=reform(lon_3d_cfr[*,*,idx_cfr]) height_var=z_km_cfr[idx_cfr] r=where(data_var ne -9999 and data_var lt 1e6,c) dmax=max(data_var[r]) dmin=min(data_var[r]) data_image=bytscl(data_var,top=top_color,max=dmax,min=dmin) r=where(data_var eq -9999,c) if c gt 0 then data_image[r]=253 r=where(data_var eq 0,c) if c gt 0 then data_image[r]=255 p0=image(data_image,x_km_cfr,y_km_cfr,/current,$ image_dimensions=[isx,isy],$ dimensions=[pxdim,pydim],position=pos[pnum,*],$ font_size=fsm,rgb_table=mytable) c0=contour(data_var,x_km_cfr,y_km_cfr,/current,position=pos[pnum,*],/nodata,$ /xstyle,/ystyle,xtickdir=1,ytickdir=1) c0=colorbar(rgb_table=mycbtable,range=[dmin,dmax],orientation=0,$ position=reform(cbpos[pnum,*]),font_size=10,tickdir=1,border_on=0,title='K') t1=text(pos[pnum,0]+0*dx,pos[pnum,3]+1*dy,$ 'HID '+string(height_var,format='(F5.2)')+'km',font_size=fs1) ;t1=text(pos[pnum,0]+0*dx,pos[pnum,3]+3*dy,cfile,font_size=10) ;*** PLOT zdr-cfr pnum=4 data_var=reform(zdr_cfr[*,*,idx_cfr]) lat_var=reform(lat_3d_cfr[*,*,idx_cfr]) lon_var=reform(lon_3d_cfr[*,*,idx_cfr]) height_var=z_km_cfr[idx_cfr] r=where(data_var ne -9999 and data_var lt 1e6,c) dmax=max(data_var[r]) dmin=min(data_var[r]) data_image=bytscl(data_var,top=top_color,max=dmax,min=dmin) r=where(data_var eq -9999,c) if c gt 0 then data_image[r]=253 r=where(finite(data_var) eq 0 ,c) if c gt 0 then data_image[r]=255 p0=image(data_image,x_km_cfr,y_km_cfr,/current,$ image_dimensions=[isx,isy],$ dimensions=[pxdim,pydim],position=pos[pnum,*],$ font_size=fsm,rgb_table=mytable) c0=contour(data_var,x_km_cfr,y_km_cfr,/current,position=pos[pnum,*],/nodata,$ /xstyle,/ystyle,xtickdir=1,ytickdir=1) c0=colorbar(rgb_table=mycbtable,range=[dmin,dmax],orientation=0,$ position=reform(cbpos[pnum,*]),font_size=10,tickdir=1,border_on=0,title='dB') t1=text(pos[pnum,0]+0*dx,pos[pnum,3]+1*dy,$ 'ZDR '+string(height_var,format='(F5.2)')+'km',font_size=fs1) ;t1=text(pos[pnum,0]+0*dx,pos[pnum,3]+3*dy,cfile,font_size=10) t1=text(0.1,0.98,cfile,font_size=fs1) hstr=string(height_var,format='(F4.2)')+'km' p0.save,'cfrad'+hdate+'.'+htime+'.'+hstr+'.4panel.png',height=pydim stop end