;
;  Plot a monthly summary
;
pro plot_average_2d,site,searchdate,avg_int
;site='sgp'
;searchdate='200003'
;avg_int='1hr'
;
;  common containing the monthly averaged data
;
common monthly_averaged_data, base_time,time_offset,height,reflectivity,$
  reflectp, reflectc, cbase, ctop, velocity, clutter,maxclutter,spectrum,$
  lwp,vceil,numtimes,numheights,allclutter,goodclutter,temp,press,$
  vap,dir,dsol,sfpres,numlayers,twsflux,zenith,solarratio,sfcalbedo,$
  column_cfrac
;
;  Read in the monthly averaged data
;
get_average_data,site,searchdate,avg_int
;
;  Convert height from meters to km
;
height=height/1000.  ;convert from meters to km
;
;  Select how many plots to make
;
num_plots=9
;
;  Select which data to plot
;
dataname=strarr(num_plots)
dataname=['Radar Reflectivity','Doppler Velocity','Spectrum Width',$
          'Cloud Fraction','Total Points',$
          'Cloudy Points','Clutter Flag','Temperature','Pressure']
datatag=strarr(num_plots)
datatag=[ 'reflectivity',      'velocity','spectrum',$
          'cfrac',         'reflectp',$
          'reflectc',   'clutter','temp','press']
cb_title=strarr(num_plots)
cb_title=['DbZe','m/s','m/s','%','number','number','flag','kelvin',$
          'mb']
;
;  Loop through the variables to plot
;
for datanum=0,num_plots-1 do begin  ;loop through the different variables
;for datanum=6,6 do begin

print, datatag[datanum]

;******************
;  Pick up the data
;*******************

;***Reflectivity***
if datanum eq 0 then begin
  data=reflectivity

;  Manipulate the data if it is Reflectivity
  data=float(data)
  result=where(data gt -8888, count)
  if count gt 0 then data[result]=data[result]/100.

;  Filter with clutter for sgp
;  Assign all bad clutter reflectivity data -9999
  if site eq 'sgp' and n_elements(clutter) ne 0 then begin
   result=where(clutter eq 0, count)
   if count gt 0 then data[result]=-9999
  endif

;  Turn missing data to gray
  result=where(reflectp eq 0,count)
  if count gt 0 then data[result]=-8888
endif
  
;***Velocity***  
if datanum eq 1 then begin
  data=velocity

;  Manipulate the data if it is Velocity or Spectrum Width
  data=float(data)
  result=where(data gt -8888,count)
  if count gt 0 then data[result]=data[result]/1000.

;  Turn missing data to gray
  result=where(reflectp eq 0,count)
  if count gt 0 then data[result]=-8888
endif  
  
;***Spectrum***    
if datanum eq 2 then begin
  data=spectrum

;  Manipulate the data if it is Velocity or Spectrum Width
  data=float(data)
  result=where(data gt -8888,count)
  if count gt 0 then data[result]=data[result]/1000.

;  Turn missing data to gray
  result=where(reflectp eq 0,count)
  if count gt 0 then data[result]=-8888
endif  
  
;**************************
;  Calculate cloud fraction
;**************************
;
;  First if statement for twp,nsa
;
if datanum eq 3 and $
  (site eq 'nsa' or site eq 'twpc1' or site eq 'twpc2') then begin

;  Percentage of points that are cloudy
  data=float(reflectc)/float(reflectp)*100.

;  Where possible points is 0, there is missing data.
;  Where possible points is -8888, there is missing data
;  Give these points the missing data value,-8888,gray
  result=where(reflectp le 0.0 or reflectp eq -8888,count)
  print,'reflectp is 0 for ',count,' points'
  if count ne 0 then data[result]=-8888
;***
;**Old way to filter with clutter**
;  Filter with clutter if sgp
;  Where clutter indicates not a cloud, give 0 cloud fraction
;
;  if n_elements(clutter) ne 0 and site eq 'sgp' then begin
;    result=where(clutter eq 0, count)
;    print,'clutter is 0 for ',count,' points'
;    if count ne 0 then data[result]=0
;  endif
;***
;
;  Second if statement for sgp
;
endif else if datanum eq 3 and $
  (site eq 'sgp' or site eq 'crystal') then begin

;  Percentage of points that are cloudy
  data=float(goodclutter)/float(allclutter)*100.

;  Where possible points is 0, there is missing data.
;  Where possible points is -8888, there is missing data.
;  Give these points the missing data value,-8888,gray
  result=where(allclutter eq 0.0 or allclutter eq -8888,count)
  print,'allclutter is 0 for ',count,' points'
  if count ne 0 then data[result]=-8888

endif
;**************************
;  End of cloud fraction
;**************************

;***Reflectp***
if datanum eq 4 then begin
  if site eq 'sgp' then data=allclutter else data=reflectp
endif

;***Reflectc***
if datanum eq 5 then begin
  if site eq 'sgp' then data=goodclutter else data=reflectc
endif

;***Clutter***
if datanum eq 6 then data=clutter

;***Temperature***
if datanum eq 7 then begin
  data=temp
; Set bad temp values to -9999
  result=where(data lt 0,count)
  if count gt 0 then data[result]=-9999
endif

;***Pressure***  
if datanum eq 8 then begin
  data=press
;  Set values less than zero to 0
  result=where(data lt 0 and data gt -8888,count)
  if count ge 0 then data[result]=0
;  Convert pressure from pascals to millibars
  result=where(data ne -9999 and data ne -8888,count)
  if count ne 0 then data[result]=data[result]*0.01
endif
;
;  Find the max and min data values
;
dmax=max(data)
result=where(data gt -8888, count)
if count gt 0 then dmin=min(data[result]) else dmin=min(data)
if datanum eq 0 then begin  ;reflectivity
  dmin=-65.0
  dmax=25.0
endif
print, 'dmax',dmax,'dmin',dmin
;
;  Calculate the file starting and ending date and time
;
standard_date_time,base_time+time_offset(0),syear,smonth,sday,shour,smin,ssec
print,'file start',syear,smonth,sday,shour,smin,ssec
standard_date_time,base_time+time_offset(numtimes-1),eyear,emonth,eday,ehour,emin,esec
print,'file end',eyear,emonth,eday,ehour,emin,esec
;
;  Calculate the standard date arrays
;
ijday_ihrfrac_fm_numsec,base_time,time_offset,hrfrac,jday,yy,mm,dd,hh,mi,ss
;
;  Calculate IDL Julian day for plotting
;
julian_day=make_array(n_elements(jday),/double,value=-9999)
julian_day=julday(mm,dd,yy,hh,mi,ss)
;
;  Color the data
;
image=bytscl(data,max=dmax,min=dmin,top=253)
;
;  Put all no cloud values to white
;
result=where(data eq -9999,count)
if count gt 0 then image[result]=255
;
;  Put all missing values to gray
;
result=where(data eq -8888,count)
if count gt 0 then image[result]=254

;*****************************
;  Set up the plot
;*****************************

;  Set up colors
loadct,39 ;5  ;load the color table
!p.background=!d.n_colors-1  ;make the background white
common colors,r_orig,g_orig,b_orig,r_curr,g_curr,b_curr
;  make gray for missing data in color 253
r_curr[254]=200
g_curr[254]=200
b_curr[254]=200

;  Set up the rest of the plot
!y.omargin=[2,5]  ;set aside some room at the top of the page for the title
!x.margin=[8,15]  ;change to accomodate the colorbar
!y.margin=[2,2]  ;change to move the images closer together
!p.multi=[0,1,5]  ;put 5 plots on a page
!p.charsize=2  ;increase the character size
!x.ticklen=-0.02
xwin=800  ;xsize of the plot window
ywin=800  ;ysize of the plot window
window,0,xsize=xwin,ysize=ywin,title='average',/pixmap
dummy=label_date(date_format=[' %D %M !C %H:%I '])

;********
;  Plot 1
;********
plot1='yes'
if plot1 eq 'yes' then begin

week1=num_sec_ep(syear,smonth,1,0,0,0)
week2=num_sec_ep(syear,smonth,8,0,0,0)
result=where(base_time+time_offset ge week1 and $
             base_time+time_offset lt week2,count)
if count eq 0 then goto,skip_plot_1
week1=result[0]
week2=result[count-1]

week=image[week1:week2,*]
weektime=julian_day[week1:week2]
contour,week,weektime,height,xstyle=1,ystyle=1,/nodata
px=!x.window*!d.x_vsize
py=!y.window*!d.y_vsize
sx=px[1]-px[0]+1
sy=py[1]-py[0]+1
tv,congrid(week,sx,sy),px[0],py[0]
axis,weektime[0],height[0],xaxis=0,color=0,xstyle=1,xtickunits='time',$
  xtickformat='label_date',xrange=[weektime[0],weektime[n_elements(weektime)-1]]
axis,weektime[0],height[numheights-1],xaxis=1,color=0,xstyle=1,$
  xtickunits='time',xtickformat='(A1)',$
  xrange=[weektime[0],weektime[n_elements(weektime)-1]]
axis,weektime[0],height[0],yaxis=0,color=0,ystyle=1,$
  yrange=[height[0],height[numheights-1]],ytitle='Height (km)'
axis,weektime[n_elements(weektime)-1],height[0],yaxis=1,color=0,ystyle=1,$
  yrange=[height[0],height[numheights-1]],ytickformat='(A1)'
colorbarm,dataname[datanum],data,dmax,dmin,cb_title[datanum],px,py,sx,sy
skip_plot_1:
endif
;********
;  Plot 2
;********
plot2='yes'
if plot2 eq 'yes' then begin

week1=num_sec_ep(syear,smonth,8,0,0,0)
week2=num_sec_ep(syear,smonth,15,0,0,0)
result=where(base_time+time_offset ge week1 and $
             base_time+time_offset lt week2,count)
if count eq 0 then goto,skip_plot_2
week1=result[0]
week2=result[count-1]

week=image[week1:week2,*]
weektime=julian_day[week1:week2]
contour,week,weektime,height,xstyle=1,ystyle=1,/nodata
px=!x.window*!d.x_vsize
py=!y.window*!d.y_vsize
sx=px[1]-px[0]+1
sy=py[1]-py[0]+1
tv,congrid(week,sx,sy),px[0],py[0]
axis,weektime[0],height[0],xaxis=0,color=0,xstyle=1,xtickunits='time',$
  xtickformat='label_date',xrange=[weektime[0],weektime[n_elements(weektime)-1]]
axis,weektime[0],height[numheights-1],xaxis=1,color=0,xstyle=1,$
  xtickunits='time',xtickformat='(A1)',$
  xrange=[weektime[0],weektime[n_elements(weektime)-1]]
axis,weektime[0],height[0],yaxis=0,color=0,ystyle=1,$
  yrange=[height[0],height[numheights-1]],ytitle='Height (km)'
axis,weektime[n_elements(weektime)-1],height[0],yaxis=1,color=0,ystyle=1,$
  yrange=[height[0],height[numheights-1]],ytickformat='(A1)'
colorbarm,dataname[datanum],data,dmax,dmin,cb_title[datanum],px,py,sx,sy
skip_plot_2:
endif
;********
;  Plot 3
;********
plot3='yes'
if plot3 eq 'yes' then begin

week1=num_sec_ep(syear,smonth,15,0,0,0)
week2=num_sec_ep(syear,smonth,22,0,0,0)
result=where(base_time+time_offset ge week1 and $
             base_time+time_offset lt week2,count)
if count eq 0 then goto,skip_plot_3
week1=result[0]
week2=result[count-1]

week=image[week1:week2,*]
weektime=julian_day[week1:week2]
contour,week,weektime,height,xstyle=1,ystyle=1,/nodata
px=!x.window*!d.x_vsize
py=!y.window*!d.y_vsize
sx=px[1]-px[0]+1
sy=py[1]-py[0]+1
tv,congrid(week,sx,sy),px[0],py[0]
axis,weektime[0],height[0],xaxis=0,color=0,xstyle=1,xtickunits='time',$
  xtickformat='label_date',xrange=[weektime[0],weektime[n_elements(weektime)-1]]
axis,weektime[0],height[numheights-1],xaxis=1,color=0,xstyle=1,$
  xtickunits='time',xtickformat='(A1)',$
  xrange=[weektime[0],weektime[n_elements(weektime)-1]]
axis,weektime[0],height[0],yaxis=0,color=0,ystyle=1,$
  yrange=[height[0],height[numheights-1]],ytitle='Height (km)'
axis,weektime[n_elements(weektime)-1],height[0],yaxis=1,color=0,ystyle=1,$
  yrange=[height[0],height[numheights-1]],ytickformat='(A1)'
colorbarm,dataname[datanum],data,dmax,dmin,cb_title[datanum],px,py,sx,sy
skip_plot_3:
endif
;********
;  Plot 4
;********
plot4='yes'
if plot4 eq 'yes' then begin

week1=num_sec_ep(syear,smonth,22,0,0,0)
week2=num_sec_ep(syear,smonth,29,0,0,0)
result=where(base_time+time_offset ge week1 and $
             base_time+time_offset lt week2,count)
if count eq 0 then goto,skip_plot_4
week1=result[0]
week2=result[count-1]

week=image[week1:week2,*]
weektime=julian_day[week1:week2]
contour,week,weektime,height,xstyle=1,ystyle=1,/nodata
px=!x.window*!d.x_vsize
py=!y.window*!d.y_vsize
sx=px[1]-px[0]+1
sy=py[1]-py[0]+1
tv,congrid(week,sx,sy),px[0],py[0]
axis,weektime[0],height[0],xaxis=0,color=0,xstyle=1,xtickunits='time',$
  xtickformat='label_date',xrange=[weektime[0],weektime[n_elements(weektime)-1]]
axis,weektime[0],height[numheights-1],xaxis=1,color=0,xstyle=1,$
  xtickunits='time',xtickformat='(A1)',$
  xrange=[weektime[0],weektime[n_elements(weektime)-1]]
axis,weektime[0],height[0],yaxis=0,color=0,ystyle=1,$
  yrange=[height[0],height[numheights-1]],ytitle='Height (km)'
axis,weektime[n_elements(weektime)-1],height[0],yaxis=1,color=0,ystyle=1,$
  yrange=[height[0],height[numheights-1]],ytickformat='(A1)'
colorbarm,dataname[datanum],data,dmax,dmin,cb_title[datanum],px,py,sx,sy
skip_plot_4:
endif
;********
;  Plot 5
;********
plot5='yes'
if plot5 eq 'yes' then begin  

week1=num_sec_ep(syear,smonth,29,0,0,0)
week2=num_sec_ep(syear,smonth,31,23,59,60)
result=where(base_time+time_offset ge week1 and $
             base_time+time_offset lt week2,count)
if count eq 0 then goto,skip_plot_5
week1=result[0]
week2=result[count-1]

week=image[week1:week2,*]
weektime=julian_day[week1:week2]
contour,week,weektime,height,xstyle=1,ystyle=1,/nodata
px=!x.window*!d.x_vsize
py=!y.window*!d.y_vsize
sx=px[1]-px[0]+1
sy=py[1]-py[0]+1
if eday eq 31 then begin 
  adjust=3./7.  ;31 days
  numticks=3
endif else if eday eq 30 then begin 
  adjust=2./7.  ;30 days
  numticks=2
endif else if eday eq 29 then begin
  adjust=1./7.  ;29 days
  numticks=1
endif

tv,congrid(week,sx*adjust,sy),px[0],py[0]
sz=size(congrid(week,sx*adjust,sy))

contour,week,weektime,height,xstyle=1,ystyle=1,/nodata,/noerase,$
	position=[px[0],py[0],px[0]+sz[1]-1,py[0]+sz[2]-1],/device
tv,congrid(week,sx*adjust,sy),px[0],py[0]
axis,weektime[0],height[0],xaxis=0,color=0,xstyle=1,xtickunits='time',$
  xtickformat='label_date',$
  xrange=[weektime[0],weektime[n_elements(weektime)-1]],$
  xticks=numticks,xminor=9
axis,weektime[0],height[numheights-1],xaxis=1,color=0,xstyle=1,$
  xtickunits='time',$
  xrange=[weektime[0],weektime[n_elements(weektime)-1]],$
  xtickformat='(A1)',xticks=numticks,xminor=9
axis,weektime[0],height[0],yaxis=0,color=0,ystyle=1,$
  yrange=[height[0],height[numheights-1]],ytitle='Height (km)'
axis,weektime[n_elements(weektime)-1],height[0],yaxis=1,color=0,ystyle=1,$
  yrange=[height[0],height[numheights-1]],ytickformat='(A1)'
colorbarm,dataname[datanum],data,dmax,dmin,cb_title[datanum],px,py,sx,sy
skip_plot_5:
endif
;
;  Add a title
;
if smonth eq 1 then smon='Jan'
if smonth eq 2 then smon='Feb'
if smonth eq 3 then smon='Mar'
if smonth eq 4 then smon='Apr'
if smonth eq 5 then smon='May'
if smonth eq 6 then smon='Jun'
if smonth eq 7 then smon='Jul'
if smonth eq 8 then smon='Aug'
if smonth eq 9 then smon='Sep'
if smonth eq 10 then smon='Oct'
if smonth eq 11 then smon='Nov'
if smonth eq 12 then smon='Dec'
if site eq 'sgp' then begin
  capsite='SGP'
  datastream='average.'+avg_int
endif else if site eq 'nsa' then begin
  capsite='NSA'
  datastream='average.'+avg_int
endif else if site eq 'twpc1' then begin
  capsite='TWP  Manus'
  datastream='average.'+avg_int
endif else if site eq 'twpc2' then begin
  capsite='TWP  Nauru'
  datastream='average.'+avg_int
endif else if site eq 'crystal' then begin
  capsite='Crystal'
  datastream='average.'+avg_int
endif

syr=strcompress(string(syear),/remove_all)
if avg_int eq '1hr' then avg_int_str='Hourly'
if avg_int eq '5min' then avg_int_str='5 Minute'
xyouts,0.5,0.95,align=0.5,charsize=1.7,/normal,color=0,$
'Monthly Summary of '+ avg_int_str+' Averages for '+dataname[datanum]+'!C'+capsite+' - '+smon+' '+syr
;
;  Create gif file
;
strmon=string(smonth)
if smonth lt 10 then strmon='0'+strmon
strday=string(sday)
if sday lt 10 then strday='0'+strday
curdate=strcompress(syr+strmon+strday,/remove_all)

gifdir='/data/mace4/arm_data/temp_image/'
giffile=site+'.'+datastream+'.'+curdate+'.'+datatag[datanum]+'.gif'
print, gifdir+giffile
write_gif,gifdir+giffile,tvrd( )
;
;  Remove the .png file
;
png_file=site+'.'+datastream+'.'+curdate+'.'+datatag[datanum]+'.png'
unix_command='rm '+gifdir+png_file
spawn,unix_command
;
;  Skip to here to skip plotting if data is all missing
;
no_data_to_plot:

endfor  ;end loop through variables
end