;objective;to map the difference between 5% and 0,2.5,10,20,30
;Jan 2018
;Heba Marey
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro horcolorbar, plotcolors, n_colors, levels, bartop, fgcolor, fmtstr, units

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

  colorleft = xsize * 0.15
  colorright = xsize * 0.95
  colorlen = (colorright - colorleft) / plotcolors
  colortop = ysize * bartop
  colorbottom = ysize * (bartop - .04)
  labely = ysize * (bartop - .06 )

  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, /device

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

  endfor

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

end
;---------------------------------------------------
;SAVE,nTR,ndata,Test,lonbin,latbin,filename='grid_Aug23data.sav'
;SAVE,nTR,ndata,cdata,Test,nblon,nblat,lonbin,latbin,$,
;cl_data,retv_it,co_eror,filename='grid_Aug23TestMopRadDay.sav'

;filename='grid_Aug23data.sav'
filename='grid_Aug23dataMODIS.sav'

;filename='grid_Aug23TestMopRadDay.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=9 ;without 5
Test=['0.0','2.5','5.0','10.0','20.0','30.0','50.0','70.0','90.0','99.9']
Testd=['0.0','2.5','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,1 do begin  ;because we will not diff 5 as it is the default
  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)-ndata(ilon,ilat,2)
endfor
endfor
endfor
;;;;;;;;;;;;;;
for i=2,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+1) gt 0. then diff(ilon,ilat,i)=ndata(ilon,ilat,i+1)-ndata(ilon,ilat,2)
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, 161, 0, 0,   0,228 ,  0, 144, 255, 255, 255, 255, 255, 255])
green = byte([0,255, 179, 70,104, 204,239,255, 255, 255, 191, 102,   0,   0, 103])
blue  = byte([0,255, 223, 255,255, 255,255,0, 130,   0,   0,   0,   0, 255, 255])
tvlct,red,green,blue

!p.background=1
!p.color=0
levels=[-20,-15,-10,-5,-1,0,1,5,10,15,20,30]

fmtstr='(i3)'
  units = ''
!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

ve=Testd(itest)
epsfile='ModisDiff/map_cloudAug23diffMODIS'+ve+'.eps'
;epsfile='test.eps'

set_plot,'ps'
 print,'ps'
  ;device,/landscape,/color, file =epsfile,/encapsulated
device, file =epsfile,/bold,/color,/encapsulated,xsize=60,xoffset=0.2,/PORTRAIT,ysize=40

  !p.font=0

position=[0.15,0.15,0.95,0.9]
maptitle = 'August 23, 2000  MOPITT pixel diference number for cloud cover '+ve+'%'

cgmap_set,0,0,/continent,/Cylindrical,/advance, $
  latdel=10,londel=10,/NoErase,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')
; Write out the map title here
 x1map=0.5
 y1map=0.95
 xyouts, x1map, y1map, maptitle, CHARSIZE=3,  /NORMAL, $
  CHARTHICK=4.0,alignment=0.5
nlevels=n_elements(levels)
sf=1
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 diff[ilon,ilat,itest] gt -999. then begin

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

 dat1 = diff[ilon,ilat,itest]

 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, $
  latdel=30,londel=60,/NoErase,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=4.0
cgDrawShapes, shapefile, Thick=3, AttrName='NAME' ;to have canada boundaries

map_grid,latdel=20,londel=20,/box_axes, CHARSIZE=2.5
ncol=n_elements(levels)
horcolorbar, ncol, ncol+2, levels, 0.1, 0, fmtstr, units
print,epsfile

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