;Sally, why don’t you just resample the C-Band to match +14 dBZ sensitivity of GPM ;for this case ; Minimum dbz of 14 for gpm ; 1d histograms only histogram for cband ge 14 ; *** dates that have data in both gpm and cband pro match_cband_gpm_volumes ;date_str='20180120' date_str='20180122' ;*** ;date_str='20180124' ;date_str='20180126' ;*** ;date_str='20180128' ;date_str='20180203' ;date_str='20180204.12' ;date_str='20180204.13' ;date_str='20180205' ;date_str='20180206' ;*** ;date_str='20180209' ;*** ;date_str='20180210' ;date_str='20180215' ;date_str='20180217' ;*** ;date_str='20180220' ; Path prefix ;path_prefix='/Volumes/' path_prefix='/uufs/chpc.utah.edu/common/home/' ; Hid - fine grid - not using anymore ;hdir='mace-group5/gpm/hid/update_20210728/HID_SS_20180209_12-00z/' ;hfile='OPOL-SS_rnkHID_F_20180209_140448.lat_lon.nc' ; Cfradial pyart cartesian interpolation cdir='mace-group5/field_programs/investigator/in2018_v01/cband/'+$ strmid(date_str,0,8)+'/' if date_str eq '20180120' then $ ;cfile='9776HUB-PPIVol-20180120-140234-0000.cfradial.pyart.81.401.nc' cfile='9776HUB-PPIVol-20180120-140910-0000.cfradial.pyart.81.401.nc' if date_str eq '20180122' then $ cfile='9776HUB-PPIVol-20180122-135650-0000.cfradial.pyart.81.401.nc' ;cfile='9776HUB-PPIVol-20180122-140326-0000.cfradial.pyart.81.401.nc' if date_str eq '20180124' then $ ;cfile='9776HUB-PPIVol-20180124-134507-0000.cfradial.pyart.81.401.nc' cfile='9776HUB-PPIVol-20180124-135144-0000.cfradial.pyart.81.401.nc' if date_str eq '20180126' then $ cfile='9776HUB-PPIVol-20180126-133919-0000.cfradial.pyart.81.401.nc' ;cfile='9776HUB-PPIVol-20180126-134555-0000.cfradial.pyart.81.401.nc' if date_str eq '20180128' then $ ;cfile='9776HUB-PPIVol-20180128-132659-0000.cfradial.pyart.81.401.nc' ;cfile='9776HUB-PPIVol-20180128-133335-0000.cfradial.pyart.81.401.nc' cfile='9776HUB-PPIVol-20180128-134011-0000.cfradial.pyart.81.401.nc' if date_str eq '20180203' then $ cfile='9776HUB-PPIVol-20180203-130325-0000.cfradial.pyart.81.401.nc' ;cfile='9776HUB-PPIVol-20180203-131002-0000.cfradial.pyart.81.401.nc' if date_str eq '20180204.12' then $ cfile='9776HUB-PPIVol-20180204-121054-0000.cfradial.pyart.81.401.nc' ;cfile='9776HUB-PPIVol-20180204-121730-0000.cfradial.pyart.81.401.nc' if date_str eq '20180204.13' then $ cfile='9776HUB-PPIVol-20180204-134959-0000.cfradial.pyart.81.401.nc' ;cfile='9776HUB-PPIVol-20180204-135639-0000.cfradial.pyart.81.401.nc' if date_str eq '20180205' then $ cfile='9776HUB-PPIVol-20180205-111831-0000.cfradial.pyart.81.401.nc' ;cfile='9776HUB-PPIVol-20180205-112507-0000.cfradial.pyart.81.401.nc' if date_str eq '20180206' then $ cfile='9776HUB-PPIVol-20180206-102552-0000.cfradial.pyart.81.401.nc' ;cfile='9776HUB-PPIVol-20180206-103228-0000.cfradial.pyart.81.401.nc' if date_str eq '20180209' then $ cfile='9776HUB-PPIVol-20180209-140504-0000.cfradial.pyart.81.401.nc' ;cfile='9776HUB-PPIVol-20180209-141140-0000.cfradial.pyart.81.401.nc' if date_str eq '20180210' then $ cfile='9776HUB-PPIVol-20180210-100732-0000.cfradial.pyart.81.401.nc' ;cfile='9776HUB-PPIVol-20180210-101409-0000.cfradial.pyart.81.401.nc' if date_str eq '20180215' then $ cfile='9776HUB-PPIVol-20180215-134057-0000.cfradial.pyart.81.401.nc' ;cfile='9776HUB-PPIVol-20180215-134733-0000.cfradial.pyart.81.401.nc' if date_str eq '20180217' then $ ;******* cfile='9776HUB-PPIVol-20180217-071136-0000.cfradial.pyart.81.401.nc' ;cfile='9776HUB-PPIVol-20180217-071811-0000.cfradial.pyart.81.401.nc' if date_str eq '20180220' then $ ;cfile='9776HUB-PPIVol-20180220-123005-0000.cfradial.pyart.81.401.nc' cfile='9776HUB-PPIVol-20180220-123644-0000.cfradial.pyart.81.401.nc' ; raw cfradial file - not used ;rdir='mace-group4/Capricorn2/CBand/cfradial/' ;rfile='9776HUB-PPIVol-20180209-140504-0000.cfradial.nc' ; gpm gdir='mace-group5/gpm/GPM_L2/GPM_2AKu.06/' ;if date_str eq '20180111' then $ ;no cband ; gfile='2A.GPM.Ku.V8-20180723.20180111-S002115-E015349.021988.V06A.SUB.HDF5' ;if date_str eq '20180113' then $ ;no cband ; gfile='2A.GPM.Ku.V8-20180723.20180113-S140450-E153724.022028.V06A.SUB.HDF5' if date_str eq '20180120' then $ gfile='2A.GPM.Ku.V8-20180723.20180120-S124419-E141653.022136.V06A.SUB.HDF5' if date_str eq '20180122' then $ gfile='2A.GPM.Ku.V8-20180723.20180122-S123430-E140704.022167.V06A.SUB.HDF5' if date_str eq '20180124' then $ gfile='2A.GPM.Ku.V8-20180723.20180124-S122439-E135713.022198.V06A.SUB.HDF5' if date_str eq '20180126' then $ gfile='2A.GPM.Ku.V8-20180723.20180126-S121448-E134722.022229.V06A.SUB.HDF5' if date_str eq '20180128' then $ gfile='2A.GPM.Ku.V8-20180723.20180128-S120455-E133730.022260.V06A.SUB.HDF5' if date_str eq '20180203' then $ gfile='2A.GPM.Ku.V8-20180723.20180203-S113514-E130748.022353.V06A.SUB.HDF5' if date_str eq '20180204.12' then $ gfile='2A.GPM.Ku.V8-20180723.20180204-S104359-E121633.022368.V06A.SUB.HDF5' if date_str eq '20180204.13' then $ gfile='2A.GPM.Ku.V8-20180723.20180204-S134909-E152143.022370.V06A.SUB.HDF5' if date_str eq '20180205' then $ gfile='2A.GPM.Ku.V8-20180723.20180205-S095244-E112518.022383.V06A.SUB.HDF5' if date_str eq '20180206' then $ ;******* gfile='2A.GPM.Ku.V8-20180723.20180206-S090128-E103402.022398.V06A.HDF5' if date_str eq '20180209' then $ ;******* gfile='2A.GPM.Ku.V8-20180723.20180209-S141035-E154309.022448.V06A.HDF5' if date_str eq '20180210' then $ ;******* gfile='2A.GPM.Ku.V8-20180723.20180210-S084134-E101408.022460.V06A.SUB.HDF5' if date_str eq '20180215' then $ gfile='2A.GPM.Ku.V8-20180723.20180215-S134018-E151250.022541.V06A.SUB.HDF5' if date_str eq '20180217' then $ ;******* gfile='2A.GPM.Ku.V8-20180723.20180217-S054649-E071922.022567.V06A.SUB.HDF5' if date_str eq '20180220' then $ gfile='2A.GPM.Ku.V8-20180723.20180220-S122714-E135947.022618.V06A.SUB.HDF5' ; Get the date and time from the filename parts=strsplit(cfile,'-',/extract) hdate=parts[2] htime=parts[3] hhour=strmid(htime,0,2) hmin=strmid(htime,2,2) hsec=strmid(htime,4,2) ; Read the hid file - no longer used do_hid='no' if do_hid eq 'yes' then begin 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 endif ; 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 time to julian_day day1=julday(1,1,1970,0,0,0) sec_per_day=long(24L*60L*60L) julian_day_cfr=double(day1+((base_time_cfr+time_offset_cfr)/sec_per_day)) caldat,julian_day_cfr,mm,dd,yy,hh,mi,ss for i=0,n_elements(julian_day_cfr)-1 do begin print,yy[i],mm[i],dd[i],hh[i],mi[i],ss[i],'cfr time' endfor ; Dimensions of cband s=size(alt_3d_cfr,/dimensions) numx_cfr=s[0] numy_cfr=s[1] numz_cfr=s[2] ; Read the gpm hdf5 file file_id=h5f_open(path_prefix+gdir+gfile) ; Group NS ns_id=h5g_open(file_id,'NS') dset_id=h5d_open(ns_id,'Latitude') lat_gpm=h5d_read(dset_id) h5d_close,dset_id dset_id=h5d_open(ns_id,'Longitude') lon_gpm=h5d_read(dset_id) h5d_close,dset_id ; Group NS/SLV slv_id=h5g_open(ns_id,'SLV') dset_id=h5d_open(slv_id,'zFactorCorrected') dbz_gpm=h5d_read(dset_id) h5d_close,dset_id h5g_close,slv_id ; Group NS/ScanTime time_id=h5g_open(ns_id,'ScanTime') dset_id=h5d_open(time_id,'DayOfYear') doy=h5d_read(dset_id) h5d_close,dset_id dset_id=h5d_open(time_id,'Year') year=h5d_read(dset_id) h5d_close,dset_id dset_id=h5d_open(time_id,'Month') month=h5d_read(dset_id) h5d_close,dset_id dset_id=h5d_open(time_id,'DayOfMonth') day=h5d_read(dset_id) h5d_close,dset_id dset_id=h5d_open(time_id,'Hour') hour=h5d_read(dset_id) h5d_close,dset_id dset_id=h5d_open(time_id,'Minute') minute=h5d_read(dset_id) h5d_close,dset_id dset_id=h5d_open(time_id,'Second') second=h5d_read(dset_id) h5d_close,dset_id dset_id=h5d_open(time_id,'MilliSecond') msecond=h5d_read(dset_id) h5d_close,dset_id h5g_close,time_id ; Group NS/PRE pre_id=h5g_open(ns_id,'PRE') dset_id=h5d_open(pre_id,'elevation') elevation_gpm=h5d_read(dset_id) ;meters dset_id=h5d_open(pre_id,'zFactorMeasured') dbz_measured_gpm=h5d_read(dset_id) h5d_close,dset_id h5g_close,pre_id ; Close HDF5 file h5g_close,ns_id h5f_close,file_id ; Non transpose way s=size(dbz_gpm,/dimensions) ;print,s nbin=s[0] ;176 nray=s[1] ;49 nscan=s[2] ;7937 ; Caculate colongitude because the pacific crosses the 180 line colon_gpm=lon_gpm result=where(lon_gpm lt 0 and lon_gpm ne -9999.90,count) if count gt 0 then colon_gpm[result]=colon_gpm[result]+360.0 ; Calculate time msecond=msecond*1e-3 second=second+msecond julian_day_gpm=julday(month,day,year,hour,minute,second) caldat,julian_day_gpm,mm,dd,yy,hh,mi,ss print,yy[0],mm[0],dd[0],hh[0],mi[0],ss[0],'gpm start time' print,yy[-1],mm[-1],dd[-1],hh[-1],mi[-1],ss[-1],'gpm end time' ;check to see if there are any dbz in this data if 1 eq 0 then begin for i=0,nbin-1 do begin for j=0,nray-1 do begin if max(dbz_gpm[i,j,*]) ne -9999.90 then print, max(dbz_gpm[i,j,*]) endfor endfor endif ; Print the limits of the data print,min(lat_3d_cfr),max(lat_3d_cfr),'lat cfr' print,min(lon_3d_cfr),max(lon_3d_cfr),'lon cfr' print,min(alt_3d_cfr),max(alt_3d_cfr),'alt cfr' r=where(lat_gpm ne -9999.90) print,min(lat_gpm[r]),max(lat_gpm[r]),'lat gpm' r=where(lon_gpm ne -9999.90) print,min(lon_gpm[r]),max(lon_gpm[r]),'lon gpm' ; Limits of cband data lat_min=min(lat_3d_cfr) lat_max=max(lat_3d_cfr) lon_min=min(lon_3d_cfr) lon_max=max(lon_3d_cfr) alt_min=min(alt_3d_cfr) alt_max=max(alt_3d_cfr) ; Subset gpm time to cband area lat1d=reform(lat_gpm[0,*]) lon1d=reform(lon_gpm[0,*]) r=where(lat1d ge lat_min and lat1d le lat_max and $ lon1d ge lon_min and lon1d le lon_max,c) julian_day_gpm=julian_day_gpm[r] caldat,julian_day_gpm,mm,dd,yy,hh,mi,ss print,yy[0],mm[0],dd[0],hh[0],mi[0],ss[0],'subset gpm start time' print,yy[-1],mm[-1],dd[-1],hh[-1],mi[-1],ss[-1],'subset gpm end time' ; Subset gpm to be in cband lat-lon range r=where(lat_gpm ge lat_min and lat_gpm le lat_max and $ lon_gpm ge lon_min and lon_gpm le lon_max,c) sub_lat_gpm=lat_gpm[r] sub_lon_gpm=lon_gpm[r] sub_dbz_gpm=make_array(c,nbin,/float,value=-8888) sub_dbz_measured_gpm=make_array(c,nbin,/float,value=-8888) sub_height_gpm=make_array(c,nbin,/float,value=-8888) ; gpm height array ;m ku_deltah=250 ;m for NS height_ku=(findgen(nbin)*ku_deltah) height_ku_reverse=reverse(height_ku) ; Loop through the gpm heights to subset each level for i=0,nbin-1 do begin dbz_gpm1=reform(dbz_gpm[i,*,*]) sub_dbz_gpm[*,i]=dbz_gpm1[r] dbz_measured_gpm1=reform(dbz_measured_gpm[i,*,*]) sub_dbz_measured_gpm[*,i]=dbz_measured_gpm1[r] sub_height_gpm[*,i]=height_ku_reverse[i] ;print,max(dbz_gpm1),max(dbz_measured_gpm1) endfor ;check to see if there are any dbz in this data if 1 eq 0 then begin for i=0,nbin-1 do begin if max(sub_dbz_gpm[*,i]) ne -9999.90 then print, max(sub_dbz_gpm[*,i]) if max(sub_dbz_measured_gpm[*,i]) ne -9999.90 then print, max(sub_dbz_measured_gpm[*,i]) endfor endif ; Now put gpm on cband lat-lon-alt grid sub_dbz_gpm_match=make_array(numx_cfr,numy_cfr,numz_cfr,/float,value=-8888) sub_dbz_measured_gpm_match=make_array(numx_cfr,numy_cfr,numz_cfr,/float,value=-8888) ;sub_dbz_gpm_match=make_array(numx_cfr,numy_cfr,n_elements(height_ku),/float,value=-8888) ;sub_dbz_measured_gpm_match=make_array(numx_cfr,numy_cfr,n_elements(height_ku),/float,value=-8888) for i=0,numx_cfr-1 do begin for j=0,numy_cfr-1 do begin lat1=lat_3d_cfr[i,j,0] lon1=lon_3d_cfr[i,j,0] ; check to see if there is any cband ;dbz1=dbz_cfr[i,j,*] ;print,lon1,lat1,max(dbz1) ; Lat-lon match r=where(sub_lat_gpm ge lat1-0.1 and sub_lat_gpm le lat1+0.1 and $ sub_lon_gpm ge lon1-0.1 and sub_lon_gpm le lon1+0.1,c) close_dist=1e6 close_k=1e6 for k=0,c-1 do begin dist_point=map_2points(lon1,lat1,sub_lon_gpm[r[k]],sub_lat_gpm[r[k]],/meters) if dist_point lt close_dist then begin close_dist=dist_point close_k=k endif ;print, k,c,dist_point,sub_lon_gpm[r[k]],sub_lat_gpm[r[k]] endfor ; puts in the full gpm profile ; if close_k ne 1e6 then begin ; ;print, close_k,c,close_dist,sub_lon_gpm[r[close_k]],sub_lat_gpm[r[close_k]] ; sub_dbz_gpm_match[i,j,*]=sub_dbz_gpm[r[close_k],*] ; sub_dbz_measured_gpm_match[i,j,*]=sub_dbz_measured_gpm[r[close_k],*] ; endif ; height match if close_k ne 1e6 and close_dist lt 10000 then begin ;print, close_k,c,close_dist,sub_lon_gpm[r[close_k]],sub_lat_gpm[r[close_k]] ;p2=plot(dbz_cfr[i,j,*],alt_3d_cfr[i,j,*],position=pos[pnum,*],/current) dbz2=sub_dbz_gpm[r[close_k],*] dbz3=sub_dbz_measured_gpm[r[close_k],*] for l=0,numz_cfr-1 do begin closest=min(abs(height_ku_reverse-alt_3d_cfr[i,j,l]),cidx) ;print,closest ;;print,dbz2[cidx] ;****screening for minimum****** ;if dbz2[cidx] gt 14 then begin ;14 dbz is minimum gpm sensitivity (found a min val of 13.72) if closest lt 25.0 then begin ;height match sub_dbz_gpm_match[i,j,l]=dbz2[cidx] sub_dbz_measured_gpm_match[i,j,l]=dbz3[cidx] endif ;endif endfor endif endfor endfor ; Check to see if there is any dbz in this data if 1 eq 0 then begin for i=0,numz_cfr-1 do begin ;for i=0,n_elements(height_ku)-1 do begin print, max(sub_dbz_gpm_match[*,*,i]),max(sub_dbz_measured_gpm_match[*,*,i]) endfor endif ; Quick check of gpm-matched data ;sub_dbz_gpm_match ;sub_dbz_gpm_measured_match r=where(sub_dbz_measured_gpm_match ne -28888.0 and $ ;gpm fill sub_dbz_measured_gpm_match ne -8888.0 and $ ;my fill sub_dbz_measured_gpm_match ge -20 and $ ;min data value sub_dbz_measured_gpm_match lt 25) start_bin=min(sub_dbz_measured_gpm_match[r]) end_bin=max(sub_dbz_measured_gpm_match[r]) dbin=(end_bin-start_bin)/150.0 data_hist=sub_dbz_measured_gpm_match hist_generic,data_hist,start_bin,end_bin,dbin,bins,data_freq,data_counts ;p0=plot(bins,data_freq) ;******************************************* ; Color Table ;******************************************* ; Top is the last color to scale 256 colors, 0-255 top_color=253 ; Colortable 0-253s mytable=colortable(43,ncolors=254) ;254=white ;gray=255 mytable=[mytable,transpose([255,255,255]),transpose([220,220,220])] mycbtable=mytable[0:top_color,*] ;******************************************* ; Start plotting ;******************************************* ; Set up plot size pxdim=1000 & pydim=900 ; Position the plots xl=0.1 & xr=0.90 yb=0.2 & yt=0.97 sx=0.15 sy=0.18 numplots_x=2 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.08 cbpos[*,2]=cbpos[*,0]+0.007 ;cbpos[*,1]=cbpos[*,1]+0.02 ;cbpos[*,3]=cbpos[*,3]-0.02 dx=0.01 & dy=0.01 fs1=12 ;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]) ;*** Dbz limits min_dbz=0.0 max_dbz=35.0 dbz_floor=14.0 ;*** Volume plot of CBAND dbz (all values) pnum=0 if 1 eq 0 then begin data_image=bytscl(dbz_cfr,min=min_dbz,max=max_dbz) p0=volume(data_image,rgb_table0=43,position=pos[pnum,*],/buffer,$ dimensions=[pxdim,pydim],axis_style=2,$ volume_dimensions=[(lon_max-lon_min),(lat_max-lat_min),(alt_max-alt_min)/1000.0],$ volume_location=[lon_min,lat_min,alt_min/1000.0],$ xtitle='Longitude',ytitle='Latitude',ztitle='km') ax=p0.axes ;ax[6].hide=1 ;removes the closest top cube edge ;ax[0].xrange=[lon_min,lon_max] c0=colorbar(range=[min_dbz,max_dbz],$ orientation=1,position=reform(cbpos[pnum,*]),font_size=fs1,rgb_table=43,$ tickdir=1,border_on=1,title='dBZ') endif else begin ;*** Image sllice of cband dbz (all values) alt1=reform(alt_3d_cfr[200,200,*]) r=where(alt1 ge 2000,c) data_var=reform(dbz_cfr[*,*,r[0]]) data_image=bytscl(data_var,min=min_dbz,max=max_dbz,top=top_color) r2=where(data_var le 0,c2) if c gt 0 then data_image[r2]=255 p0=image(data_image,$ rgb_table=mytable,position=pos[pnum,*],/buffer,$ dimensions=[pxdim,pydim],image_dimensions=[isx,isy]) c0=contour(data_image,lon_3d_cfr[*,0,0],lat_3d_cfr[0,*,0],$ xtitle='longitude',ytitle='latitiude',xstyle=1,ystyle=1,/nodata,$ /current,position=pos[pnum,*],font_size=fs1,xtickdir=1,ytickdir=1,$ axis_style=2) t1=text(pos[pnum,0]+20*dx,pos[pnum,3]+2*dy,$ string(alt1[r[0]]/1000.0,format='(F5.2)')+' km',font_size=fs2) endelse t1=text(pos[pnum,0]+2*dx,pos[pnum,3]+2*dy,'CBAND Z (dBZ)',font_size=fs2) ;Time of gpm caldat,julian_day_cfr[0],mm,dd,yy,hh,mi,ss cdate_str=string(yy,format='(I4)')+string(mm,format='(I02)')+string(dd,format='(I02)') ctime_str=string(hh,format='(I02)')+':'+string(mi,format='(I02)')+':'+string(ss,format='(I02)') ctime_str2=string(hh,format='(I02)')+string(mi,format='(I02)')+string(ss,format='(I02)') t1=text(pos[pnum,0]+2*dx,pos[pnum,3]+4*dy,cdate_str+' '+ctime_str,font_size=fs2) ;*** Volume plot of GPM dbz (all values) pnum=1 if 1 eq 0 then begin r=where(sub_dbz_gpm_match gt -8888,c) ;min_dbz=min(sub_dbz_gpm_match[r]) ;max_dbz=max(sub_dbz_gpm_match[r]) data_image=bytscl(sub_dbz_gpm_match,min=min_dbz,max=max_dbz) p1=volume(data_image,rgb_table0=43,position=pos[pnum,*],/current,$ dimensions=[pxdim,pydim],axis_style=2,$ volume_dimensions=[(lon_max-lon_min),(lat_max-lat_min),(alt_max-alt_min)/1000.0],$ volume_location=[lon_min,lat_min,alt_min/1000.0],$ xtitle='Longitude',ytitle='Latitude',ztitle='km') ax1=p1.axes ;ax1[6].hide=1 ;removes the closest top cube edge endif else begin ;*** Image slice of GPM dbz (all values) alt1=reform(alt_3d_cfr[200,200,*]) r=where(alt1 ge 2000,c) data_var=reform(sub_dbz_gpm_match[*,*,r[0]]) data_image=bytscl(data_var,min=min_dbz,max=max_dbz,top=top_color) r2=where(data_var le 0,c2) if c gt 0 then data_image[r2]=255 p0=image(data_image,$ rgb_table=mytable,position=pos[pnum,*],/current,$ dimensions=[pxdim,pydim],image_dimensions=[isx,isy]) c0=contour(data_image,lon_3d_cfr[*,0,0],lat_3d_cfr[0,*,0],$ xtitle='longitude',ytitle='latitiude',xstyle=1,ystyle=1,/nodata,$ /current,position=pos[pnum,*],font_size=fs1,xtickdir=1,ytickdir=1,$ axis_style=2) t1=text(pos[pnum,0]+20*dx,pos[pnum,3]+2*dy,$ string(alt1[r[0]]/1000.0,format='(F5.2)')+' km',font_size=fs2) endelse t1=text(pos[pnum,0]+2*dx,pos[pnum,3]+2*dy,'GPM Z (dBZ)',font_size=fs2) c0=colorbar(range=[min_dbz,max_dbz],$ orientation=1,position=reform(cbpos[pnum,*]),font_size=fs1,rgb_table=43,$ tickdir=1,border_on=1,title='dBZ') ;Time of gpm caldat,julian_day_gpm[0],mm,dd,yy,hh,mi,ss gdate_str=string(yy,format='(I4)')+string(mm,format='(I02)')+string(dd,format='(I02)') gtime_str=string(hh,format='(I02)')+':'+string(mi,format='(I02)')+':'+string(ss,format='(I02)') gtime_str2=string(hh,format='(I02)')+string(mi,format='(I02)')+string(ss,format='(I02)') t1=text(pos[pnum,0]+2*dx,pos[pnum,3]+4*dy,gdate_str+' '+gtime_str,font_size=fs2) ; 1D histogram of DBZ ge gpm noise floor start_bin=-20 end_bin=40 dbin=0.5 ; histogram cband data gt than gpm noise floor (14?) rgnoise=where(dbz_cfr ge dbz_floor,cgnoise) data_hist=dbz_cfr[rgnoise] ; histogram all cband data ;data_hist=dbz_cfr hist_generic,data_hist,start_bin,end_bin,dbin,bins,data_freq,data_counts pnum=2 p0=plot(bins,data_freq,position=pos[pnum,*],/current,dimensions=[pxdim,pydim],/histogram,$ xtitle='dBZ',ytitle='Frequency',color='red',$ ;xrange=[0,0.8], ;yrange=[0,0.45], font_size=fs1,xstyle=1) r=where(sub_dbz_gpm_match ge dbz_floor,c) if c gt 0 then begin data_hist=sub_dbz_gpm_match[r] hist_generic,data_hist,start_bin,end_bin,dbin,bins,data_freq,data_counts p2=plot(bins,data_freq,/overplot,/histogram,color='purple') endif t1=text(pos[pnum,0]+2*dx,pos[pnum,3]-4*dy,'CBAND',color='red',font_size=fs2) t1=text(pos[pnum,0]+2*dx,pos[pnum,3]-8*dy,'GPM',color='purple',font_size=fs2) ; 1D histogram of HEIGHT of DBZ above noise floor start_bin=0 end_bin=15000 dbin=125 ; All Height ;r=where(dbz_cfr gt -9999,c) ;data_hist=alt_3d_cfr[r] ; histogram cband data gt noise floor (14) data_hist=alt_3d_cfr[rgnoise] hist_generic,data_hist,start_bin,end_bin,dbin,bins,data_freq,data_counts pnum=3 p0=plot(bins/1000.0,data_freq,position=pos[pnum,*],/current,dimensions=[pxdim,pydim],/histogram,$ xtitle='Height (m)',ytitle='Frequency',color='red',$ ;xrange=[0,0.8], ;yrange=[0,0.45], ;xtext_orientation=90,$ clip=0,$ font_size=fs1,xstyle=1) r=where(sub_dbz_gpm_match ge dbz_floor,c) if c gt 0 then begin data_hist=alt_3d_cfr[r] hist_generic,data_hist,start_bin,end_bin,dbin,bins,data_freq,data_counts p2=plot(bins/1000.0,data_freq,/overplot,/histogram,color='purple') endif ;p0.save,'cband_gpm_match.'+gdate_str+'.'+gtime_str2+'.g14dbz.png',height=pydim ;cband file time imagename='cband_gpm_match.'+hdate+'.'+htime+'.g14dbz.xc.png' p0.save,imagename,height=pydim print,imagename ;****************************************** ; Same plot but with gpm measured (not masked) dbz ;****************************************** ;*** Volume plot of CBAND dbz (all values) pnum=0 if 1 eq 0 then begin r=where(dbz_cfr ne -9999,c) ;min_dbz=min(dbz_cfr[r]) ;max_dbz=max(dbz_cfr[r]) data_image=bytscl(dbz_cfr,min=min_dbz,max=max_dbz) p0=volume(data_image,rgb_table0=43,position=pos[pnum,*],/buffer,$ dimensions=[pxdim,pydim],axis_style=2,$ volume_dimensions=[(lon_max-lon_min),(lat_max-lat_min),(alt_max-alt_min)/1000.0],$ volume_location=[lon_min,lat_min,alt_min/1000.0],$ xtitle='Longitude',ytitle='Latitude',ztitle='km') ax=p0.axes ;ax[6].hide=1 ;removes the closest top cube edge ;ax[0].xrange=[lon_min,lon_max] c0=colorbar(range=[min_dbz,max_dbz],$ orientation=1,position=reform(cbpos[pnum,*]),font_size=fs1,rgb_table=43,$ tickdir=1,border_on=1,title='dBZ') endif else begin ;*** Image sllice of cband dbz (all values) alt1=reform(alt_3d_cfr[200,200,*]) r=where(alt1 ge 2000,c) data_var=reform(dbz_cfr[*,*,r[0]]) data_image=bytscl(data_var,min=min_dbz,max=max_dbz,top=top_color) r2=where(data_var le 0,c2) if c gt 0 then data_image[r2]=255 p0=image(data_image,$ rgb_table=mytable,position=pos[pnum,*],/buffer,$ dimensions=[pxdim,pydim],image_dimensions=[isx,isy]) c0=contour(data_image,lon_3d_cfr[*,0,0],lat_3d_cfr[0,*,0],$ xtitle='longitude',ytitle='latitiude',xstyle=1,ystyle=1,/nodata,$ /current,position=pos[pnum,*],font_size=fs1,xtickdir=1,ytickdir=1,$ axis_style=2) t1=text(pos[pnum,0]+20*dx,pos[pnum,3]+2*dy,$ string(alt1[r[0]]/1000.0,format='(F5.2)')+' km',font_size=fs2) endelse t1=text(pos[pnum,0]+2*dx,pos[pnum,3]+2*dy,'CBAND Z (dBZ)',font_size=fs2) ;Time of cband caldat,julian_day_cfr[0],mm,dd,yy,hh,mi,ss cdate_str=string(yy,format='(I4)')+string(mm,format='(I02)')+string(dd,format='(I02)') ctime_str=string(hh,format='(I02)')+':'+string(mi,format='(I02)')+':'+string(ss,format='(I02)') ctime_str2=string(hh,format='(I02)')+string(mi,format='(I02)')+string(ss,format='(I02)') t1=text(pos[pnum,0]+2*dx,pos[pnum,3]+4*dy,cdate_str+' '+ctime_str,font_size=fs2) ;*** Volume plot of GPM dbz (above noise floor) pnum=1 if 1 eq 0 then begin data_var=sub_dbz_measured_gpm_match ;r=where(data_var lt -20,c) ;data min r=where(data_var lt dbz_floor,c) ;noise floor min if c gt 0 then data_var[r]=-9999 ;min_dbz=min(sub_dbz_gpm_match[r]) ;max_dbz=max(sub_dbz_gpm_match[r]) data_image=bytscl(data_var,min=min_dbz,max=max_dbz) p1=volume(data_image,rgb_table0=43,position=pos[pnum,*],/current,$ dimensions=[pxdim,pydim],axis_style=2,$ volume_dimensions=[(lon_max-lon_min),(lat_max-lat_min),(alt_max-alt_min)/1000.0],$ volume_location=[lon_min,lat_min,alt_min/1000.0],$ xtitle='Longitude',ytitle='Latitude',ztitle='km') ax1=p1.axes ;ax1[6].hide=1 ;removes the closest top cube edge endif else begin ;*** Image slice of GPM dbz (all values) alt1=reform(alt_3d_cfr[200,200,*]) r=where(alt1 ge 2000,c) data_var=reform(sub_dbz_measured_gpm_match[*,*,r[0]]) data_image=bytscl(data_var,min=min_dbz,max=max_dbz,top=top_color) r2=where(data_var le 0,c2) if c gt 0 then data_image[r2]=255 p0=image(data_image,$ rgb_table=mytable,position=pos[pnum,*],/current,$ dimensions=[pxdim,pydim],image_dimensions=[isx,isy]) c0=contour(data_image,lon_3d_cfr[*,0,0],lat_3d_cfr[0,*,0],$ xtitle='longitude',ytitle='latitiude',xstyle=1,ystyle=1,/nodata,$ /current,position=pos[pnum,*],font_size=fs1,xtickdir=1,ytickdir=1,$ axis_style=2) t1=text(pos[pnum,0]+20*dx,pos[pnum,3]+2*dy,$ string(alt1[r[0]]/1000.0,format='(F5.2)')+' km',font_size=fs2) endelse t1=text(pos[pnum,0]+2*dx,pos[pnum,3]+2*dy,'GPM Z (dBZ)',font_size=fs2) c0=colorbar(range=[min_dbz,max_dbz],$ orientation=1,position=reform(cbpos[pnum,*]),font_size=fs1,rgb_table=43,$ tickdir=1,border_on=1,title='dBZ') ;Time of gpm caldat,julian_day_gpm[0],mm,dd,yy,hh,mi,ss gdate_str=string(yy,format='(I4)')+string(mm,format='(I02)')+string(dd,format='(I02)') gtime_str=string(hh,format='(I02)')+':'+string(mi,format='(I02)')+':'+string(ss,format='(I02)') gtime_str2=string(hh,format='(I02)')+string(mi,format='(I02)')+string(ss,format='(I02)') t1=text(pos[pnum,0]+2*dx,pos[pnum,3]+4*dy,gdate_str+' '+gtime_str,font_size=fs2) ; 1D histogram of DBZ start_bin=-20 end_bin=40 dbin=0.5 ; histogram cband data gt 14 rgnoise=where(dbz_cfr ge dbz_floor,cgnoise) data_hist=dbz_cfr[rgnoise] ; histogram all cband data ;data_hist=dbz_cfr hist_generic,data_hist,start_bin,end_bin,dbin,bins,data_freq,data_counts pnum=2 p0=plot(bins,data_freq,position=pos[pnum,*],/current,dimensions=[pxdim,pydim],/histogram,$ xtitle='dBZ',ytitle='Frequency',color='red',$ ;xrange=[0,0.8], ;yrange=[0,0.45], font_size=fs1,xstyle=1) r=where(sub_dbz_measured_gpm_match ge dbz_floor,c) if c gt 0 then begin data_hist=sub_dbz_measured_gpm_match[r] hist_generic,data_hist,start_bin,end_bin,dbin,bins,data_freq,data_counts endif p2=plot(bins,data_freq,/overplot,/histogram,color='purple') t1=text(pos[pnum,0]+2*dx,pos[pnum,3]-4*dy,'CBAND',color='red',font_size=fs2) t1=text(pos[pnum,0]+2*dx,pos[pnum,3]-8*dy,'GPM',color='purple',font_size=fs2) ; 1D histogram of HEIGHT of DBZ gt noise floor start_bin=0 end_bin=15000 dbin=125 ; All Height ;r=where(dbz_cfr gt -9999,c) ;data_hist=alt_3d_cfr[r] ; histogram cband data gt 14 data_hist=alt_3d_cfr[rgnoise] hist_generic,data_hist,start_bin,end_bin,dbin,bins,data_freq,data_counts pnum=3 p0=plot(bins/1000.0,data_freq,position=pos[pnum,*],/current,dimensions=[pxdim,pydim],/histogram,$ xtitle='Height (km)',ytitle='Frequency',color='red',$ ;xrange=[0,0.8], ;yrange=[0,0.45], ;xtext_orientation=90,$ clip=0,$ font_size=fs1,xstyle=1) r=where(sub_dbz_measured_gpm_match ge dbz_floor and sub_dbz_measured_gpm_match lt 20,c) ;should be -20 if c gt 0 then begin data_hist=alt_3d_cfr[r] hist_generic,data_hist,start_bin,end_bin,dbin,bins,data_freq,data_counts p2=plot(bins/1000.0,data_freq,/overplot,/histogram,color='purple') endif ;p0.save,'cband_gpm_measured_match.'+gdate_str+'.'+gtime_str2+'.g14dbz.png',height=pydim ; cband time p0.save,'cband_gpm_measured_match.'+hdate+'.'+htime+'.g14dbz.xc.png',height=pydim stop ; DBZ start_bin=-20 end_bin=50 dbin=0.5 data_hist=dbz_cfr hist_generic,data_hist,start_bin,end_bin,dbin,bins,data_freq,data_counts pnum=0 p0=plot(bins,data_freq,position=pos[pnum,*],/buffer,dimensions=[pxdim,pydim],/histogram,$ xtitle='cfradial dbz',ytitle='Frequency',color='blue',$ ;xrange=[0,0.8], ;yrange=[0,0.45], font_size=fs1,xstyle=1) r=where(sub_height_gpm le max(alt_3d_cfr),c) data_hist=sub_dbz_measured_gpm[r] hist_generic,data_hist,start_bin,end_bin,dbin,bins,data_freq,data_counts p2=plot(bins,data_freq,/overplot,/histogram,color='green') ; HEIGHT start_bin=0 end_bin=15000 dbin=10 r=where(dbz_cfr gt -9999 and alt_3d_cfr lt 10000,c) data_hist=alt_3d_cfr[r] hist_generic,data_hist,start_bin,end_bin,dbin,bins,data_freq,data_counts pnum=1 p0=plot(bins,data_freq,position=pos[pnum,*],/current,dimensions=[pxdim,pydim],/histogram,$ xtitle='cfradial height',ytitle='Frequency',color='blue',$ ;xrange=[0,0.8], ;yrange=[0,0.45], xtext_orientation=90,$ font_size=fs1,xstyle=1) r=where(sub_dbz_measured_gpm gt -50 and sub_height_gpm lt 10000,c) data_hist=sub_height_gpm[r] hist_generic,data_hist,start_bin,end_bin,dbin,bins,data_freq,data_counts p2=plot(bins,data_freq,/overplot,/histogram,color='green') ; loop through cband and plot closest gpm profile ; Position the plots xl=0.1 & xr=0.90 yb=0.35 & yt=0.95 sx=0.15 sy=0.13 numplots_x=4 numplots_y=2 position_plots,xl,xr,yb,yt,sx,sy,numplots_x,numplots_y,pos pnum=-1 for i=0,numx_cfr-1 do begin for j=0,numy_cfr-1 do begin lat1=lat_3d_cfr[i,j,0] lon1=lon_3d_cfr[i,j,0] ; check to see if there is any cband dbz1=dbz_cfr[i,j,*] print,lat1,lon1,max(dbz1) rr=where(dbz1 gt -50,cc) if cc gt 0 then begin ;print,lat1,lon1,max(dbz1) r=where(sub_lat_gpm ge lat1-0.1 and sub_lat_gpm le lat1+0.1 and $ sub_lon_gpm ge lon1-0.2 and sub_lat_gpm le lon1+0.2,c) close_dist=1e6 close_k=1e6 for k=0,c-1 do begin dist_point=map_2points(lon1,lat1,sub_lon_gpm[r[k]],sub_lat_gpm[r[k]],/meters) if dist_point lt close_dist then begin close_dist=dist_point close_k=k endif ;print, k,c,close_dist,close_k,sub_lon_gpm[r[k]],sub_lat_gpm[r[k]] endfor if close_k ne 1e6 then begin pnum=pnum+1 if pnum eq 8 then goto,finish_plot p2=plot(dbz_cfr[i,j,*],alt_3d_cfr[i,j,*],position=pos[pnum,*],/current) p3=plot(sub_dbz_measured_gpm[r[close_k],*],sub_height_gpm[r[close_k],*],/overplot,color='red') endif endif ;found some cband dbz endfor endfor finish_plot: p0.save,'cband_gpm_histograms2.png',height=pydim stop ; dbz_gpm=dbz_gpm[r] ; dbz_measured_gpm=dbz_measured_gpm[r] ; elevation_gpm=elevation_gpm[r] stop ;; 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