;****************************************************
;  Reads in a feature volume and plots cross-sections
;****************************************************
PRO ProgramDefineROI_plot_ret

;***  Read in the data

;  Read in the feature volume cut from the radar pyart interpolated file
fname='9776HUB-PPIVol-20180209-105330_roi.00.data.cdf'
fid=ncdf_open(fname)
vid=ncdf_varid(fid,'cropImg_dbz') & ncdf_varget,fid,vid,img_dbz
vid=ncdf_varid(fid,'cropImg_zdr') & ncdf_varget,fid,vid,img_zdr
vid=ncdf_varid(fid,'cropImg_kdp') & ncdf_varget,fid,vid,img_kdp
vid=ncdf_varid(fid,'cropImg_rho') & ncdf_varget,fid,vid,img_rho
vid=ncdf_varid(fid,'volImg_dbz') & ncdf_varget,fid,vid,vol_dbz
vid=ncdf_varid(fid,'volImg_zdr') & ncdf_varget,fid,vid,vol_zdr
vid=ncdf_varid(fid,'volImg_kdp') & ncdf_varget,fid,vid,vol_kdp
vid=ncdf_varid(fid,'volImg_rho') & ncdf_varget,fid,vid,vol_rho
ncdf_close,fid
;  Strings from the filename
fstr=strmid(fname,0,37)
parts=strsplit(fname,'-',/extract)
date_str=parts[2]
time_str=strmid(parts[3],0,6)

;  Read in the radar retrieval file
fname2='/Volumes/u0079358/cband/Citation_Inverse_CBand_Count_File_1ht2.cdf'
fid=ncdf_open(fname2)
vid=ncdf_varid(fid,'dbz_bins') & ncdf_varget,fid,vid,dbz_bins
vid=ncdf_varid(fid,'kdp_bins') & ncdf_varget,fid,vid,kdp_bins
vid=ncdf_varid(fid,'zdr_bins') & ncdf_varget,fid,vid,zdr_bins
vid=ncdf_varid(fid,'ice_water_content') & ncdf_varget,fid,vid,iwc_ret
vid=ncdf_varid(fid,'ice_water_content_uncertainty') & ncdf_varget,fid,vid,iwc_unc
vid=ncdf_varid(fid,'Dmean') & ncdf_varget,fid,vid,dmean_ret
vid=ncdf_varid(fid,'Dmean_uncertainty') & ncdf_varget,fid,vid,dmean_unc
vid=ncdf_varid(fid,'Aspect_Ratio') & ncdf_varget,fid,vid,ar_ret
vid=ncdf_varid(fid,'Aspect_Ratio_uncertainty') & ncdf_varget,fid,vid,ar_unc
vid=ncdf_varid(fid,'bm') & ncdf_varget,fid,vid,bm_ret
vid=ncdf_varid(fid,'bm_uncertainty') & ncdf_varget,fid,vid,bm_unc
vid=ncdf_varid(fid,'dbz_w') & ncdf_varget,fid,vid,dbz_w_ret
vid=ncdf_varid(fid,'dbz_ka') & ncdf_varget,fid,vid,dbz_ka_ret
vid=ncdf_varid(fid,'dbz_ku') & ncdf_varget,fid,vid,dbz_ku_ret
vid=ncdf_varid(fid,'precip_rate') & ncdf_varget,fid,vid,precip_rate_ret
ncdf_close,fid

;*** Size of the radar volume arrays
s=size(vol_dbz,/dimensions)
numx=s[0]
numy=s[1]
numz=s[2]
x_arr=indgen(s[0])
y_arr=indgen(s[1])
z_arr=indgen(s[2])

;*** Vertical cross section
;  I give each line segment a letter
;  Define the 2 endpoints of the line OR an array of points to polyfit to
;  top_pix is the vertical height of the slice.  
;  Line segment
seg_str='c'
if date_str eq '20180209' and time_str eq '105330' and seg_str eq 'a' then begin
  x_a=60.0 & y_a=0.0
  x_b=30.0 & y_b=40.0
  top_pix=25
endif else if date_str eq '20180209' and time_str eq '105330' and seg_str eq 'b' then begin
  x_a=3.0 & y_a=98.0  
  x_b=25.0 & y_b=81.0
  top_pix=15
endif else if date_str eq '20180209' and time_str eq '105330' and seg_str eq 'c' then begin
  x_a=3.0 & y_a=99.0
  x_b=45.0 & y_b=54.0
  top_pix=20
endif else if date_str eq '20180209' and time_str eq '105330' and seg_str eq 'd' then begin
  x_a=[  4,15,20,28,35,44]
  y_a=[ 90,85,80,70,60,55]
  top_pix=20
endif

if n_elements(x_a) eq 1 then begin
  ;  Line fit between 2 points
  slope_m=(y_a-y_b)/(x_a-x_b)
  intercept=y_a-(slope_m*x_a)

  num_hires=100.0
  dxx=(max([x_a,x_b])-min([x_a,x_b]))/num_hires
  hiresx=(findgen(num_hires)*dxx)+min([x_a,x_b])

  yy=(slope_m*hiresx)+intercept
  r=where(yy ge min([y_a,y_b]) and yy le max([y_a,y_b]))  ;line between my x,y endpoints
  yy=yy[r]
  xx=hiresx[r]
endif else begin
  ;  Polyfit between many points
  coeff=poly_fit(x_a,y_a,2)
  
  num_hires=100.0
  dxx=(max(x_a)-min(x_a))/num_hires
  hiresx=(findgen(num_hires)*dxx)+min(x_a)
  
  yy=coeff[0]+coeff[1]*hiresx+coeff[2]*hiresx^2;+coeff[3]*hiresx^3
  r=where(yy ge min(y_a) and yy le max(y_a))
  yy=yy[r]
  xx=hiresx[r]
endelse

;  Extract the data along the line
dbz=make_array(n_elements(xx),numz)
zdr=make_array(n_elements(xx),numz)
kdp=make_array(n_elements(xx),numz)
rho=make_array(n_elements(xx),numz)
for ii=0,n_elements(yy)-1 do begin
  closest=min(abs(yy[ii]-y_arr),idx_y)
  closest=min(abs(xx[ii]-x_arr),idx_x)
  dbz[ii,*]=vol_dbz[idx_x,idx_y,*]
  zdr[ii,*]=vol_zdr[idx_x,idx_y,*]
  kdp[ii,*]=vol_kdp[idx_x,idx_y,*]
  rho[ii,*]=vol_rho[idx_x,idx_y,*]
endfor

;  Cut height at top_pix
dbz=dbz[*,0:top_pix]
zdr=zdr[*,0:top_pix]
kdp=kdp[*,0:top_pix]
rho=rho[*,0:top_pix]

;  Get the retrieval for the cross section 
;  Size of the cross section
s=size(dbz,/dimensions)
numx=s[0]
numy=s[1]
x_arr=indgen(s[0])
y_arr=indgen(s[1])

;  arrays to hold the retrieval for the cross section
iwc=make_array(numx,numy,/float,value=-9999)
dmean=make_array(numx,numy,/float,value=-9999)
bm=make_array(numx,numy,/float,value=-9999)
aspect_ratio=make_array(numx,numy,/float,value=-9999)

iwc_uncertainty=make_array(numx,numy,/float,value=-9999)
dmean_uncertainty=make_array(numx,numy,/float,value=-9999)
bm_uncertainty=make_array(numx,numy,/float,value=-9999)
aspect_ratio_uncertainty=make_array(numx,numy,/float,value=-9999)

;  loop through each point on the cross section and choose the
;  retrieval based on dbz,zdr,kdp
for i=0,numx-1 do begin
  for j=0,numy-1 do begin  
    if dbz[i,j] ne 0 then begin
      closest1=min(abs(dbz[i,j]-dbz_bins),dbz_idx)
      closest2=min(abs(zdr[i,j]-zdr_bins),zdr_idx)
      closest3=min(abs(kdp[i,j]-kdp_bins),kdp_idx)
      iwc[i,j]=iwc_ret[dbz_idx,zdr_idx,kdp_idx]
      dmean[i,j]=dmean_ret[dbz_idx,zdr_idx,kdp_idx]
      bm[i,j]=bm_ret[dbz_idx,zdr_idx,kdp_idx]
      aspect_ratio[i,j]=ar_ret[dbz_idx,zdr_idx,kdp_idx]
      
      iwc_uncertainty[i,j]=iwc_unc[dbz_idx,zdr_idx,kdp_idx]
      dmean_uncertainty[i,j]=dmean_unc[dbz_idx,zdr_idx,kdp_idx]
      bm_uncertainty[i,j]=bm_unc[dbz_idx,zdr_idx,kdp_idx]
      aspect_ratio_uncertainty[i,j]=ar_unc[dbz_idx,zdr_idx,kdp_idx]     
      ;print,dbz[i,j],dbz_bins[dbz_idx]
      ;print,zdr[i,j],zdr_bins[zdr_idx]
      ;print,kdp[i,j],kdp_bins[kdp_idx]
      ;print,'iwc',iwc[i,j]
      ;stop
    endif
  endfor
endfor

;**************************************************
; 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,*]

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 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=3
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])

;****************************
;  Plot a retrieval image
;****************************

;  *** dbz obs plan
pnum=0
data_var=img_dbz
r=where(data_var ne -9999)
dmin_dbz=-12.;min(data_var[r])
dmax_dbz=37.;max(data_var[r])
data_image=bytscl(data_var,min=dmin_dbz,max=dmax_dbz,top=top_color)
r=where(data_var eq 0,c)
if c gt 0 then data_image[r]=253  ;white
s=size(img_dbz,/dimensions)
x_arr_img=indgen(s[0])
y_arr_img=indgen(s[1])
p8=image(data_image,x_arr_img,y_arr_img,/current,position=pos[pnum,*],$
  rgb_table=mytable,image_dimensions=[isx,isy],irregular=0,dimensions=[pxdim,pydim])
p9=contour(data_var,x_arr_img,y_arr_img,/current,position=pos[pnum,*],$
  axis_style=2,font_size=fs1,yrange=[y_arr_img[0],y_arr_img[-1]],$
  xrange=[x_arr_img[0],x_arr_img[-1]],$
  xstyle=1,ystyle=1,xtickdir=1,ytickdir=1,/nodata)
;  overplot cross section line in black
p13=plot(xx,yy,/overplot,color='black',thick=1,linestyle=0)
t1=text(pos[pnum,0],pos[pnum,3]+2*dy,'DBZ',font_size=fs2)
c0=colorbar(range=[dmin_dbz,dmax_dbz],$
  orientation=1,position=reform(cbpos[pnum,*]),font_size=fs1,rgb_table=mycbtable,$
  tickdir=1,border_on=1,title='dBZ')

;  *** dbz obs xsec
pnum=1
data_var=dbz
r=where(data_var ne -9999)
dmin_dbz=-12.;min(data_var[r])
dmax_dbz=37.;max(data_var[r])
data_image=bytscl(data_var,min=dmin_dbz,max=dmax_dbz,top=top_color)
r=where(data_var eq -9999,c)
if c gt 0 then data_image[r]=253  ;white
p8=image(data_image,x_arr,y_arr,/current,position=pos[pnum,*],$
  rgb_table=mytable,image_dimensions=[isx,isy],irregular=0,dimensions=[pxdim,pydim])
p9=contour(data_var,x_arr,y_arr,/current,position=pos[pnum,*],$
  axis_style=2,font_size=fs1,yrange=[y_arr[0],y_arr[-1]],$
  xrange=[x_arr[0],x_arr[-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=[dmin_dbz,dmax_dbz],$
  orientation=1,position=reform(cbpos[pnum,*]),font_size=fs1,rgb_table=mycbtable,$
  tickdir=1,border_on=1,title='dBZ')

;  *** zdr obs xsec
pnum=2
data_var=zdr
dmin_zdr=-1.0;min(zdr)
dmax_zdr=1.0;max(zdr)
data_image=bytscl(data_var,min=dmin_zdr,max=dmax_zdr,top=top_color)
r=where(data_var eq 0,c)
if c gt 0 then data_image[r]=253  ;white
r=where(finite(data_var,/NAN),c)
if c gt 0 then data_image[r]=253  ;white
p8=image(data_image,x_arr,y_arr,/current,position=pos[pnum,*],$
  rgb_table=mytable2,image_dimensions=[isx,isy],irregular=0)
p9=contour(data_var,x_arr,y_arr,/current,position=pos[pnum,*],$
  axis_style=2,font_size=fs1,yrange=[y_arr[0],y_arr[-1]],$
  xrange=[x_arr[0],x_arr[-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_zdr,dmax_zdr],$
  orientation=1,position=reform(cbpos[pnum,*]),font_size=fs1,rgb_table=mycbtable2,$
  tickdir=1,border_on=1,title='ZDR')

;  *** kdp obs xsec
pnum=3
rho_filter=1
data_var=kdp
data_rho=rho
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 1 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,x_arr,y_arr,/current,position=pos[pnum,*],$
  rgb_table=mytable3,image_dimensions=[isx,isy],irregular=0)
p9=contour(data_var,x_arr,y_arr,/current,position=pos[pnum,*],$
  axis_style=2,font_size=fs1,yrange=[y_arr[0],y_arr[-1]],$
  xrange=[x_arr[0],x_arr[-1]],$
  xstyle=1,ystyle=1,xtickdir=1,ytickdir=1,/nodata)
  t1=text(pos[pnum,0],pos[pnum,3]+2*dy,'KDP',font_size=fs2)
if rho_filter eq 1 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,$
  orientation=1,tickname=kdp_int_str,textpos=0,major=9,$
  position=reform(cbpos[pnum,*]),$
  font_size=10,$
  tickdir=1,border_on=1)

;  *** RHO obs xsec
pnum=4
dmin_rho=0.0;min(rho[r])
dmax_rho=1.0;max(rho[r])
data_var=rho
data_image=bytscl(data_var,top=top_color,min=dmin_rho,max=dmax_rho)
r=where(finite(data_var,/NAN) or data_var eq 0,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,x_arr,y_arr,/current,position=pos[pnum,*],$
  rgb_table=mytable,image_dimensions=[isx,isy],irregular=0)
p9=contour(data_var,x_arr,y_arr,/current,position=pos[pnum,*],$
  axis_style=2,font_size=fs1,yrange=[y_arr[0],y_arr[-1]],$
  xrange=[x_arr[0],x_arr[-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)
c0=colorbar(range=[dmin_rho,dmax_rho],$
  orientation=1,position=reform(cbpos[pnum,*]),font_size=fs1,rgb_table=mycbtable,$
  tickdir=1,border_on=1);,title='dBZ')

;  *** iwc xsec
pnum=5
data_var=iwc
r=where(data_var gt 0,c)
if c gt 0 then begin
  data_var[r]=data_var[r]*1e6
  data_var[r]=alog10(data_var[r])
endif
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 0,c)
if c gt 0 then data_image[r]=253
p8=image(data_image,x_arr,y_arr,/current,position=pos[pnum,*],$
  rgb_table=mytable,image_dimensions=[isx,isy],irregular=0)
p9=contour(data_var,x_arr,y_arr,/current,position=pos[pnum,*],$
  axis_style=2,font_size=fs1,yrange=[y_arr[0],y_arr[-1]],$
  xrange=[x_arr[0],x_arr[-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='log(g/m3)')
  
;  *** dmean xsec  
pnum=6
data_var=dmean
;r=where(data_var gt 0,c)
;if c gt 0 then data_var[r]=alog10(data_var[r])
r=where(data_var ne -9999 and data_var ne 0,c)
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 0,c)
if c gt 0 then data_image[r]=253
p8=image(data_image,x_arr,y_arr,/current,position=pos[pnum,*],$
  rgb_table=mytable,image_dimensions=[isx,isy],irregular=0)
p9=contour(data_var,x_arr,y_arr,/current,position=pos[pnum,*],$
  axis_style=2,font_size=fs1,yrange=[y_arr[0],y_arr[-1]],$
  xrange=[x_arr[0],x_arr[-1]],$
  xstyle=1,ystyle=1,xtickdir=1,ytickdir=1,/nodata)
t1=text(pos[pnum,0],pos[pnum,3]+2*dy,'Dmean 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='cm')
  
;  *** bm xsec
pnum=7
data_var=bm
r=where(data_var gt 0,c)
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 0,c)
if c gt 0 then data_image[r]=253
p8=image(data_image,x_arr,y_arr,/current,position=pos[pnum,*],$
  rgb_table=mytable,image_dimensions=[isx,isy],irregular=0)
p9=contour(data_var,x_arr,y_arr,/current,position=pos[pnum,*],$
  axis_style=2,font_size=fs1,yrange=[y_arr[0],y_arr[-1]],$
  xrange=[x_arr[0],x_arr[-1]],$
  xstyle=1,ystyle=1,xtickdir=1,ytickdir=1,/nodata)
t1=text(pos[pnum,0],pos[pnum,3]+2*dy,'Bm 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='')  
    
;  *** aspect ratio xsec
pnum=8
data_var=aspect_ratio
r=where(data_var gt 0,c)
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 0,c)
if c gt 0 then data_image[r]=253
p8=image(data_image,x_arr,y_arr,/current,position=pos[pnum,*],$
  rgb_table=mytable,image_dimensions=[isx,isy],irregular=0)
p9=contour(data_var,x_arr,y_arr,/current,position=pos[pnum,*],$
  axis_style=2,font_size=fs1,yrange=[y_arr[0],y_arr[-1]],$
  xrange=[x_arr[0],x_arr[-1]],$
  xstyle=1,ystyle=1,xtickdir=1,ytickdir=1,/nodata)
t1=text(pos[pnum,0],pos[pnum,3]+2*dy,'Aspect Ratio 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='')   
t1=text(xl,0.97,fname,font_size=fs1)

;  Save the plot    
p0.save,fstr+'.'+seg_str+'_xsec2.png';,height=pydim

;********************************
;  Start a new plot with more variables
;******************************** 
 
;  *** iwc xsec uncertainty
pnum=0
data_var=iwc_uncertainty
r=where(data_var gt 0,c)
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 0,c)
if c gt 0 then data_image[r]=253
p8=image(data_image,x_arr,y_arr,/buffer,position=pos[pnum,*],dimensions=[pxdim,pydim],$
  rgb_table=mytable,image_dimensions=[isx,isy],irregular=0)
p9=contour(data_var,x_arr,y_arr,/current,position=pos[pnum,*],$
  axis_style=2,font_size=fs1,yrange=[y_arr[0],y_arr[-1]],$
  xrange=[x_arr[0],x_arr[-1]],$
  xstyle=1,ystyle=1,xtickdir=1,ytickdir=1,/nodata)
t1=text(pos[pnum,0],pos[pnum,3]+2*dy,'IWC Uncertainty',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='Fractional')

;  *** dmean uncertainty
pnum=1
data_var=dmean_uncertainty
r=where(data_var ne -9999 and data_var ne 0,c)
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 0,c)
if c gt 0 then data_image[r]=253  ;white
p8=image(data_image,x_arr,y_arr,/current,position=pos[pnum,*],$
  rgb_table=mytable,image_dimensions=[isx,isy],irregular=0)
p9=contour(data_var,x_arr,y_arr,/current,position=pos[pnum,*],$
  axis_style=2,font_size=fs1,yrange=[y_arr[0],y_arr[-1]],$
  xrange=[x_arr[0],x_arr[-1]],$
  xstyle=1,ystyle=1,xtickdir=1,ytickdir=1,/nodata)
t1=text(pos[pnum,0],pos[pnum,3]+2*dy,'Dmean Uncertainty',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='fractional')

;  *** bm xsec
pnum=2
data_var=bm_uncertainty
r=where(data_var gt 0,c)
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 0,c)
if c gt 0 then data_image[r]=253
p8=image(data_image,x_arr,y_arr,/current,position=pos[pnum,*],$
  rgb_table=mytable,image_dimensions=[isx,isy],irregular=0)
p9=contour(data_var,x_arr,y_arr,/current,position=pos[pnum,*],$
  axis_style=2,font_size=fs1,yrange=[y_arr[0],y_arr[-1]],$
  xrange=[x_arr[0],x_arr[-1]],$
  xstyle=1,ystyle=1,xtickdir=1,ytickdir=1,/nodata)
t1=text(pos[pnum,0],pos[pnum,3]+2*dy,'Bm Uncertainty',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='')

;  *** aspect ratio xsec
pnum=3
data_var=aspect_ratio_uncertainty
r=where(data_var gt 0,c)
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 0,c)
if c gt 0 then data_image[r]=253
p8=image(data_image,x_arr,y_arr,/current,position=pos[pnum,*],$
  rgb_table=mytable,image_dimensions=[isx,isy],irregular=0)
p9=contour(data_var,x_arr,y_arr,/current,position=pos[pnum,*],$
  axis_style=2,font_size=fs1,yrange=[y_arr[0],y_arr[-1]],$
  xrange=[x_arr[0],x_arr[-1]],$
  xstyle=1,ystyle=1,xtickdir=1,ytickdir=1,/nodata)
t1=text(pos[pnum,0],pos[pnum,3]+2*dy,'Aspect Ratio Uncertainty',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='')
    
t1=text(xl,0.97,fname,font_size=fs1)

;  Save image
p8.save,fstr+'.'+seg_str+'_unc2.png';,height=pydim


stop
  
 end