;
; Plots calibration statistics for the MXD against temperature.
; Also infills the few grid boxes without calibration coefficients.
;
restore,filename='calibmxd1.idlsave'
; Gets:  g,mxdyear,mxdnyr,fdcorr,fdalph,fdbeta,fdvexp,fdcalib,mxdfd2,$
;        fdrver,fdvver,timey
;
; Prepare for plotting
;
loadct,39
multi_plot,nrow=2
if !d.name eq 'X' then begin
  window,ysize=850
  !p.font=-1
endif else begin
  !p.font=0
  device,/helvetica,/bold,font_size=9
endelse
;
outfd=total(finite(mxdfd2),3)
list4=where(outfd eq 0)
list2=where(outfd gt 0)
outfd(list4)=4.
outfd(list2)=2.
;
map=def_map(/npolar) & map.limit(0)=25.
labels=def_labels(/off)
;
inter_boxfd,fdcorr,g.x,g.y,$
  labels=labels,map=map,$
  levels=[0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9],/scale,$
  title='Correlation over calibration period'
fdin={ fd: outfd, x: g.x, y: g.y, nx: g.nx, ny: g.ny }
whizz_fd,fdin=fdin,fdout=f,limit=map.limit
boxplot,f.fd,f.x,f.y,/overplot,/overmap,highlight=f.fd,thick=1.5
;
inter_boxfd,fdrver,g.x,g.y,$
  labels=labels,map=map,$
  levels=[0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9],/scale,$
  title='Correlation over verification period'
boxplot,f.fd,f.x,f.y,/overplot,/overmap,highlight=f.fd,thick=1.5
;
pause
inter_boxfd,fdalph,g.x,g.y,$
  labels=labels,map=map,$
  levels=0.1*[-1.1,-0.4,-0.2,-0.1,0.,0.1,0.2,0.4,1.1],/scale,$
  title='Regression constant'
boxplot,f.fd,f.x,f.y,/overplot,/overmap,highlight=f.fd,thick=1.5
;
inter_boxfd,fdbeta,g.x,g.y,$
  labels=labels,map=map,$
  levels=[-1.,0.2,0.4,0.6,0.8,1.0,1.2,1.4,1.6],/scale,$
  title='Regression slope'
boxplot,f.fd,f.x,f.y,/overplot,/overmap,highlight=f.fd,thick=1.5
;
pause
inter_boxfd,fdvexp,g.x,g.y,$
  labels=labels,map=map,$
  levels=[-1,0,0.1,0.2,0.3,0.4,0.5,0.6,0.7],/scale,$
  title='% variance explained (calibration)'
boxplot,f.fd,f.x,f.y,/overplot,/overmap,highlight=f.fd,thick=1.5
;
inter_boxfd,fdvver,g.x,g.y,$
  labels=labels,map=map,$
  levels=[-1.1,0,0.1,0.2,0.3,0.4,0.5,0.6,0.7],/scale,$
  title='% variance explained (verification)'
boxplot,f.fd,f.x,f.y,/overplot,/overmap,highlight=f.fd,thick=1.5
;
; Now need to infill the 6 boxes without regression coefficients.  This is
; done using a gaussian-weighted mean of nearby coefficients.
;
print,'Infilling coefficients'
allx=fltarr(g.nx,g.ny)
for iy = 0 , g.ny-1 do allx(*,iy)=g.x(*)
ally=fltarr(g.nx,g.ny)
for ix = 0 , g.nx-1 do ally(ix,*)=g.y(*)
statx=allx(list2)
staty=ally(list2)
;
statval=fdalph(list2)
misslist=where(finite(statval) eq 0,nmiss)
if nmiss gt 0 then begin
  fd_extend,statval,statx,staty,search=1000.,wavelen=400.
  fdalph(list2(misslist))=statval(misslist)
endif
;
statval=fdbeta(list2)
misslist=where(finite(statval) eq 0,nmiss)
if nmiss gt 0 then begin
  fd_extend,statval,statx,staty,search=1000.,wavelen=400.
  fdbeta(list2(misslist))=statval(misslist)
endif
;
; Now plot the infilled coefficients
;
pause
inter_boxfd,fdalph,g.x,g.y,$
  labels=labels,map=map,$
  levels=0.1*[-1.1,-0.4,-0.2,-0.1,0.,0.1,0.2,0.4,1.1],/scale,$
  title='Regression constant (infilled)'
boxplot,f.fd,f.x,f.y,/overplot,/overmap,highlight=f.fd,thick=1.5
;
inter_boxfd,fdbeta,g.x,g.y,$
  labels=labels,map=map,$
  levels=[-1.,0.2,0.4,0.6,0.8,1.0,1.2,1.4,1.6],/scale,$
  title='Regression slope (infilled)'
boxplot,f.fd,f.x,f.y,/overplot,/overmap,highlight=f.fd,thick=1.5
;
; Now compute the calibrated values now that all boxes have coefficients
;
for iyr = 0 , mxdnyr-1 do begin
  fdcalib(*,*,iyr)=fdalph(*,*)+fdbeta(*,*)*mxdfd2(*,*,iyr)
endfor
;
; Now save the data for later analysis
;
save,filename='calibmxd2.idlsave',$
  g,mxdyear,mxdnyr,fdcorr,fdalph,fdbeta,fdvexp,fdcalib,mxdfd2,$
  fdrver,fdvver,timey,fdseas
;
end