;********************************* ; 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_xsec ;******************************* ; Read in the data ;******************************* fdir='/Volumes/mace-group4/Capricorn2/CBand/cfradial/' fname='9776HUB-PPIVol-20180209-140504-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 xid=ncdf_varid(fid,'corrected_reflectivity') & ncdf_varget,fid,xid,dbz 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 sidx=where(z_km ge 2.0) sidx=sidx[0] ; horizontal slice at sidx height level img=dbz[*,*,sidx] dims = SIZE(img, /DIMENSIONS) ;************************************************** ; Create a mask that identifies the pixels of interest ; 1=pixel gt -15 dbz, 0=pixel le -15 dbz ;************************************************** threshImg=(img 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 ;************************************************** ; Set up the plot ;************************************************** ; 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,*] ; Set up the positions pxdim=1000 & pydim=800 xl=0.06 & xr=0.90 yb=0.05 & yt=0.95 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, 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 = 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)] 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 = maskImg[x1:x2,y1:y2] cropx=x[x1:x2] cropy=y[y1:y2] cropDims = SIZE(cropImg, /DIMENSIONS) ; 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 ; 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 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 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 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 ; Diagonal line ;x_a=min(cropx) ;x_b=max(cropx) ;y_a=min(cropy) ;y_b=max(cropy) ; Ellipse major axis x_a=ellipse_x_center x_b=xprime y_a=ellipse_y_center y_b=yprime 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]) yy=yy[r] ;xx=cropx[r] xx=hiresx[r] data_var=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) data_var[ii,*]=volImg[idx_lon,idx_lat,*] endfor ; Choose height range to plot 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 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) ; ; 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') ; Image name fstr=strmid(fname,0,25) istr=string(i,format='(I02)') p0.save,fstr+'_roi.'+istr+'.png',height=pydim ; Remove each unneeded object reference after ; displaying it. OBJ_DESTROY, oROI print,'next plot' ; 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