load "./shapefile_utils.ncl"
;----------------------------------------------------------;

begin

prefix  = get_script_prefix_name()

;-- Output file
fnameo  = "Gridded_daily_ratio.nc"
system("rm " + fnameo)                          ; remove any pre-existing file
out     = addfile(fnameo,"c")                   ; open netCDF file

;-- Set grid size
nlat = 1800
nlon = 3600

lon = fspan(-179.9,180.0,nlon)
lon@units = "degrees_east"
lat = fspan(-89.95,89.95,nlat)
lat@units = "degrees_north"

;-- Grid for masking
dgrid = new((/nlat,nlon/),integer)
dgrid!0 = "lat"
dgrid!1 = "lon"
dgrid&lat = lat
dgrid&lon = lon
dgrid@long_name = "Region mask"
dgrid@_FillValue = -9999

dgrid   = 1

;-- Select smaller area to make the code faster
min_lat = 0.
max_lat = 25.
min_lon =-20.
max_lon = 0.

dgrid_fs = dgrid({min_lat:max_lat},{min_lon:max_lon})  
printVarSummary(dgrid_fs)
print("-------------------------------------------------")

;-- Read ascii file containing daily trends for the different regions in Senegal 
trend_dir  = "./"
ftrend     = "trend_by_region_Senegal.csv"

delimiter  = ";"                        ; field delimiter

;-- Read in data as strings.
hdata      = asciiread(trend_dir+ftrend,-1,"string")

header     = hdata(0)                   ; Header. Use for variable names.
read_hdata = hdata(1:)                  ; Get rid of first line which is a header.
nrows      = dimsizes(read_hdata)       ; Number of rows.

nfields    = str_fields_count(read_hdata(0),delimiter)
print("nfields = " + nfields)
print("nrows = " + nrows)

;-- Region names
reg_names  = new(nfields-1,string)
do i=1,nfields-1
   reg_names(i-1) = str_get_field(header,i+1,delimiter)
end do
;print(reg_names)

;-- Time dimension
date_str = str_get_field(read_hdata,1,delimiter)
ntime    = dimsizes(date_str)
newtime  = new(ntime,integer)
do i=0,ntime-1
   day = toint(str_get_field(date_str(i), 1, "/"))
   mon = toint(str_get_field(date_str(i), 2, "/"))
   year= 2000+toint(str_get_field(date_str(i), 3, "/"))

   ;convert to YYYYMMDD. Use "%02i" to force a leading '0' for day and month
   newtime(i) = toint(year + sprinti("%02i",mon) + sprinti("%02i",day))
end do 
;print(date_str + "= "+newtime)

;-- Daily trends at individual region 
vtrend    = new((/nrows,nfields-1/),float)
do j=1,nfields-1
   vtrend(:,j-1) = stringtofloat(str_get_field(read_hdata,j+1,delimiter))
end do
   
vtrend!0    = "time"
vtrend&time = newtime
vtrend!1    = "city_name" 
vtrend&city_name  = reg_names
vtrend@_FillValue = -9999
printVarSummary(vtrend)
print("-------------------------------------------------")

;-- Open shapefile 
dirShapf = "./"
shp_name = "gadm36_1.shp"
shp_path = dirShapf+shp_name
print_shapefile_info(shp_path)
;plot_shapefile(shp_path)

;-- Set masking options for shapefile_mask_data function
; and loop over the region
do ir=0,dimsizes(reg_names)-1
   
   vname = reg_names(ir)
   print("region = "+vname)

   opt             = True     
   opt@shape_var   = "NAME_1" 
   opt@shape_names = (/vname/) 
   opt@return_mask = True    
   opt@debug       = True    
  
  ;-- Mask the data
   data            = (/shapefile_mask_data(dgrid_fs,shp_path,opt)/)
   copy_VarMeta(dgrid_fs, data) 
   printVarSummary(data)
   ;printMinMax(data,True)

  ;-- Put trend values to the areas that we don't want masked
  ;get the size of the lat/lon grid of variable data
   iNumlat = dimsizes(data&lat)
   iNumlon = dimsizes(data&lon)
   
   ; Create a new variable (float) to put in the daily trends for each region 
   ; and fill the other grids with the value 1.
   data_mask = new((/iNumlat,iNumlon/),float)
   copy_VarMeta(data, data_mask)

   ;data_mask = 0

   do j=0,iNumlat-1
      do i=0,iNumlon-1
        if ((.not.ismissing(data(j,i))).and.(data(j,i) .eq. 1)) then
           data_mask(j,i) = vtrend(0,ir)/100.
        end if
      end do
   end do

end do   ; end looping over the region
printVarSummary(data_mask)
printMinMax(data_mask, True)
print("-------------------------------------------------")

;-------------------- Plot-----------------------------------------------------
out_dir = "./"

map_name = prefix                              
outfile  = out_dir + map_name

wks_type = "ps"

plot = new(2,graphic)

wks = gsn_open_wks(wks_type, outfile)         

colmaps = (/"WhBlGrYeRe","ViBlGrWhYeOrRe"/)                           
;gsn_define_colormap(wks, colmap)               ; définition de la table de couleur

mres                       = True                 ;-- plot mods desired
mres@gsnDraw               = False                ;-- don't draw the plot yet
mres@gsnFrame              = False                ;-- don't advance
mres@gsnAddCyclic          = False                ;-- add cyclic point
mres@gsnMaximize           = True  

mres@cnFillPalette         = colmaps(0)
mres@cnFillOn              = True                 ;-- color Fill
mres@cnFillMode            = "RasterFill"         ;-- faster
mres@cnLinesOn             = False                ;-- Turn off contour lines
mres@pmTickMarkDisplayMode = "Always"         
mres@cnLevelSelectionMode  = "ManualLevels"       

mres@cnMinLevelValF       =  0.
mres@cnMaxLevelValF       =  1.
mres@cnLevelSpacingF      =  0.1
mres@cnConstFEnableFill   = True
mres@cnConstFLabelOn      = False

;-- A décommanter pour ajouter des grilles 
mres@mpPerimOn             = True                 ;-- encadrer la figure

;-- We want county outlines.
mres@mpDataSetName         = "Earth..4"
mres@mpDataBaseVersion     = "MediumRes"         ;-- set map data base
mres@mpOutlineBoundarySets = "AllBoundaries"
mres@mpNationalLineThicknessF= 0.6
mres@mpOutlineOn           = True
;mres@mpOutlineDrawOrder    = "PostDraw"

;-- Focus on west Africa
mres@mpMinLatF             = min_lat               ;-- lat min. value
mres@mpMaxLatF             = max_lat               ;-- lat max. value
mres@mpMinLonF             = min_lon               ;-- lon min. value
mres@mpMaxLonF             = max_lon              ;-- lon max. value

;-- Ajout de la barre de légende
;mres@lbBoxLinesOn          = False
mres@lbLabelBarOn          = True                ; don't draw a labelbar below each plot
mres@lbOrientation         = "Vertical"          ; vertical label bar
mres@lbLabelFontHeightF    = 0.012               ; 
;mres@lbBoxMinorExtentF     = 0.16                ; decrease height of labelbar
mres@pmLabelBarWidthF      = 0.04                ; augmenter/diminuer la longueur label box
mres@pmLabelBarHeightF     = 0.5
mres@lbPerimOn             = False 

mres@tmXBLabelFontHeightF  = 0.012               ; change maj lat tm spacing
mres@tmYLLabelFontHeightF  = 0.012               ; change maj lon tm spacing
mres@tmXBMajorLengthF      = 0.01                ; change the tickmark length

mres@gsnLeftString       = "Mask"
mres@gsnStringFontHeightF  = 0.022
;mres@gsnLeftStringOrthogonalPosF = .010  

plot(0) = gsn_csm_contour_map(wks, data, mres)

mres@cnFillPalette        = colmaps(1)
mres@cnMinLevelValF       =  -1.
mres@cnMaxLevelValF       =   1.
mres@cnLevelSpacingF      =  0.2
mres@gsnLeftString        = "Gridded daily weights"

plot(1) = gsn_csm_contour_map(wks, data_mask, mres)   

;---Panel both plots
pres                    = True
pres@gsnMaximize        = True
pres@gsnPanelLabelBar   = False

gsn_panel(wks,(/plot(0),plot(1)/),(/2,1/),pres)

;---------------------------------------------------------------------
; output
;---------------------------------------------------------------------
;-- Save the variable in a netcdf file

  delete(data)
  delete(data_mask)

end