;********************************************
;  This plots a plan view (horizontal slice on a map)
;  Constant Latitude slice, constant longitude slice, and line segment slice
;  4 plot.
;**************************************************

pro cband_stats_observables_1hr

;rtag='oceanic'
;rtag='mountain'


;  Make a 1hr time series
julian_day_1hr=timegen(start=julday(02,21,2018,0,0,0),$
  final=julday(02,22,2018,0,0,0),$
  ;final=julday(01,15,2016,0,0,0),$  ;whole data set
  units='Hours',step_size=1)
caldat,julian_day_1hr,mm,dd,yy,hh,mi,ss
;for i=0,n_elements(julian_day_1hr)-1 do begin
;  print,yy[i],mm[0],dd[i],hh[i],mi[i],ss[i]
;endfor
;print,yy[0],mm[0],dd[0],hh[0],mi[0],ss[0]
;print,yy[-1],mm[-1],dd[-1],hh[-1],mi[-1],ss[-1]  

;  Number of hours in the time series
numhours=n_elements(julian_day_1hr)

;  chpc or imac
;path_prefix='/Volumes/'  
path_prefix='/uufs/chpc.utah.edu/common/home/'

for k=0,numhours-1 do begin
  
  ;  Get year month day strings
  syy=string(yy[k],format='(I4)')
  smm=string(mm[k],format='(I02)')
  sdd=string(dd[k],format='(I02)')
  shh=string(hh[k],format='(I02)')

  ;  Choose date and time of cband 
  ;date_str='20180209' & time_str='23*' 
  ;date_str='20180209' & time_str='140504'
  ;date_str='20151204' & time_str='163710'  ;olympex

  ;  Get year month day strings
  ;syy=strmid(date_str,0,4)
  ;smm=strmid(date_str,4,2)
  ;sdd=strmid(date_str,6,2)

  ;shh=strmid(time_str,0,2)
  ;smi=strmid(time_str,2,2)
  ;sss=strmid(time_str,4,2)
  date_str=syy+smm+sdd;+shh+smi+sss
  ;ymd_str=syy+smm+sdd

  ;  directory
  if syy eq '2018' then begin
    ;filedir=path_prefix+'mace-group4/Capricorn2/CBand/'
    ;filedir=path_prefix+'mace-group4/Capricorn2/CBand/cfradial/'
    filedir=path_prefix+'mace-group5/field_programs/investigator/in2018_v01/cband/'+date_str+'/'
    filestr='9776HUB-PPIVol-'+syy+smm+sdd+'-'+shh+'*stats_obs_histo.rho_filter.cdf'
    radar_name='cband'
  endif else if syy eq '2019' then begin
    filedir=path_prefix+'mace-group5/field_programs/investigator/in2019_v01_gpm_passes/out_vol/'
    filestr='9776HUB-PPIVol-'+syy+smm+sdd+'-'+time_str+'*pyart*.nc'
  endif else if syy eq '2015' or syy eq '2016' then begin
    filedir=path_prefix+'mace-group6/field_programs/olympex/sband/'+date_str+'/ppi/'
    filestr='olympex_NPOL1_'+syy+smm+sdd+'_'+shh+'*pyart_stats_obs_histo.rho_filter.'+rtag+'.cdf'
    radar_name='sband'
    ;  Get the flight path
    ;read_nav_olympex,path_prefix,date_str,julian_day_plane,temp,tas,theta,pres,ppmvw,lat_plane,lon_plane,$
    ;  alt_plane,wwind,uwind,vwind,cwc_kng,twc_nev,lwc_nev,iwc_nev,tot_s,cwc_s,rh,cwc_csi,tot_c,tot_p
  endif

  ;  file name
  filenames=file_search(filedir+filestr,count=numfiles)
  print,'numfiles',numfiles,' date ',date_str+shh
  obs_counts4d_1hr=!NULL  
  
  ;  Only continue if there are files for this hour
  if numfiles gt 0 then begin
  
    ;  Loop through the cband files
    for f=0,numfiles-1 do begin

      filename=filenames[f]
      print,filename
 
      ;print,'start reading data'  
      fid=ncdf_open(filename)
      vid=ncdf_varid(fid,'bins_dbz') & ncdf_varget,fid,vid,bins_dbz
      vid=ncdf_varid(fid,'bins_zdr') & ncdf_varget,fid,vid,bins_zdr 
      vid=ncdf_varid(fid,'bins_kdp') & ncdf_varget,fid,vid,bins_kdp
      vid=ncdf_varid(fid,'bins_hgt') & ncdf_varget,fid,vid,bins_hgt 
      vid=ncdf_varid(fid,'counts_array') & ncdf_varget,fid,vid,obs_counts4d
      ncdf_close,fid
      ;print,'end reading data'
  
      s=size(obs_counts4d,/dimensions)
      num_dbz=s[0]
      num_zdr=s[1]
      num_kdp=s[2]
      num_hgt=s[3]
  
      ;  Filetag
      file_tag='counts.'+radar_name+'.'+syy+smm+sdd+'.'+shh  
  
      ;  Add up the counts in the 1hr
      if obs_counts4d_1hr eq !NULL then begin
        obs_counts4d_1hr=obs_counts4d
      endif else begin
        obs_counts4d_1hr=obs_counts4d_1hr+obs_counts4d
      endelse
    endfor  ;end of loop through cband files for this hour

    ;  Write to a netcdf file
    if 1 eq 1 then begin
      if rtag eq !NULL then begin
        output_filename=filedir+'stats_obs_histo.rho_filter.'+syy+smm+sdd+'.'+shh+'.cdf'
      endif else begin
        output_filename=filedir+'stats_obs_histo.rho_filter.'+rtag+'.'+syy+smm+sdd+'.'+shh+'.cdf'
      endelse
      ;output_filename=filedir+file_tag+'_stats_obs_histo.cdf'
      fid=ncdf_create(output_filename,/clobber)
      print, output_filename,fid,'start write'
      bins_dbz_did=ncdf_dimdef(fid,'num_bins_dbz',n_elements(bins_dbz))
      bins_zdr_did=ncdf_dimdef(fid,'num_bins_zdr',n_elements(bins_zdr))
      bins_kdp_did=ncdf_dimdef(fid,'num_bins_kdp',n_elements(bins_kdp))
      bins_hgt_did=ncdf_dimdef(fid,'num_bins_hgt',n_elements(bins_hgt))
      bins_dbz_id=ncdf_vardef(fid,'bins_dbz',bins_dbz_did,/float)
      ncdf_attput,fid,bins_dbz_id,'long_name','dbz bins'
      ncdf_attput,fid,bins_dbz_id,'units','dbz'
      bins_zdr_id=ncdf_vardef(fid,'bins_zdr',bins_zdr_did,/float)
      ncdf_attput,fid,bins_zdr_id,'long_name','zdr bins'
      ncdf_attput,fid,bins_zdr_id,'units','zdr'
      bins_kdp_id=ncdf_vardef(fid,'bins_kdp',bins_kdp_did,/float)
      ncdf_attput,fid,bins_kdp_id,'long_name','kdp bins'
      ncdf_attput,fid,bins_kdp_id,'units','kdp'
      bins_hgt_id=ncdf_vardef(fid,'bins_hgt',bins_hgt_did,/float)
      ncdf_attput,fid,bins_hgt_id,'long_name','height bins'
      ncdf_attput,fid,bins_hgt_id,'units','km'
      hist_id=ncdf_vardef(fid,'counts_array',[bins_dbz_did,bins_zdr_did,bins_kdp_did,bins_hgt_did],/long)
      ncdf_attput,fid,hist_id,'long_name','counts of obs'
      ncdf_control,fid,/endef
      ncdf_varput,fid,bins_dbz_id,bins_dbz
      ncdf_varput,fid,bins_zdr_id,bins_zdr
      ncdf_varput,fid,bins_kdp_id,bins_kdp
      ncdf_varput,fid,bins_hgt_id,bins_hgt
      ncdf_varput,fid,hist_id,obs_counts4d_1hr
      ncdf_close,fid
      print,'end write'

    endif ;  turn of write

;**************************************
;  Make the plot
;**************************************
  if 1 eq 1 then begin
  ;  Set up plot size
  pxdim=1000 & pydim=1000
  ;  Position the plots
  xl=0.1 & xr=0.95
  yb=0.06 & yt=0.85
  sx=0.11
  sy=0.15
  numplots_x=2
  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.08  ;vertical
  ;cbpos[*,2]=cbpos[*,0]+0.01
  cbpos[*,1]=cbpos[*,3]+0.05 ;horizontal
  cbpos[*,3]=cbpos[*,1]+0.01

  ;  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(43,ncolors=253)
  ;254=hot pink               ;gray=255
  mytable=[mytable,transpose([255,255,255]),transpose([238,18,137]),transpose([230,230,230])]
  mycbtable=mytable[0:top_color,*]

  ;top_color=254
  ;mytable=colortable(39,ncolors=255)
  ;mytable[0,*]=[250,250,250];  first color =gray
  ;mycbtable=mytable[1:-1,*]

  ;  Strings
  dx=0.01 & dy=0.01
  fs1=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])

  ;obs_counts4d=lonarr(bins_dbz,bins_zdr,bins_kdp,bins_hgt)

  ;**************************
  ;  hold dbz constant
  ;**************************

  pnum=0
  ;  Constant dbz
  constant_dbz=20
  r_sub=where(bins_dbz eq constant_dbz,c_sub)
  data_var3d=reform(obs_counts4d_1hr[r_sub[0],*,*,*])
  data_var=make_array(num_zdr,num_kdp)
  for i=0,num_zdr-1 do begin
    for j=0,num_kdp-1 do begin
      data_var[i,j]=total(data_var3d[i,j,*])
    endfor
  endfor
  dmin=min(data_var)
  dmax=max(data_var)
  data_image=bytscl(data_var,top=top_color,min=dmin,max=dmax)
  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
  p4=image(data_image,/current,image_dimensions=[isx,isy],$
    ;min_value=dmin_cband,max_value=dmax_cband,$
    position=pos[pnum,*],rgb_table=mytable)
  c1=contour(data_var,bins_zdr,bins_kdp,/nodata,/current,$
    ytitle='kdp',xtitle='zdr',irregular=0,xtickdir=1,ytickdir=1,$
    xstyle=1,ystyle=1,position=pos[pnum,*],font_size=fs1,axis_style=2)
  cb0=colorbar(rgb_table=mycbtable,range=[dmin,dmax],orientation=0,$
    position=cbpos[pnum,*],font_size=9,tickdir=1,border_on=1,title='counts')
  t1=text(pos[pnum,0]+dx,cbpos[pnum,3]+dy,'Constant DBZ '+string(constant_dbz,format='(F6.2)'),font_size=fs1)

  pnum=1
  ;  Constant dbz
  constant_dbz=22
  r_sub=where(bins_dbz eq constant_dbz,c_sub)
  data_var3d=reform(obs_counts4d_1hr[r_sub[0],*,*,*])
  data_var=make_array(num_zdr,num_kdp)
  for i=0,num_zdr-1 do begin
    for j=0,num_kdp-1 do begin
      data_var[i,j]=total(data_var3d[i,j,*])
    endfor
  endfor
  dmin=min(data_var)
  dmax=max(data_var)
  data_image=bytscl(data_var,top=top_color,min=dmin,max=dmax)
  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
  p4=image(data_image,/current,image_dimensions=[isx,isy],$
    ;min_value=dmin_cband,max_value=dmax_cband,$
    position=pos[pnum,*],rgb_table=mytable)
  c1=contour(data_var,bins_zdr,bins_kdp,/nodata,/current,$
    ytitle='kdp',xtitle='zdr',irregular=0,xtickdir=1,ytickdir=1,$
    xstyle=1,ystyle=1,position=pos[pnum,*],font_size=fs1,axis_style=2)
  cb0=colorbar(rgb_table=mycbtable,range=[dmin,dmax],orientation=0,$
    position=cbpos[pnum,*],font_size=9,tickdir=1,border_on=1,title='counts')
  t1=text(pos[pnum,0]+dx,cbpos[pnum,3]+dy,'Constant DBZ '+string(constant_dbz,format='(F6.2)'),font_size=fs1)

  pnum=2
  ;  Constant dbz
  constant_dbz=24
  r_sub=where(bins_dbz eq constant_dbz,c_sub)
  data_var3d=reform(obs_counts4d_1hr[r_sub[0],*,*,*])
  data_var=make_array(num_zdr,num_kdp)
  for i=0,num_zdr-1 do begin
    for j=0,num_kdp-1 do begin
      data_var[i,j]=total(data_var3d[i,j,*])
    endfor
  endfor
  dmin=min(data_var)
  dmax=max(data_var)
  data_image=bytscl(data_var,top=top_color,min=dmin,max=dmax)
  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
  p4=image(data_image,/current,image_dimensions=[isx,isy],$
    ;min_value=dmin_cband,max_value=dmax_cband,$
    position=pos[pnum,*],rgb_table=mytable)
  c1=contour(data_var,bins_zdr,bins_kdp,/nodata,/current,$
    ytitle='kdp',xtitle='zdr',irregular=0,xtickdir=1,ytickdir=1,$
    xstyle=1,ystyle=1,position=pos[pnum,*],font_size=fs1,axis_style=2)
  cb0=colorbar(rgb_table=mycbtable,range=[dmin,dmax],orientation=0,$
    position=cbpos[pnum,*],font_size=9,tickdir=1,border_on=1,title='counts')
  t1=text(pos[pnum,0]+dx,cbpos[pnum,3]+dy,'Constant DBZ '+string(constant_dbz,format='(F6.2)'),font_size=fs1)

  pnum=3
  ;  Constant dbz
  constant_dbz=28
  r_sub=where(bins_dbz eq constant_dbz,c_sub)
  data_var3d=reform(obs_counts4d_1hr[r_sub[0],*,*,*])
  data_var=make_array(num_zdr,num_kdp)
  for i=0,num_zdr-1 do begin
    for j=0,num_kdp-1 do begin
      data_var[i,j]=total(data_var3d[i,j,*])
    endfor
  endfor
  dmin=min(data_var)
  dmax=max(data_var)
  data_image=bytscl(data_var,top=top_color,min=dmin,max=dmax)
  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
  p4=image(data_image,/current,image_dimensions=[isx,isy],$
    ;min_value=dmin_cband,max_value=dmax_cband,$
    position=pos[pnum,*],rgb_table=mytable)
  c1=contour(data_var,bins_zdr,bins_kdp,/nodata,/current,$
    ytitle='kdp',xtitle='zdr',irregular=0,xtickdir=1,ytickdir=1,$
    xstyle=1,ystyle=1,position=pos[pnum,*],font_size=fs1,axis_style=2)
  cb0=colorbar(rgb_table=mycbtable,range=[dmin,dmax],orientation=0,$;tickformat='(E7.1)',$
    position=cbpos[pnum,*],font_size=9,tickdir=1,border_on=1,title='counts')
  t1=text(pos[pnum,0]+dx,cbpos[pnum,3]+dy,'Constant DBZ '+string(constant_dbz,format='(F6.2)'),font_size=fs1)

  t1=text(0.1,0.97,radar_name+' '+file_basename(output_filename),font_size=fs1)

  if rtag eq !NULL then begin
    imagename='stats_obs_histo.rho_filter.'+syy+smm+sdd+'.'+shh+'.png'
  endif else begin
    imagename='stats_obs_histo.rho_filter.'+rtag+'.'+syy+smm+sdd+'.'+shh+'.png'
  endelse
  p0.save,filedir+imagename
  
endif  ;end of turn off plot
endif  ;found files for the hour
endfor  ;end of loop through hours


stop


end