;************************************* ; This program reads in the raw cfradial format cband data. ; It plots the cloudy rays and then 1 plan view. ;************************************ ;***************************** ; 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_cfradial,julian_day_ship,grid_var,grid_colon,grid_lat,$ nbins,rstart,rscale,julian_day,elangle,ulat,llat,llon,rlon,clon,clat ;*** INPUT *** ; RV Investigator ;date_str='20180209-113309' ;date_str='20180125-234004' ;datadir='mace-group5/field_programs/investigator/in2018_v01/cband/' ;CSU Sea-Pol date_str='20190924_055005' datadir='mace-group6/cpol_dolan/MACE_SEAPOL_PPI_SEPT24/' 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,*] ; directory filedir=path_prefix+datadir ; Date-time strings ;caldat,julian_day_ship,smm,sdd,syy,shh,smi,sss fyy=strmid(date_str,0,4) fmm=strmid(date_str,4,2) fdd=strmid(date_str,6,2) fhh=strmid(date_str,9,2) fmi=strmid(date_str,11,2) fss=strmid(date_str,13,2) ; File search string if date_str eq '20190924_055005' then begin filestr='SEA'+fyy+fmm+fdd+'_'+fhh+fmi+'*ppi.nc' endif else begin filestr=fyy+fmm+fdd+'/9776HUB-PPIVol-'+fyy+fmm+fdd+'-'+fhh+fmi+'*.cfradial.nc' endelse ; file name filenames=file_search(filedir+filestr,count=numfiles) ; 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 ; Continue if you found a file if numfiles gt 0 then begin ; Read in the data from the files for kk=0,numfiles-1 do begin ; Filename filename=filenames[kk] print,filename ;file_tag=strmid(file_basename(filename),15,20) ; Open cfradial netcdf file fid=ncdf_open(filename) vid=ncdf_varid(fid,'longitude') & ncdf_varget,fid,vid,clon vid=ncdf_varid(fid,'latitude') & ncdf_varget,fid,vid,clat vid=ncdf_varid(fid,'altitude') & ncdf_varget,fid,vid,calt & calt=float(calt) vid=ncdf_varid(fid,'time_coverage_start') & ncdf_varget,fid,vid,time_str time_str=string(time_str) syy=strmid(time_str,0,4) smm=strmid(time_str,5,2) sdd=strmid(time_str,8,2) shh=strmid(time_str,11,2) smi=strmid(time_str,14,2) sss=strmid(time_str,17,2) print, 'lon/lat/alt = ', clon, clat, calt print, 'time_str = ', syy,smm,sdd,shh,smi,sss ; Range to center of each bin (m) n=960 vid=ncdf_varid(fid,'range') & ncdf_varget,fid,vid,range ; Time at the center of each ray (s) n=3978 vid=ncdf_varid(fid,'time') & ncdf_varget,fid,vid,time vid=ncdf_varid(fid,'azimuth') & ncdf_varget,fid,vid,azimuth vid=ncdf_varid(fid,'elevation') & ncdf_varget,fid,vid,elevation ;julian_day=julday(smm,sdd,syy,shh,smi,sss)+double(time)/(60D*60D*24D) ; (time,range) vid=ncdf_varid(fid,'temperature') if vid ne -1 then ncdf_varget,fid,vid,temp vid=ncdf_varid(fid,'air_echo_classification') if vid ne -1 then ncdf_varget,fid,vid,air vid=ncdf_varid(fid,'radar_echo_classification') if vid ne -1 then ncdf_varget,fid,vid,hclas vid=ncdf_varid(fid,'radar_estimated_rain_rate') if vid ne -1 then ncdf_varget,fid,vid,rain vid=ncdf_varid(fid,'velocity') if vid ne -1 then ncdf_varget,fid,vid,dvel vid=ncdf_varid(fid,'corrected_velocity') if vid ne -1 then ncdf_varget,fid,vid,dvel_cor vid=ncdf_varid(fid,'total_power') if vid ne -1 then ncdf_varget,fid,vid,power ; DZ or CZ vid=ncdf_varid(fid,'corrected_reflectivity') if vid ne -1 then ncdf_varget,fid,vid,dbz if vid eq -1 then begin vid=ncdf_varid(fid,'DZ') ;reflectivity dBZ ncdf_varget,fid,vid,dbz dbz=float(dbz) r=where(dbz eq -32768,c) if c gt 0 then dbz[r]=-9999 r=where(dbz ne -9999,c) if c gt 0 then dbz[r]=dbz[r]*0.01 endif ; RH vid=ncdf_varid(fid,'cross_correlation_ratio') if vid ne -1 then ncdf_varget,fid,vid,rhohv ;correlation H-to-V if vid eq -1 then begin vid=ncdf_varid(fid,'RH') ncdf_varget,fid,vid,rhohv rhohv=float(rhohv) r=where(rhohv eq -32768,c) if c gt 0 then rhohv[r]=-9999 r=where(rhohv ne -9999,c) if c gt 0 then rhohv[r]=rhohv[r]*0.001 endif ; ZD vid=ncdf_varid(fid,'differential_reflectivity') if vid ne -1 then ncdf_varget,fid,vid,zdr if vid eq -1 then begin vid=ncdf_varid(fid,'ZD') ;differential reflectivity dB ncdf_varget,fid,vid,zdr zdr=float(zdr) r=where(zdr eq -32768,c) if c gt 0 then zdr[r]=-9999 r=where(zdr ne -9999,c) if c gt 0 then zdr[r]=zdr[r]*0.01 endif vid=ncdf_varid(fid,'corrected_differential_reflectivity') if vid ne -1 then ncdf_varget,fid,vid,zdr_cor ; PH vid=ncdf_varid(fid,'differential_phase') if vid ne -1 then ncdf_varget,fid,vid,phidp ;differential phase (deg) if vid eq -1 then begin vid=ncdf_varid(fid,'PH') ncdf_varget,fid,vid,phidp phidp=float(phidp) r=where(phidp eq -32768,c) if c gt 0 then phidp[r]=-9999 r=where(phidp ne -9999,c) if c gt 0 then phidp[r]=phidp[r]*0.01 endif vid=ncdf_varid(fid,'corrected_differential_phase') if vid ne -1 then ncdf_varget,fid,vid,phidp_cor ; KD vid=ncdf_varid(fid,'corrected_specific_differential_phase') if vid ne -1 then ncdf_varget,fid,vid,kdp if vid eq -1 then begin vid=ncdf_varid(fid,'KD') ;differential propagation of phase (deg/km) ncdf_varget,fid,vid,kdp kdp=float(kdp) r=where(kdp eq -32768,c) if c gt 0 then kdp[r]=-9999 r=where(kdp ne -9999,c) if c gt 0 then kdp[r]=kdp[r]*0.01 endif ; HD is all 0 vid=ncdf_varid(fid,'HD') ncdf_varget,fid,vid,hd r=where(hd eq -32768,c) if c gt 0 then hd[r]=-9999 ; (sweep) vid=ncdf_varid(fid,'sweep_number') & ncdf_varget,fid,vid,sweep_number vid=ncdf_varid(fid,'sweep_start_ray_index') & ncdf_varget,fid,vid,sweep_start_ray_index vid=ncdf_varid(fid,'sweep_end_ray_index') & ncdf_varget,fid,vid,sweep_end_ray_index ncdf_close,fid nrays_sweep=make_array(n_elements(sweep_number),value=-9999,/long) nbins_sweep=make_array(n_elements(sweep_number),value=-9999,/long) ; Loop through the sweeps ;for i=0,n_elements(sweep_number)-1 do begin for i=0,1 do begin ; Pull out the sweep data azimuth_sweep=azimuth[sweep_start_ray_index[i]:sweep_end_ray_index[i]] elevation_sweep=elevation[sweep_start_ray_index[i]:sweep_end_ray_index[i]] time_sweep=time[sweep_start_ray_index[i]:sweep_end_ray_index[i]] if temp ne !NULL then begin temp_sweep=temp[*,sweep_start_ray_index[i]:sweep_end_ray_index[i]] air_sweep=air[*,sweep_start_ray_index[i]:sweep_end_ray_index[i]] hclas_sweep=hclas[*,sweep_start_ray_index[i]:sweep_end_ray_index[i]] rain_sweep=rain[*,sweep_start_ray_index[i]:sweep_end_ray_index[i]] dvel_sweep=dvel[*,sweep_start_ray_index[i]:sweep_end_ray_index[i]] dvel_cor_sweep=dvel_cor[*,sweep_start_ray_index[i]:sweep_end_ray_index[i]] power_sweep=power[*,sweep_start_ray_index[i]:sweep_end_ray_index[i]] zdr_cor_sweep=zdr_cor[*,sweep_start_ray_index[i]:sweep_end_ray_index[i]] phidp_cor_sweep=phidp_cor[*,sweep_start_ray_index[i]:sweep_end_ray_index[i]] endif dbz_sweep=dbz[*,sweep_start_ray_index[i]:sweep_end_ray_index[i]] rhohv_sweep=rhohv[*,sweep_start_ray_index[i]:sweep_end_ray_index[i]] zdr_sweep=zdr[*,sweep_start_ray_index[i]:sweep_end_ray_index[i]] phidp_sweep=phidp[*,sweep_start_ray_index[i]:sweep_end_ray_index[i]] kdp_sweep=kdp[*,sweep_start_ray_index[i]:sweep_end_ray_index[i]] hd_sweep=hd[*,sweep_start_ray_index[i]:sweep_end_ray_index[i]] ; Size of sweep array s=size(dbz_sweep,/dimensions) nrays=s[1] nbins=s[0] nrays_sweep[i]=nrays nbins_sweep[i]=nbins ; Time per gate in the sweep if i lt n_elements(sweep_number)-2 then begin delta_t=(time[sweep_start_ray_index[i+1]]-time[sweep_start_ray_index[i]])/float(nrays)/float(nbins) endif print,delta_t,'seconds per gate' ; Arrays to hold the sweep data n=0L ;number of points in sweep numcount=long(nrays)*long(nbins) x=fltarr(numcount) y=fltarr(numcount) z=fltarr(numcount) if temp ne !NULL then begin temp_sweep_array=fltarr(numcount) air_sweep_array=fltarr(numcount) hclas_sweep_array=fltarr(numcount) rain_sweep_array=fltarr(numcount) dvel_sweep_array=fltarr(numcount) dvel_cor_sweep_array=fltarr(numcount) power_sweep_array=fltarr(numcount) zdr_cor_sweep_array=fltarr(numcount) phidp_cor_sweep_array=fltarr(numcount) endif dbz_sweep_array=fltarr(numcount) rhohv_sweep_array=fltarr(numcount) zdr_sweep_array=fltarr(numcount) phidp_sweep_array=fltarr(numcount) kdp_sweep_array=fltarr(numcount) hd_sweep_array=fltarr(numcount) azimuth_sweep_array=fltarr(numcount) elevation_sweep_array=fltarr(numcount) time_sweep_array=dblarr(numcount) lat_sweep_array=fltarr(numcount) lon_sweep_array=fltarr(numcount) print,nrays,'nrays',nbins,'nbins' for j=0,nrays-1 do begin for k=0,nbins-1 do begin ;bins ; azimuth, 0degNorth and positive clockwise azimuth_sweep_array[n]=azimuth_sweep[j] ; elevation, positive above the reference plane elevation_sweep_array[n]=elevation_sweep[j] ; pseudo radius of earth pse_km=(4./3.)*6374.0 ;km pse_m=pse_km*1000.0 ;m ; height r=range, lambda=azimuth phi=elevation term1=range[k]^2+pse_m^2+2.0*range[k]*pse_m*sin(elevation_sweep[j]*(!pi/180.)) height=sqrt(term1)-pse_m+calt ; msl z[n]=height ; distance x[n]=range[k]*cos(elevation_sweep[j]*!pi/180.)*sin(azimuth_sweep[j]*!pi/180.) y[n]=range[k]*cos(elevation_sweep[j]*!pi/180.)*cos(azimuth_sweep[j]*!pi/180.) ; radar variables if temp ne !null then begin temp_sweep_array[n]=temp_sweep[k,j] air_sweep_array[n]=air_sweep[k,j] hclas_sweep_array[n]=hclas_sweep[k,j] rain_sweep_array[n]=rain_sweep[k,j] dvel_sweep_array[n]=dvel_sweep[k,j] dvel_cor_sweep_array[n]=dvel_cor_sweep[k,j] power_sweep_array[n]=power_sweep[k,j] zdr_cor_sweep_array[n]=zdr_cor_sweep[k,j] phidp_cor_sweep_array[n]=phidp_cor_sweep[k,j] endif dbz_sweep_array[n]=dbz_sweep[k,j] rhohv_sweep_array[n]=rhohv_sweep[k,j] zdr_sweep_array[n]=zdr_sweep[k,j] phidp_sweep_array[n]=phidp_sweep[k,j] kdp_sweep_array[n]=kdp_sweep[k,j] hd_sweep_array[n]=hd_sweep[k,j] ; Time seconds time_sweep_array[n]=time_sweep[j]+double(k)*delta_t ; Lat and Lon arc_dist=range[k]/(6371D*1000D) ;radius of earth in meters result=ll_arc_distance([clon,clat],arc_dist,azimuth_sweep[j],/degrees) lon_sweep_array[n]=result[0] lat_sweep_array[n]=result[1] n=n+1l endfor ;end of loop through range gates endfor ;end of loop through rays, finished sweep ; Calculate colongitude because the pacific crosses the 180 line colon_sweep_array=lon_sweep_array result=where(lon_sweep_array lt 0,count) if count gt 0 then colon[result]=colon[result]+360.0 ; Calculate julian day julian_day_sweep_array=julday(smm,sdd,syy,shh,smi,sss)+double(time_sweep_array)/(60D*60D*24D) ; Save the sweep data if x_all eq !NULL then begin x_all=x y_all=y z_all=z if temp ne !NULL then begin temp_sweep_array_all=temp_sweep_array air_sweep_array_all=air_sweep_array hclas_sweep_array_all=hclas_sweep_array rain_sweep_array_all=rain_sweep_array dvel_sweep_array_all=dvel_sweep_array dvel_cor_sweep_array_all=dvel_cor_sweep_array power_sweep_array_all=power_sweep_array zdr_cor_sweep_array_all=zdr_cor_sweep_array phidp_cor_sweep_array_all=phidp_cor_sweep_array endif dbz_sweep_array_all=dbz_sweep_array rhohv_sweep_array_all=rhohv_sweep_array zdr_sweep_array_all=zdr_sweep_array phidp_sweep_array_all=phidp_sweep_array kdp_sweep_array_all=kdp_sweep_array hd_sweep_array_all=hd_sweep_array azimuth_sweep_array_all=azimuth_sweep_array elevation_sweep_array_all=elevation_sweep_array time_sweep_array_all=time_sweep_array lat_sweep_array_all=lat_sweep_array lon_sweep_array_all=lon_sweep_array colon_sweep_array_all=colon_sweep_array julian_day_sweep_array_all=julian_day_sweep_array endif else begin x_all=[x_all,x] y_all=[y_all,y] z_all=[z_all,z] if temp ne !NULL then begin temp_sweep_array_all=[temp_sweep_array_all,temp_sweep_array] air_sweep_array_all=[air_sweep_array_all,air_sweep_array] hclas_sweep_array_all=[hclas_sweep_array_all,hclas_sweep_array] rain_sweep_array_all=[rain_sweep_array_all,rain_sweep_array] dvel_sweep_array_all=[dvel_sweep_array_all,dvel_sweep_array] dvel_cor_sweep_array_all=[dvel_cor_sweep_array_all,dvel_cor_sweep_array] power_sweep_array_all=[power_sweep_array_all,power_sweep_array] zdr_cor_sweep_array_all=[zdr_cor_sweep_array_all,zdr_cor_sweep_array] phidp_cor_sweep_array_all=[phidp_cor_sweep_array_all,phidp_cor_sweep_array] endif dbz_sweep_array_all=[dbz_sweep_array_all,dbz_sweep_array] rhohv_sweep_array_all=[rhohv_sweep_array_all,rhohv_sweep_array] zdr_sweep_array_all=[zdr_sweep_array_all,zdr_sweep_array] phidp_sweep_array_all=[phidp_sweep_array_all,phidp_sweep_array] kdp_sweep_array_all=[kdp_sweep_array_all,kdp_sweep_array] hd_sweep_array_all=[hd_sweep_array_all,hd_sweep_array] azimuth_sweep_array_all=[azimuth_sweep_array_all,azimuth_sweep_array] elevation_sweep_array_all=[elevation_sweep_array_all,elevation_sweep_array] time_sweep_array_all=[time_sweep_array_all,time_sweep_array] lat_sweep_array_all=[lat_sweep_array_all,lat_sweep_array] lon_sweep_array_all=[lon_sweep_array_all,lon_sweep_array] colon_sweep_array_all=[colon_sweep_array_all,colon_sweep_array] julian_day_sweep_array_all=[julian_day_sweep_array_all,julian_day_sweep_array] endelse ;*************************************************** ; Plot the cloudy rays if 1 eq 0 then begin ; Create plot pxdim=900 & pydim=1000 xl=0.1 & xr=0.95 yb=0.03 & yt=0.95 sx=0.07 sy=0.03 numplots_x=1 numplots_y=9 position_plots,xl,xr,yb,yt,sx,sy,numplots_x,numplots_y,pos dx=0.01 & dy=0.01 & fs1=12 ; Min and max in the sweep ridx=where(dbz_sweep_array ne -9999,cidx) dbz_min=min(dbz_sweep_array[ridx]) dbz_max=max(dbz_sweep_array[ridx]) rhohv_min=min(rhohv_sweep_array[ridx]) rhohv_max=max(rhohv_sweep_array[ridx]) zdr_min=-5.0;min([zdr_sweep_array,zdr_cor_sweep_array]) zdr_max=5.0;max([zdr_sweep_array,zdr_cor_sweep_array]) air_sweep_array=make_array(numcount,value=3) kdp_min=min(kdp_sweep_array) kdp_max=max(kdp_sweep_array) hclas_min=0.5;min(air_sweep_array) hclas_max=10.5;max(air_sweep_array) pidxs=0 pidxe=600 sidx=0 for j=1,nrays do begin ;loop through rays sub_dbz_sweep_array=dbz_sweep_array[sidx:sidx+nbins-1] r=where(sub_dbz_sweep_array ne -9999,c) print,sidx,sidx+nbins-1 if c gt 10 then begin ;if max(air_sweep_array[sidx+400:sidx+nbins-1]) eq 3 then begin pnum=0 p0=plot(sub_dbz_sweep_array[pidxs:pidxe],/buffer,dimensions=[pxdim,pydim],/xstyle,$ yrange=[dbz_min,dbz_max],position=pos[pnum,*],ytitle='dbz',font_size=fs1) t0=text(pos[pnum,0]+1*dx,pos[pnum,3]+1*dy,'corrected reflectivity (dbz)',font_size=fs1) pnum=1 sub_rhohv_sweep_array=rhohv_sweep_array[sidx:sidx+nbins-1] p1=plot(sub_rhohv_sweep_array[pidxs:pidxe],/current,dimensions=[pxdim,pydim],/xstyle,$ yrange=[rhohv_min,rhohv_max],position=pos[pnum,*],ytitle='rhohv',font_size=fs1,xtickformat='(A1)') t1=text(pos[pnum,0]+1*dx,pos[pnum,3]+1*dy,'cross correlation ratio',font_size=fs1) if air ne !NULL then begin pnum=2 sub_air_sweep_array=air_sweep_array[sidx:sidx+nbins-1] p2=plot(sub_air_sweep_array[pidxs:pidxe],/current,dimensions=[pxdim,pydim],/xstyle,$ yrange=[-0.5,3.5],position=pos[pnum,*],ytitle='air',font_size=fs1,xtickformat='(A1)') t2=text(pos[pnum,0]+1*dx,pos[pnum,3]+1*dy,'air echo classification 0: N/A, 1: Clutter, 2: Clear Air, 3: Meteorological echoes',font_size=fs1) endif pnum=3 sub_zdr_sweep_array=zdr_sweep_array[sidx:sidx+nbins-1] p0=plot(sub_zdr_sweep_array[pidxs:pidxe],/current,dimensions=[pxdim,pydim],/xstyle,thick=1,$ yrange=[zdr_min,zdr_max],position=pos[pnum,*],ytitle='zdr',font_size=fs1,xtickformat='(A1)') t1=text(pos[pnum,0]+1*dx,pos[pnum,3]+1*dy,'differential_reflectivity (db)',font_size=fs1) if zdr_cor ne !NULL then begin sub_zdr_cor_sweep_array=zdr_cor_sweep_array[sidx:sidx+nbins-1] p3=plot(sub_zdr_cor_sweep_array[pidxs:pidxe],/overplot,color='red',thick=2) t1=text(pos[pnum,0]+30*dx,pos[pnum,3]+1*dy,'corrected_differential_reflectivity (db)',font_size=fs1,color='red') endif pnum=4 if phidp_cor ne !NULL then begin sub_phidp=phidp_cor_sweep_array[sidx:sidx+nbins-1] sub_air=air_sweep_array[sidx:sidx+nbins-1] r=where(sub_air eq 3,c) phidp_min=min(sub_phidp[r]) phidp_max=max(sub_phidp[r])+1.0 p4=plot(sub_phidp[pidxs:pidxe],/current,dimensions=[pxdim,pydim],/xstyle,$ yrange=[phidp_min,phidp_max],position=pos[pnum,*],ytitle='phidp',font_size=fs1,xtickformat='(A1)') t1=text(pos[pnum,0]+1*dx,pos[pnum,3]+1*dy,'corrected_differential_phase (degrees)',font_size=fs1) endif else begin sub_phidp=phidp_sweep_array[sidx:sidx+nbins-1] r=where(sub_phidp ne -9999) phidp_min=min(sub_phidp[r]) phidp_max=max(sub_phidp[r]) p4=plot(sub_phidp[pidxs:pidxe],/current,dimensions=[pxdim,pydim],/xstyle,$ yrange=[phidp_min,phidp_max],position=pos[pnum,*],ytitle='phidp',font_size=fs1,xtickformat='(A1)') t1=text(pos[pnum,0]+1*dx,pos[pnum,3]+1*dy,'differential_phase (degrees)',font_size=fs1) endelse pnum=5 sub_kdp=kdp_sweep_array[sidx:sidx+nbins-1] sub_air=air_sweep_array[sidx:sidx+nbins-1] r=where(sub_air eq 3 and sub_kdp ne -9999,c) kdp_min=min(sub_kdp[r]) kdp_max=max(sub_kdp[r]) p5=plot(sub_kdp[pidxs:pidxe],/current,dimensions=[pxdim,pydim],/xstyle,$ yrange=[kdp_min,kdp_max],position=pos[pnum,*],ytitle='kdp',font_size=fs1,xtickformat='(A1)') t1=text(pos[pnum,0]+1*dx,pos[pnum,3]+1*dy,'corrected_specific_differential_phase (degrees/km)',font_size=fs1) pnum=6 if hcla ne !NULL then begin p5=plot(hclas_sweep_array[sidx:sidx+nbins-1],/current,dimensions=[pxdim,pydim],/xstyle,$ yrange=[hclas_min,hclas_max],position=pos[pnum,*],ytitle='hclas',font_size=fs1,xtickformat='(A1)') ctable=['red','orange','lime green','green','dodger blue','blue','purple','blue violet','sienna','magenta'] hydro_str=['1=Drizzle','2=Rain','3=Ice Crystals','4=Aggregates','5=Wet Snow','6=Vertical Ice','7=LD Graupel','8=HD Graupel','9=Hail','10=Big Drops'] hydro_x=[0,8,14,25,36,46,57,68,79,84] sub_array=hclas_sweep_array[sidx:sidx+nbins-1] sub_x=findgen(n_elements(sub_array)) for k=0,9 do begin r=where(sub_array eq k+1,c) if c gt 0 then begin p5=plot(sub_x[r],sub_array[r],/overplot,symbol='o',sym_size=0.3,/sym_filled,color=ctable[k],linestyle=6) endif t6=text(pos[pnum,0]+hydro_x[k]*dx,pos[pnum,3]-2*dy,hydro_str[k],font_size=10,color=ctable[k]) endfor t1=text(pos[pnum,0]+1*dx,pos[pnum,3]-4*dy,'radar_echo_classification',font_size=fs1) endif pnum=7 sub_x=x[sidx:sidx+nbins-1] sub_y=y[sidx:sidx+nbins-1] sub_z=z[sidx:sidx+nbins-1] p7=plot( sqrt((sub_x[pidxs:pidxe])^2+(sub_y[pidxs:pidxe])^2)/1000.0,/current,/xstyle,$ position=pos[pnum,*],font_size=fs1,ytitle='distance (km)',xtickformat='(A1)');,$ pnum=8 p8=plot(sub_z[pidxs:pidxe]/1000.0,/current,/xstyle,xtickformat='(A1)',$ position=pos[pnum,*],font_size=fs1,ytitle='Z (km)',yrange=[min(z)/1000.0,max(z)/1000.0]) caldat,julian_day_sweep_array[sidx],jmm,jdd,jyy,jhr,jmi,jss sweep_str=string(jyy,format='(I4)')+string(jmm,format='(I02)')+string(jdd,format='(I02)')+' '+$ string(jhr,format='(I02)')+':'+string(jmi,format='(I02)')+':'+string(jss,format='(I02)') sweep_time_str=string(jyy,format='(I4)')+string(jmm,format='(I02)')+string(jdd,format='(I02)')+'.'+$ string(jhr,format='(I02)')+string(jmi,format='(I02)')+string(jss,format='(I02)') elevation_str=string(elevation_sweep_array[sidx],format='(F05.2)') azimuth_str=string(azimuth_sweep_array[sidx],format='(F07.2)') if azimuth_sweep_array[sidx] gt 177.0 then print,azimuth_str t1=text(0.1,0.97,'Cband '+sweep_str+' Elev='+elevation_str+' Azm='+azimuth_str,font_size=fs1) imagename='ray.'+sweep_time_str+'.el'+elevation_str+'.az'+azimuth_str print,imagename p0.save,imagename+'.'+string(pidxe,format='(I03)')+'.png',height=pydim endif ;cloudy, plot sidx=sidx+nbins endfor ;end of loop through rays endif ;do the plot endfor ;end of loop through sweeps endfor ;end of loop through files endif ;end of found files x=x_all y=y_all z=z_all if temp ne !NULL then begin temp_sweep_array=temp_sweep_array_all air_sweep_array=air_sweep_array_all hclas_sweep_array=hclas_sweep_array_all rain_sweep_array=rain_sweep_array_all dvel_sweep_array=dvel_sweep_array_all dvel_cor_sweep_array=dvel_cor_sweep_array_all power_sweep_array=power_sweep_array_all zdr_cor_sweep_array=zdr_cor_sweep_array_all phidp_cor_sweep_array=phidp_cor_sweep_array_all endif dbz_sweep_array=dbz_sweep_array_all rhohv_sweep_array=rhohv_sweep_array_all zdr_sweep_array=zdr_sweep_array_all phidp_sweep_array=phidp_sweep_array_all kdp_sweep_array=kdp_sweep_array_all azimuth_sweep_array=azimuth_sweep_array_all elevation_sweep_array=elevation_sweep_array_all time_sweep_array=time_sweep_array_all lat_sweep_array=lat_sweep_array_all lon_sweep_array=lon_sweep_array_all colon_sweep_array=colon_sweep_array_all julian_day_sweep_array=julian_day_sweep_array_all ;******************* ; Image a sweep ;******************* ; Number of elements in each sweep sweep_num=nrays_sweep*nbins_sweep sidx=0 & eidx=sweep_num[0] ; Subset out the sweep lon_sweep_array=lon_sweep_array[sidx:eidx] colon_sweep_array=colon_sweep_array[sidx:eidx] lat_sweep_array=lat_sweep_array[sidx:eidx] dbz_sweep_array=dbz_sweep_array[sidx:eidx] elevation_sweep_array=elevation_sweep_array[sidx:eidx] julian_day_sweep_array=julian_day_sweep_array[sidx:eidx] ; Lat,lon limits llon=min(lon_sweep_array) rlon=max(lon_sweep_array) llat=min(lat_sweep_array) ulat=max(lat_sweep_array) lcolon=llon if lcolon lt 0 then lcolon=lcolon+360.0 rcolon=rlon if rcolon lt 0 then rcolon=rcolon+360.0 ; 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 ; ; First grid cloud - these array only hold points with cloud r=where(dbz_sweep_array gt -9999,c) dbz_cloud=dbz_sweep_array[r] colon_cloud=colon_sweep_array[r] lat_cloud=lat_sweep_array[r] ;grid_input,colon_cloud,lat_cloud,dbz_cloud,xyz,newdata_var,$ ; /degree,/sphere,duplicates='Max',epsilon=0.01 ;lon = !radeg * atan(xyz[1,*],xyz[0,*]) & lat = !radeg * asin(xyz[2,*]) ;qhull,lon,lat,triangles,/delaunay,sphere=s qhull,colon_cloud,lat_cloud,triangles,/delaunay,sphere=s ;grid_cloud=griddata(lon,lat,newdata_var,$ grid_cloud=griddata(colon_cloud,lat_cloud,dbz_cloud,$ xout=grid_colon,yout=grid_lat,$ /grid,/degrees,/sphere,triangles=triangles,$ ;/KRIGING,min_points=10,sectors=8,empty_sectors=4,$ /KRIGING,min_points=20,sectors=8,empty_sectors=4,$ ;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) ; Second create mask grid_input,colon_sweep_array,lat_sweep_array,dbz_sweep_array,xyz,newdata_var,$ /degree,/sphere,duplicates='Max',epsilon=0.001 lon = !radeg * atan(xyz[1,*],xyz[0,*]) & lat = !radeg * asin(xyz[2,*]) qhull,lon,lat,triangles,/delaunay,sphere=s ;qhull,colon_sweep_array,lat_sweep_array,triangles,/delaunay,sphere=s grid_mask=griddata(lon,lat,newdata_var,$ ;grid_mask=griddata(colon_sweep_array,lat_sweep_array,dbz_sweep_array,$ xout=grid_colon,yout=grid_lat,$ /grid,/degrees,/sphere,triangles=triangles,$ ;/KRIGING,min_points=10,sectors=8,empty_sectors=4,$ /KRIGING,min_points=16,sectors=8,empty_sectors=4,$ ;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=-8888) ; Put the two data sets together grid_var=grid_cloud r=where(grid_mask eq -8888,c) if c gt 0 then grid_var[r]=-8888 ; Create plot pxdim=900 & pydim=900 xl=0.1 & xr=0.95 yb=0.05 & yt=0.5 sx=0.07 sy=0.05 numplots_x=3 numplots_y=3 position_plots,xl,xr,yb,yt,sx,sy,numplots_x,numplots_y,pos cbpos=pos cbpos[*,0]=pos[*,2]+0.14 cbpos[*,2]=cbpos[*,0]+0.02 dx=0.01 & dy=0.01 & fs1=12 ; Original data pnum=0 data_var=dbz_sweep_array r=where(data_var ne -9999.0,c) start_bin=min(data_var[r]) & start_bin_max=start_bin end_bin=max(data_var) & end_bin_max=end_bin dbin=(end_bin-start_bin)/100.0 r=where(data_var eq -9999.0,c) hist_generic,data_var,start_bin,end_bin,dbin,bins,data_freq,data_counts p0=plot(bins,data_freq,/hist,/buffer,position=pos[pnum,*],ytitle='dbz_sweep_array',$ dimensions=[pxdim,pydim],font_size=fs1,xrange=[start_bin_max,end_bin_max],xtitle='dbz') t1=text(pos[pnum,0]+1*dx,pos[pnum,3]-3*dy,'-9999='+strcompress(c,/remove_all),font_size=fs1) pnum=1 data_var=colon_sweep_array r=where(data_var ne -9999.0,c) start_bin=min(data_var[r]) end_bin=max(data_var) dbin=(end_bin-start_bin)/100.0 r=where(data_var eq -9999.0,c) hist_generic,data_var,start_bin,end_bin,dbin,bins,data_freq,data_counts p0=plot(bins,data_freq,/hist,/current,position=pos[pnum,*],font_size=fs1,xtitle='lon') t1=text(pos[pnum,0]+1*dx,pos[pnum,3]-3*dy,'-9999='+strcompress(c,/remove_all),font_size=fs1) pnum=2 data_var=lat_sweep_array r=where(data_var ne -9999.0,c) start_bin=min(data_var[r]) end_bin=max(data_var) dbin=(end_bin-start_bin)/100.0 r=where(data_var eq -9999.0,c) hist_generic,data_var,start_bin,end_bin,dbin,bins,data_freq,data_counts p0=plot(bins,data_freq,/hist,/current,position=pos[pnum,*],font_size=fs1,xtitle='lat') t1=text(pos[pnum,0]+1*dx,pos[pnum,3]-3*dy,'-9999='+strcompress(c,/remove_all),font_size=fs1) ; cloudy points only if 1 eq 1 then begin pnum=3 data_var=dbz_cloud;newdata_var r=where(data_var ne -9999.0,c) start_bin=min(data_var[r]) end_bin=max(data_var) dbin=(end_bin-start_bin)/100.0 r=where(data_var eq -9999.0,c) hist_generic,data_var,start_bin,end_bin,dbin,bins,data_freq,data_counts p0=plot(bins,data_freq,/hist,/current,position=pos[pnum,*],ytitle='dbz_cloud',font_size=fs1,$ xrange=[start_bin_max,end_bin_max]) t1=text(pos[pnum,0]+1*dx,pos[pnum,3]-3*dy,'-9999='+strcompress(c,/remove_all),font_size=fs1) pnum=4 data_var=colon_cloud r=where(data_var ne -9999.0,c) start_bin=min(data_var[r]) end_bin=max(data_var) dbin=(end_bin-start_bin)/100.0 r=where(data_var eq -9999.0,c) hist_generic,data_var,start_bin,end_bin,dbin,bins,data_freq,data_counts p0=plot(bins,data_freq,/hist,/current,position=pos[pnum,*],font_size=fs1) t1=text(pos[pnum,0]+1*dx,pos[pnum,3]-3*dy,'-9999='+strcompress(c,/remove_all),font_size=fs1) pnum=5 data_var=lat_cloud r=where(data_var ne -9999.0,c) start_bin=min(data_var[r]) end_bin=max(data_var) dbin=(end_bin-start_bin)/100.0 r=where(data_var eq -9999.0,c) hist_generic,data_var,start_bin,end_bin,dbin,bins,data_freq,data_counts p0=plot(bins,data_freq,/hist,/current,position=pos[pnum,*],font_size=fs1) t1=text(pos[pnum,0]+1*dx,pos[pnum,3]-3*dy,'-9999='+strcompress(c,/remove_all),font_size=fs1) endif ; run through griddata if 1 eq 1 then begin pnum=6 data_var=grid_var r=where(data_var ne -9999 and data_var ne -8888,c) start_bin=min(data_var[r]) end_bin=max(data_var) dbin=(end_bin-start_bin)/100.0 r=where(data_var eq -9999.0,c) hist_generic,data_var,start_bin,end_bin,dbin,bins,data_freq,data_counts p0=plot(bins,data_freq,/hist,/current,position=pos[pnum,*],ytitle='grid_var',font_size=fs1,$ xrange=[start_bin_max,end_bin_max]) r=where(data_var eq -9999,c) t1=text(pos[pnum,0]+1*dx,pos[pnum,3]-3*dy,'-9999='+strcompress(c,/remove_all),font_size=fs1) r=where(data_var eq -8888,c) t1=text(pos[pnum,0]+1*dx,pos[pnum,3]-6*dy,'-8888='+strcompress(c,/remove_all),font_size=fs1) pnum=7 data_var=grid_colon r=where(data_var ne -9999 and data_var ne -8888,c) start_bin=min(data_var[r]) end_bin=max(data_var) dbin=(end_bin-start_bin)/100.0 r=where(data_var eq -9999.0,c) hist_generic,data_var,start_bin,end_bin,dbin,bins,data_freq,data_counts p0=plot(bins,data_freq,/hist,/current,position=pos[pnum,*],font_size=fs1) t1=text(pos[pnum,0]+1*dx,pos[pnum,3]-5*dy,'-9999='+strcompress(c,/remove_all),font_size=fs1) pnum=8 data_var=grid_lat r=where(data_var ne -9999 and data_var ne -8888,c) start_bin=min(data_var[r]) end_bin=max(data_var) dbin=(end_bin-start_bin)/100.0 r=where(data_var eq -9999.0,c) hist_generic,data_var,start_bin,end_bin,dbin,bins,data_freq,data_counts p0=plot(bins,data_freq,/hist,/current,position=pos[pnum,*],font_size=fs1) t1=text(pos[pnum,0]+1*dx,pos[pnum,3]-5*dy,'-9999='+strcompress(c,/remove_all),font_size=fs1) xl=0.1 & xr=0.95 yb=0.6 & yt=0.97 sx=0.07 sy=0.05 numplots_x=1 numplots_y=1 position_plots,xl,xr,yb,yt,sx,sy,numplots_x,numplots_y,pos cbpos=pos cbpos[*,0]=pos[*,2]-0.14 cbpos[*,2]=cbpos[*,0]+0.02 ; 253=white 255=gray ; Plot ppi image pnum=0 r=where(grid_var gt -8888) dmin=start_bin_max dmax=end_bin_max data_image=bytscl(grid_var,top=top_color,min=dmin,max=dmax) result=where(grid_var eq -8888,count) if count gt 0 then data_image[result]=255 ;gray result=where(grid_var eq -9999,count) if count gt 0 then data_image[result]=253 ;white p9=image(data_image,grid_colon,grid_lat,limit=[llat,lcolon,ulat,rcolon],$ /current,dimensions=[pxdim,pydim],grid_units='degrees',label_position=0,linestyle='dot',$ center_longitude=(rcolon-lcolon)/2.0,$ /box_axes,rgb_table=mytable,$;,min_value=dmin,max_value=dmax,$ map_projection='Mercator',label_format='MapGrid_Labels',$ position=pos[pnum,*],font_size=fs1) p1=mapcontinents(color='black') caldat,julian_day_sweep_array[0],jmm,jdd,jyy,jhr,jmi,jss sweep_str=string(jyy,format='(I4)')+string(jmm,format='(I02)')+string(jdd,format='(I02)')+' '+$ string(jhr,format='(I02)')+':'+string(jmi,format='(I02)')+':'+string(jss,format='(I02)') sweep_time_str=string(jyy,format='(I4)')+string(jmm,format='(I02)')+string(jdd,format='(I02)')+'.'+$ string(jhr,format='(I02)')+string(jmi,format='(I02)')+string(jss,format='(I02)') t1=text(pos[pnum,0],pos[pnum,3],'Cband '+sweep_str,font_size=fs1) t1=text(pos[pnum,0],pos[pnum,3]-3*dy,'elevation='+string(elevation_sweep_array[0],format='(F6.2)'),font_size=fs1) c0=colorbar($ rgb_table=cbtable,$ range=[dmin,dmax],$ orientation=1,position= reform(cbpos[pnum,*]),font_size=fs1,$ tickdir=1,border_on=1,title='DbZ') endif p0.save,'test_grid_dolan.png',height=pydim stop ; ; ; ; ; ; ; ; dmin=10 ; dmax=35 ; data_image=bytscl(grid_var,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]) ; ; p0=image(data_image,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(elevation_sweep[0],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') ; ; p0.save,'test.png',height=pydim ; stop ; ; ; ; 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+'.2.png' ; ;p0.save,imagename,height=pydim ; ;stop ; ; ; ; ; ; endfor ;end of loop through sweeps ;endfor ;end of loop through files ;endif ;end of found files ; stop ; ; 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 ; ;******************* ; ; ;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+'.2.png' ;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