pro plot_hid_cfradial

;  Path prefix
path_prefix='/Volumes/'
;path_prefix='/uufs/chpc.utah.edu/common/home/'

;  Hid - fine grid
hdir='mace-group5/gpm/hid/update_20210728/HID_SS_20180209_12-00z/'
hfile='OPOL-SS_rnkHID_F_20180209_140448.lat_lon.nc'

;  Cfradial
cdir='mace-group4/Capricorn2/CBand/cfradial/'
cfile='9776HUB-PPIVol-20180209-140504-0000.cfradial.pyart.81.401.nc'

;  Get the date and time from the filename
parts=strsplit(hfile,'_',/extract)
hdate=parts[3]
htime=parts[4]
hhour=strmid(htime,0,2)
hmin=strmid(htime,2,2)
hsec=strmid(htime,4,2)

;  Read the hid file
fid=ncdf_open(path_prefix+hdir+hfile)
vid=ncdf_varid(fid,'time') & ncdf_varget,fid,vid,time_hid
xid=ncdf_varid(fid,'x') & ncdf_varget,fid,xid,x_m_hid
xid=ncdf_varid(fid,'y') & ncdf_varget,fid,xid,y_m_hid
xid=ncdf_varid(fid,'z') & ncdf_varget,fid,xid,z_m_hid
xid=ncdf_varid(fid,'origin_latitude') & ncdf_varget,fid,xid,olat_hid
xid=ncdf_varid(fid,'origin_longitude') & ncdf_varget,fid,xid,olon_hid
xid=ncdf_varid(fid,'origin_altitude') & ncdf_varget,fid,xid,oalt_hid
xid=ncdf_varid(fid,'point_latitude') & if xid ne -1 then ncdf_varget,fid,xid,lat_3d_hid
xid=ncdf_varid(fid,'point_longitude') & if xid ne -1 then ncdf_varget,fid,xid,lon_3d_hid
xid=ncdf_varid(fid,'point_altitude') & if xid ne -1 then ncdf_varget,fid,xid,alt_3d_hid
xid=ncdf_varid(fid,'corrected_specific_differential_phase') & ncdf_varget,fid,xid,kdp_hid
xid=ncdf_varid(fid,'cross_correlation_ratio') & ncdf_varget,fid,xid,rhohv_hid
xid=ncdf_varid(fid,'log_reflectivity') & ncdf_varget,fid,xid,dbz_hid
xid=ncdf_varid(fid,'log_zdr') & ncdf_varget,fid,xid,zdr_hid
xid=ncdf_varid(fid,'HID_class_01') & ncdf_varget,fid,xid,hid
xid=ncdf_varid(fid,'hid_score_DZ') & ncdf_varget,fid,xid,hid_dz ;drizzle
xid=ncdf_varid(fid,'hid_score_RN') & ncdf_varget,fid,xid,hid_rn ;rain
xid=ncdf_varid(fid,'hid_score_CR') & ncdf_varget,fid,xid,hid_cr ;ice crystals
xid=ncdf_varid(fid,'hid_score_AG') & ncdf_varget,fid,xid,hid_ag ;aggregates
xid=ncdf_varid(fid,'hid_score_WS') & ncdf_varget,fid,xid,hid_ws ;wet snow
xid=ncdf_varid(fid,'hid_score_VI') & ncdf_varget,fid,xid,hid_vi ;vertical ice
xid=ncdf_varid(fid,'hid_score_GL') & ncdf_varget,fid,xid,hid_gl ;graupel
xid=ncdf_varid(fid,'hid_score_GH') & ncdf_varget,fid,xid,hid_gh ;graupel HD
xid=ncdf_varid(fid,'hid_score_HA') & ncdf_varget,fid,xid,hid_ha ;hail
xid=ncdf_varid(fid,'hid_score_HA') & ncdf_varget,fid,xid,hid_ha ;hail
xid=ncdf_varid(fid,'hid_score_BD') & ncdf_varget,fid,xid,hid_bd ;big drops/melting hail
ncdf_close,fid

fid=ncdf_open(path_prefix+cdir+cfile)
vid=ncdf_varid(fid,'base_time') & ncdf_varget,fid,vid,base_time_cfr
vid=ncdf_varid(fid,'time_offset') & ncdf_varget,fid,vid,time_offset_cfr
xid=ncdf_varid(fid,'x') & ncdf_varget,fid,xid,x_m_cfr
xid=ncdf_varid(fid,'y') & ncdf_varget,fid,xid,y_m_cfr
xid=ncdf_varid(fid,'z') & ncdf_varget,fid,xid,z_m_cfr
vid=ncdf_varid(fid,'point_latitude') & ncdf_varget,fid,vid,lat_3d_cfr
vid=ncdf_varid(fid,'point_longitude') & ncdf_varget,fid,vid,lon_3d_cfr
vid=ncdf_varid(fid,'point_altitude') & ncdf_varget,fid,vid,alt_3d_cfr
vid=ncdf_varid(fid,'corrected_specific_differential_phase') & ncdf_varget,fid,vid,kdp_cfr
vid=ncdf_varid(fid,'corrected_reflectivity') & ncdf_varget,fid,vid,dbz_cfr
ncdf_close,fid

;  Convert distances to km
x_km_hid=x_m_hid/1000.0
y_km_hid=y_m_hid/1000.0
z_km_hid=z_m_hid/1000.0
x_km_cfr=x_m_cfr/1000.0
y_km_cfr=y_m_cfr/1000.0
z_km_cfr=z_m_cfr/1000.0

;  Choose the 2km height level
r=where(z_km_hid le 2.0,c)
idx_hid=r[-1]
r=where(z_km_cfr le 2.0,c)
idx_cfr=r[-1]

r=where(x_km_hid lt 0 and x_km_hid ge min(x_km_hid),c)
lat_3d_hid=lat_3d_hid[r,*,*]
lon_3d_hid=lon_3d_hid[r,*,*]
kdp_hid=kdp_hid[r,*,*]
;hid=hid[r,*,*]
;rhohv=rhohv[r,*,*]
dbz_hid=dbz_hid[r,*,*]
;zdr=zdr[r,*,*]
x_km_hid=x_km_hid[r]

r=where(x_km_cfr lt 0 and x_km_cfr ge min(x_km_hid),c)
lat_3d_cfr=lat_3d_cfr[r,*,*]
lon_3d_cfr=lon_3d_cfr[r,*,*]
kdp_cfr=kdp_cfr[r,*,*]
;hid=hid[r,*,*]
;rhohv=rhohv[r,*,*]
dbz_cfr=dbz_cfr[r,*,*]
;zdr=zdr[r,*,*]
x_km_cfr=x_km_cfr[r]
r=where(y_km_cfr lt max(y_km_hid) and y_km_cfr ge min(y_km_hid),c)
lat_3d_cfr=lat_3d_cfr[*,r,*]
lon_3d_cfr=lon_3d_cfr[*,r,*]
kdp_cfr=kdp_cfr[*,r,*]
;hid=hid[*,r,*]
;rhohv=rhohv[*,r,*]
dbz_cfr=dbz_cfr[*,r,*]
;zdr=zdr[*,r,*]
y_km_cfr=y_km_cfr[r]


;  Find the lat,lon range of the plot area
ulat_hid=max(lat_3d_hid)
llat_hid=min(lat_3d_hid)
llon_hid=min(lon_3d_hid)
rlon_hid=max(lon_3d_hid)
if llon_hid lt 0 then lcolon_hid=360.0+llon_hid else lcolon_hid=llon_hid
if rlon_hid lt 0 then rcolon_hid=360.0+rlon_hid else rcolon_hid=rlon_hid

;*************************
;  Temp image
;*************************
;  Set up plot size
pxdim=900 & pydim=900
;  Position the plots
xl=0.1 & xr=0.90
yb=0.05 & yt=0.85
sx=0.1
sy=0.1
numplots_x=2
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.02  ;vertical
cbpos[*,2]=cbpos[*,0]+0.01
;cbpos[*,1]=pos[*,3]+0.05  ;horizontal
;cbpos[*,3]=cbpos[*,1]+0.01
dx=0.01 & dy=0.01
;Clear plot 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])

;  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([230,230,230])]
mycbtable=mytable[0:top_color,*]
;                  drizzle    rain          ice crystals  aggregates   wet snow
hidct=['white smoke','sky blue','medium blue','orange',     'pink',      'aqua',$
  ;vertical ice  LD Graupel   HD Graupel  Hail  Big Drops/Melting Hail
  'dark gray',  'lime green', 'yellow',  'red','magenta']

;*** PLOT KDP
pnum=0
data_var=reform(kdp_hid[*,*,idx_hid])
lat_var=reform(lat_3d_hid[*,*,idx_hid])
lon_var=reform(lon_3d_hid[*,*,idx_hid])
height_var=z_km_hid[idx_hid]
dmax=max(data_var)
r=where(data_var ne -9999,c)
dmin=min(data_var[r])
data_image=bytscl(data_var,top=top_color,max=dmax,min=dmin)
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
p0=image(data_image,x_km_hid,y_km_hid,/current,$
    dimensions=[pxdim,pydim],position=pos[pnum,*],$
    font_size=fsm,rgb_table=mytable)
c0=contour(data_var,x_km_hid,y_km_hid,/current,position=pos[pnum,*],/nodata,$
    /xstyle,/ystyle,xtickdir=1,ytickdir=1)
c0=colorbar(rgb_table=mycbtable,range=[dmin,dmax],orientation=1,$
    position=reform(cbpos[pnum,*]),font_size=10,tickdir=1,border_on=1,title='deg/km')
t1=text(pos[pnum,0]+0*dx,pos[pnum,3]+1*dy,$
    'KDP '+string(height_var,format='(F5.2)')+'km',font_size=fs1)
t1=text(pos[pnum,0]+0*dx,pos[pnum,3]+3*dy,$
  hfile,font_size=10)

pnum=1
data_var=reform(kdp_cfr[*,*,idx_cfr])
lat_var=reform(lat_3d_cfr[*,*,idx_cfr])
lon_var=reform(lon_3d_cfr[*,*,idx_cfr])
height_var=z_km_cfr[idx_cfr]
;dmax=max(data_var)
;r=where(data_var ne -9999,c)
;dmin=min(data_var[r])
data_image=bytscl(data_var,top=top_color,max=dmax,min=dmin)
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
p0=image(data_image,x_km_cfr,y_km_cfr,/current,$
  dimensions=[pxdim,pydim],position=pos[pnum,*],$
  font_size=fsm,rgb_table=mytable)
c0=contour(data_var,x_km_cfr,y_km_cfr,/current,position=pos[pnum,*],/nodata,$
  /xstyle,/ystyle,xtickdir=1,ytickdir=1)
c0=colorbar(rgb_table=mycbtable,range=[dmin,dmax],orientation=1,$
  position=reform(cbpos[pnum,*]),font_size=10,tickdir=1,border_on=1,title='deg/km')
t1=text(pos[pnum,0]+0*dx,pos[pnum,3]+1*dy,$
  'KDP '+string(height_var,format='(F5.2)')+'km',font_size=fs1)
t1=text(pos[pnum,0]+0*dx,pos[pnum,3]+3*dy,$
  cfile,font_size=10)








hstr=string(height_var,format='(F4.2)')+'km'
p0.save,'hid.cfrad'+hdate+'.'+htime+'.'+hstr+'.kdp.c5.png',height=pydim

stop

;*** PLOT HID
do_hid='no'
if do_hid eq 'yes' then begin
pnum=0
data_var=reform(hid[*,*,sidx])
data_lat=lat_hid_3d[*,*,sidx]
data_lon=lon_hid_3d[*,*,sidx]
height_var=z_km[sidx]
data_image=byte(data_var)
r=where(data_var eq -9999,c)
if c gt 0 then data_image[r]=0
p0=image(data_image,/buffer,$
  dimensions=[pxdim,pydim],position=pos[pnum,*],$
  rgb_table=hidct,font_size=fsm)
c0=contour(data_var,x_km,y_km,/current,position=pos[pnum,*],/nodata,$
  /xstyle,/ystyle,xtickdir=1,ytickdir=1)
t1=text(pos[pnum,0]+0*dx,pos[pnum,3]+1*dy,$
  'HID '+string(height_var,format='(F5.2)')+'km',font_size=fs1)
endif

;*** PLOT KDP
pnum=1
data_var=reform(kdp[*,*,sidx])
data_lat=lat_hid_3d[*,*,sidx]
data_lon=lon_hid_3d[*,*,sidx]
height_var=z_km[sidx]
kmax=max(data_var)
r=where(data_var ne -9999,c)
kmin=min(data_var[r])
data_image=bytscl(data_var,max=kmax,min=kmin,top=254)
r=where(data_var eq -9999,c)
if c gt 0 then data_image[r]=255
p0=image(data_image,/current,$
  dimensions=[pxdim,pydim],position=pos[pnum,*],$
  font_size=fsm,rgb_table=39)
c0=contour(data_var,x_km,y_km,/current,position=pos[pnum,*],/nodata,$
  /xstyle,/ystyle,xtickdir=1,ytickdir=1)
c0=colorbar(rgb_table=39,range=[kmin,kmax],orientation=1,$
    position=reform(cbpos[pnum,*]),font_size=10,tickdir=1,border_on=1,title='deg/km')
t1=text(pos[pnum,0]+0*dx,pos[pnum,3]+1*dy,$
  'KDP '+string(height_var,format='(F5.2)')+'km',font_size=fs1)
  
;*** PLOT DBZ
pnum=2
data_var=reform(dbz[*,*,sidx])
data_lat=lat_hid_3d[*,*,sidx]
data_lon=lon_hid_3d[*,*,sidx]
height_var=z_km[sidx]
kmax=max(data_var)
r=where(data_var ne -9999,c)
kmin=min(data_var[r])
data_image=bytscl(data_var,max=kmax,min=kmin,top=254)
r=where(data_var eq -9999,c)
if c gt 0 then data_image[r]=255
p0=image(data_image,x_km,y_km,/current,$
  dimensions=[pxdim,pydim],position=pos[pnum,*],$
  font_size=fsm,rgb_table=39)
c0=contour(data_var,x_km,y_km,/current,position=pos[pnum,*],/nodata,$
  /xstyle,/ystyle,xtickdir=1,ytickdir=1)
c0=colorbar(rgb_table=39,range=[kmin,kmax],orientation=1,$
  position=reform(cbpos[pnum,*]),font_size=10,tickdir=1,border_on=1,title='dBZ')
t1=text(pos[pnum,0]+0*dx,pos[pnum,3]+1*dy,$
  'Reflectivity '+string(height_var,format='(F5.2)')+'km',font_size=fs1)


;*** PLOT ZDR
pnum=0
data_var=reform(dbz_zdr[*,*,sidx])
data_lat=lat_hid_3d[*,*,sidx]
data_lon=lon_hid_3d[*,*,sidx]
height_var=z_km[sidx]
kmax=max(data_var)
r=where(data_var ne -9999,c)
kmin=min(data_var[r])
data_image=bytscl(data_var,max=kmax,min=kmin,top=254)
r=where(data_var eq -9999,c)
if c gt 0 then data_image[r]=255
p0=image(data_image,x_km,y_km,/current,$
  dimensions=[pxdim,pydim],position=pos[pnum,*],$
  font_size=fsm,rgb_table=39)
c0=contour(data_var,x_km,y_km,/current,position=pos[pnum,*],/nodata,$
  /xstyle,/ystyle,xtickdir=1,ytickdir=1)
c0=colorbar(rgb_table=39,range=[kmin,kmax],orientation=1,$
  position=reform(cbpos[pnum,*]),font_size=10,tickdir=1,border_on=1,title='dBZ')
t1=text(pos[pnum,0]+0*dx,pos[pnum,3]+1*dy,$
  'ZDR '+string(height_var,format='(F5.2)')+'km',font_size=fs1)
  
;*** PLOT rhohv
pnum=3
data_var=reform(rhohv[*,*,sidx])
data_lat=lat_hid_3d[*,*,sidx]
data_lon=lon_hid_3d[*,*,sidx]
height_var=z_km[sidx]
kmax=max(data_var)
r=where(data_var ne -9999,c)
kmin=min(data_var[r])
data_image=bytscl(data_var,max=kmax,min=kmin,top=254)
r=where(data_var eq -9999,c)
if c gt 0 then data_image[r]=255
p0=image(data_image,/current,$
  dimensions=[pxdim,pydim],position=pos[pnum,*],$
  font_size=fsm,rgb_table=39)
c0=contour(data_var,x_km,y_km,/current,position=pos[pnum,*],/nodata,$
  /xstyle,/ystyle,xtickdir=1,ytickdir=1)
c0=colorbar(rgb_table=39,range=[kmin,kmax],orientation=1,$
  position=reform(cbpos[pnum,*]),font_size=10,tickdir=1,border_on=1,title='ratio')
t1=text(pos[pnum,0]+0*dx,pos[pnum,3]+1*dy,$
  'RHOHV '+string(height_var,format='(F5.2)')+'km',font_size=fs1)

;*** TITLES
t1=text(0.1,0.98,'CBAND '+hdate+' '+$
  hhour+':'+hmin+':'+hsec,font_size=fs1)
;t1=text(pos[pnum,0]+0*dx,pos[pnum,3]+6*dy,$
;  'HID '+string(height_var,format='(F5.2)')+'km',font_size=fs1)
;t1=text(pos[pnum,0]+0*dx,pos[pnum,3]+4*dy,$
;  hid_tag+' grid',font_size=fs1)
if do_hid eq 'yes' then begin
xh=0.3  & yh=0.91    
;t1=text(pos[pnum,0],pos[pnum,3]+2*dy,'HID',font_size=fs2)
t1=text(xh+15*dx,yh+6.0*dy,'Drizzle',font_size=fs2,color=hidct[1])
t1=text(xh+15*dx,yh+4.0*dy,'Rain',font_size=fs2,color=hidct[2])
t1=text(xh+15*dx,yh+2.0*dy,'Ice Crystals',font_size=fs2,color=hidct[3])
t1=text(xh+29*dx,yh+6.0*dy,'Aggregates',font_size=fs2,color=hidct[4])
t1=text(xh+29*dx,yh+4.0*dy,'Wet Snow',font_size=fs2,color=hidct[5])
t1=text(xh+29*dx,yh+2.0*dy,'Vertical Ice',font_size=fs2,color=hidct[6])
t1=text(xh+42*dx,yh+6.0*dy,'LD Graupel',font_size=fs2,color=hidct[7])
t1=text(xh+42*dx,yh+4.0*dy,'HD Graupel',font_size=fs2,color=hidct[8])
t1=text(xh+42*dx,yh+2.0*dy,'Hail',font_size=fs2,color=hidct[9])
t1=text(xh+55*dx,yh+6.0*dy,'Big Drops',font_size=fs2,color=hidct[10])
t1=text(xh+55*dx,yh+4.0*dy,'Melting Hail',font_size=fs2,color=hidct[10])
endif



hstr=string(height_var,format='(F4.2)')+'km'
p0.save,'hid.'+hdate+'.'+hhour+hmin+hsec+'.'+hstr+'.'+hid_tag+'.v4.png',height=pydim
stop
end