;******************************************** ; This plots a plan view (horizontal slice on a map) ; Constant Latitude slice, constant longitude slice, and line segment slice ; 4 plot. ;************************************************** ;function MapGrid_Labels, orientation, location, fractional, defaultlabel ; ;if (location eq 0) then $ ; ; return, orientation ? 'Equator' : 'Prime Meridian' ; degree = '!M' + STRING(176b) ; Use the Math symbol ; label = STRTRIM(ROUND(ABS(location)),2) + degree ; label=locaton ; suffix = orientation ? ((location lt 0) ? 'S' : 'N') : $ ; ((location lt 0) ? 'W' : 'E') ;return, label + suffix ;***************************** ; Function to plot a circle ;***************************** ;function circle, xcenter,ycenter,radius ; points = (2 * !PI / 99.0) * FINDGEN(100) ; x = xcenter + radius * COS(points ) ; y = ycenter + radius * SIN(points ) ; RETURN, TRANSPOSE([[x],[y]]) ;end ;***************************** ; Program to plot cband radar ;***************************** pro cband_12plot ; Choose date and time of cband ;date_str='20190126' & time_str='063231' ;date_str='20190128' & time_str='045515' ;date_str='20190129' & time_str='035625' ;date_str='20180128' & time_str='0*' ;date_str='20180209' & time_str='105330' ;date_str='20151214' & time_str='*';'163711' ;olympex ;date_str='20200313' & time_str='0*' date_str='20180209' & time_str='121923' rho_filter=0 ;no ;rho_filter=1 ;yes ; 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='/uufs/chpc.utah.edu/common/home/' path_prefix='/Volumes/' ; directory 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*.nc' 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' then begin filedir=path_prefix+'mace-group6/field_programs/olympex/sband/'+strmid(date_str,4,4)+'/ppi/' filestr='olympex_NPOL1_'+date_str+'_'+time_str+'*pyart*.nc' ; Get the flight path nav_date_str='20151213' read_nav_olympex,path_prefix,nav_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 else if syy eq '2020' then begin filedir=path_prefix+'mace-group5/field_programs/comble/cband/'+date_str+'/' filestr=syy+smm+sdd+time_str+'*pyart*.nc' endif ; File search string ;filestr='9776HUB-PPIVol-'+syy+smm+sdd+'-'+shh+'*pyart.41.201.nc' ;filestr='9776HUB-PPIVol-'+syy+smm+sdd+'*81.401.nc' ;filestr='9776HUB-PPIVol-*81.401.nc' ;filestr='9776HUB-PPIVol-'+syy+smm+sdd+'*.pyart.*.nc' ;filestr='9776HUB-PPIVol-'+syy+smm+sdd+'-'+time_str+'*pyart*.nc' ; file name filenames=file_search(filedir+filestr,count=numfiles) print,'numfiles',numfiles ; Loop through the cband files for idx=0,numfiles-1 do begin filename=filenames[idx] print,filename if syy eq '2020' then begin shh=strmid(file_basename(filename),8,2) smi=strmid(file_basename(filename),10,2) sss=strmid(file_basename(filename),12,2) endif else 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) endelse time_str=shh+smi+sss 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_time') & ncdf_varget,fid,vid,radar_time ;ncdf_attget,fid,vid,'units',astr & astr=string(astr) 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') & if vid ne -1 then 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') if vid ne -1 then ncdf_varget,fid,vid,kdp vid=ncdf_varid(fid,'corrected_differential_reflectivity') & if vid ne -1 then ncdf_varget,fid,vid,zdr vid=ncdf_varid(fid,'corrected_reflectivity') if vid eq -1 then vid=ncdf_varid(fid,'DBZH_gmm') if vid eq -1 then vid=ncdf_varid(fid,'reflectivity') if vid ne -1 then ncdf_varget,fid,vid,dbz ncdf_close,fid ; 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 ;time_offset=0L 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 ; Zoom in on radar, dist_km amount (max distance is 123) dist_km=150.0 ;75 x1=point_x[*,0,0] xr=where(x1 ge -1*dist_km and x1 le dist_km,xc) y1=point_y[0,*,0] yr=where(y1 ge -1*dist_km and y1 le dist_km,yc) if dbz ne !NULL then begin dbz=dbz[xr,*,*] dbz=dbz[*,yr,*] lat=lat[xr,*,*] lon=lon[xr,*,*] lat=lat[*,yr,*] lon=lon[*,yr,*] endif if kdp ne !NULL then begin kdp=kdp[xr,*,*] zdr=zdr[xr,*,*] rho=rho[xr,*,*] kdp=kdp[*,yr,*] zdr=zdr[*,yr,*] rho=rho[*,yr,*] endif ; Height array height_cband=alt[0,0,*]/1000.0 ;************************************* ; Find flight track at this time ;************************************* if julian_day_plane ne !NULL then begin xtime=5D/(60D*24D) r=where(julian_day_plane ge julian_day-xtime and julian_day_plane le julian_day+xtime and $ lat_plane ge min(lat) and lat_plane le max(lat) and $ lon_plane ge min(lon) and lon_plane le max(lon) ,c) if c gt 1 then begin sub_lat_plane=lat_plane[r] sub_lon_plane=lon_plane[r] sub_alt_plane=alt_plane[r]/1000.0 print,'found track' endif else begin print,'no flight track' endelse endif ; Filetag parts=strsplit(file_basename(filename),'.',/extract) file_tag=parts[0] ;file_tag=strmid(file_basename(filename),0,19) ; Set up plot size pxdim=1000 & pydim=1000 ; Position the plots xl=0.07 & xr=0.93 yb=0.08 & yt=0.92 sx=0.00 sy=0.00 numplots_x=4 numplots_y=4 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.005 ;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,*] mytable2=colortable(72,ncolors=253) ;254=hot pink ;gray=255 mytable2=[mytable2,transpose([255,255,255]),transpose([238,18,137]),transpose([230,230,230])] mycbtable2=mytable2[0:top_color,*] ;mytable3=colortable(72,ncolors=10) ;white11 ;gray12 ;mytable3=[mytable3,transpose([255,255,255]),transpose([230,230,230])] ; light or blue yel mytable3=[[165,0,38],[241,102,64],[254,212,133],[235,247,228],[175,219,234],[108,164,204],[61,90,167],$ [100,100,100],[255,255,255],[230,230,230]] mycbtable3=mytable3[*,0:7] ;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]) ; radar min and max ;dmin_cband=-17 & dmax_cband=31 ;dmin_cband=-5 & dmax_cband=25 ;r=where(kdp ne -9999,c) ;dmin_cband=min(kdp[r]) ;dmax_cband=max(kdp[r]) ;r=where(reflect ne -9999,c) ;dmin_cband=min(reflect[r]) ;dmax_cband=max(reflect[r]) ;plot_var=reflect ;var_name='reflect' ;radar_name='SBAND' ;*********************************** ; Row 1 is dbz ;*********************************** ;r=where(dbz ne -9999,c) dmin_cband=5;min(dbz[r]) dmax_cband=50;max(dbz[r]) ; 0.5 km pnum=0 r=where(height_cband ge 0.5,c) sidx=r[0] data_var=reform(dbz[*,*,sidx]) lat_var=reform(lat[*,*,sidx]) lon_var=reform(lon[*,*,sidx]) if kdp ne !NULL then kdp_var=reform(kdp[*,*,sidx]) llat=min(lat) & ulat=max(lat) llon=min(lon) & rlon=max(lon) data_image=bytscl(data_var,top=top_color,min=dmin_cband,max=dmax_cband) r=where(data_var eq -9999,c) if c gt 0 then data_image[r]=253 ;white if kdp ne !NULL then begin r=where(kdp_var eq -9999,c) if c gt 0 then data_image[r]=255 ;gray endif p1=image(data_image,lon_var,lat_var,limit=[llat,llon,ulat,rlon],$ /current,dimensions=[pxdim,pydim],grid_units='degrees',$ rgb_table=mytable,linestyle='dot',$ label_format='MapGrid_Labels',/box_axes,label_position=0,$ label_show=1,irregular=0,$ map_projection='Mercator',position=pos[pnum,*],font_size=9) p2=mapcontinents(color='black',/hires) if sub_lon_plane ne !NULL then p7=plot(sub_lon_plane,sub_lat_plane,color='blue',thick=2,/overplot) ; calculate range rings around ship location (25km) lon_cir=fltarr(4,360) lat_cir=fltarr(4,360) for r=1,4 do begin arc_dist=(float(r)*25000.0)/(6371D*1000D) ;m for s=0,359 do begin result=ll_arc_distance([lon_radar,lat_radar],arc_dist,s,/degrees) lon_cir[r-1,s]=result[0] lat_cir[r-1,s]=result[1] endfor endfor for r=0,3 do begin p7=plot(lon_cir[r,*],lat_cir[r,*],/current,/overplot,color='gray') endfor ;t1=text(pos[pnum,0]-5*dx,pos[pnum,3]+2*dy,$ t1=text(pos[12,0]+3*dx,pos[12,3]+3*dy,$ string(height_cband[sidx],format='(F4.2)')+'km',font_size=fs1) ; 1.0 km pnum=1 r=where(height_cband ge 1.0,c) sidx=r[0] data_var=reform(dbz[*,*,sidx]) lat_var=reform(lat[*,*,sidx]) lon_var=reform(lon[*,*,sidx]) if kdp ne !NULL then kdp_var=reform(kdp[*,*,sidx]) llat=min(lat) & ulat=max(lat) llon=min(lon) & rlon=max(lon) data_image=bytscl(data_var,top=top_color,min=dmin_cband,max=dmax_cband) r=where(data_var eq -9999,c) if c gt 0 then data_image[r]=253 ;white if kdp ne !NULL then begin r=where(kdp_var eq -9999,c) if c gt 0 then data_image[r]=255 ;gray endif p1=image(data_image,lon_var,lat_var,limit=[llat,llon,ulat,rlon],$ /current,dimensions=[pxdim,pydim],grid_units='degrees',$ rgb_table=mytable,linestyle='dot',$ label_format='MapGrid_Labels',/box_axes,label_position=0,$ label_show=0,irregular=0,$ map_projection='Mercator',position=pos[pnum,*],font_size=fs1) p2=mapcontinents(color='black',/hires) if sub_lon_plane ne !NULL then p7=plot(sub_lon_plane,sub_lat_plane,color='blue',thick=2,/overplot) ;t1=text(pos[pnum,0]-5*dx,pos[pnum,3]+2*dy,$ t1=text(pos[13,0]+3*dx,pos[13,3]+3*dy,$ string(height_cband[sidx],format='(F4.2)')+'km',font_size=fs1) ; 3.0 km pnum=2 r=where(height_cband ge 3.0,c) sidx=r[0] data_var=reform(dbz[*,*,sidx]) lat_var=reform(lat[*,*,sidx]) lon_var=reform(lon[*,*,sidx]) if kdp ne !NULL then kdp_var=reform(kdp[*,*,sidx]) llat=min(lat) & ulat=max(lat) llon=min(lon) & rlon=max(lon) data_image=bytscl(data_var,top=top_color,min=dmin_cband,max=dmax_cband) r=where(data_var eq -9999,c) if c gt 0 then data_image[r]=253 ;white if kdp ne !NULL then begin r=where(kdp_var eq -9999,c) if c gt 0 then data_image[r]=255 ;gray endif p1=image(data_image,lon_var,lat_var,limit=[llat,llon,ulat,rlon],$ /current,dimensions=[pxdim,pydim],grid_units='degrees',$ rgb_table=mytable,linestyle='dot',$ label_format='MapGrid_Labels',/box_axes,label_position=0,$ label_show=0,irregular=0,$ map_projection='Mercator',position=pos[pnum,*],font_size=fs1) p2=mapcontinents(color='black',/hires) if sub_lon_plane ne !NULL then p7=plot(sub_lon_plane,sub_lat_plane,color='blue',thick=2,/overplot) ;t1=text(pos[pnum,0]-5*dx,pos[pnum,3]+2*dy,$ t1=text(pos[14,0]+3*dx,pos[14,3]+3*dy,$ string(height_cband[sidx],format='(F4.2)')+'km',font_size=fs1) ; 5.0 km pnum=3 r=where(height_cband ge 5.0,c) sidx=r[0] data_var=reform(dbz[*,*,sidx]) lat_var=reform(lat[*,*,sidx]) lon_var=reform(lon[*,*,sidx]) llat=min(lat) & ulat=max(lat) llon=min(lon) & rlon=max(lon) if kdp ne !NULL then kdp_var=reform(kdp[*,*,sidx]) data_image=bytscl(data_var,top=top_color,min=dmin_cband,max=dmax_cband) r=where(data_var eq -9999,c) if c gt 0 then data_image[r]=253 ;white if kdp ne !NULL then begin r=where(kdp_var eq -9999,c) if c gt 0 then data_image[r]=255 ;gray endif p1=image(data_image,lon_var,lat_var,limit=[llat,llon,ulat,rlon],$ /current,dimensions=[pxdim,pydim],grid_units='degrees',$ rgb_table=mytable,linestyle='dot',$ label_format='MapGrid_Labels',/box_axes,label_position=0,$ label_show=0,irregular=0,$ map_projection='Mercator',position=pos[pnum,*],font_size=fs1) p2=mapcontinents(color='black',/hires) if sub_lon_plane ne !NULL then p7=plot(sub_lon_plane,sub_lat_plane,color='blue',thick=2,/overplot) ;t1=text(pos[pnum,0]-5*dx,pos[pnum,3]+2*dy,$ t1=text(pos[15,0]+3*dx,pos[15,3]+3*dy,$ string(height_cband[sidx],format='(F4.2)')+'km',font_size=fs1) ; Colorbar - horizontal ;cb0=colorbar(rgb_table=mycbtable,range=[dmin_cband,dmax_cband],$ ; orientation=0,position=[pos[1,0],pos[1,3]+0.5*dy,pos[2,2],pos[2,3]+1.5*dy],font_size=fs1,$ ; tickdir=1,border_on=1);,title='dBZ') ;t1=text(pos[1,0]-5*dx,pos[1,3]-1*dy,'dBZ',font_size=fs1) ; Colorbar - vertical cb0=colorbar(rgb_table=mycbtable,range=[dmin_cband,dmax_cband],$ orientation=1,$ position=[pos[3,2]+4.5*dx,pos[3,1],pos[3,2]+5.5*dx,pos[3,3]-3*dy],font_size=fs1,$ tickdir=1,border_on=1);,title='dBZ') t1=text(pos[3,2]+3*dx,pos[3,3]-2*dy,'dBZ',font_size=fs1) ;*********************************** ; Row 2 is zdr ;*********************************** if zdr ne !NULL then begin ;r=where(dbz ne -9999,c) dmin_zdr=-1;min(dbz[r]) dmax_zdr=1;max(dbz[r]) ; 0.5 km pnum=4 r=where(height_cband ge 0.5,c) sidx=r[0] data_var=reform(zdr[*,*,sidx]) lat_var=reform(lat[*,*,sidx]) lon_var=reform(lon[*,*,sidx]) kdp_var=reform(kdp[*,*,sidx]) llat=min(lat) & ulat=max(lat) llon=min(lon) & rlon=max(lon) data_image=bytscl(data_var,top=top_color,min=dmin_zdr,max=dmax_zdr) r=where(finite(data_var,/NAN),c) if c gt 0 then data_image[r]=253 ;white r=where(kdp_var eq -9999,c) if c gt 0 then data_image[r]=255 ;gray p1=image(data_image,lon_var,lat_var,limit=[llat,llon,ulat,rlon],$ /current,dimensions=[pxdim,pydim],grid_units='degrees',$ rgb_table=mytable2,linestyle='dot',$ label_format='MapGrid_Labels',/box_axes,label_position=0,$ label_show=0,irregular=0,$ map_projection='Mercator',position=pos[pnum,*],font_size=fs1) p2=mapcontinents(color='black',/hires) if sub_lon_plane ne !NULL then p7=plot(sub_lon_plane,sub_lat_plane,color='blue',thick=2,/overplot) ;t1=text(pos[pnum,0]-5*dx,pos[pnum,3]+2*dy,string(height_cband[sidx],format='(F4.2)')+'km',$ ; font_size=fs1) ; 1.0 km pnum=5 r=where(height_cband ge 1.0,c) sidx=r[0] data_var=reform(zdr[*,*,sidx]) lat_var=reform(lat[*,*,sidx]) lon_var=reform(lon[*,*,sidx]) kdp_var=reform(kdp[*,*,sidx]) llat=min(lat) & ulat=max(lat) llon=min(lon) & rlon=max(lon) data_image=bytscl(data_var,top=top_color,min=dmin_zdr,max=dmax_zdr) r=where(finite(data_var,/NAN),c) if c gt 0 then data_image[r]=253 ;white r=where(kdp_var eq -9999,c) if c gt 0 then data_image[r]=255 ;gray p1=image(data_image,lon_var,lat_var,limit=[llat,llon,ulat,rlon],$ /current,dimensions=[pxdim,pydim],grid_units='degrees',$ rgb_table=mytable2,linestyle='dot',$ label_format='MapGrid_Labels',/box_axes,label_position=0,$ label_show=0,irregular=0,$ map_projection='Mercator',position=pos[pnum,*],font_size=fs1) p2=mapcontinents(color='black',/hires) if sub_lon_plane ne !NULL then p7=plot(sub_lon_plane,sub_lat_plane,color='blue',thick=2,/overplot) ;t1=text(pos[pnum,0]-5*dx,pos[pnum,3]+2*dy,string(height_cband[sidx],format='(F4.2)')+'km',$ ; font_size=fs1) ; 3.0 km pnum=6 r=where(height_cband ge 3.0,c) sidx=r[0] data_var=reform(zdr[*,*,sidx]) lat_var=reform(lat[*,*,sidx]) lon_var=reform(lon[*,*,sidx]) kdp_var=reform(kdp[*,*,sidx]) llat=min(lat) & ulat=max(lat) llon=min(lon) & rlon=max(lon) data_image=bytscl(data_var,top=top_color,min=dmin_zdr,max=dmax_zdr) r=where(finite(data_var,/NAN),c) if c gt 0 then data_image[r]=253 ;white r=where(kdp_var eq -9999,c) if c gt 0 then data_image[r]=255 ;gray p1=image(data_image,lon_var,lat_var,limit=[llat,llon,ulat,rlon],$ /current,dimensions=[pxdim,pydim],grid_units='degrees',$ rgb_table=mytable2,linestyle='dot',$ label_format='MapGrid_Labels',/box_axes,label_position=0,$ label_show=0,irregular=0,$ map_projection='Mercator',position=pos[pnum,*],font_size=fs1) p2=mapcontinents(color='black',/hires) if sub_lon_plane ne !NULL then p7=plot(sub_lon_plane,sub_lat_plane,color='blue',thick=2,/overplot) ;t1=text(pos[pnum,0]-5*dx,pos[pnum,3]+2*dy,string(height_cband[sidx],format='(F4.2)')+'km',$ ; font_size=fs1) ; 5.0 km pnum=7 r=where(height_cband ge 5.0,c) sidx=r[0] data_var=reform(zdr[*,*,sidx]) lat_var=reform(lat[*,*,sidx]) lon_var=reform(lon[*,*,sidx]) kdp_var=reform(kdp[*,*,sidx]) llat=min(lat) & ulat=max(lat) llon=min(lon) & rlon=max(lon) data_image=bytscl(data_var,top=top_color,min=dmin_zdr,max=dmax_zdr) r=where(finite(data_var,/NAN),c) if c gt 0 then data_image[r]=253 ;white r=where(kdp_var eq -9999,c) if c gt 0 then data_image[r]=255 ;gray p1=image(data_image,lon_var,lat_var,limit=[llat,llon,ulat,rlon],$ /current,dimensions=[pxdim,pydim],grid_units='degrees',$ rgb_table=mytable2,linestyle='dot',$ label_format='MapGrid_Labels',/box_axes,label_position=0,$ label_show=0,irregular=0,$ map_projection='Mercator',position=pos[pnum,*],font_size=fs1) p2=mapcontinents(color='black',/hires) if sub_lon_plane ne !NULL then p7=plot(sub_lon_plane,sub_lat_plane,color='blue',thick=2,/overplot) ;t1=text(pos[pnum,0]-5*dx,pos[pnum,3]+2*dy,string(height_cband[sidx],format='(F4.2)')+'km',$ ; font_size=fs1) ; Colorbar ;cb0=colorbar(rgb_table=mycbtable2,range=[dmin_zdr,dmax_zdr],$ ; orientation=0,position=[pos[5,0],pos[5,3]+0.5*dy,pos[6,2],pos[6,3]+1.5*dy],font_size=fs1,$ ; tickdir=1,border_on=1);,title='dBZ') ;t1=text(pos[5,0]-5*dx,pos[5,3]-1*dy,'ZDR',font_size=fs1) ; Colorbar - vertical cb0=colorbar(rgb_table=mycbtable2,range=[dmin_zdr,dmax_zdr],$ orientation=1,$ position=[pos[7,2]+4.5*dx,pos[7,1],pos[7,2]+5.5*dx,pos[7,3]-3*dy],font_size=fs1,$ tickdir=1,border_on=1);,title='dBZ') t1=text(pos[7,2]+3*dx,pos[7,3]-2*dy,'ZDR',font_size=fs1) endif ;*********************************** ; Row 3 is kdp ;*********************************** if kdp ne !NULL then begin ; light or blue yel ;mytable3=[[165,0,38],[241,102,64],[254,212,133],[235,247,228],[175,219,234],[108,164,204],[61,90,167],$ ; [255,255,255],[230,230,230]] kdp_int=[-100,-0.1,-0.05,0,0.05,0.1,1.0,100.0] lrho=0.95 dmin_kdp=-0.2;min(dbz[r]) dmax_kdp=0.2;max(dbz[r]) ; 0.5 km pnum=8 r=where(height_cband ge 0.5,c) sidx=r[0] data_var=reform(kdp[*,*,sidx]) lat_var=reform(lat[*,*,sidx]) lon_var=reform(lon[*,*,sidx]) data_rho=reform(rho[*,*,sidx]) llat=min(lat) & ulat=max(lat) llon=min(lon) & rlon=max(lon) ;data_image=bytscl(data_var,top=top_color,min=dmin_kdp,max=dmax_kdp) data_image=byte(data_var) & data_image[*]=0B for i=0,6 do begin r=where(data_var ge kdp_int[i] and data_var lt kdp_int[i+1],c) if c gt 0 then data_image[r]=i print,kdp_int[i],kdp_int[i+1],c endfor r=where(data_var eq 0,c) if c gt 0 then data_image[r]=8;253 ;white r=where(data_var eq -9999,c) if c gt 0 then data_image[r]=9;255 ;gray if rho_filter eq 1 then begin r=where(data_rho lt lrho and data_var ne 0 and data_var ne -9999,c) if c gt 0 then data_image[r]=7;254 ;pink endif p1=image(data_image,lon_var,lat_var,limit=[llat,llon,ulat,rlon],$ /current,dimensions=[pxdim,pydim],grid_units='degrees',$ rgb_table=mytable3,linestyle='dot',$ label_format='MapGrid_Labels',/box_axes,label_position=0,$ label_show=0,irregular=0,$ map_projection='Mercator',position=pos[pnum,*],font_size=fs1) p2=mapcontinents(color='black',/hires) if sub_lon_plane ne !NULL then p7=plot(sub_lon_plane,sub_lat_plane,color='blue',thick=2,/overplot) ;t1=text(pos[pnum,0]-5*dx,pos[pnum,3]+2*dy,string(height_cband[sidx],format='(F4.2)')+'km',$ ; font_size=fs1) ; 1.0 km pnum=9 r=where(height_cband ge 1.0,c) sidx=r[0] data_var=reform(kdp[*,*,sidx]) lat_var=reform(lat[*,*,sidx]) lon_var=reform(lon[*,*,sidx]) data_rho=reform(rho[*,*,sidx]) llat=min(lat) & ulat=max(lat) llon=min(lon) & rlon=max(lon) ;data_image=bytscl(data_var,top=top_color,min=dmin_kdp,max=dmax_kdp) data_image=byte(data_var) & data_image[*]=0B for i=0,6 do begin r=where(data_var ge kdp_int[i] and data_var lt kdp_int[i+1],c) if c gt 0 then data_image[r]=i print,kdp_int[i],kdp_int[i+1],c endfor r=where(data_var eq 0,c) if c gt 0 then data_image[r]=8;253 ;white r=where(data_var eq -9999,c) if c gt 0 then data_image[r]=9;255 ;gray if rho_filter eq 1 then begin r=where(data_rho lt lrho and data_var ne 0 and data_var ne -9999,c) if c gt 0 then data_image[r]=7;254 ;pink endif p1=image(data_image,lon_var,lat_var,limit=[llat,llon,ulat,rlon],$ /current,dimensions=[pxdim,pydim],grid_units='degrees',$ rgb_table=mytable3,linestyle='dot',$ label_format='MapGrid_Labels',/box_axes,label_position=0,$ label_show=0,irregular=0,$ map_projection='Mercator',position=pos[pnum,*],font_size=fs1) p2=mapcontinents(color='black',/hires) if sub_lon_plane ne !NULL then p7=plot(sub_lon_plane,sub_lat_plane,color='blue',thick=2,/overplot) ;t1=text(pos[pnum,0]-5*dx,pos[pnum,3]+2*dy,string(height_cband[sidx],format='(F4.2)')+'km',$ ; font_size=fs1) ; 3.0 km pnum=10 r=where(height_cband ge 3.0,c) sidx=r[0] data_var=reform(kdp[*,*,sidx]) lat_var=reform(lat[*,*,sidx]) lon_var=reform(lon[*,*,sidx]) data_rho=reform(rho[*,*,sidx]) llat=min(lat) & ulat=max(lat) llon=min(lon) & rlon=max(lon) ;data_image=bytscl(data_var,top=top_color,min=dmin_kdp,max=dmax_kdp) data_image=byte(data_var) & data_image[*]=0B for i=0,6 do begin r=where(data_var ge kdp_int[i] and data_var lt kdp_int[i+1],c) if c gt 0 then data_image[r]=i print,kdp_int[i],kdp_int[i+1],c endfor r=where(data_var eq 0,c) if c gt 0 then data_image[r]=8;253 ;white r=where(data_var eq -9999,c) if c gt 0 then data_image[r]=9;255 ;gray if rho_filter eq 1 then begin r=where(data_rho lt lrho and data_var ne 0 and data_var ne -9999,c) if c gt 0 then data_image[r]=7;254 ;pink endif p1=image(data_image,lon_var,lat_var,limit=[llat,llon,ulat,rlon],$ /current,dimensions=[pxdim,pydim],grid_units='degrees',$ rgb_table=mytable3,linestyle='dot',$ label_format='MapGrid_Labels',/box_axes,label_position=0,$ label_show=0,irregular=0,$ map_projection='Mercator',position=pos[pnum,*],font_size=fs1) p2=mapcontinents(color='black',/hires) if sub_lon_plane ne !NULL then p7=plot(sub_lon_plane,sub_lat_plane,color='blue',thick=2,/overplot) ;t1=text(pos[pnum,0]-5*dx,pos[pnum,3]+2*dy,string(height_cband[sidx],format='(F4.2)')+'km',$ ; font_size=fs1) ; 5.0 km pnum=11 r=where(height_cband ge 5.0,c) sidx=r[0] data_var=reform(kdp[*,*,sidx]) lat_var=reform(lat[*,*,sidx]) lon_var=reform(lon[*,*,sidx]) data_rho=reform(rho[*,*,sidx]) llat=min(lat) & ulat=max(lat) llon=min(lon) & rlon=max(lon) ;data_image=bytscl(data_var,top=top_color,min=dmin_kdp,max=dmax_kdp) data_image=byte(data_var) & data_image[*]=0B for i=0,6 do begin r=where(data_var ge kdp_int[i] and data_var lt kdp_int[i+1],c) if c gt 0 then data_image[r]=i print,kdp_int[i],kdp_int[i+1],c endfor r=where(data_var eq 0,c) if c gt 0 then data_image[r]=8;253 ;white r=where(data_var eq -9999,c) if c gt 0 then data_image[r]=9;255 ;gray if rho_filter eq 1 then begin r=where(data_rho lt lrho and data_var ne 0 and data_var ne -9999,c) if c gt 0 then data_image[r]=7;254 ;pink endif p1=image(data_image,lon_var,lat_var,limit=[llat,llon,ulat,rlon],$ /current,dimensions=[pxdim,pydim],grid_units='degrees',$ rgb_table=mytable3,linestyle='dot',$ label_format='MapGrid_Labels',/box_axes,label_position=0,$ label_show=0,irregular=0,$ map_projection='Mercator',position=pos[pnum,*],font_size=fs1) p2=mapcontinents(color='black',/hires) if sub_lon_plane ne !NULL then p7=plot(sub_lon_plane,sub_lat_plane,color='blue',thick=2,/overplot) ;t1=text(pos[pnum,0]-5*dx,pos[pnum,3]+2*dy,string(height_cband[sidx],format='(F4.2)')+'km',$ ; font_size=fs1) ; Colorbar ;cb0=colorbar(rgb_table=mycbtable,range=[dmin_kdp,dmax_kdp],$ ; orientation=0,position=[pos[9,0],pos[9,3]+0.5*dy,pos[10,2],pos[10,3]+1.5*dy],font_size=fs1,$ ; tickdir=1,border_on=1);,title='dBZ') ;t1=text(pos[9,0]-5*dx,pos[9,3]-1*dy,'KDP',font_size=fs1) ; Colorbar - vertical if rho_filter eq 1 then begin kdp_int_str=['<-0.1','-0.1','-0.05','0','0.05','0.1','1','>1','rho<0.95'] endif else begin kdp_int_str=['<-0.1','-0.1','-0.05','0','0.05','0.1','1','>1',''] endelse cb0=colorbar(rgb_table=mycbtable3,$;range=[dmin_kdp,dmax_kdp],$ ;mytable2 orientation=1,tickname=kdp_int_str,textpos=0,major=9,$ position=[pos[11,2]+4.5*dx,pos[11,1],pos[11,2]+5.5*dx,pos[11,3]-3*dy],font_size=10,$ tickdir=1,border_on=1);,title='dBZ') t1=text(pos[11,2]+3*dx,pos[11,3]-2*dy,'KDP',font_size=fs1) ;p0.save,'test5.png' endif ;*********************************** ; Row 4 is rhohv ;*********************************** if rho ne !NULL then begin r=where(rho ne -9999,c) dmin_rho=0.0;min(rho[r]) dmax_rho=1.0;max(rho[r]) ; 0.5 km pnum=12 r=where(height_cband ge 0.5,c) sidx=r[0] data_var=reform(rho[*,*,sidx]) lat_var=reform(lat[*,*,sidx]) lon_var=reform(lon[*,*,sidx]) llat=min(lat) & ulat=max(lat) llon=min(lon) & rlon=max(lon) data_image=bytscl(data_var,top=top_color,min=dmin_rho,max=dmax_rho) r=where(finite(data_var,/NAN),c) if c gt 0 then data_image[r]=253 r=where(data_var eq -9999,c) if c gt 0 then data_image[r]=255 p1=image(data_image,lon_var,lat_var,limit=[llat,llon,ulat,rlon],$ /current,dimensions=[pxdim,pydim],grid_units='degrees',$ rgb_table=mytable,linestyle='dot',$ label_format='MapGrid_Labels',/box_axes,label_position=0,$ label_show=0,irregular=0,$ map_projection='Mercator',position=pos[pnum,*],font_size=fs1) p2=mapcontinents(color='black',/hires) if sub_lon_plane ne !NULL then p7=plot(sub_lon_plane,sub_lat_plane,color='blue',thick=2,/overplot) ;t1=text(pos[pnum,0]-5*dx,pos[pnum,3]+2*dy,string(height_cband[sidx],format='(F4.2)')+'km',$ ; font_size=fs1) ; 1.0 km pnum=13 r=where(height_cband ge 1.0,c) sidx=r[0] data_var=reform(rho[*,*,sidx]) lat_var=reform(lat[*,*,sidx]) lon_var=reform(lon[*,*,sidx]) llat=min(lat) & ulat=max(lat) llon=min(lon) & rlon=max(lon) data_image=bytscl(data_var,top=top_color,min=dmin_rho,max=dmax_rho) r=where(finite(data_var,/NAN),c) if c gt 0 then data_image[r]=253 r=where(data_var eq -9999,c) if c gt 0 then data_image[r]=255 p1=image(data_image,lon_var,lat_var,limit=[llat,llon,ulat,rlon],$ /current,dimensions=[pxdim,pydim],grid_units='degrees',$ rgb_table=mytable,linestyle='dot',$ label_format='MapGrid_Labels',/box_axes,label_position=0,$ label_show=0,irregular=0,$ map_projection='Mercator',position=pos[pnum,*],font_size=fs1) p2=mapcontinents(color='black',/hires) if sub_lon_plane ne !NULL then p7=plot(sub_lon_plane,sub_lat_plane,color='blue',thick=2,/overplot) ;t1=text(pos[pnum,0]-5*dx,pos[pnum,3]+2*dy,string(height_cband[sidx],format='(F4.2)')+'km',$ ; font_size=fs1) ; 3.0 km pnum=14 r=where(height_cband ge 3.0,c) sidx=r[0] data_var=reform(rho[*,*,sidx]) lat_var=reform(lat[*,*,sidx]) lon_var=reform(lon[*,*,sidx]) llat=min(lat) & ulat=max(lat) llon=min(lon) & rlon=max(lon) data_image=bytscl(data_var,top=top_color,min=dmin_rho,max=dmax_rho) r=where(finite(data_var,/NAN),c) if c gt 0 then data_image[r]=253 r=where(data_var eq -9999,c) if c gt 0 then data_image[r]=255 p1=image(data_image,lon_var,lat_var,limit=[llat,llon,ulat,rlon],$ /current,dimensions=[pxdim,pydim],grid_units='degrees',$ rgb_table=mytable,linestyle='dot',$ label_format='MapGrid_Labels',/box_axes,label_position=0,$ label_show=0,irregular=0,$ map_projection='Mercator',position=pos[pnum,*],font_size=fs1) p2=mapcontinents(color='black',/hires) if sub_lon_plane ne !NULL then p7=plot(sub_lon_plane,sub_lat_plane,color='blue',thick=2,/overplot) ;t1=text(pos[pnum,0]-5*dx,pos[pnum,3]+2*dy,string(height_cband[sidx],format='(F4.2)')+'km',$ ; font_size=fs1) ; 5.0 km pnum=15 r=where(height_cband ge 5.0,c) sidx=r[0] data_var=reform(rho[*,*,sidx]) lat_var=reform(lat[*,*,sidx]) lon_var=reform(lon[*,*,sidx]) llat=min(lat) & ulat=max(lat) llon=min(lon) & rlon=max(lon) data_image=bytscl(data_var,top=top_color,min=dmin_rho,max=dmax_rho) r=where(finite(data_var,/NAN),c) if c gt 0 then data_image[r]=253 r=where(data_var eq -9999,c) if c gt 0 then data_image[r]=255 p1=image(data_image,lon_var,lat_var,limit=[llat,llon,ulat,rlon],$ /current,dimensions=[pxdim,pydim],grid_units='degrees',$ rgb_table=mytable,linestyle='dot',$ label_format='MapGrid_Labels',/box_axes,label_position=0,$ label_show=0,irregular=0,$ map_projection='Mercator',position=pos[pnum,*],font_size=fs1) p2=mapcontinents(color='black',/hires) if sub_lon_plane ne !NULL then p7=plot(sub_lon_plane,sub_lat_plane,color='blue',thick=2,/overplot) ;t1=text(pos[pnum,0]-5*dx,pos[pnum,3]+2*dy,string(height_cband[sidx],format='(F4.2)')+'km',$ ; font_size=fs1) ; Colorbar - vertical cb0=colorbar(rgb_table=mycbtable,range=[dmin_rho,dmax_rho],$ orientation=1,$ position=[pos[15,2]+4.5*dx,pos[15,1],pos[15,2]+5.5*dx,pos[15,3]-3*dy],font_size=fs1,$ tickdir=1,border_on=1);,title='dBZ') t1=text(pos[15,2]+1*dx,pos[15,3]-2*dy,'RHOHV',font_size=fs1) endif t1=text(0.4,0.97,file_basename(filename),font_size=fs1) t1=text(xl,0.97,'25 km range rings',font_size=fs1) ;p0.save,'test2.png' if rho_filter eq 1 then begin imagename=file_tag+'.12plot.rho_filter.png' endif else begin imagename=file_tag+'.12plot.png' endelse print,imagename ;p0.save,filedir+imagename;,height=pydim p0.save,imagename endfor ;if 1 eq 0 then begin ; This plots a bunch of image panels ; dmin=-17 ; dmax=31 ; ; for i=0,49,2 do begin ; var=reform(reflect[*,i,*]) ; var_image=bytscl(var,min=dmin,max=dmax,top=top_color) ; r=where(var eq -9999,c) ; if c gt 0 then var_image[r]=255 ; p5=image(var_image,/current,image_dimensions=[isx,isy],$ ; rgb_table=mytable,position=pos[pnum,*]) ; pnum=pnum+1 ; endfor ; p0.save,'test.png',height=pydim ; pxdim=1000 & pydim=1000 ; ; Position the plots ; xl=0.1 & xr=0.90 ; yb=0.1 & yt=0.90 ;; sx=0.05 ; sy=0.08 ; numplots_x=3 ; numplots_y=3 ; position_plots,xl,xr,yb,yt,sx,sy,numplots_x,numplots_y,pos ; pnum=0 ; alt1=alt[0,0,*]/1000.0 ; closest=min(abs(alt1-3.0),idx) ; data_var=reform(reflect[*,*,idx]) ; lat_image=reform(lat[*,*,idx]) ; lon_image=reform(lon[*,*,idx]) ; llat=min(lat) & ulat=max(lat) ; llon=min(lon) & rlon=max(lon) ; p0=image(data_var,lon_image,lat_image,limit=[llat,llon,ulat,rlon],$ ; /current,dimensions=[pxdim,pydim],grid_units='degrees',$ ; ; ;center_longitude=(rcolon_ovp-lcolon_ovp)/2.0, ; rgb_table=mytable,min_value=dmin,max_value=dmax,$ ; label_format='MapGrid_Labels',/box_axes,label_position=0,$ ; label_show=0,$ ; map_projection='Mercator',position=pos[pnum,*],font_size=fs1) ; t1=text(pos[pnum,0],pos[pnum,3]+2*dy,string(alt[idx]/1000.0,format='(F4.2)')+'km',$ ; font_size=12) ; p1=mapcontinents(color='black') ; lat1=lat[0,*,0] ; endif ; if 1 eq 0 then begin ; ;***************************** ; ; Plot map - Ku or Ka - dbz ; ;***************************** ; pnum=0 ; ; Calculate regular grid array ; ; llat=-65;min(lat) & ; ulat=-57;max(lat) ; llon=min(lon) & rlon=max(lon) ; ; delta=0.01 ; ynum=fix((max(lat)-min(lat))/delta) ; grid_lat=findgen(ynum)*delta+min(lat) ; xnum=fix((max(lon)-min(lon))/delta) ; grid_lon=findgen(xnum)*delta+min(lon) ; ; ; pull out surface ; var_lat=reform(lat[*,*,0]) ; var_lon=reform(lon[*,*,0]) ; var=reform(reflect[*,*,0]) ; result=where(var le -9999,count) ;-9999 is missing value ; var[result]=-7777 ; ; ;grid_input,colon,lat,var,xyz,newvar,/degree,/sphere,epsilon = 0.5 ; ;newcolon = !radeg * atan(xyz[1,*],xyz[0,*]) ; ;newlat = !radeg * asin(xyz[2,*]) ; ;qhull,newcolon,newlat,triangles,/delaunay,sphere=s ; ; Triangulate the data ; qhull,var_lon,var_lat,triangles,/delaunay,sphere=s ; ;triangulate,var_lon,var_lat,triangles,/degrees,sphere=s ; ; grid_var=griddata(var_lon,var_lat,var,xout=grid_lon,yout=grid_lat,$ ; ;/KRIGING,min_points=10,sectors=8,empty_sectors=4,$ ; ;/KRIGING,min_points=10,sectors=8,empty_sectors=4,$ ; ;method='NearestNeighbor',$ ; /inverse_distance,min_points=10,sectors=8,empty_sectors=3,$;max_per_sector=10,$ ; ;/inverse_distance,min_points=1,sectors=1,empty_sectors=1,max_per_sector=1,$ ; missing=-9999,triangles=triangles,/grid,/degrees,/sphere) ; ; ; Bytscal the data ; ;result=where(grid_var gt 0) ; dmin_radar=10;min(grid_var[result]) ; dmax_radar=35;max(grid_var[result]) ; data_image=bytscl(grid_var,top=top_color,min=dmin_radar,max=dmax_radar) ; ;result=where(grid_var gt -8000 and grid_var le -7000,count) ; result=where(grid_var gt -8000 and grid_var lt -10,count) ; if count gt 0 then data_image[result]=255 ;gray ; result=where(grid_var gt -10000 and grid_var le -9000,count) ; if count gt 0 then data_image[result]=253 ;white ; ; p0=image(data_image,grid_lon,grid_lat,limit=[llat,llon,ulat,rlon],$ ; /current,dimensions=[pxdim,pydim],grid_units='degrees',$ ; ;center_longitude=(rcolon_ovp-lcolon_ovp)/2.0, ; rgb_table=mytable,$ ; map_projection='Mercator',position=pos[pnum,*],font_size=fs1) ; ; p1=mapcontinents(color='black') ; p0.mapgrid.box_axes=1 ; p0.mapgrid.box_color='black' ; p0.mapgrid.color='black' ; p0.mapgrid.label_color='black' ; p0.mapgrid.linestyle='dot' ; p0.mapgrid.label_position=0 ; ; p0.save,'test_cband_pyart.png',height=pydim ; endif ; endif ; endfor ;endif ;scan within 2 hours end