;******************************************** ; I want to expand this to the Olympex NPOL data. ; For that data set, I want to create separate histograms ; for oceanic versus mountains. The oceanic would be ; everything west of 124.5 West and the mountains would ; be a sector north of 47 N and east of 124.75. ;************************************************** pro cband_stats_observables,filenames,rtag ;filenames='/uufs/chpc.utah.edu/common/home/mace-group5/field_programs/investigator/in2018_v01/cband/9776HUB-PPIVol-20180209-000559-0000.cfradial.pyart.81.401.nc' ;filenames='/Volumes/mace-group5/field_programs/investigator/in2018_v01/cband/9776HUB-PPIVol-20180209-000559-0000.cfradial.pyart.81.401.nc' ;filenames='/Volumes/mace-group5/field_programs/investigator/in2018_v01/cband/20180116/9776HUB-PPIVol-20180116-235315-0000.cfradial.nc' if filenames eq !NULL then begin ; Choose date and time of cband ;date_str='20180209' & time_str='*' ;capricorn2 ;date_str='20180209' & time_str='20*' ;;date_str='20180209' & time_str='164335' ;can't read this pyart file. may need to remake it date_str='20180209' & time_str='2' ;date_str='20151204' & time_str='163710' ;olympex ;date_str='20180221' & time_str='*';'145702'; 101624';163706' ;olympex ;date_str='20200328' & time_str='*';'085004';'145702'; 101624';163706' ;olympex ; Get year month day strings syy=strmid(date_str,0,4) smm=strmid(date_str,4,2) sdd=strmid(date_str,6,2) shh=strmid(time_str,0,2) smi=strmid(time_str,2,2) sss=strmid(time_str,4,2) ;date_str=syy+smm+sdd+shh+smi+sss ;ymd_str=syy+smm+sdd ; chpc or imac ;path_prefix='/Volumes/' path_prefix='/uufs/chpc.utah.edu/common/home/' ; directory if syy eq '2020' then begin ;comble filedir=path_prefix+'mace-group5/field_programs/comble/cband_somna/' filestr='nosmn.ss.vol-dbz*noblock.'+syy+smm+sdd+'T'+time_str+'*pyart.41.201.nc' radar_name='nosmn_cband' endif else if syy eq '2018' then begin ;filedir=path_prefix+'mace-group4/Capricorn2/CBand/' ;filedir=path_prefix+'mace-group4/Capricorn2/CBand/cfradial/' filedir=path_prefix+'mace-group5/field_programs/investigator/in2018_v01/cband/'+date_str+'/' filestr='9776HUB-PPIVol-'+syy+smm+sdd+'-'+time_str+'*pyart.81.401.nc' radar_name='cband' endif else if syy eq '2019' then begin filedir=path_prefix+'mace-group5/field_programs/investigator/in2019_v01_gpm_passes/out_vol/' filestr='9776HUB-PPIVol-'+syy+smm+sdd+'-'+time_str+'*pyart*.nc' endif else if syy eq '2015' or syy eq '2016' then begin filedir=path_prefix+'mace-group6/field_programs/olympex/sband/'+date_str+'/ppi/' filestr='olympex_NPOL1_'+date_str+'_'+time_str+'*pyart*.nc' radar_name='sband' rtag='oceanic' ;rtag='mountain' ; Get the flight path ;read_nav_olympex,path_prefix,date_str,julian_day_plane,temp,tas,theta,pres,ppmvw,lat_plane,lon_plane,$ ; alt_plane,wwind,uwind,vwind,cwc_kng,twc_nev,lwc_nev,iwc_nev,tot_s,cwc_s,rh,cwc_csi,tot_c,tot_p endif ; file name filenames=file_search(filedir+filestr,count=numfiles) print,'numfiles',numfiles endif else begin ;filenames='/uufs/chpc.utah.edu/common/home/mace-group6/field_programs/olympex/sband/20151130/olympex_NPOL1_20151130*.pyart.81.401.nc" filedir=file_dirname(filenames)+'/' ; Get year month day strings if strmid(file_basename(filenames),0,7) eq 'olympex' then begin parts=strsplit(file_basename(filenames),'_',/extract) endif else begin parts=strsplit(file_basename(filenames),'-',/extract) endelse syy=strmid(parts[2],0,4) smm=strmid(parts[2],4,2) sdd=strmid(parts[2],6,2) shh=strmid(parts[3],0,2) smi=strmid(parts[3],2,2) sss=strmid(parts[3],4,2) date_str=syy+smm+sdd if syy eq '2018' then radar_name='cband' numfiles=1 endelse if syy eq '2015' or syy eq '2016' then begin radar_name='sband' ; choose one rtag if rtag eq 'oceanic' then olon=-124.5 if rtag eq 'mountain' then begin mlon=-124.75 & mlat=47 endif endif else begin rtag='' endelse ; Loop through the cband files for f=0,numfiles-1 do begin filename=filenames[f] print,filename if syy eq '2015' then begin parts=strsplit(file_basename(filename),'_',/extract) shh=strmid(parts[3],0,2) smi=strmid(parts[3],2,2) sss=strmid(parts[3],4,2) time_str=shh+smi+sss endif else if syy eq '2018' then begin parts=strsplit(file_basename(filename),'-',/extract) shh=strmid(parts[3],0,2) smi=strmid(parts[3],2,2) sss=strmid(parts[3],4,2) time_str=shh+smi+sss endif else if syy eq '2020' then begin parts=strsplit(file_basename(filename),'.',/extract) shh=strmid(parts[4],9,2) smi=strmid(parts[4],11,2) sss=strmid(parts[4],13,2) time_str=shh+smi+sss endif print,'start reading data' if syy eq '2020' then begin read_nosmn_cband,filename,lat,lon,point_z,rho,kdp,zdr,dbz,$ point_x,point_y,base_time,time_offset endif else begin ; all data in one file fid=ncdf_open(filename) ;vid=ncdf_varid(fid,'base_time') & ncdf_varget,fid,vid,base_time ;vid=ncdf_varid(fid,'time_offset') & ncdf_varget,fid,vid,time_offset ;vid=ncdf_varid(fid,'radar_latitude') & ncdf_varget,fid,vid,lat_radar ;vid=ncdf_varid(fid,'radar_longitude') & ncdf_varget,fid,vid,lon_radar vid=ncdf_varid(fid,'point_latitude') & ncdf_varget,fid,vid,lat vid=ncdf_varid(fid,'point_longitude') & ncdf_varget,fid,vid,lon ;vid=ncdf_varid(fid,'point_altitude') & ncdf_varget,fid,vid,alt ;vid=ncdf_varid(fid,'point_x') & ncdf_varget,fid,vid,point_x ;vid=ncdf_varid(fid,'point_y') & ncdf_varget,fid,vid,point_y vid=ncdf_varid(fid,'point_z') & ncdf_varget,fid,vid,point_z vid=ncdf_varid(fid,'cross_correlation_ratio') & ncdf_varget,fid,vid,rho vid=ncdf_varid(fid,'specific_differential_phase') if vid eq -1 then vid=ncdf_varid(fid,'corrected_specific_differential_phase') ncdf_varget,fid,vid,kdp vid=ncdf_varid(fid,'corrected_differential_reflectivity') & ncdf_varget,fid,vid,zdr vid=ncdf_varid(fid,'corrected_reflectivity') if vid eq -1 then vid=ncdf_varid(fid,'DBZH_gmm') ncdf_varget,fid,vid,dbz ncdf_close,fid endelse print,'end reading data' ; convert point_x,y,z from m to km ;point_x=point_x/1000.0 ;point_y=point_y/1000.0 point_z=point_z/1000.0 ; convert time to julian_day ;day1=julday(1,1,1970,0,0,0) ;sec_per_day=long(24L*60L*60L) ;julian_day=double(day1+((base_time+time_offset)/sec_per_day)) ;caldat,julian_day,mm,dd,yy,hh,mi,ss ;print,yy,mm,dd,hh,mi,ss ; Calculate colongitude ;colon=lon ;result=where(lon lt 0,count) ;if count gt 0 then colon[result]=colon[result]+360.0 ; if oceanic for sband only keep data over ocean if rtag eq 'oceanic' and radar_name eq 'sband' then begin print,'histogram over ocean' r=where(lon le olon,c) lat=lat[r] lon=lon[r] point_z=point_z[r] rho=rho[r] kdp=kdp[r] zdr=zdr[r] dbz=dbz[r] print,mean(lon) endif ; if mountain for sband only keep data over mountains if rtag eq 'mountain' and radar_name eq 'sband' then begin print,'histogram over mountain' r=where(lon ge mlon and lat ge mlat,c) lat=lat[r] lon=lon[r] point_z=point_z[r] rho=rho[r] kdp=kdp[r] zdr=zdr[r] dbz=dbz[r] print,mean(lat),mean(lon) endif ; Filter out where rho is less than 0.95 r=where(rho lt 0.95,c) if c gt 0 then begin kdp[r]=-9999 zdr[r]=-9999 dbz[r]=-9999 endif ; Dbz Vector ;start_bin=-10.0 ;end_bin=50.0 ;dbin_dbz=1.0 start_bin=0.0 end_bin=40.0 ;dbin_dbz=0.5/2.0 ;regular files dbin_dbz=0.25/2.0 ;hires files bins_dbz=start_bin while bins_dbz[n_elements(bins_dbz)-1] lt end_bin $ do bins_dbz=[bins_dbz,bins_dbz[n_elements(bins_dbz)-1]+(2.*dbin_dbz)] hist_generic,dbz,start_bin,end_bin,dbin_dbz,bins,freq_dbz,counts_dbz ; Zdr Vector start_bin=-2.0 end_bin=5.0 ;dbin_zdr=0.1/2.0 ;regular files dbin_zdr=0.05/2.0 ;hires files ;r=where(zdr ne -9999) ;start_bin=floor(min(zdr[r])) ;end_bin=ceil(max(zdr[r])) ;dbin_zdr=(end_bin-start_bin)/50.0 bins_zdr=start_bin while bins_zdr[n_elements(bins_zdr)-1] lt end_bin $ do bins_zdr=[bins_zdr,bins_zdr[n_elements(bins_zdr)-1]+(2.*dbin_zdr)] hist_generic,zdr,start_bin,end_bin,dbin_zdr,bins,zdr_freq,counts_zdr ; Kdp vector ;start_bin=-0.5 ;end_bin=2.5 ;dbin_kdp=0.02/2.0 start_bin=-0.2 end_bin=1.0 ;dbin_kdp=0.02/2.0 ;regular files dbin_kdp=0.01/2.0 ;hires files bins_kdp=start_bin while bins_kdp[n_elements(bins_kdp)-1] lt end_bin $ do bins_kdp=[bins_kdp,bins_kdp[n_elements(bins_kdp)-1]+(2.*dbin_kdp)] hist_generic,kdp,start_bin,end_bin,dbin_kdp,bins,kdp_freq,counts_kdp ; height vector ;start_bin=0 ;end_bin=20.0 ;dbin_hgt=0.25/2.0 start_bin=0 end_bin=5.0 dbin_hgt=0.25/2.0 bins_hgt=start_bin while bins_hgt[n_elements(bins_hgt)-1] lt end_bin $ do bins_hgt=[bins_hgt,bins_hgt[n_elements(bins_hgt)-1]+(2.*dbin_hgt)] hist_generic,point_z,start_bin,end_bin,dbin_hgt,bins,hgt_freq,counts_hgt ;Old way that takes too long ;obs_counts=lonarr(n_elements(bins_dbz),n_elements(bins_zdr),n_elements(bins_kdp)) ;for i=0,n_elements(bins_dbz)-1 do begin ; print,'i',i,bins_dbz[i]-dbin_dbz, bins_dbz[i]+dbin_dbz,' dbz' ; for j=0,n_elements(bins_zdr)-1 do begin ; for k=0,n_elements(bins_kdp)-1 do begin ; r=where(dbz ge bins_dbz[i]-dbin_dbz and dbz lt bins_dbz[i]+dbin_dbz and $ ; zdr ge bins_zdr[j]-dbin_zdr and zdr lt bins_zdr[j]+dbin_zdr and $ ; kdp ge bins_kdp[k]-dbin_kdp and kdp lt bins_kdp[k]+dbin_kdp,count) ; obs_counts[i,j,k]=count ; ;if count gt 0 then begin ; ; print,bins_dbz[i]-dbin_dbz, bins_dbz[i]+dbin_dbz,' dbz' ; ; print,dbz[r] ; ; print,bins_zdr[j]-dbin_zdr, bins_zdr[j]+dbin_zdr,' zdr' ; ; print,zdr[r] ; ; print,bins_kdp[k]-dbin_kdp, bins_kdp[k]+dbin_kdp,' kdp' ; ; print,kdp[r] ; ; print, count ; ;endif ; endfor ; endfor ;endfor print,'start histo' r_valid=where(dbz ne -9999 and zdr ne -9999 and kdp ne -9999,count_valid) obs_counts3d=lonarr(n_elements(bins_dbz),n_elements(bins_zdr),n_elements(bins_kdp)) obs_counts4d=lonarr(n_elements(bins_dbz),n_elements(bins_zdr),n_elements(bins_kdp),n_elements(bins_hgt)) for i=0,count_valid-1 do begin if i mod 10000 eq 0 then print,'i',i ddbz=min(abs(dbz[r_valid[i]]-bins_dbz),dbz_idx) ;print,dbz[r_valid[i]],bins_dbz[dbz_idx],ddbz if radar_name eq 'nosmn_cband' then begin ;can you rerun the somna counts files and add 1db to Zdr? dzdr=min(abs((zdr[r_valid[i]]+1)-bins_zdr),zdr_idx) endif else begin dzdr=min(abs(zdr[r_valid[i]]-bins_zdr),zdr_idx) endelse dkdp=min(abs(kdp[r_valid[i]]-bins_kdp),kdp_idx) ;if dkdp ge 0.5 then print,kdp[r_valid[i]],bins_kdp[kdp_idx],dkdp close=min(abs(point_z[r_valid[i]]-bins_hgt),hgt_idx) if dkdp lt 0.1 then begin obs_counts3d[dbz_idx,zdr_idx,kdp_idx]=obs_counts3d[dbz_idx,zdr_idx,kdp_idx]+1L obs_counts4d[dbz_idx,zdr_idx,kdp_idx,hgt_idx]=obs_counts4d[dbz_idx,zdr_idx,kdp_idx,hgt_idx]+1L endif endfor print,'end histo' ; Filetag if radar_name eq 'sband' then file_tag=strmid(file_basename(filename),0,35) if radar_name eq 'cband' then file_tag=strmid(file_basename(filename),0,44) if radar_name eq 'nosmn_cband' then file_tag=strmid(file_basename(filename),0,75) ; Write to a netcdf file if 1 eq 1 then begin output_filename=filedir+file_tag+'_stats_obs_histo.rho_filter.hires' if rtag ne '' then output_filename=output_filename+'.'+rtag output_filename=output_filename+'.cdf' ;output_filename=filedir+file_tag+'_stats_obs_histo.cdf' fid=ncdf_create(output_filename,/clobber) print, output_filename,fid,'start write' bins_dbz_did=ncdf_dimdef(fid,'num_bins_dbz',n_elements(bins_dbz)) bins_zdr_did=ncdf_dimdef(fid,'num_bins_zdr',n_elements(bins_zdr)) bins_kdp_did=ncdf_dimdef(fid,'num_bins_kdp',n_elements(bins_kdp)) bins_hgt_did=ncdf_dimdef(fid,'num_bins_hgt',n_elements(bins_hgt)) bins_dbz_id=ncdf_vardef(fid,'bins_dbz',bins_dbz_did,/float) ncdf_attput,fid,bins_dbz_id,'long_name','dbz bins' ncdf_attput,fid,bins_dbz_id,'units','dbz' bins_zdr_id=ncdf_vardef(fid,'bins_zdr',bins_zdr_did,/float) ncdf_attput,fid,bins_zdr_id,'long_name','zdr bins' ncdf_attput,fid,bins_zdr_id,'units','zdr' bins_kdp_id=ncdf_vardef(fid,'bins_kdp',bins_kdp_did,/float) ncdf_attput,fid,bins_kdp_id,'long_name','kdp bins' ncdf_attput,fid,bins_kdp_id,'units','kdp' bins_hgt_id=ncdf_vardef(fid,'bins_hgt',bins_hgt_did,/float) ncdf_attput,fid,bins_hgt_id,'long_name','height bins' ncdf_attput,fid,bins_hgt_id,'units','km' hist_id=ncdf_vardef(fid,'counts_array',[bins_dbz_did,bins_zdr_did,bins_kdp_did,bins_hgt_did],/long) ncdf_attput,fid,hist_id,'long_name','counts of obs' ncdf_control,fid,/endef ncdf_varput,fid,bins_dbz_id,bins_dbz ncdf_varput,fid,bins_zdr_id,bins_zdr ncdf_varput,fid,bins_kdp_id,bins_kdp ncdf_varput,fid,bins_hgt_id,bins_hgt ncdf_varput,fid,hist_id,obs_counts4d ncdf_close,fid print,'end write' endif ; turn of write ;************************************** ; Make the plot ;************************************** if 1 eq 1 then begin ; Set up plot size pxdim=1000 & pydim=1000 ; Position the plots xl=0.1 & xr=0.95 yb=0.06 & yt=0.85 sx=0.11 sy=0.15 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 ;vertical ;cbpos[*,2]=cbpos[*,0]+0.01 cbpos[*,1]=cbpos[*,3]+0.05 ;horizontal cbpos[*,3]=cbpos[*,1]+0.01 ; Top is the last color to scale 256 colors, 0-255 top_color=252 ; Colortable 0-252 253=white ;mytable=colortable(39,ncolors=254) mytable=colortable(43,ncolors=253) ;254=hot pink ;gray=255 mytable=[mytable,transpose([255,255,255]),transpose([238,18,137]),transpose([230,230,230])] mycbtable=mytable[0:top_color,*] ;top_color=254 ;mytable=colortable(39,ncolors=255) ;mytable[0,*]=[250,250,250]; first color =gray ;mycbtable=mytable[1:-1,*] ; Strings 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]) ;obs_counts4d=lonarr(bins_dbz,bins_zdr,bins_kdp,bins_hgt) ;************************** ; hold dbz constant ;************************** pnum=0 ;r_sub=where(counts_dbz gt 0,c_sub) ;1d vector ;r_sub=where(bins_dbz eq 10,c_sub) ;10 dbz maxval=max(counts_dbz,r_sub) ;dbz with max counts data_var=reform(obs_counts3d[r_sub,*,*]) dmin=min(data_var) dmax=max(data_var) data_image=bytscl(data_var,top=top_color,min=dmin,max=dmax) 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 p4=image(data_image,/current,image_dimensions=[isx,isy],$ ;min_value=dmin_cband,max_value=dmax_cband,$ position=pos[pnum,*],rgb_table=mytable) c1=contour(data_var,bins_zdr,bins_kdp,/nodata,/current,$ ytitle='kdp',xtitle='zdr',irregular=0,xtickdir=1,ytickdir=1,$ xstyle=1,ystyle=1,position=pos[pnum,*],font_size=fs1,axis_style=2) cb0=colorbar(rgb_table=mycbtable,range=[dmin,dmax],orientation=0,$ position=cbpos[pnum,*],font_size=fs1,tickdir=1,border_on=1,title='counts') t1=text(pos[pnum,0]+dx,cbpos[pnum,3]+dy,'Constant DBZ '+string(bins_dbz[r_sub],format='(F6.2)'),font_size=fs1) pnum=1 ;ridx=where(bins_dbz eq 20,c_sub) ;20 dbz if r_sub+1 ge n_elements(bins_dbz) then ridx=r_sub-1 else ridx=r_sub+1 data_var=reform(obs_counts3d[ridx,*,*]) dmin=min(data_var) dmax=max(data_var) data_image=bytscl(data_var,top=top_color,min=dmin,max=dmax) 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 p4=image(data_image,/current,image_dimensions=[isx,isy],$ ;min_value=dmin_cband,max_value=dmax_cband,$ position=pos[pnum,*],rgb_table=mytable) c1=contour(data_var,bins_zdr,bins_kdp,/nodata,/current,$ ytitle='kdp',xtitle='zdr',irregular=0,xtickdir=1,ytickdir=1,$ xstyle=1,ystyle=1,position=pos[pnum,*],font_size=fs1,axis_style=2) cb0=colorbar(rgb_table=mycbtable,range=[dmin,dmax],orientation=0,$ position=cbpos[pnum,*],font_size=fs1,tickdir=1,border_on=1,title='counts') t1=text(pos[pnum,0]+dx,cbpos[pnum,3]+dy,'Constant DBZ '+string(bins_dbz[ridx],format='(F6.2)'),font_size=fs1) pnum=2 ;ridx=where(bins_dbz eq 25,c_sub) ;25 dbz if r_sub+2 ge n_elements(bins_dbz) then ridx=r_sub-2 else ridx=r_sub+2 data_var=reform(obs_counts3d[ridx,*,*]) dmin=min(data_var) dmax=max(data_var) data_image=bytscl(data_var,top=top_color,min=dmin,max=dmax) 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 p4=image(data_image,/current,image_dimensions=[isx,isy],$ ;min_value=dmin_cband,max_value=dmax_cband,$ position=pos[pnum,*],rgb_table=mytable) c1=contour(data_var,bins_zdr,bins_kdp,/nodata,/current,$ ytitle='kdp',xtitle='zdr',irregular=0,xtickdir=1,ytickdir=1,$ xstyle=1,ystyle=1,position=pos[pnum,*],font_size=fs1,axis_style=2) cb0=colorbar(rgb_table=mycbtable,range=[dmin,dmax],orientation=0,$ position=cbpos[pnum,*],font_size=fs1,tickdir=1,border_on=1,title='counts') t1=text(pos[pnum,0]+dx,cbpos[pnum,3]+dy,'Constant DBZ '+string(bins_dbz[ridx],format='(F6.2)'),font_size=fs1) pnum=3 ;maxval=max(counts_dbz,ridx) ;dbz with max counts if r_sub+4 ge n_elements(bins_dbz) then ridx=r_sub-4 else ridx=r_sub+4 data_var=reform(obs_counts3d[ridx,*,*]) dmin=min(data_var) dmax=max(data_var) data_image=bytscl(data_var,top=top_color,min=dmin,max=dmax) 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 p4=image(data_image,/current,image_dimensions=[isx,isy],$ ;min_value=dmin_cband,max_value=dmax_cband,$ position=pos[pnum,*],rgb_table=mytable) c1=contour(data_var,bins_zdr,bins_kdp,/nodata,/current,$ ytitle='kdp',xtitle='zdr',irregular=0,xtickdir=1,ytickdir=1,$ xstyle=1,ystyle=1,position=pos[pnum,*],font_size=fs1,axis_style=2) cb0=colorbar(rgb_table=mycbtable,range=[dmin,dmax],orientation=0,$ position=cbpos[pnum,*],font_size=fs1,tickdir=1,border_on=1,title='counts') t1=text(pos[pnum,0]+dx,cbpos[pnum,3]+dy,'Constant DBZ '+string(bins_dbz[ridx],format='(F6.2)'),font_size=fs1) t1=text(0.1,0.97,radar_name+' '+date_str+' '+shh+':'+smi+':'+sss,font_size=fs1) ;imagename=filedir+file_tag+'_stats_dbz.rho_filter.10dbz' imagename=filedir+file_tag+'_stats_dbz.rho_filter.hires' if rtag ne '' then imagename=imagename+'.'+rtag print,imagename ;p0.save,'test.png' p0.save,imagename+'.png' endif ;end of turn off plot ;************************************************* ; Do some height slicing ;************************************************* if 1 eq 1 then begin pnum=0 const_hgt=1;2 h=where(bins_hgt ge const_hgt) h=h[0] r_sub=where(bins_dbz eq 10,c_sub) ;10dbz ;maxval=max(counts_dbz,r_sub) data_var=reform(obs_counts4d[r_sub,*,*,h]) dmin=min(data_var) dmax=max(data_var) data_image=bytscl(data_var,top=top_color,min=dmin,max=dmax) 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 p4=image(data_image,/buffer,image_dimensions=[isx,isy],$ position=pos[pnum,*],rgb_table=mytable,dimensions=[pxdim,pydim]) c1=contour(data_var,bins_zdr,bins_kdp,/nodata,/current,$ ytitle='kdp',xtitle='zdr',irregular=0,xtickdir=1,ytickdir=1,$ xstyle=1,ystyle=1,position=pos[pnum,*],font_size=fs1,axis_style=2) cb0=colorbar(rgb_table=mycbtable,range=[dmin,dmax],orientation=0,$ position=cbpos[pnum,*],font_size=fs1,tickdir=1,border_on=1,title='counts') t1=text(pos[pnum,0]+dx,cbpos[pnum,3]+dy,'Constant DBZ '+string(bins_dbz[r_sub],format='(F6.2)')+$ ' Constant Height '+string(bins_hgt[h],format='(F6.2)')+'km',font_size=fs1) pnum=1 const_hgt=1;3 h=where(bins_hgt ge const_hgt) h=h[0] r_sub=where(bins_dbz eq 20,c_sub) data_var=reform(obs_counts4d[r_sub,*,*,h]) dmin=min(data_var) dmax=max(data_var) data_image=bytscl(data_var,top=top_color,min=dmin,max=dmax) 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 p4=image(data_image,/current,image_dimensions=[isx,isy],$ position=pos[pnum,*],rgb_table=mytable) c1=contour(data_var,bins_zdr,bins_kdp,/nodata,/current,$ ytitle='kdp',xtitle='zdr',irregular=0,xtickdir=1,ytickdir=1,$ xstyle=1,ystyle=1,position=pos[pnum,*],font_size=fs1,axis_style=2) cb0=colorbar(rgb_table=mycbtable,range=[dmin,dmax],orientation=0,$ position=cbpos[pnum,*],font_size=fs1,tickdir=1,border_on=1,title='counts') t1=text(pos[pnum,0]+dx,cbpos[pnum,3]+dy,'Constant DBZ '+string(bins_dbz[r_sub],format='(F6.2)')+$ ' Constant Height '+string(bins_hgt[h],format='(F6.2)')+'km',font_size=fs1) pnum=2 const_hgt=1;4 h=where(bins_hgt ge const_hgt) h=h[0] r_sub=where(bins_dbz eq 25,c_sub) data_var=reform(obs_counts4d[r_sub,*,*,h]) dmin=min(data_var) dmax=max(data_var) data_image=bytscl(data_var,top=top_color,min=dmin,max=dmax) 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 p4=image(data_image,/current,image_dimensions=[isx,isy],$ position=pos[pnum,*],rgb_table=mytable) c1=contour(data_var,bins_zdr,bins_kdp,/nodata,/current,$ ytitle='kdp',xtitle='zdr',irregular=0,xtickdir=1,ytickdir=1,$ xstyle=1,ystyle=1,position=pos[pnum,*],font_size=fs1,axis_style=2) cb0=colorbar(rgb_table=mycbtable,range=[dmin,dmax],orientation=0,$ position=cbpos[pnum,*],font_size=fs1,tickdir=1,border_on=1,title='counts') t1=text(pos[pnum,0]+dx,cbpos[pnum,3]+dy,'Constant DBZ '+string(bins_dbz[r_sub],format='(F6.2)')+$ ' Constant Height '+string(bins_hgt[h],format='(F6.2)')+'km',font_size=fs1) pnum=3 const_hgt=1;5 h=where(bins_hgt ge const_hgt) h=h[0] maxval=max(counts_dbz,r_sub) data_var=reform(obs_counts4d[r_sub,*,*,h]) dmin=min(data_var) dmax=max(data_var) data_image=bytscl(data_var,top=top_color,min=dmin,max=dmax) 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 p4=image(data_image,/current,image_dimensions=[isx,isy],$ position=pos[pnum,*],rgb_table=mytable) c1=contour(data_var,bins_zdr,bins_kdp,/nodata,/current,$ ytitle='kdp',xtitle='zdr',irregular=0,xtickdir=1,ytickdir=1,$ xstyle=1,ystyle=1,position=pos[pnum,*],font_size=fs1,axis_style=2) cb0=colorbar(rgb_table=mycbtable,range=[dmin,dmax],orientation=0,$ position=cbpos[pnum,*],font_size=fs1,tickdir=1,border_on=1,title='counts') t1=text(pos[pnum,0]+dx,cbpos[pnum,3]+dy,'Constant DBZ '+string(bins_dbz[r_sub],format='(F6.2)')+$ ' Constant Height '+string(bins_hgt[h],format='(F6.2)')+'km',font_size=fs1) t1=text(0.1,0.97,radar_name+' '+date_str+' '+shh+':'+smi+':'+sss,font_size=fs1) imagename=filedir+file_tag+'_stats_height.rho_filter.'+string(const_hgt,format='(I1)')+'km.hires' if rtag ne '' then imagename=imagename+'.'+rtag print,imagename ;p4.save,'test.png' p4.save,imagename+'.png' endif endfor ;end of loop through cband files end