;********************************* ; Probably start entering the retrieval code at line 206 ; I only loop through the 3 largest features. This can be ; changed at line 122 PRO ProgramDefineROI_convcells_ret do_ret='no' ;do the retrieval equations rho_filter='yes' ;filter kdp with rho < 0.95 do_ellipse='no' ;Fit the feature to an ellipse ;******************************* ; Read in the data ;******************************* ;fdir='/Volumes/mace-group4/Capricorn2/CBand/cfradial/' ;fname='9776HUB-PPIVol-20180209-140504-0000.cfradial.pyart.81.401.nc' fdir='/Volumes/mace-group5/field_programs/investigator/in2018_v01/cband/' fname='9776HUB-PPIVol-20180209-105330-0000.cfradial.pyart.81.401.nc' fid=ncdf_open(fdir+fname) xid=ncdf_varid(fid,'x') & ncdf_varget,fid,xid,x_m xid=ncdf_varid(fid,'y') & ncdf_varget,fid,xid,y_m xid=ncdf_varid(fid,'z') & ncdf_varget,fid,xid,z_m vid=ncdf_varid(fid,'cross_correlation_ratio') & ncdf_varget,fid,vid,rho vid=ncdf_varid(fid,'specific_differential_phase') if vid eq -1 then vid=ncdf_varid(fid,'corrected_specific_differential_phase') ncdf_varget,fid,vid,kdp vid=ncdf_varid(fid,'corrected_differential_reflectivity') & ncdf_varget,fid,vid,zdr vid=ncdf_varid(fid,'corrected_reflectivity') & ncdf_varget,fid,vid,dbz vid=ncdf_varid(fid,'velocity') & ncdf_varget,fid,vid,vcy ncdf_close,fid ; Convert m to km x_km=x_m/1000.0 y_km=y_m/1000.0 z_km=z_m/1000.0 ; Number of elements in each dimension numx=n_elements(x_m) numy=n_elements(y_m) numz=n_elements(z_m) ; Arrays of index numbers in each direction x=indgen(numx) y=indgen(numy) z=indgen(numz) ;************************************************** ; 2d array horizontal slice for identifying features ;************************************************** ; 2km height level - changed this to 1km sidx=where(z_km ge 1.0) sidx=sidx[0] ; horizontal slice at sidx height level img_dbz=dbz[*,*,sidx] img_kdp=kdp[*,*,sidx] img_zdr=zdr[*,*,sidx] img_vcy=vcy[*,*,sidx] img_rho=rho[*,*,sidx] dims = SIZE(img_dbz, /DIMENSIONS) ;************************************************** ; Create a mask that identifies the pixels of interest ; 1=pixel gt -15 dbz, 0=pixel le -15 dbz ;************************************************** threshImg=(img_dbz gt -15 ) ;************************************************** ; Get rid of gaps, applying a 3x3 element to the image ; using the erosion and dilation morphological ; operators. ;************************************************** ;strucElem = REPLICATE(1, 3, 3) ;threshImg = ERODE(DILATE(TEMPORARY(threshImg), $ ; strucElem), strucElem) ;************************************************* ; Extract the contours of the thresholded image. ; path_info=path information about the contours ; path_xy=the xy coordinates of the contours ;************************************************* CONTOUR, threshImg, $ ;x,y,$ ;can't use this because compute mask only uses pixels LEVEL = 1, $ XMARGIN = [0, 0], YMARGIN = [0, 0], $ /NOERASE, PATH_INFO = pathInfo, PATH_XY = pathXY, $ XSTYLE = 5, YSTYLE = 5, /PATH_DATA_COORDS ;************************************************** ; Color tables ;************************************************** ; Top is the last color to scale 256 colors, 0-255 top_color=252 ; Colortable 0-253 254=white mytable=colortable(43,ncolors=253) ; 253=white 254=hot pink ;gray=255 mytable=[mytable,transpose([255,255,255]),transpose([238,18,137]),transpose([180,180,180])] 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,*] ; 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] ;************************************************** ; Set up the plot ;************************************************** ; ; Set up the positions pxdim=1000 & pydim=850 xl=0.06 & xr=0.90 yb=0.05 & yt=0.92 sx=0.13 sy=0.1 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.07 cbpos[*,2]=cbpos[*,0]+0.007 cbpos[*,1]=cbpos[*,1]+0.02 cbpos[*,3]=cbpos[*,3]-0.02 ; For text positioning dx=0.01 & dy=0.01 ; Image 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]) ;************************************************** ; For each region, feed the contours into an IDLgrROI ; object for display with DRAW_ROI. ;************************************************** ; Loop through first 3 contours FOR I = 0, 3 do begin ; Loop through every contour ;FOR I = 0, (N_ELEMENTS(pathInfo) - 1 ) DO BEGIN ; Initialize the IDLgrROI object with the contour ; information of the current region with the FOR ; loop. line = [LINDGEN(pathInfo(I).N), 0] ;oROI = OBJ_NEW('IDLanROI', $ ;an roi and gr roi are the same oROI = OBJ_NEW('IDLgrROI', $ (pathXY(*,pathInfo(I).OFFSET + line))[0, *], $ (pathXY(*,pathInfo(I).OFFSET + line))[1, *]) ; Get out the x and y vectors of the ROI oROI -> GetProperty, DATA = ROIdata x_roi = reform(ROIdata[0,*]) y_roi = reform(ROIdata[1,*]) ; Use ComputeMask in conjunction with ; IMAGE_STATISTICS to obtain the number of pixels ; covered by the regions when displayed. maskResult = oROI -> ComputeMask(DIMENSIONS = dims) IMAGE_STATISTICS, img_dbz, MASK = maskResult, $ COUNT = maskArea, MEAN=maskMean ; Use ComputeGeometry to obtain the geometric area ; and perimeter of each region where 1 pixel = ; 1.2 x 1.2 mm. ROIStats = oROI -> ComputeGeometry($ AREA = geomArea, PERIMETER = perimeter,$ SPATIAL_SCALE = [500.,500.,500.]) ;From the ROI mask, create a binary mask, consisting of only zeros and ones. ;Multiply the binary mask times the original image to retain only the portion ;of the image that was defined in the original ROI: mask = (maskResult GT 0) maskImg_dbz = img_dbz * mask maskImg_kdp = img_kdp * mask maskImg_zdr = img_zdr * mask maskImg_vcy = img_vcy * mask maskImg_rho = img_rho * mask ;Using the minimum and maximum values of the ROI array, create a cropped array, ;cropImg, and get its dimensions: ;cropImg = maskImg[min(x):max(x), min(y): max(y)] x1=where(x le min(x_roi),c) & x1=x1[-1] x2=where(x ge max(x_roi),c) & x2=x2[0] y1=where(y le min(y_roi),c) & y1=y1[-1] y2=where(y ge max(y_roi),c) & y2=y2[0] cropImg_dbz = maskImg_dbz[x1:x2,y1:y2] cropImg_zdr = maskImg_zdr[x1:x2,y1:y2] cropImg_kdp = maskImg_kdp[x1:x2,y1:y2] cropImg_vcy = maskImg_vcy[x1:x2,y1:y2] cropImg_rho = maskImg_rho[x1:x2,y1:y2] cropx=x[x1:x2] cropy=y[y1:y2] cropDims = SIZE(cropImg_dbz, /DIMENSIONS) ; Get the 3d volume of the cropimage volImg_dbz=dbz[x1:x2,y1:y2,*] volImg_zdr=zdr[x1:x2,y1:y2,*] volImg_kdp=kdp[x1:x2,y1:y2,*] volImg_vcy=vcy[x1:x2,y1:y2,*] volImg_rho=rho[x1:x2,y1:y2,*] ;################################## ; Run retrieval on cropImg ;################################## if do_ret eq 'yes' then begin common in_situ_snow_properties, Lx,Nx,alpha,mass_agg,dm_agg,dbz_cband_agg,$ mass_rimed,dm_rimed,dbz_cband_rimed,mass_graup,dm_graup,dbz_cband_graup common observables_snow_response, dBZe_h_resp,zdr_resp,kdp_resp,rhohv_resp,$ precip_rate_resp,iwc_resp,alpha_resp,bm_resp,aspect_ratio_resp,Dm_resp,$ temperature_resp,sigma_resp ; following read fills in the observables_snow_response common fid=ncdf_open('C-Band_Observable_Response_to_Snow_Parameters2.cdf') xid=ncdf_varid(fid,'dBZe_h') & ncdf_varget,fid,xid, dBZe_h_resp xid=ncdf_varid(fid,'zdr') & ncdf_varget,fid,xid, zdr_resp xid=ncdf_varid(fid,'kdp') & ncdf_varget,fid,xid, kdp_resp xid=ncdf_varid(fid,'rhohv') & ncdf_varget,fid,xid, rhohv_resp xid=ncdf_varid(fid,'precip_rate') & ncdf_varget,fid,xid, precip_rate_resp xid=ncdf_varid(fid,'iwc') & ncdf_varget,fid,xid, iwc_resp xid=ncdf_varid(fid,'alpha') & ncdf_varget,fid,xid, alpha_resp xid=ncdf_varid(fid,'bm') & ncdf_varget,fid,xid, bm_resp xid=ncdf_varid(fid,'aspect_ratio') & ncdf_varget,fid,xid, aspect_ratio_resp xid=ncdf_varid(fid,'Dm') & ncdf_varget,fid,xid, Dm_resp xid=ncdf_varid(fid,'temperature') & ncdf_varget,fid,xid, temperature_resp xid=ncdf_varid(fid,'sigma') & ncdf_varget,fid,xid, sigma_resp ncdf_close,fid ; get snow stats. Presently reads in data from Socrates and Olympex. ; am/bm are assumed as aggregates, rimed, and graupel. ; populate the in_situ_snow_properties common snow_path='' ;snow_path='/Users/u0029340/Documents/data/Capricorn/misc_graphics/precip_properties/' get_insitu_snow_properties,snow_path ; set these zh_unc=1. ;db zdr_unc=0.1 ; db canting_angle_assumed=0. freq_ghz=6. alpha_guess=2. r_guess=0.6 bm_guess=2.4 am_regress=10.^((bm_guess-3.7038)/0.75) temperature=260. am_uncert=0.1 alpha_uncert=1. canting_sdv=10. iwc_ret_xy=make_array(cropDims[0],cropDims[1],/float); iwc retrieved iwc_ret_unc_xy=make_array(cropDims[0],cropDims[1],/float) ; iwc retrieval fractional uncertainty dm_ret_xy=make_array(cropDims[0],cropDims[1],/float) ; mass mean size retrieved dm_ret_unc_xy=make_array(cropDims[0],cropDims[1],/float) ; mass mean size retrieved r_ret_xy=make_array(cropDims[0],cropDims[1],/float) ; aspect ratio retrieved r_ret_unc_xy=make_array(cropDims[0],cropDims[1],/float) ; aspect ratio retrieved bm_ret_xy=make_array(cropDims[0],cropDims[1],/float) ; bm size retrieved bm_ret_unc_xy=make_array(cropDims[0],cropDims[1],/float) ; bm size retrieved dbz_sim=make_array(cropDims[0],cropDims[1],/float) zdr_sim=make_array(cropDims[0],cropDims[1],/float) kdp_sim=make_array(cropDims[0],cropDims[1],/float) for j=0,cropDims[0]-1 do begin for k=0,cropDims[1]-1 do begin if j mod 10 eq 0 then print,j,k,cropImg_dbz[j,k] if cropImg_dbz[j,k] ne 0.0 then begin ; -10., is a temprature estiamte. From Hogan et al. 2006 for a 3 Ghz radar: log10(IWC) 0.060Z 0.0197T 1.70 iwc_guess=10.^((0.060*((cropImg_dbz[j,k])))-(0.0197*(-10.))-1.7) kdp_unc=max([0.25*abs(cropImg_kdp[j,k]),0.025]) ; need to iterate on the mass mean size, dm, for a given dbz: dm=0.005 & dbz_iter=-12. while dbz_iter lt cropImg_dbz[j,k] do begin ;zz=zhv_forward_model_rayleigh_oblate(aspect_ratio_viv,am_regress,2.2,iwc_viv*1.e-6,Dm,alpha_viv,6.,temperature) ; This one errored out because only 8 subscripts allowed. ;zz=zhv_forward_model_rayleigh_oblate_canting(r_guess,am_regress,bm_guess,iwc_guess*1.e-6,Dm,$ ; alpha_guess,canting_angle_assumed,cropImg_vcy[j,k],freq_ghz,temperature) dbz_iter=10.*alog10(zz[0]) dm=dm*1.05 endwhile first_guess_vector=[iwc_guess*1.e-6,dm/1.05,r_guess,bm_guess] ; this runs the inversion invert_zh_zdr_kdp_for_iwc_dm_r_bm_canting_test,cropImg_dbz[j,k],cropImg_zdr[j,k],cropImg_kdp[j,k],$ ret_vector_r,unc_vector_r,first_guess_vector,zh_unc,zdr_unc,kdp_unc,$ alpha_guess,temperature,fx_vector,fx_vector_uncertainty,am_uncert,alpha_uncert,canting_sdv ; store/plot the following into the output variables - here assumed to be in 3d j,k,l iwc_ret_xy[j,k]=ret_vector_r[0] ; iwc retrieved iwc_ret_unc_xy[j,k]=unc_vector_r[0] ; iwc retrieval fractional uncertainty dm_ret_xy[j,k]=ret_vector_r[1] ; mass mean size retrieved dm_ret_unc_xy[j,k]=unc_vector_r[1] ; mass mean size retrieved r_ret_xy[j,k]=ret_vector_r[2] ; aspect ratio retrieved r_ret_unc_xy[j,k]=unc_vector_r[2] ; aspect ratio retrieved bm_ret_xy[j,k]=ret_vector_r[3] ; bm size retrieved bm_ret_unc_xy[j,k]=ret_vector_r[3] ; bm size retrieved ; plot the difference between the simulated observations and the actual observations dbz_sim[j,k]=fx_vector[0] zdr_sim[j,k]=fx_vector[1] kdp_sim[j,k]=fx_vector[2] endif endfor endfor endif ;turn off retrieval ;********* ; Write out the cband data for the feature to a file ;********* fstr=strmid(fname,0,30) istr=string(i,format='(I02)') ; Write to a netcdf file if 1 eq 1 then begin output_filename=fstr+'_roi.'+istr+'.data.cdf' fid=ncdf_create(output_filename,/clobber) print, output_filename,fid,'start write' x_did=ncdf_dimdef(fid,'x',cropDims[0]) y_did=ncdf_dimdef(fid,'y',cropDims[1]) z_did=ncdf_dimdef(fid,'z',numz) dbz_id=ncdf_vardef(fid,'cropImg_dbz',[x_did,y_did],/float) ncdf_attput,fid,dbz_id,'long_name','1km high plan observations dbz' ncdf_attput,fid,dbz_id,'comment','1km high plan used to identify feature' zdr_id=ncdf_vardef(fid,'cropImg_zdr',[x_did,y_did],/float) ncdf_attput,fid,zdr_id,'long_name','1km high plan observations zdr' kdp_id=ncdf_vardef(fid,'cropImg_kdp',[x_did,y_did],/float) ncdf_attput,fid,kdp_id,'long_name','1km high plan observations kdp' vcy_id=ncdf_vardef(fid,'cropImg_vcy',[x_did,y_did],/float) ncdf_attput,fid,vcy_id,'long_name','1km high plan observations velocity' rho_id=ncdf_vardef(fid,'cropImg_rho',[x_did,y_did],/float) ncdf_attput,fid,rho_id,'long_name','1km high plan observations rhohv' dbzv_id=ncdf_vardef(fid,'volImg_dbz',[x_did,y_did,z_did],/float) ncdf_attput,fid,dbzv_id,'long_name','volume observations dbz' zdrv_id=ncdf_vardef(fid,'volImg_zdr',[x_did,y_did,z_did],/float) ncdf_attput,fid,zdrv_id,'long_name','volume observations zdr' kdpv_id=ncdf_vardef(fid,'volImg_kdp',[x_did,y_did,z_did],/float) ncdf_attput,fid,kdpv_id,'long_name','volume observations kdp' vcyv_id=ncdf_vardef(fid,'volImg_vcy',[x_did,y_did,z_did],/float) ncdf_attput,fid,vcyv_id,'long_name','volume observations velocity' rhov_id=ncdf_vardef(fid,'volImg_rho',[x_did,y_did,z_did],/float) ncdf_attput,fid,rhov_id,'long_name','volume observations rhohv' if do_ret eq 'yes' then begin dbz_sim_id=ncdf_vardef(fid,'dbz_sim',[x_did,y_did],/float) ncdf_attput,fid,dbz_sim_id,'long_name','simulated dbz' zdr_sim_id=ncdf_vardef(fid,'zdr_sim',[x_did,y_did],/float) ncdf_attput,fid,zdr_sim_id,'long_name','simulated zdr' kdp_sim_id=ncdf_vardef(fid,'kdp_sim',[x_did,y_did],/float) ncdf_attput,fid,kdp_sim_id,'long_name','simulated kdp' iwc_ret_id=ncdf_vardef(fid,'iwc_ret',[x_did,y_did],/float) ncdf_attput,fid,iwc_ret_id,'long_name','iwc retrieved' iwc_ret_unc_id=ncdf_vardef(fid,'iwc_ret_unc',[x_did,y_did],/float) ncdf_attput,fid,iwc_ret_unc_id,'long_name','iwc retrieval fractional uncertainty' dm_ret_id=ncdf_vardef(fid,'dm_ret',[x_did,y_did],/float) ncdf_attput,fid,dm_ret_id,'long_name','mass mean size retrieved' dm_ret_unc_id=ncdf_vardef(fid,'dm_ret_unc',[x_did,y_did],/float) ncdf_attput,fid,dm_ret_unc_id,'long_name','mass mean size retrieved uncertainty' r_ret_id=ncdf_vardef(fid,'r_ret',[x_did,y_did],/float) ncdf_attput,fid,r_ret_id,'long_name','aspect ratio retrieved' r_ret_unc_id=ncdf_vardef(fid,'r_ret_unc',[x_did,y_did],/float) ncdf_attput,fid,r_ret_unc_id,'long_name','aspect ratio retrieved uncertainty' bm_ret_id=ncdf_vardef(fid,'bm_ret',[x_did,y_did],/float) ncdf_attput,fid,bm_ret_id,'long_name','bm size retrieved' bm_ret_unc_id=ncdf_vardef(fid,'bm_ret_unc',[x_did,y_did],/float) ncdf_attput,fid,bm_ret_unc_id,'long_name','bm size retrieved uncertainty' endif ncdf_control,fid,/endef ncdf_varput,fid,dbz_id,cropImg_dbz ncdf_varput,fid,zdr_id,cropImg_zdr ncdf_varput,fid,kdp_id,cropImg_kdp ncdf_varput,fid,vcy_id,cropImg_vcy ncdf_varput,fid,rho_id,cropImg_rho ncdf_varput,fid,dbzv_id,volImg_dbz ncdf_varput,fid,zdrv_id,volImg_zdr ncdf_varput,fid,kdpv_id,volImg_kdp ncdf_varput,fid,vcyv_id,volImg_vcy ncdf_varput,fid,rhov_id,volImg_rho if do_ret eq 'yes' then begin ncdf_varput,fid,dbz_sim_id,dbz_sim ncdf_varput,fid,zdr_sim_id,zdr_sim ncdf_varput,fid,kdp_sim_id,kdp_sim ncdf_varput,fid,iwc_ret_id,iwc_ret_xy ncdf_varput,fid,iwc_ret_unc_id,iwc_ret_unc_xy ncdf_varput,fid,dm_ret_id,dm_ret_xy ncdf_varput,fid,dm_ret_unc_id,dm_ret_unc_xy ncdf_varput,fid,r_ret_id,r_ret_xy ncdf_varput,fid,r_ret_unc_id,r_ret_unc_xy ncdf_varput,fid,bm_ret_id,bm_ret_xy ncdf_varput,fid,bm_ret_unc_id,bm_ret_unc_xy endif ncdf_close,fid print,'end write' endif ; turn of write ;################################## ; Work with ellipse ;################################## if do_ellipse eq 'yes' then begin ; Get the 3d volume of the cropimage volImg=dbz[x1:x2,y1:y2,*] ; Fit an ellipse to the cropimage mask_crop=mask[x1:x2,y1:y2] indices=where(mask_crop gt 0) s=size(mask_crop,/dimensions) ellipsePts=Fit_Ellipse(indices,xsize=s[0],ysize=s[1],$ center=center,orientation=orientation,axes=axes,semiaxes=semiaxes) ellipse_x=ellipsePts[0,*]+x1 ellipse_y=ellipsePts[1,*]+y1 ellipse_x_center=center[0]+x1 ellipse_y_center=center[1]+y1 semimajor=semiaxes[0] semiminor=semiaxes[1] phi=0 xcm=center[0]+x1 ycm=center[1]+y1 ; Position angle in radians. ang = orientation / !RADEG ; Sin and cos of angle. cosang = Cos(ang) sinang = Sin(ang) ; Parameterized equation of ellipse. xx = semimajor * Cos(phi) yy = semiminor * Sin(phi) ; Rotate to desired position angle. xprime = xcm + (xx * cosang) - (yy * sinang) yprime = ycm + (xx * sinang) + (yy * cosang) print, xprime,yprime endif ;end of do ellipse ;********************** ; Start the image ;********************** ; Define the image 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]) ; ***Plot dbz image pnum=0 data_var=img_dbz r=where(data_var ne -9999) dbz_min=min(data_var[r]) dbz_max=max(data_var[r]) data_image=bytscl(data_var,min=dbz_min,max=dbz_max,top=top_color) r=where(data_var eq -9999,c) if c gt 0 then data_image[r]=253 ;put -9999 to white p1=image(data_image,x,y,/current,position=pos[pnum,*],$ rgb_table=mytable,image_dimensions=[isx,isy],irregular=0) p2=contour(data_var,x,y,/current,position=pos[pnum,*],$ axis_style=2,font_size=fs1,$ xstyle=1,ystyle=1,xtickdir=1,ytickdir=1,/nodata) t1=text(pos[pnum,0],pos[pnum,3]+2*dy,'Reflectivity (dBZ)',font_size=fs2) c0=colorbar(range=[dbz_min,dbz_max],$ orientation=1,position=reform(cbpos[pnum,*]),font_size=fs1,rgb_table=mycbtable,$ tickdir=1,border_on=1,title='dBZ') ; Overplot the ROI outline p3=plot(x_roi,y_roi,/overplot,color='black',$ xrange=[x[0],x[-1]],yrange=[y[0],y[-1]],thick=2) ; ; ***Plot mask from compute mask ; pnum=1 ; data_var=maskResult ; r=where(data_var ne -9999) ; dmin=min(data_var[r]) ; dmax=max(data_var[r]) ; data_image=bytscl(data_var,min=dmin,max=dmax,top=top_color) ; r=where(data_var eq -9999,c) ; if c gt 0 then data_image[r]=253 ; p1=image(data_image,x,y,/current,position=pos[pnum,*],$ ; rgb_table=mytable,image_dimensions=[isx,isy],irregular=0) ; p2=contour(data_var,x,y,/current,position=pos[pnum,*],$ ; axis_style=2,font_size=fs1,$ ; xstyle=1,ystyle=1,xtickdir=1,ytickdir=1,/nodata) ; t1=text(pos[pnum,0],pos[pnum,3]+2*dy,'mask result from ComputeMask',font_size=fs2) ; c0=colorbar(range=[dmin,dmax],$ ; orientation=1,position=reform(cbpos[pnum,*]),font_size=fs1,rgb_table=mycbtable,$ ; tickdir=1,border_on=1,title='') ; *** Plot binary mask of 0 and 1 pnum=1 data_var=mask r=where(data_var ne -9999) dmin=min(data_var[r]) dmax=max(data_var[r]) data_image=bytscl(data_var,min=dmin,max=dmax,top=top_color) r=where(data_var eq -9999,c) if c gt 0 then data_image[r]=253 p4=image(data_image,x,y,/current,position=pos[pnum,*],$ rgb_table=mytable,image_dimensions=[isx,isy],irregular=0) p5=contour(data_var,x,y,/current,position=pos[pnum,*],$ axis_style=2,font_size=fs1,$ xstyle=1,ystyle=1,xtickdir=1,ytickdir=1,/nodata) t1=text(pos[pnum,0],pos[pnum,3]+2*dy,'mask of binary 0 and 1',font_size=fs2) c0=colorbar(range=[dmin,dmax],$ ;the colorbar is to match the image labeled p1 orientation=1,position=reform(cbpos[pnum,*]),font_size=fs1,rgb_table=mycbtable,$ tickdir=1,border_on=1,title='') ; *** Plot mask*image pnum=3 data_var=maskImg_dbz data_image=bytscl(data_var,min=dbz_min,max=dbz_max,top=top_color) r=where(data_var eq 0,c) if c gt 0 then data_image[r]=253 p6=image(data_image,x,y,/current,position=pos[pnum,*],$ rgb_table=mytable,image_dimensions=[isx,isy],irregular=0) p7=contour(data_var,x,y,/current,position=pos[pnum,*],$ axis_style=2,font_size=fs1,$ xstyle=1,ystyle=1,xtickdir=1,ytickdir=1,/nodata) t1=text(pos[pnum,0],pos[pnum,3]+2*dy,'mask*image',font_size=fs2) c0=colorbar(range=[dbz_min,dbz_max],$ orientation=1,position=reform(cbpos[pnum,*]),font_size=fs1,rgb_table=mycbtable,$ tickdir=1,border_on=1,title='dBZ') ; *** Plot crop image pnum=4 data_var=cropImg_dbz data_image=bytscl(data_var,min=dbz_min,max=dbz_max,top=top_color) r=where(data_var eq 0,c) if c gt 0 then data_image[r]=253 p8=image(data_image,cropx,cropy,/current,position=pos[pnum,*],$ rgb_table=mytable,image_dimensions=[isx,isy],irregular=0) p9=contour(data_var,cropx,cropy,/current,position=pos[pnum,*],$ axis_style=2,font_size=fs1,yrange=[cropy[0],cropy[-1]],$ xrange=[cropx[0],cropx[-1]],$ xstyle=1,ystyle=1,xtickdir=1,ytickdir=1,/nodata) t1=text(pos[pnum,0],pos[pnum,3]+2*dy,'crop image',font_size=fs2) c0=colorbar(range=[dbz_min,dbz_max],$ orientation=1,position=reform(cbpos[pnum,*]),font_size=fs1,rgb_table=mycbtable,$ tickdir=1,border_on=1,title='dBZ') ; *** Plot the image slice if 1 eq 1 then begin pnum=5 if do_ellipse eq 'yes' then begin ; Ellipse major axis x_a=ellipse_x_center x_b=xprime y_a=ellipse_y_center y_b=yprime endif else begin ; Diagonal line ;x_a=min(cropx) ;x_b=max(cropx) ;y_a=min(cropy) ;y_b=max(cropy) ; Line segment x_a=290 x_b=260 y_a=min(cropy) y_b=150 endelse slope_m=(y_a-y_b)/(x_a-x_b) intercept=y_a-(slope_m*x_a) num_hires=100.0 dxx=(cropx[-1]-cropx[0])/num_hires hiresx=(findgen(num_hires)*dxx)+cropx[0] ;y=(slope_m*x)+intercept ;yy=(slope_m*cropx)+intercept yy=(slope_m*hiresx)+intercept ;r=where(cropx ge x_a and cropx le x_b and $ ; yy ge y_a and yy le y_b,c) ;yy=yy[r] ;lat ;xx=cropx[r] ;r=where(yy ge cropy[0] and yy le cropy[-1]) r=where(yy ge y_a and yy le y_b) ;line between my x,y endpoints yy=yy[r] ;xx=cropx[r] xx=hiresx[r] xsec_dbz=make_array(n_elements(xx),numz) xsec_zdr=make_array(n_elements(xx),numz) xsec_kdp=make_array(n_elements(xx),numz) for ii=0,n_elements(yy)-1 do begin closest=min(abs(yy[ii]-cropy),idx_lat) closest=min(abs(xx[ii]-cropx),idx_lon) xsec_dbz[ii,*]=volImg_dbz[idx_lon,idx_lat,*] xsec_zdr[ii,*]=volImg_zdr[idx_lon,idx_lat,*] xsec_kdp[ii,*]=volImg_kdp[idx_lon,idx_lat,*] endfor ; Choose height range to plot data_var=xsec_dbz hidx=where(z_km lt 10) data_image=bytscl(data_var,top=top_color,min=dbz_min,max=dbz_max) r=where(data_var eq -9999,c) if c gt 0 then data_image[r]=253 r=where(data_var eq 0,c) if c gt 0 then data_image[r]=255 p10=image(data_image[*,hidx],/current,image_dimensions=[isx,isy],$ position=pos[pnum,*],rgb_table=mytable) p11=contour(data_var[*,hidx],xx,z_km[hidx],$ ytitle='height (km)',xtitle='X',irregular=0,xtickdir=1,ytickdir=1,$ xstyle=1,ystyle=1,/nodata,/current,position=pos[pnum,*],font_size=fs1,$ axis_style=2,$ yrange=[z_km[hidx[0]],z_km[hidx[-1]]]);,xrange=[lat_var[0,0],lat_var[-1,0]]) ;xaxis=axis('X',location=-2,title='X',target=p11,$ ; coord_transform=[intercept,slope_m],tickfont_size=fs1) p12=plot(xx,make_array(n_elements(xx),value=z_km[sidx]),$ /overplot,color='dim gray',thick=3,linestyle=4) ; Plot the diagonal line on the crop plane image ;p13=plot(xx,yy,overplot=p8,color='gray',thick=3,linestyle=4) p13=plot(xx,yy,/current,position=pos[4,*],$ font_size=fs1,yrange=[cropy[0],cropy[-1]],xrange=[cropx[0],cropx[-1]],$ xstyle=1,ystyle=1,xtickdir=1,ytickdir=1,color='dim gray',thick=3,linestyle=4) ;t1=text(pos[pnum,0]+dx,pos[pnum,3]+2*dy,'Constant Latitude '+string(xlat,format='(F6.1)'),font_size=fs1) t1=text(pos[pnum,0]+dx,pos[pnum,3]+2*dy,'cross-section',font_size=fs1) ; Colorbar cb0=colorbar(rgb_table=mycbtable,range=[dbz_min,dbz_max],$ orientation=1,position=reform(cbpos[pnum,*]),font_size=fs1,$ tickdir=1,border_on=1,title='dBZ') endif ;*** Plot the ellipse if do_ellipse eq 'yes' then begin pnum=2 ;p14=plot(ellipsepts,/current,position=pos[pnum,*]) p14=plot(ellipse_x,ellipse_y,/current,position=pos[pnum,*],$ yrange=[cropy[0],cropy[-1]],xrange=[cropx[0],cropx[-1]]) p15=symbol(ellipse_x_center,ellipse_y_center,symbol='o',sym_color='red',/data) p15=symbol(xprime,yprime,symbol='o',sym_color='blue',/data) p16=plot(xx,yy,/overplot,color='dim gray',thick=3,linestyle=4) endif ; ; Print the statistics of each ROI when it is ; pnum=2 ; t2=text(pos[pnum,0],pos[pnum,3]-6*dy,'ROI mask area = '+strcompress(maskArea,/remove_all)+' pixels') ; t2=text(pos[pnum,0],pos[pnum,3]-8*dy,'ROI mask mean = '+strcompress(maskMean,/remove_all)+' dBZ') ; t2=text(pos[pnum,0],pos[pnum,3]-10*dy,'ROI geometric area = '+strcompress(geomArea,/remove_all)+' m') ; t2=text(pos[pnum,0],pos[pnum,3]-12*dy,'ROI perimeter = '+strcompress(perimeter,/remove_all)+' m') t1=text(xl,0.97,fname,font_size=fs1) ; Image name fstr=strmid(fname,0,30) istr=string(i,format='(I02)') p0.save,fstr+'_roi.'+istr+'.line.png',height=pydim ;**************************** ; Plot a retrieval image ;**************************** ; *** Plot crop image pnum=0 data_var=cropImg_dbz data_image=bytscl(data_var,min=dbz_min,max=dbz_max,top=top_color) r=where(data_var eq 0,c) if c gt 0 then data_image[r]=253 p8=image(data_image,cropx,cropy,/buffer,position=pos[pnum,*],$ rgb_table=mytable,image_dimensions=[isx,isy],irregular=0,dimensions=[pxdim,pydim]) p9=contour(data_var,cropx,cropy,/current,position=pos[pnum,*],$ axis_style=2,font_size=fs1,yrange=[cropy[0],cropy[-1]],$ xrange=[cropx[0],cropx[-1]],$ xstyle=1,ystyle=1,xtickdir=1,ytickdir=1,/nodata) t1=text(pos[pnum,0],pos[pnum,3]+2*dy,'DBZ',font_size=fs2) c0=colorbar(range=[dbz_min,dbz_max],$ orientation=1,position=reform(cbpos[pnum,*]),font_size=fs1,rgb_table=mycbtable,$ tickdir=1,border_on=1,title='dBZ') pnum=1 data_var=cropImg_zdr dmin=-1.0;min(cropImg_zdr) dmax=1.0;max(cropImg_zdr) data_image=bytscl(data_var,min=dmin,max=dmax,top=top_color) r=where(data_var eq 0,c) if c gt 0 then data_image[r]=253 r=where(finite(data_var,/NAN),c) if c gt 0 then data_image[r]=253 p8=image(data_image,cropx,cropy,/current,position=pos[pnum,*],$ rgb_table=mytable2,image_dimensions=[isx,isy],irregular=0) p9=contour(data_var,cropx,cropy,/current,position=pos[pnum,*],$ axis_style=2,font_size=fs1,yrange=[cropy[0],cropy[-1]],$ xrange=[cropx[0],cropx[-1]],$ xstyle=1,ystyle=1,xtickdir=1,ytickdir=1,/nodata) t1=text(pos[pnum,0],pos[pnum,3]+2*dy,'ZDR',font_size=fs2) c0=colorbar(range=[dmin,dmax],$ orientation=1,position=reform(cbpos[pnum,*]),font_size=fs1,rgb_table=mycbtable2,$ tickdir=1,border_on=1,title='ZDR') pnum=2 data_var=cropImg_kdp data_rho=cropImg_rho ;dmin=min(cropImg_kdp) ;dmax=max(cropImg_kdp) ;data_image=bytscl(data_var,min=dmin,max=dmax,top=top_color) ;r=where(data_var eq 0,c) ;if c gt 0 then data_image[r]=253 kdp_int=[-100,-0.1,-0.05,0,0.05,0.1,1.0,100.0] lrho=0.95 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 'yes' 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 p8=image(data_image,cropx,cropy,/current,position=pos[pnum,*],$ rgb_table=mytable3,image_dimensions=[isx,isy],irregular=0) p9=contour(data_var,cropx,cropy,/current,position=pos[pnum,*],$ axis_style=2,font_size=fs1,yrange=[cropy[0],cropy[-1]],$ xrange=[cropx[0],cropx[-1]],$ xstyle=1,ystyle=1,xtickdir=1,ytickdir=1,/nodata) t1=text(pos[pnum,0],pos[pnum,3]+2*dy,'KDP',font_size=fs2) ;c0=colorbar(range=[dmin,dmax],$ ; orientation=1,position=reform(cbpos[pnum,*]),font_size=fs1,rgb_table=mycbtable3,$ ; tickdir=1,border_on=1,title='KDP') if rho_filter eq 'yes' 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],$ position=reform(cbpos[pnum,*]),$ font_size=10,$ tickdir=1,border_on=1);,title='dBZ') dmin_rho=0.0;min(rho[r]) dmax_rho=1.0;max(rho[r]) ; 0.5 km pnum=3 data_var=cropImg_rho 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 p8=image(data_image,cropx,cropy,/current,position=pos[pnum,*],$ rgb_table=mytable,image_dimensions=[isx,isy],irregular=0) p9=contour(data_var,cropx,cropy,/current,position=pos[pnum,*],$ axis_style=2,font_size=fs1,yrange=[cropy[0],cropy[-1]],$ xrange=[cropx[0],cropx[-1]],$ xstyle=1,ystyle=1,xtickdir=1,ytickdir=1,/nodata) t1=text(pos[pnum,0],pos[pnum,3]+2*dy,'RHO HV',font_size=fs2) if do_ret eq 'yes' then begin pnum=3 data_var=iwc_ret_xy dmin=min(data_var) dmax=max(data_var) data_image=bytscl(data_var,min=dmin,max=dmax,top=top_color) r=where(data_var eq 0,c) if c gt 0 then data_image[r]=253 p8=image(data_image,cropx,cropy,/current,position=pos[pnum,*],$ rgb_table=mytable,image_dimensions=[isx,isy],irregular=0) p9=contour(data_var,cropx,cropy,/current,position=pos[pnum,*],$ axis_style=2,font_size=fs1,yrange=[cropy[0],cropy[-1]],$ xrange=[cropx[0],cropx[-1]],$ xstyle=1,ystyle=1,xtickdir=1,ytickdir=1,/nodata) t1=text(pos[pnum,0],pos[pnum,3]+2*dy,'IWC Retrieval',font_size=fs2) c0=colorbar(range=[dmin,dmax],$ orientation=1,position=reform(cbpos[pnum,*]),font_size=fs1,rgb_table=mycbtable,$ tickdir=1,border_on=1,title='') endif t1=text(xl,0.97,fname,font_size=fs1) p8.save,fstr+'_roi.'+istr+'.retrieval.png',height=pydim ; Remove each unneeded object reference after ; displaying it. OBJ_DESTROY, oROI print,'next plot' stop ; End the FOR loop. ENDFOR stop ;; Initialize the IDLgrROI object with the contour ;; information of the current region with the FOR ;; loop. ;I=0 ;line = [LINDGEN(pathInfo(I).N), 0] ;;oROI = OBJ_NEW('IDLanROI', $ ;oROI = OBJ_NEW('IDLgrROI', $ ; (pathXY(*,pathInfo(I).OFFSET + line))[0, *], $ ; (pathXY(*,pathInfo(I).OFFSET + line))[1, *]) ;;assign the roi data to the arrays x,y ;oROI -> GetProperty, DATA = ROIdata ;x = ROIdata[0,*] ;y = ROIdata[1,*] ;;Initialize an IDLgrImage object containing the original image data: ;oImg = OBJ_NEW('IDLgrImage', img,DIMENSIONS = dims) ;;Create a window in which to display the image and the ROI: ;oWindow = OBJ_NEW('IDLgrWindow', DIMENSIONS = dims, RETAIN = 2, TITLE = 'Selected ROI') ;;Create the view plane and initialize the view: ;viewRect = [0, 0, dims[0], dims[1]] ;oView = OBJ_NEW('IDLgrView', VIEWPLANE_RECT = viewRect) ;;Initialize a model object and add the image and ROI to the model. ;;Add the model to the view and draw the view in the window to display ;;the ROI overlaid onto the original image: ;oModel = OBJ_NEW('IDLgrModel') ;oModel -> Add, oImg ;oModel -> Add, oROI ;ROIout ;oView -> Add, oModel ;oWindow -> Draw, oView ; ;;Use the IDLanROI::ComputeMask function to create a 2D mask of the region. ;;Pixels that fall outside of the ROI will be assigned a value of 0: ;;maskResult = ROIout -> ComputeMask(DIMENSIONS = dims) ;maskResult = oROI -> ComputeMask(DIMENSIONS = dims) ; ;;Use the IMAGE_STATISTICS procedure to compute the area of the mask, ;;inputting maskResult as the MASK value. Print count to view the ;;number of pixels occurring within the masked region: ;IMAGE_STATISTICS, img, MASK = MaskResult, COUNT = count ;PRINT, 'area of mask = ', count,' pixels' ; ;;From the ROI mask, create a binary mask, consisting of only zeros and ones. ;;Multiply the binary mask times the original image to retain only the portion ;;of the image that was defined in the original ROI: ;mask = (maskResult GT 0) ;maskImg = img * mask ;;Using the minimum and maximum values of the ROI array, create a cropped array, ;;cropImg, and get its dimensions: ;cropImg = maskImg[min(x):max(x), min(y): max(y)] ;cropDims = SIZE(cropImg, /DIMENSIONS) ; ;;Initialize an image object with the cropped region data: ;oMaskImg = OBJ_NEW('IDLgrImage', cropImg, DIMENSIONS = dims) ;;Using the cropped region dimensions, create an offset window. ;;Multiply the x and y dimensions times the value by which you wish to magnify the ROI: ;oMaskWindow = OBJ_NEW('IDLgrWindow', DIMENSIONS = 2 * cropDims, RETAIN = 2, $ ; TITLE = 'Magnified ROI', LOCATION = dims) ;;Create the display objects and display the cropped and magnified ROI: ;oMaskView = OBJ_NEW('IDLgrView', VIEWPLANE_RECT = viewRect) ;oMaskModel = OBJ_NEW('IDLgrModel') ;oMaskModel -> Add, oMaskImg ;oMaskView -> Add, oMaskModel ;OMaskWindow -> Draw, oMaskView ;;Clean up object references that are not destroyed by the window manager ;;when you close the Object Graphics displays: ;;OBJ_DESTROY, [oView, oMaskView, ROIout] ;OBJ_DESTROY, [oView, oMaskView, oROI] END