pro plot_gpm_2adpr

;******************
;  Input
;******************

;  Change this if you are chpc or imac
path_prefix='/Volumes/'  ;imac
;path_prefix='/uufs/chpc.utah.edu/common/home/'  ;chpc

;  File directory and name
fdir=path_prefix+'mace-group5/gpm/GPM_L2/GPM_2ADPR.06/'
fname='2A.GPM.DPR.V8-20180723.20180209-S141035-E154309.022448.V06A.HDF5'

;radar='ku'
;radar='ka'
;if radar eq 'ku' then begin
;  radarname='KuPR'
;  radarfreq='13.6 GHz'
;  radarswath='245 km'
;endif else if radar eq 'ka' then begin
;  radarname='KaPR'
;  radarfreq='35.5 GHz'
;  radarswath='120 km'
;endif

;  Set up the map bounds 
;  ulat=upper latitude of box
;  llat=lower latitude of box
;  llon=left longitude of box
;  rlon=right longitude of box
ulat=-50.0 & llat=-70.0
llon=120 & rlon=170
;  colongitude
if llon lt 0 then lcolon=360.0+llon else lcolon=llon
if rlon lt 0 then rcolon=360.0+rlon else rcolon=rlon

;******************
;  Read data from hdf5 file
;******************
;
;  This is the idl hdf5 ncdump
;presult=h5_parse(fdir+fname)

;  Open an existing hdf5 file
file_id=h5f_open(fdir+fname)
ns_id=h5g_open(file_id,'NS')  ;MS
;ns_id=h5g_open(file_id,'HS')
dset_id=h5d_open(ns_id,'Latitude')
lat=h5d_read(dset_id)
h5d_close,dset_id

dset_id=h5d_open(ns_id,'Longitude')
lon=h5d_read(dset_id)
h5d_close,dset_id

slv_id=h5g_open(ns_id,'SLV')
dset_id=h5d_open(slv_id,'zFactorCorrected')
dbz=h5d_read(dset_id)
h5d_close,dset_id
h5g_close,slv_id

time_id=h5g_open(ns_id,'ScanTime')
dset_id=h5d_open(time_id,'DayOfYear')
doy=h5d_read(dset_id)
h5d_close,dset_id
dset_id=h5d_open(time_id,'Year')
year=h5d_read(dset_id)
h5d_close,dset_id
dset_id=h5d_open(time_id,'Month')
month=h5d_read(dset_id)
h5d_close,dset_id
dset_id=h5d_open(time_id,'DayOfMonth')
day=h5d_read(dset_id)
h5d_close,dset_id
dset_id=h5d_open(time_id,'Hour')
hour=h5d_read(dset_id)
h5d_close,dset_id
dset_id=h5d_open(time_id,'Minute')
minute=h5d_read(dset_id)
h5d_close,dset_id
dset_id=h5d_open(time_id,'Second')
second=h5d_read(dset_id)
h5d_close,dset_id
dset_id=h5d_open(time_id,'MilliSecond')
msecond=h5d_read(dset_id)
h5d_close,dset_id
h5g_close,time_id

h5g_close,ns_id
h5f_close,file_id

;****************
;  nbin,nray,nscan
;****************

s=size(dbz,/dimensions)
nbin=s[0]  ;176
nray=s[1]  ;49
nscan=s[2] ;7937

;*****************
;  Calculate colongitude
;*****************

colon=lon
r=where(lon lt 0,c)
if c gt 0 then colon[r]=colon[r]+360.0

;*****************
;  Calculate time
;*****************

msecond=msecond*1e-3
second=second+msecond
julian_day=julday(month,day,year,hour,minute,second)

;  loop through the range bins to see where the data is
;for i=0,nbin-1 do begin
;  dbz1=reform(dbz[i,*,*])
;  result=where(dbz1 gt -9999,count)
;  if count le 0 then begin
;    print,i,count
;  endif else if count gt 0 then begin
;    print,i,count,min(dbz1[result]),max(dbz1[result])
;  endif
;endfor

;*****************************
;  Subset the data in the region of interest box          
;*****************************

if (llon lt rlon and lcolon gt rcolon) or $
  (llon lt rlon and rcolon gt lcolon) then begin
  print,'stradeling 0 degree'
  idx=where(lon ge llon and lon le rlon and $
            lat ge llat and lat le ulat,cidx)
  julian_day=julian_day[idx]
  lat=lat[idx]
  lon=lon[idx]
  colon=colon[idx]
  ;  Create a regular grid to interpolated to.
  dlat=0.05
  dlon=0.1
  grid_lat=floor(llat) & while max(grid_lat) lt ulat do grid_lat=[grid_lat,max(grid_lat)+dlat]
  grid_lon=floor(llon) & while max(grid_lon) lt rlon do grid_lon=[grid_lon,max(grid_lon)+dlon]
  grid_colon=grid_lon
  result=where(grid_lon lt 0,count)
  if count gt 0 then grid_colon[result]=grid_colon[result]+360.0
endif else begin
  ;  Subset data within colon bounds
  idx=where(colon ge lcolon and colon le rcolon and $
    lat ge llat and lat le ulat,cidx)
  julian_day=julian_day[idx]
  lat=lat[idx]
  lon=lon[idx]
  colon=colon[idx]
  ;  Create a regular grid that the seviri data will be interpolated to.
  dlat=0.05
  dlon=0.1
  grid_lat=floor(llat) & while max(grid_lat) lt ulat do grid_lat=[grid_lat,max(grid_lat)+dlat]
  grid_colon=floor(lcolon) & while max(grid_colon) lt rcolon do grid_colon=[grid_colon,max(grid_colon)+dlon]
  grid_lon=grid_colon
  result=where(grid_lon gt 180,count)
  if count gt 0 then grid_lon[result]=grid_lon[result]-360.0
endelse
print,llon,rlon
print,lcolon,rcolon

;  Subset the dbz
dbz_sub=make_array(nbin,cidx)
;  loop through the range bins 
for i=0,nbin-1 do begin
  ;  Get the data at a range bin
  dbz1=reform(dbz[i,*,*])
  dbz_sub[i,*]=dbz1[idx]
  dbz2=dbz_sub[i,*]
  result=where(dbz2 gt -9999,count)
  if count le 0 then begin
    print,i,count
  endif else if count gt 0 then begin
    print,i,count,min(dbz2[result]),max(dbz2[result])
  endif       
endfor  

;*************************************************
;  Plot lat-lon track and timeseries of variables
;*************************************************

do_plot='yes'
if do_plot eq 'yes' then begin
  pxdim=900
  pydim=700
  dummy=label_date(date_format=['%Z/%M/%D!C%H:%I'])
  ;  Set up the map coordinates
  p0=map('cylindrical',dimensions=[pxdim,pydim],/buffer,$
    position=[0.1,0.4,0.9,0.9])
  m1=mapcontinents(fill_color='pale goldenrod')
  p1=plot(lon[*,0],lat[*,0],color='blue',thick=3,/overplot)  ;geoprof
  var=dbz_sub[158,*]
  result=where(var gt -9999.90,count)  
  dmin=min(var[result])
  dmax=max(var[result])
  p2=plot(julian_day,var,/current,$
    position=[0.1,0.1,0.9,0.4],xtickformat='label_date',$
    /xstyle,/ystyle,color='red',yrange=[dmin,dmax])
  imagename='map_2adpr.png'
  p0.save,imagename,height=pydim
endif  ;end of plot to check

;***************************
;  Make a plot of this subset
;***************************

;  Set up plot size
pxdim=900 & pydim=600
;  Position the plots
xl=0.07 & xr=0.80
yb=0.10 & yt=0.85
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
;  Colorbar position
cbpos=pos
cbpos[*,0]=pos[*,2]+0.1
cbpos[*,2]=cbpos[*,0]+0.01

;********************************
;  My 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([200,200,200])]
mycbtable=mytable[0:top_color,*]

;  Strings
dx=0.01 & dy=0.01
font_size=12

;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])

;  Choose the data to plot
hidx=158
data_var=dbz_sub[hidx,*]
data_lat=lat
data_lon=lon
;  mask interpolation
grid_input,data_lon,data_lat,data_var,xyz,newdata_var,/degree,/sphere,duplicates='Avg',epsilon=0.01
lon = !radeg * atan(xyz[1,*],xyz[0,*]) &   lat = !radeg * asin(xyz[2,*])
;  Triangulate the data
qhull,lon,lat,triangles,/delaunay,sphere=s
;  Put satellite data on regular grid
grid_mask=griddata(lon,lat,newdata_var,xout=grid_lon,yout=grid_lat,$
  /grid,/degrees,/sphere,triangles=triangles,$
  /kriging,min_points=16,sectors=8,empty_sectors=3,missing=-8888)
;  Restrict the data to valid values
r=where(data_var ne -9999.9,c)
good_var=data_var[r]
good_lat=data_lat[r]
good_lon=data_lon[r]
;  Pass valid data through grid_input to get rid of duplicates
print,'grid input'
grid_input,good_lon,good_lat,good_var,xyz,newdata_var,/degree,/sphere,duplicates='Avg',epsilon=0.01
lon = !radeg * atan(xyz[1,*],xyz[0,*]) &   lat = !radeg * asin(xyz[2,*])
print,'qull'
;  Triangulate the data
qhull,lon,lat,triangles,/delaunay,sphere=s
;  Put satellite data on regular grid
print,'griddata'
grid_var=griddata(lon,lat,newdata_var,xout=grid_lon,yout=grid_lat,$
  /grid,/degrees,/sphere,triangles=triangles,$
  /kriging,min_points=16,sectors=8,empty_sectors=3,missing=-9999)
;  Find the gridded data min and max
result=where(grid_var gt -7777)
dmin=min(grid_var[result])
dmax=max(grid_var[result])
;  Bytscal the data
data_image=bytscl(grid_var,top=top_color,min=dmin,max=dmax)
;  Put missing data to gray color
result=where(grid_var eq -9999,count)
if count gt 0 then data_image[result]=253
result=where(grid_mask eq -8888,count)
if count gt 0 then data_image[result]=255
;  Plot the mapped image of the data
p0=image(data_image,grid_lon,grid_lat,limit=[llat,llon,ulat,rlon],$
  /buffer,dimensions=[pxdim,pydim],grid_units='degrees',$
  center_longitude=((rlon-llon)/2.0)+llon,/box_axes,$
  label_position=0,font_size=12,linestyle=2,$
  rgb_table=mytable,$
  map_projection='Mercator',$
  position=pos[pnum,*])
;  Add continent lines
p1=mapcontinents(color='white',/hires,thick=2)
;  Put on a colorbar
c0=colorbar(range=[dmin,dmax],$
  rgb_table=mycbtable,$
  orientation=1,position=reform(cbpos[pnum,*]),font_size=10,$
  tickdir=1,border_on=1,title='Kelvin')



;  Bytscal the data
;result=where(grid_var gt 0)
;dmin=min(grid_var[result])
;dmax=max(grid_var[result])
;data_image=bytscl(grid_var,top=top_color,min=dmin,max=dmax)
;;result=where(grid_var gt -8000 and grid_var le -7000,count)
;if count gt 0 then data_image[result]=255
;result=where(grid_var gt -10000 and grid_var le -9000,count)
;if count gt 0 then data_image[result]=253

;p0=image(data_image,grid_colon,grid_lat,limit=[llat,lcolon,ulat,rcolon],$  
;  /buffer,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,*])
;  
;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
;
;  c0=colorbar($
;    rgb_table=mycbtable,$
;    range=[dmin,dmax],$
;    orientation=1,position= reform(cbpos[pnum,*]),font_size=10,$
;    tickdir=1,border_on=1,title='DbZ')
print,dmin,dmax 
stime=min(julian_day)
etime=max(julian_day)
caldat,stime,smm,sdd,syy,shh,smi,sss
caldat,etime,emm,edd,eyy,ehh,emi,ess

sdate_str=string(syy,format='(I4)')+string(smm,format='(I02)')+string(sdd,format='(I02)')
stime_str=string(shh,format='(I02)')+string(smi,format='(I02)')+string(sss,format='(I02)')  
edate_str=string(eyy,format='(I4)')+string(emm,format='(I02)')+string(edd,format='(I02)')
etime_str=string(ehh,format='(I02)')+string(emi,format='(I02)')+string(ess,format='(I02)')  

stime_cstr=string(shh,format='(I02)')+':'+string(smi,format='(I02)')+':'+string(sss,format='(I02)') 
etime_cstr=string(ehh,format='(I02)')+':'+string(emi,format='(I02)')+':'+string(ess,format='(I02)')
datastream='GPM_2ADPR.05'
;t3=text(0.1,0.90,'Channel: '+channels[chan])
t1=text(0.1,0.96,datastream+' '+sdate_str+' '+stime_cstr+' - '+edate_str+' '+etime_cstr)
t2=text(0.1,0.93,'Dual-Frequency Precipitation Radar (DPR)')
t3=text(0.1,0.90,'Reflectivity Factor with attenuation correction (dBZ) at height bin '+string(hidx,format='(I03)') )
imagename='map.'+datastream+'.dbz.'+$
  sdate_str+stime_str+'.'+edate_str+etime_str+'.NS1.png'
    
p0.save,imagename,height=pydim
stop
end