;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_pyart_shallowcumulus,julian_day_ship,ulat,llat,llon,rlon,$ reflect,lat,lon,alt,julian_day,colon ; Date to test ;julian_day_ship=julday(1,25,2018,19,22,29) ;overcast ;julian_day_ship=julday(1,28,2018,16,44,36) ;julian_day_ship=julday(2,10,2018,14,58,13) ;caldat,julian_day_ship,smm,sdd,syy,shh,smi,sss ;syy=string(syy,format='(I4)') ;smm=string(smm,format='(I02)') ;sdd=string(sdd,format='(I02)') ;shh=string(shh,format='(I02)') ;smi=string(smi,format='(I02)') ;sss=string(sss,format='(I02)') ;date_str=syy+smm+sdd+shh+smi+sss ;ymd_str=syy+smm+sdd ; directory filedir='/uufs/chpc.utah.edu/common/home/mace-group4/Capricorn2/CBand/' ; File search string ;filestr='9776HUB-PPIVol-'+syy+smm+sdd+'-'+shh+'*nc' ;filestr='9776HUB-PPIVol-'+syy+smm+sdd+'*81.401.nc' filestr='9776HUB-PPIVol-*81.401.nc' ; file name filenames=file_search(filedir+filestr,count=numfiles) ; Get the times of the file start and find the closest to ship time ;time_str=strmid(file_basename(filenames),24,6) ;fhh=strmid(time_str,0,2) ;fmi=strmid(time_str,2,2) ;fss=strmid(time_str,4,2) ;julian_day_files=julday(smm,sdd,syy,fhh,fmi,fss) ;closest=min(abs(julian_day_files-julian_day_ship),idx) ; Closest must be in 2 hours ;if closest le 2D/24D then begin ; caldat,julian_day_files[idx],xmm,xdd,xyy,xhh,xmi,xss ; print,xyy,xmm,xdd,xhh,xmi,xss ; numfiles=1 ; kk1=idx for idx=0,numfiles-1 do begin ; for kk=kk1,kk1 do begin filename=filenames[idx] print,filename 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,'reflectivity') & ncdf_varget,fid,vid,reflect 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 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)) ; Calculate colongitude colon=lon result=where(lon lt 0,count) if count gt 0 then colon[result]=colon[result]+360.0 ; Zoom in on 75 km x1=point_x[*,0,0] xr=where(x1 ge -75 and x1 le 75,xc) y1=point_y[0,*,0] yr=where(y1 ge -75 and y1 le 75,yc) ; reflect=reflect[xr,*,*] lat=lat[xr,*,*] lon=lon[xr,*,*] reflect=reflect[*,yr,*] lat=lat[*,yr,*] lon=lon[*,yr,*] ;do_plot='yes' ;if do_plot eq 'yes' then begin file_tag=strmid(file_basename(filename),16,19) ; Set up plot size 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 ; Colorbar position cbpos=pos ;cbpos[*,0]=pos[*,2]+0.08 ;vertical ;cbpos[*,2]=cbpos[*,0]+0.01 cbpos[*,1]=pos[*,3]+0.005 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) ;254=hot pink ;gray=255 ;mytable=[mytable,transpose([238,18,137]),transpose([150,150,150])] ;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,*] dmin=-17 & dmax=31 ; Strings dx=0.01 & dy=0.01 fs1=9 ;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]) ;REFLECT FLOAT = Array[201, 201, 41] for i=0,17,2 do begin data_var=reform(reflect[*,*,i]) lat_image=reform(lat[*,*,i]) lon_image=reform(lon[*,*,i]) 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=1,irregular=0,$ map_projection='Mercator',position=pos[pnum,*],font_size=fs1) t1=text(pos[pnum,0]-5*dx,pos[pnum,3]+2*dy,string(alt[0,0,i]/1000.0,format='(F4.2)')+'km',$ font_size=13) p1=mapcontinents(color='black') ;p0.mapgrid.box_axes=1 ;p0.mapgrid.box_color='black' ;0.mapgrid.color='black' ;p0.mapgrid.label_color='black' p0.mapgrid.linestyle='dot' ;p0.mapgrid.label_position=0 pnum=pnum+1 endfor t1=text(pos[0,0],0.95,file_basename(filename),font_size=12) c0=colorbar(rgb_table=mycbtable,range=[dmin,dmax],$ orientation=0,position=[0.75,0.97,0.9,0.98],font_size=10,$ tickdir=1,border_on=1,title='dBZ') p0.save,file_basename(filename)+'.zoom.png',height=pydim 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 endfor stop ; 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