;objective;to map the difference between 5% and 0,2.5,10,20,30
;Jan 2018
;Heba Marey
;modified March 2018 to add chi2
;also modified to map all tests in one figure
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro horcolorbar, plotcolors, n_colors, levels, bartop, fgcolor, fmtstr, units,cr,cl,s

  xsize = !D.x_size
  ysize = !D.y_size

  colorleft = xsize * cl
  colorright = xsize * cr
  colorlen = (colorright - colorleft) / plotcolors
  colortop = ysize * bartop
  colorbottom = ysize * (bartop - 0.04) ;width of the bars
  labely = ysize * (bartop - 0.03)

  for i = 2, n_colors-1 do begin
    colleft = colorleft + colorlen*(i-2)
    colright = colleft + colorlen

    polyfill, [colleft,colright,colright,colleft], $
              [colorbottom,colorbottom,colortop,colortop], $
              color = i+s-2, /device

    ;levels goes from 0 to n_colors - 3, so adjust accordingly (subtract 2)
    xyouts, colleft, labely, string(levels[i-2], format=fmtstr), $
          charsize = 1, alignment = 0.0, color = fgcolor, /device

  endfor

  xyouts, colorright, labely, units, $
          charsize = 1, align=0.5, color = fgcolor, /device

end
;---------------------------------------------------
pro map_all,sf, levels,lonbin,latbin,data,position,title,fmtstr,bt,cr,cl,units 
z=1.0
cgmap_set,0,0,/continent,/Cylindrical,/advance, title=title+'!C',CHARSIZE=z-0.4,$
 latdel=10,londel=10,/NoErase,position=position;[0.15,0.15,0.95,0.9]

;map_set,0,0.5*(lonmin+lonmax),/continent,/CYLINDRICAL,/advance, $
 ; latdel=10,londel=10,/NoErase,limit = [latmin,lonmin,latmax,lonmax],position=position;,/isotropi

shapefile = Filepath(subdir=['resource','maps', 'shape'], 'canadaprovince.shp')

lobin = lonbin[1]-lonbin[0]
labin = latbin[1]-latbin[0]

s=2
for ilon=0,n_elements(lonbin)-1 do begin
 for ilat=0,n_elements(latbin)-1 do begin
  if data[ilon,ilat] gt -9999. then begin

   lat1 = latbin[ilat]
   lon1 = lonbin[ilon]
   lat2 = lat1 + labin
   lon2 = lon1 + lobin

 dat1 = data[ilon,ilat]

 icol = Max(where(levels*sf lt dat1)) + s > s
;print,dat1,' icol',icol
   polyfill,[lon1,lon1,lon2,lon2],[lat1,lat2,lat2,lat1],col=icol,noclip=0
   endif
  endfor
endfor

cgmap_set,0,0,/continent,/Cylindrical,/advance,CHARSIZE=z-0.4, $
 /NoErase,position=position;[0.15,0.15,0.95,0.9]
;map_set,0,0.5*(lonmin+lonmax),/continent,/CYLINDRICAL,/advance, $
 ; latdel=10,londel=10,/NoErase,limit = [latmin,lonmin,latmax,lonmax],position=position;,/isotropic


cgMap_Continents, /countries, /usa,CHARTHICK=1
cgDrawShapes, shapefile, Thick=1, AttrName='NAME' ;to have canada boundaries

map_grid,latdel=60,londel=60,/box_axes, CHARSIZE=z-0.3
ncol=n_elements(levels)
horcolorbar, ncol, ncol+2, levels, bt, 0, fmtstr, units,cr,cl,s

end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;SAVE,nTR,ndata,Test,lonbin,latbin,filename='grid_Aug23data.sav'
;filename='grid_Aug23data.sav'
;restore,filename
;stdata=ndata(*,*,2)  ;data of 5% and 0.95

;SAVE,nTR,ndata,cdata,Test,nblon,nblat,lonbin,latbin,$,
;dgFr,retv_it,co_eror,qaiSeq,filename='grid_Aug23TestMopRadDay.sav'
;filename='grid_Aug23TestMopRad.sav'

filename='grid_Aug23dataMODIS.sav'
;filename='grid_Aug23TestMopRadNight.sav'
restore,filename
ndata95=ndata(*,*,2)
co_eror95=co_eror(*,*,2)

filename='grid_Aug23TestMopRad.sav'

restore,filename

a=4
;;;;;;;;;;;;;;;;;;;;;
if a eq 0 then begin
area='AF'
latmin =-40
latmax = 40
lonmin = -20.
lonmax = 40.0
endif
;;;;;;;;;;;;;;;;;
;;;;;;;;;
if a eq 2 then begin

area='EU'
latmin =30
latmax =  70
lonmin = -20.
lonmax = 60.0
end

;;;;;;;;;;;;;;;;;;;;;
if a eq 1 then begin
area="NA"

latmin =10
latmax =  65
lonmin = -170
lonmax = -40
end

;;;;;;;;;;;;;;;;;;;;;
if a eq 3 then begin

area='AS'
latmin =5.0
latmax =70
lonmin = 30.
lonmax = 150.0
end

;;;;;;;;;;;;;;;;;;;;;
if a eq 4 then begin

area='AU'
latmin =-45
latmax =0 
lonmin = 110.
lonmax = 160
end

;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;
if a eq 5 then begin

area='SA'
latmin =-60
latmax =20
lonmin = -80.
lonmax = -30
end

;;;;;;;;;;;;;;;;;;;;;



ntest=11
nt=12 ;with 0.95

;Testd=['0.5','0.7','0.9','1.0','1.2','0.5_99.9']
;Testd=['0.01-modis-99.9','0.01','0.1','0.5','0.5-modis-99.9','0.7','0.8','0.8-1.2','0.9','1.0','1.2']
 Testd=['0.01-99.9%','0.01','0.1','0.5','0.5-99.9%','0.7','0.8','0.8-1.2','0.9','1.0','1.2']
Test=['0.01-modis-99.9','0.01','0.1','0.5','0.5-modis-99.9','0.7','0.8','0.8-1.2','0.9','0.95','1.0','1.2']

diff=fltarr(nblon,nblat,ntest)-99999.
diff95=fltarr(nblon,nblat)-99999.
diffall=fltarr(nblon,nblat,nt)-99999.
 ;there will be 2-2 so we will 

;to calculate average of 0% then subtract it from average for each test each pixel
co_eror0=co_eror(*,*,9)
in0=where(co_eror0 gt 0)

co_eror0=co_eror0(in0)
ReMean0=float(Mean(co_eror0))

;;;;;;;;;;;;;;
for i=0,ntest-1 do begin
  for ilon=0,n_elements(lonbin)-1 do begin
 for ilat=0,n_elements(latbin)-1 do begin

co_erorr=co_eror(ilon,ilat,i)
in=where(co_erorr gt 0)

co_ero=co_erorr(in)
ReMean=float(Mean(co_ero))
  if ndata(ilon,ilat,i) gt 0. then diff(ilon,ilat,i)=((ReMean-ReMean0)/ReMean0)*100
;if ndata(ilon,ilat,i) gt 0. then diff(ilon,ilat,i)=(chi2Mean-chi2Mean0)

endfor
endfor
endfor
;***************to add 0.95 to the total array*********

for ilon=0,n_elements(lonbin)-1 do begin
 for ilat=0,n_elements(latbin)-1 do begin
c95=co_eror95(ilon,ilat)
in95=where(c95 gt 0)

co_ero95=c95(in95)
ReMean95=float(Mean(co_ero95))
  if ndata95(ilon,ilat) gt 0. then diffall(ilon,ilat,9)=((ReMean95-ReMean0)/ReMean0)*100
;if ndata(ilon,ilat,i) gt 0. then diff(ilon,ilat,i)=(chi2Mean-chi2Mean0)

endfor
endfor

   diffall(*,*,0:8)=diff(*,*,0:8)
  diffall(*,*,10:11)=diff(*,*,9:10)

;**************

;***********************************
;red   = byte([0,255, 228, 161,   0,   0,   0, 144, 255, 255, 255, 255, 255, 255])
;green = byte([0,255, 239, 179, 104, 204, 255, 255, 255, 191, 102,   0,   0, 103])
;blue  = byte([0,255, 255, 223, 255, 255,   0, 130,   0,   0,   0,   0, 255, 255])

;red   = byte([0,255,228, 161,  0,   0,0,  0, 144, 255, 255, 255, 255, 255, 255])
;green = byte([0,255, 239,179,104, 204,0,255, 255, 255, 191, 102,   0,   0, 103])
;blue  = byte([0,255, 255,223, 255, 255,255,0, 130,   0,   0,   0,   0, 255, 255])

red   = byte([0,255,228, 161, 0, 0,   0,  0, 144, 255, 255, 255, 255, 255, 255])
green = byte([0,255, 239,179, 70,104, 204,255, 255, 255, 191, 102,   0,   0, 103])
blue  = byte([0,255, 255,223, 255,255, 255,0, 130,   0,   0,   0,   0, 255, 255])
tvlct,red,green,blue

!p.background=1
!p.color=0

!p.charsize=1.2
;;;;;;;;;;;;;;;;;;;;;;mapping;;;;;;;;;;;;;;;;;;;;
;!P.Multi = [0, 3, 2]
;pos = cgLayout([2, 5], OXMargin=[50, 12], OYMargin=[12, 12], XGap=30, YGap=15)

epsfile='ParamMopRad/map_Aug23_diffMopRad_Re.eps'
;epsfile='test_ch.eps'
set_plot,'ps'
 print,epsfile
  ;device,/landscape,/color, file =epsfile,/encapsulated
;device, file =epsfile,/bold,/color,/encapsulated,xsize=70,xoffset=0.2,/PORTRAIT,ysize=60
device, file =epsfile,/bold,/color,/encapsulated,xsize=30,xoffset=0.1,/PORTRAIT,ysize=20
  !p.font=0

for itest=0,nt-1 do begin  ;to map the diference of all tests each in seperate file

pos=cgLayout([4, 3], OXMargin=[8, 8], OYMargin=[15, 6], XGap=6, YGap=8)

ve=Test(itest)

position=pos[*,itest]
;position=[0.15,0.15,0.95,0.9]
maptitle = 'MOPITT relative Re% for each MOPITT rad. ratio test'

;maptitle = 'MOPITT Retrieval error for each MODIS cloud cover test'
;levels=[5,10,20,25,30,35,40]
sf=1.e15
fmtstr='(f7.2)'
levels=[0.1,1,10,20,30,40,50,100]
;fmtstr='(i2)'
;data=co_eror(*,*,itest)
title=ve+'%'
;bt=0.68
;cr=0.9
;cl=0.1
units='1.e15'

bt=0.06
cr=0.9
cl=0.1
;map_all,sf, levels,lonbin,latbin,data,position,title,fmtstr,bt,cr,cl,units

; Write out the map title here
 x1map=0.5
 y1map=0.97

 xyouts, x1map, y1map, maptitle, CHARSIZE=1.,  /NORMAL, $
  CHARTHICK=2.0,alignment=0.5
;********************************************
;position=pos[*,1]
sf=1
fmtstr='(f5.2)'
title='Degree Freedom'
levels=[0.2,0.4,0.6,0.8,1,1.2,1.4]
;data=dgFr(*,*,itest)
bt=0.37
cr=0.9
cl=0.1
;map_all,sf, levels,lonbin,latbin,data,position,title,fmtstr,bt,cr,cl
;********************************************
;position=pos[*,2]
;data=retv_it(*,*,itest)
bt=0.06
cr=0.9
cl=0.1
fmtstr='(i2)'
levels=[1,2,3,5,6,8,10]
title='Retrieval iteration'
;map_all,sf, levels,lonbin,latbin,data,position,title,fmtstr,bt,cr,cl

;**********************************************************
position=pos[*,itest]
data=diffall(*,*,itest)
;sf=1.0e01
sf=1
bt=0.07
cr=0.9
cl=0.1
fmtstr='(f10.2)'
levels=[-100,-50,-10,0,20,40,70,100,200,300,500]
;levels=[-100,0,10]

title=ve
units = '%'
map_all,sf, levels,lonbin,latbin,data,position,title,fmtstr,bt,cr,cl,units 



endfor   ;close the loop of tests
device,/close
end