;***************************** ; 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 read_cband_latlon_allscans_volume,julian_day_ship,grid_var,grid_colon,grid_lat,$ nbins,rstart,rscale,julian_day,elangle,ulat,llat,llon,rlon,clon,clat ;*** INPUT *** ulat=-57.7 & llat=-60.0 llon=137.5 & rlon=142.0 ;julian_day_ship=julday(1,25,2018,19,22,29) julian_day_ship=julday(2,8,2018,16,56,36) ;ulat=-42.0 & llat=-66 ;llon=130.0 & rlon=152.0 path_prefix='/Volumes/' ;path_prefix='/uufs/chpc.utah.edu/common/home/' ;*** Color table *** ; 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])] cbtable=mytable[0:top_color,*] ;******************* ; Plots the first scan only, nelev=1 ; To plot other scans change nelev at line 69 ; directory filedir=path_prefix+'mace-group4/Capricorn2/CBand/' ; Date-time strings of the ship time 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)') ; File search string filestr='9776HUB-PPIVol-'+syy+smm+sdd+'-'+shh+'*hdf' ;filestr='9776HUB-PPIVol-'+syy+smm+sdd+'*hdf' ;filestr='9776HUB-PPIVol-*hdf' ; file name filenames=file_search(filedir+'/'+syy+smm+sdd+'/'+filestr,count=numfiles) if numfiles gt 0 then begin ; Only use the big files finfo=file_info(filenames) result_size=where(finfo.size gt 3e7,count_size) filenames=filenames[result_size] numfiles=count_size ; 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) caldat,julian_day_files[idx],xmm,xdd,xyy,xhh,xmi,xss print,xyy,xmm,xdd,xhh,xmi,xss numfiles=1 kk1=idx ; Read in the data from the files ;for kk=0,numfiles-1 do begin for kk=kk1,kk1 do begin filename=filenames[kk] print,filename file_tag=strmid(file_basename(filename),16,19) ; dump the hdf structure ;presult=h5_parse(filename) ; Open an existing hdf5 file file_id=h5f_open(filename) ; Root / where (lon, lat, height) group_id=h5g_open(file_id,'where') attr_id=h5a_open_name(group_id,'lon') ;154.38200 clon=h5a_read(attr_id) attr_id=h5a_open_name(group_id,'lat') ;-26.617001 clat=h5a_read(attr_id) attr_id=h5a_open_name(group_id,'height') ;22.00 height=h5a_read(attr_id) h5a_close,attr_id h5g_close,group_id ;print, 'lon/lat/ht = ', lon, lat, height ; Root / what (object, version, date, time, source) group_id=h5g_open(file_id,'what') attr_id=h5a_open_name(group_id,'source') ;WMO:00000,RAD:,PLC:CSIRO,NOD:9776HUB,ORG:0,CTY:000,CMT: csgn=h5a_read(attr_id) h5a_close,attr_id h5g_close,group_id ;print, 'station ID = ', csgn ; Root / how (wavelength, beamwidth) group_id=h5g_open(file_id,'how') attr_id=h5a_open_name(group_id,'wavelength') ;5.3399911 wavelength=h5a_read(attr_id) attr_id=h5a_open_name(group_id,'beamwH') ;1.3 beamwH=h5a_read(attr_id) h5a_close,attr_id h5g_close,group_id ;print, 'wavelength,beamwH = ', wavelength, beamwH ; Set up the arrays to hold the data ;nelev=17 ;how many scans to read, total number scans is 17 ;nelev=12 ;new file total number of scans nelev=6 ;nelev=1 ; first level scan only startdate=strarr(nelev) & starttime=strarr(nelev) & endtime=strarr(nelev) dataset_n=strarr(nelev) & elangle=dblarr(nelev) & nbins=intarr(nelev) rstart=dblarr(nelev) & rscale=dblarr(nelev) & nrays=intarr(nelev) a1gate=intarr(nelev) & ray_spacing=fltarr(nelev) & sjulian_day=dblarr(nelev) ejulian_day=dblarr(nelev) & delta_t=dblarr(nelev) ; Initialize gate_start to -0.5 gate_start=make_array(/double,nelev,value=-0.5) totdata=17 ;total number of data arrays in a dataset data_n=strarr(totdata) quantity=strarr(totdata) quality_n=strarr(totdata) gain=dblarr(totdata) offset=dblarr(totdata) nodata=dblarr(totdata) undetect=dblarr(totdata) ref=dblarr(nelev,3,960,360) ;hold the reflectivity data ; Loop through the datasets at the root level for i=0,nelev-1 do begin print,'i',i,'*****' dnum=strcompress(i+1,/remove_all) ; Root / datasetX dataset_n[i]='dataset'+dnum print,'opening ',dataset_n[i] group_id=h5g_open(file_id,dataset_n[i]) ; ... finding start/end date/time g_id=h5g_open(group_id,'what') attr_id=h5a_open_name(g_id,'startdate') ;20170822 startdate[i]=h5a_read(attr_id) ;print, 'startdate:',startdate[i] attr_id=h5a_open_name(g_id,'starttime') ;205654 starttime[i]=h5a_read(attr_id) ;print, 'starttime:',starttime[i] attr_id=h5a_open_name(g_id,'endtime') ;205730 endtime[i]=h5a_read(attr_id) ;print, 'endtime:',endtime[i] h5a_close,attr_id h5g_close,g_id h5g_close,group_id ; Constructing the time array syear=strmid(startdate[i],0,4) smonth=strmid(startdate[i],4,2) sday=strmid(startdate[i],6,2) shour=strmid(starttime[i],0,2) smin=strmid(starttime[i],2,2) ssec=strmid(starttime[i],4,2) ehour=strmid(endtime[i],0,2) emin=strmid(endtime[i],2,2) esec=strmid(endtime[i],4,2) ;date_str=syear+smonth+sday+shour+smin+ssec sjulian_day[i]=julday(smonth,sday,syear,shour,smin,ssec) ejulian_day[i]=julday(smonth,sday,syear,ehour,smin,ssec) ;print,'date and start and end time:',startdate[i],' ',starttime[i],' ',endtime[i] ; Root / datasetX / where to read (elangle, nbins, rstart, rscale, nrays, a1gate) group_id=h5g_open(file_id,dataset_n[i]) g_id=h5g_open(group_id,'where') attr_id=h5a_open_name(g_id,'elangle') ;dataset1=0.5 elangle[i]=h5a_read(attr_id) h5a_close,attr_id attr_id=h5a_open_name(g_id,'nbins') ;dataset1=960 nbins[i]=h5a_read(attr_id) h5a_close,attr_id attr_id=h5a_open_name(g_id,'rstart') ;dataset1=0 rstart[i]=h5a_read(attr_id) h5a_close,attr_id attr_id=h5a_open_name(g_id,'rscale') ;dataset1=125 rscale[i]=h5a_read(attr_id) h5a_close,attr_id attr_id=h5a_open_name(g_id,'nrays') ;dataset1=360 nrays[i]=h5a_read(attr_id) h5a_close,attr_id attr_id=h5a_open_name(g_id,'a1gate') ;dataset1=256 a1gate[i]=h5a_read(attr_id) h5a_close,attr_id delta_t[i]=(ejulian_day[i]-sjulian_day[i])/nrays[i] ;print,i,elangle[i],nbins[I],rstart[I],rscale[I],nrays[I],a1gate[I] ;print, 'dataset ', i, ' nbins,nrays = ', nbins[I],nrays[i] h5g_close,g_id ; Root / datasetX / where to read (elangle, nbins, rstart, rscale, nrays, a1gate) group_id=h5g_open(file_id,dataset_n[i]) g_id=h5g_open(group_id,'how') attr_id=h5a_open_name(g_id,'startazA') ;dataset1=0.5 startaza=h5a_read(attr_id) h5a_close,attr_id ;xstr='startaza'+dnum+'=startaza' ;this array is not always the same size ;result=execute(xstr) attr_id=h5a_open_name(g_id,'stopazA') ;dataset1=960 stopaza=h5a_read(attr_id) h5a_close,attr_id ;xstr='stopaza'+dnum+'=stopaza' ;this array is not always the same size ;result=execute(xstr) h5g_close,g_id ; loop through data arrays in dataset for j = 0,totdata-1 do begin ;print,'j',j,'&&&&&&' dnumj=strcompress(j+1,/remove_all) data_n[j]='data'+dnumj ;name of data array ;print,'OPENING data=',j,data_n[j],'dnumj',dnumj ; Root / datasetX / dataX to read the data data_id=h5g_open(group_id,data_n[j]) ; Grab the quanity value of the datsetX.dataX array ; Root / datasetX / dataX / what (gain, offset, nodata, undetect, quantity) g_id=h5g_open(data_id,'what') attr_id=h5a_open_name(g_id,'gain') ;0.5 gain[j]=h5a_read(attr_id) attr_id=h5a_open_name(g_id,'offset') ;-32 offset[j]=h5a_read(attr_id) attr_id=h5a_open_name(g_id,'nodata') ;255 nodata[j]=h5a_read(attr_id) attr_id=h5a_open_name(g_id,'undetect') ;0 undetect[j]=h5a_read(attr_id) attr_id=h5a_open_name(g_id,'quantity') quantity[j]=h5a_read(attr_id) ;print, 'quantity ',j,quantity[j] h5a_close,attr_id h5g_close,g_id ; Continue if this is a quantity we want to plot if quantity[j] eq 'DBZH' then begin print, 'quantity ',i,j,' ',quantity[j] ; read data dset_id=h5d_open(data_id,'data') dataX1=float(h5d_read(dset_id)) h5d_close,dset_id ; data h5g_close,data_id ; dataX ; apply gain and offset result=where(dataX1 eq nodata[j],count) ;no values in dataX1=missing if count gt 0 then dataX1[result] = -8888. result=where(dataX1 eq undetect[j],count) ;0 values = undetect=no cloud if count gt 0 then dataX1[result] = -9999. result=where(dataX1 gt -8888,count) if count gt 0 then dataX1[result]=dataX1[result]*gain[j]+offset[j] ; Separate the arrays by the quality string if strmid(quantity[j],0,5) eq 'VRADH' then begin xstr='dop'+dnum+'=dataX1' result=execute(xstr) endif else if strmid(quantity[j],0,4) eq 'DBZH' then begin xstr='ref'+dnum+'=dataX1' result=execute(xstr) endif else if strmid(quantity[j],0,5) eq 'CLASS' then begin xstr='radar_flag'+dnum+'=dataX1' result=execute(xstr) endif else if strmid(quantity[j],0,5) eq 'WRADH' then begin xstr='sw'+dnum+'=dataX1' result=execute(xstr) endif ;else begin endif ;end of this is a quantity I want endfor ;end j loop through data arrays in dataset h5g_close,group_id ;close dataX ;h5f_close,file_id ; ray spacing is 1 degree ray_spacing[i] = 360. / nrays[i] ;print,'ray spacing:',ray_spacing[I] ;***** ; Convert azimuth, elevation, range to x,y,z ;***** n=0L ;number of points in scan ; Define arrays at this level numcount=long(nrays[i])*long(nbins[i]) azimuth=fltarr(numcount) degree=fltarr(numcount) elevation=fltarr(numcount) range=fltarr(numcount) x=fltarr(numcount) y=fltarr(numcount) z=fltarr(numcount) dbz=fltarr(numcount) distance=fltarr(numcount) julian_day=dblarr(numcount) lat=fltarr(numcount) lon=fltarr(numcount) for j=0,nrays[i]-1 do begin ;rays ;xstr='startaza=startaza'+dnum ;result=execute(xstr) ;xstr='stopaza=stopaza'+dnum ;result=execute(xstr) if startaza[j] ge 350.0 and stopaza[j] le 2.0 then begin daza=(stopaza[j]+360.0-startaza[j])/nbins[i] ;print,'in here' endif else begin daza=(stopaza[j]-startaza[j])/nbins[i] endelse ;print, startaza[j],stopaza[j],nbins[i],daza for k=0,nbins[i]-1 do begin ;bins ; azimuth, 0degNorth ;azimuth[n] = 90.0-(float(j)*ray_spacing[i] + gate_start[i]) ;my old way ;azimuth[n] = (float(j)*ray_spacing[i] + gate_start[i]) ;degrees ;azimuth[n] = float(j)*ray_spacing[i] azimuth[n]=startaza[j]+float(k)*daza ; Convert azimuth to math coordinate degree degree[n]=450.0-azimuth[n] if degree[n] gt 360.0 then degree[n]=degree[n]-360.0 elevation[n]=elangle[i] range[n]=rstart[i]+float(k)*rscale[i] ;meters ; height height=(float(k)*rscale[i])*sin(elangle[i]*(!pi/180.))+rstart[i] ; msl z[n]=height ; distance d_hor=(float(k)*rscale[i])*cos(elangle[i]*(!pi/180.)) distance[n]=d_hor ;radians x[n]=d_hor*sin(azimuth[n]*!pi/180.) y[n]=d_hor*cos(azimuth[n]*!pi/180.) ; dbz xstr='dbz[n]=ref'+dnum+'[k,j] result=execute(xstr) ; Time julian_day[n]=sjulian_day[i]+double(j)*delta_t[i] ; Lat and Lon arc_dist=range[n]/(6371D*1000D) ;radius of earth in meters ;result=ll_arc_distance([clon,clat],arc_dist,degree[n],/degrees) result=ll_arc_distance([clon,clat],arc_dist,azimuth[n],/degrees) lon[n]=result[0] lat[n]=result[1] n=n+1l endfor ;end of loop through ranges k endfor ;end of loop through azimuths j ;******************* ; Now plot the scans ;******************* var_image=bytscl(dbz,min=-20,max=60) var_trans=fix(bytscl(dbz,min=-20,max=60,top=100)) ;p0=scatterplot3D(x,y,z,magnitude=var_image,rgb_table=39,zrange=[0,4000]) p0=plot3D(x/1000.0,y/1000.0,z,'o',/sym_filled,rgb_table=33, vert_colors=var_image,$ depth_cue=[0,2],axis_style=2,/perspective,zrange=[0,10000],/buffer,$ xy_shadow=1) ax = p0.AXES ax[2].HIDE = 1 ax[6].HIDE = 1 ax[7].HIDE = 1 p0.save,'plot_3d_elev.'+string(i,format='(I02)')+'.3.png',height=600 ;var_grid=grid3(x,y,z,dbz) ;vol=volume(var_grid) ;r=range ;theta=azm ;f=ref ; Loop throught the elevations to plot ; Caculate colongitude because the pacific crosses the 180 line colon=lon result=where(lon lt 0,count) if count gt 0 then colon[result]=colon[result]+360.0 lcolon=llon if lcolon lt 0 then lcolon=lcolon+360.0 rcolon=rlon if rcolon lt 0 then rcolon=rcolon+360.0 ; decrease size of arrays so that it doesn't take too long to plot num=2 dbz_new=dbz[0:-1:num] lon_new=lon[0:-1:num] lat_new=lat[0:-1:num] colon_new=colon[0:-1:num] ; Calculate regular grid array delta=0.005 ynum=(ulat-llat)/delta grid_lat=findgen(ynum)*delta+llat xnum=(rcolon-lcolon)/delta grid_colon=findgen(xnum)*delta+lcolon grid_lon=grid_colon result=where(grid_lon gt 180,count) if count gt 0 then grid_lon[result]=grid_lon[result]-360.0 print,'start gridding' ; Triangulate the data qhull,colon_new,lat_new,triangles,/delaunay,sphere=s var=dbz_new result=where(var eq -9999,count) var[result]=-7777 grid_var=griddata(colon_new,lat_new,var,xout=grid_colon,yout=grid_lat,$ /grid,/degrees,/sphere,$ ;/KRIGING,min_points=10,sectors=8,empty_sectors=4,$ /KRIGING,min_points=10,sectors=8,empty_sectors=5,$ ;method='NearestNeighbor',$ ;/inverse_distance,min_points=10,sectors=8,empty_sectors=3,$;max_per_sector=10,$ ;/inverse_distance,empty_sectors=1,max_per_sector=1,sectors=1,min_points=1,$ missing=-9999,$ triangles=triangles);,$ if grid_var_all eq !NULL then begin s=size(grid_var,/dimensions) grid_var_all=make_array(/float,nelev,s[0],s[1]) endif grid_var_all[i,*,*]=grid_var print,'finish griddata' endfor ; end of loop through scans to plot ;h5f_close,file_id endfor ;end of loop through files endif ;end of files were found ; Bytscal the data dmin=10 dmax=35 ;data_image=bytscl(grid_var_all,top=top_color,min=dmin,max=dmax) ; result=where(grid_var gt -8000 and grid_var lt dmin,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 ; Quick plot of dbz_grid pxdim=900 & pydim=700 ; Position the plots xl=0.05 & xr=0.95 yb=0.05 & yt=0.90 sx=0.05 sy=0.05 numplots_x=3 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.14 cbpos[*,2]=cbpos[*,0]+0.02 ; set plot area pnum=0 p0=plot([0,1],[0,1],/buffer,/nodata,axis_style=4,dimensions=[pxdim,pydim],position=pos[pnum,*]) 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]) px0=d[0,0,0] py0=d[1,0,0] for i=0,nelev-1 do begin grid1=reform(grid_var_all[i,*,*]) data_image1=bytscl(grid1,top=top_color,min=dmin,max=dmax) result=where(grid1 gt -9000 and grid1 lt dmin,count) if count gt 0 then data_image1[result]=255 ;gray result=where(grid1 gt -10000 and grid1 le -9000,count) if count gt 0 then data_image1[result]=253 ;white ;data_image1=reform(data_image[i,*,*]) p0=image(data_image1,grid_colon,grid_lat,limit=[llat,lcolon,ulat,rcolon],$ /current,dimensions=[pxdim,pydim],grid_units='degrees',$ center_longitude=(rcolon-lcolon)/2.0,$ rgb_table=mytable,$;,min_value=dmin,max_value=dmax,$ 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.mapgrid.label_show=0 t1=text(pos[pnum,0],pos[pnum,3],string(elangle[i],format='(F5.2)'),font_size=fs1) ;c0=colorbar($ ; rgb_table=cbtable,$ ; range=[dmin,dmax],$ ; orientation=1,position= reform(cbpos[pnum,*]),font_size=10,$ ; tickdir=1,border_on=1,title='DbZ') pnum=pnum+1 print, pnum endfor c0=colorbar($ rgb_table=cbtable,$ range=[dmin,dmax],$ orientation=0,position=[xl,0.95,0.5,0.98],font_size=10,$ tickdir=1,border_on=1,title='DbZ') imagename='dbz_panel_'+file_tag+'.png' print,imagename p0.save,imagename,height=pydim stop ;endif ; ;; pass out these arrays for plotting ;;dbz_grid ;;x_out/1000.0,y_out/1000.0 ;;nbins ; ;do_plot4='no' ;if do_plot4 eq 'yes' then begin ; ;*** Color table *** ; ; 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])] ; cbtable=mytable[0:top_color,*] ; ; ; Quick plot of dbz_grid ; pxdim=700 & pydim=600 ; ; Position the plots ; xl=0.15 & xr=0.80 ; yb=0.15 & yt=0.90 ; sx=0.14 ; sy=0.03 ; numplots_x=1 ; numplots_y=1 ; position_plots,xl,xr,yb,yt,sx,sy,numplots_x,numplots_y,pos ; ; Colorbar position ; cbpos=pos ; cbpos[*,0]=pos[*,2]+0.14 ; cbpos[*,2]=cbpos[*,0]+0.02 ; ; set plot area ; pnum=0 ; p0=plot([0,1],[0,1],/buffer,/nodata,axis_style=4,dimensions=[pxdim,pydim],position=pos[pnum,*]) ; 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]) ; px0=d[0,0,0] ; py0=d[1,0,0] ; ; Plot reflectivity ; pnum=0 ; ; ;data_var=dbz_mask ; ;dmin=min(data_var) ; ;dmax=max(data_var) ; ;data_image=bytscl(data_var,min=dmin,max=dmax) ; ; data_var=dbz_grid ; ;result=where(data_var ne -9999,count) ; dmin=-20.0;min(data_var[result]) ; dmax=50.0;max(data_var[result]) ; data_image=bytscl(data_var,top=top_color,min=dmin,max=dmax) ; ;result=where(data_var eq -999.,count) ; result=where(data_var eq -9999,count) ; if count gt 0 then data_image[result]=253 ; ;result=where(dbz_mask eq -8888,count) ; ;if count gt 0 then data_image[result]=255 ; p1=image(data_image,/current,image_dimensions=[isx,isy],position=pos[pnum,*],rgb_table=mytable) ; c1=contour(data_var,x_out/1000.0,y_out/1000.0,ytitle='Y (km)',xtitle='X (km)',$ ; xstyle=1,ystyle=1,font_size=14,/nodata,/current,position=pos[pnum,*],$ ; irregular=0,xtickdir=1,ytickdir=1,font_name='Helvetica',$ ; axis_style=2) ; ; plot range rings ; for r=0,nbins[i]-1,150 do begin ; ;print,'range ring',r ; ; radius in km ; range_ring=(rstart[i]+float(r)*rscale[i])/1000.0 ; result=circle(0,0,range_ring) ; p2=plot(result,/current,/overplot) ; ;stop ; endfor ; c2=colorbar(orientation=1,rgb_table=cbtable,range=[dmin,dmax],title='Z (dBZ)',$ ; position=reform(cbpos[pnum,*]),font_size=14,font_name='Helvetica',tickdir=1,border_on=1) ; caldat,julian_day1[0],mm,dd,yy,hh,mi,ss ; date_time=string(yy,format='(I04)')+string(mm,format='(I02)')+string(dd,format='(I02)')+$ ; ' '+string(hh,format='(I02)')+':'+string(mi,format='(I02)')+':'+string(ss,format='(I02)') ; scan_str=string(elangle[i],format='(F05.2)') ; t1=text(xl,0.95,'el='+scan_str,font_size=14) ; t1=text(0.5,0.95,date_time,font_size=14) ; imagename='dbz_grid_'+file_tag+'.'+dnum+'.png' ; p0.save,imagename,height=pydim ; ;p0.save,'dbz_grid.el'+scan_str+'.'+startdate[i]+'.'+starttime[i]+'.mp08.png',height=pydim ; print,imagename ;endif ; do plot stop end