;objective;to map the difference between 5% and 0,2.5,10,20,30
;Jan 2018
;Heba Marey
;modified March 2018
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
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.01)
  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 = 0.5, 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
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=6
for ilon=0,n_elements(lonbin)-1 do begin
 for ilat=0,n_elements(latbin)-1 do begin
  if data[ilon,ilat] gt 0. 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
   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
units = ''
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,filename='grid_Aug23TestMopRadDay.sav'
;filename='grid_Aug23TestMopRad.sav'

filename='grid_Aug23dataMODIS.sav'
;filename='grid_Aug23TestMopRadNight.sav'

restore,filename
;nTR Level 1 data that gives the total number of data
;ndata has all test data

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=10

;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.0','2.5','5.0','10.0','20.0','30.0','50.0','70.0','90.0','99.9']



diff=fltarr(nblon,nblat,ntest)-999.
 ;there will be 2-2 so we will 

;;;;;;;;;;;;;;
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

  ;if ndata(ilon,ilat,i) gt 0. then diff(ilon,ilat,i)=ndata(ilon,ilat,i)-stdata(ilon,ilat)
endfor
endfor
endfor

;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, 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)

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

pos=cgLayout([2, 2], OXMargin=[8, 8], OYMargin=[6, 6], XGap=6, YGap=8)

ve=Testd(itest)
epsfile='ParamMODIS/map_cloudAug23_diffMODIS'+ve+'_V2eps'
;epsfile='test.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=20,xoffset=0.1,/PORTRAIT,ysize=16
  !p.font=0

position=pos[*,0]
;position=[0.15,0.15,0.95,0.9]
maptitle = 'MODIS cloud cover % test '+ve

levels=[5,10,20,25,30,35,40]
sf=1.e16
fmtstr='(i2)'
data=co_eror(*,*,itest)
title='Retrieval error'
bt=0.68
cr=0.9
cl=0.1
map_all,sf, levels,lonbin,latbin,data,position,title,fmtstr,bt,cr,cl

; 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='(f3.1)'
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[*,3]
data=qaiSeq(*,*,itest)
sf=1.e02
bt=0.06
cr=0.9
cl=0.1
fmtstr='(i2)'
levels=[1,5,10,15,20,25,30,35,40]
title='Chi2'
map_all,sf, levels,lonbin,latbin,data,position,title,fmtstr,bt,cr,cl



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