pro plot_gpm_1bgmi

;******************
;  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_L1B/GPM_1BGMI.05/'
fname='1B.GPM.GMI.TB2016.20180209-S141035-E154309.022448.V05A.HDF5'

;  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

;Swath S1 has channels 1-9: 9channels
;10V 10H 19V 19H 23V 37V 37H 89V 89H.
;Swath S2 has channels 10-13: 4channels
;166V 166H 183+/-3V 183+/-8V.
chan=7  ;channel
channels=['10V','10H','19V','19H','23V','37V','37H','89V','89H',$
          '166V','166H','183+/-3V','183+/-8V']
          
;******************
;  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)

;  Open the top level (root)
root_id=h5g_open(file_id,'/')

; Root / S1
group_id=h5g_open(root_id,'S1')
dset_id=h5d_open(group_id,'Latitude')
lat_s1=h5d_read(dset_id)
h5d_close,dset_id
dset_id=h5d_open(group_id,'Longitude')
lon_s1=h5d_read(dset_id)
h5d_close,dset_id
dset_id=h5d_open(group_id,'Tb')
tb_s1=h5d_read(dset_id)
h5d_close,dset_id
time_id=h5g_open(group_id,'ScanTime')
dset_id=h5d_open(time_id,'DayOfYear')
doy_s1=h5d_read(dset_id)
h5d_close,dset_id
dset_id=h5d_open(time_id,'Year')
year_s1=h5d_read(dset_id)
h5d_close,dset_id
dset_id=h5d_open(time_id,'Month')
month_s1=h5d_read(dset_id)
h5d_close,dset_id
dset_id=h5d_open(time_id,'DayOfMonth')
day_s1=h5d_read(dset_id)
h5d_close,dset_id
dset_id=h5d_open(time_id,'Hour')
hour_s1=h5d_read(dset_id)
h5d_close,dset_id
dset_id=h5d_open(time_id,'Minute')
minute_s1=h5d_read(dset_id)
h5d_close,dset_id
dset_id=h5d_open(time_id,'Second')
second_s1=h5d_read(dset_id)
h5d_close,dset_id
dset_id=h5d_open(time_id,'MilliSecond')
msecond_s1=h5d_read(dset_id)
h5d_close,dset_id
h5g_close,time_id
h5g_close,group_id

; Root / S2
group_id=h5g_open(root_id,'S2')
dset_id=h5d_open(group_id,'Latitude')
lat_s2=h5d_read(dset_id)
h5d_close,dset_id
dset_id=h5d_open(group_id,'Longitude')
lon_s2=h5d_read(dset_id)
h5d_close,dset_id
dset_id=h5d_open(group_id,'Tb')
tb_s2=h5d_read(dset_id)
h5d_close,dset_id
time_id=h5g_open(group_id,'ScanTime')
dset_id=h5d_open(time_id,'DayOfYear')
doy_s2=h5d_read(dset_id)
h5d_close,dset_id
dset_id=h5d_open(time_id,'Year')
year_s2=h5d_read(dset_id)
h5d_close,dset_id
dset_id=h5d_open(time_id,'Month')
month_s2=h5d_read(dset_id)
h5d_close,dset_id
dset_id=h5d_open(time_id,'DayOfMonth')
day_s2=h5d_read(dset_id)
h5d_close,dset_id
dset_id=h5d_open(time_id,'Hour')
hour_s2=h5d_read(dset_id)
h5d_close,dset_id
dset_id=h5d_open(time_id,'Minute')
minute_s2=h5d_read(dset_id)
h5d_close,dset_id
dset_id=h5d_open(time_id,'Second')
second_s2=h5d_read(dset_id)
h5d_close,dset_id
dset_id=h5d_open(time_id,'MilliSecond')
msecond_s2=h5d_read(dset_id)
h5d_close,dset_id
h5g_close,time_id
h5g_close,group_id

h5g_close,root_id
h5f_close,file_id

;****************
;  Transpose data to be nscan, npix, nchan
;****************

lat_s1=transpose(lat_s1)
lon_s1=transpose(lon_s1)
tb_s1=transpose(tb_s1)
lat_s2=transpose(lat_s2)
lon_s2=transpose(lon_s2)
tb_s2=transpose(tb_s2)

;  Convert to colongitude
colon_s1=lon_s1
r=where(lon_s1 lt 0,c)
if c gt 0 then colon_s1[r]=colon_s1[r]+360.0
colon_s2=lon_s2
r=where(lon_s2 lt 0,c)
if c gt 0 then colon_s2[r]=colon_s2[r]+360.0

;Swath S1 has channels 1-9: 9channels
;10V 10H 19V 19H 23V 37V 37H 89V 89H.
;Swath S2 has channels 10-13: 4channels
;166V 166H 183+/-3V 183+/-8V.

s=size(tb_s1,/dimensions)
nscans=s[0]
npix=s[1]
nchan1=s[2]

s=size(tb_s2,/dimensions)
nchan2=s[2]

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

msecond_s1=msecond_s1*1e-3
second_s1=second_s1+msecond_s1
julian_day_s1=julday(month_s1,day_s1,year_s1,hour_s1,minute_s1,second_s1)

msecond_s2=msecond_s2*1e-3
second_s2=second_s2+msecond_s2
julian_day_s2=julday(month_s2,day_s2,year_s2,hour_s2,minute_s2,second_s2)

;*****************************
;  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_s1 ge llon and lon_s1 le rlon and $
            lat_s1 ge llat and lat_s1 le ulat,count)
  julian_day_s1=julian_day_s1[idx]          
  lat=lat_s1[idx]
  lon=lon_s1[idx]
  colon=colon_s1[idx]
  tb=tb_s1[*,*,chan]  ;select one channel
  tb=tb[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_s1 ge lcolon and colon_s1 le rcolon and $
            lat_s1 ge llat and lat_s1 le ulat,count)
  julian_day_s1=julian_day_s1[idx]
  lat=lat_s1[idx]
  lon=lon_s1[idx]
  colon=colon_s1[idx]
  tb=tb_s1[*,*,chan]  ;select one channel
  tb=tb[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

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

do_plot='no'
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_s1[*,0],lat_s1[*,0],color='blue',thick=3,/overplot)  
  var=tb_s1[*,0,0]
  result=where(var gt -9999.90,count)  
  dmin=min(var[result])
  dmax=max(var[result])
  p2=plot(julian_day_s1,var,/current,$
    position=[0.1,0.1,0.9,0.4],xtickformat='label_date',$
    /xstyle,/ystyle,color='red',yrange=[dmin,dmax])
  imagename='map_1bgmi.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.1 & 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)
;mytable=colortable(22,ncolors=254)
;254=hot pink               ;gray=255
mytable=[mytable,transpose([238,18,137]),transpose([200,200,200])]
mycbtable=mytable[0:top_color,*]

;  16 Color tablle
;  Top is the last color to scale 256 colors, 0-255
;top_color=252
;  Colortable  0-252  253=white
;ref_table=colortable(38,ncolors=254)
;254=hot pink               ;gray=255
;ref_table=[ref_table,transpose([238,18,137]),transpose([150,150,150])]
;ref_cbtable=ref_table[0:top_color,*]

;********Color table end********

;  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 gmi data to plot
data_var=tb
data_lat=lat
data_lon=lon
;  Restrict the data to valid values
r=where(data_var ne -9999,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 ne -9999)
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]=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=cbpos[pnum,*],font_size=10,$
    tickdir=1,border_on=1,title='Kelvin')


;  Text strings for annotation and image name 
stime=min(julian_day_s1)
etime=max(julian_day_s1)
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)')  
t3=text(0.1,0.90,'Channel: '+channels[chan])
t1=text(0.1,0.93,sdate_str+' '+stime_str+' - '+edate_str+' '+etime_str)
datastream='GPM_1BGMI.05'
t2=text(0.1,0.96,datastream)
  
imagename='map_tb.'+datastream+'.'+channels[chan]+'.'+$
  sdate_str+stime_str+'.'+edate_str+etime_str+'.png'
    
p0.save,imagename,height=pydim
stop

end