Evaluate and produce in an operational way seasonal hydrological forecasts for the Niger basin by forcing the NIGER-HYPE model with ECMWF SF system 5 data.

README.md

Evaluate and produce in an operational way seasonal hydrological forecasts for the Niger basin by forcing the NIGER-HYPE model with ECMWF SF system 5 data.

The objective is to evaluate seasonal hydrological forecasts over the 1993-2015 period and to produce seasonal hydrological forecasts on an operational basis at the scale of the Niger basin.

To use these scripts, you need to:

1. Install the software R: https://cran.r-project.org/bin/windows/base/;

2. Install the necessary packages (rgdal, raster, ncdf4);

analys_SF_syst5_Annual_time_step.R

### Script to calculate and plot annual parameters of ECMWF SF system5
### Authors: Bernard Minoungou, Jafet Andersson, Abdou Ali, Mohamed Hamatan


  # -----------
  # 1 - General
  #Remove workspace
  rm(list=ls())
  library(raster)
  library(rgdal)
  require(lattice)
  #--------------
  #Get Niger Sub-bassin
  subs<-readOGR(dsn="C:/Niger_2.22/GIS", layer="Current_Niger_onedelta_subbasins_2012-10-17")
  subs1<-subs
  #--------------
  #Magnification for plots
  cexx<-1.5
  #--------------
  #Define SF file and outputs dir
  months<-vector()
  for (i in 1:12) {
    if (i<=9) {
      months<-c(months, paste("0", i, sep=""))
    } 
    else {
      months<-c(months, paste(i))
    }
  }

  datadir<-paste("I:/ecmwf_syst5_hype_format/D03_Niger_month", months, sep="")
  outdir<-paste("I:/analysis_ecmwf_syst5/D03_Niger_month",months, sep="")
  #--------------
  #Define SF ensembles index
  ens<-c(0:24)
  enss<-vector()
  for (i in 1:25) {
    if (i<=9) {
      enss<-c(enss, paste("00", i, sep=""))
    } 
    else {
      enss<-c(enss, paste("0", i, sep=""))
    }
  }
  #--------------
  #Define variables names, variables attributes, years, plots lab, units
  var.name<-c("PObs", "Tobs", "TMAXobs", "TMINobs")
  var.long.name<-c("Seasonal Precipitation", "Seasonal Mean Temperature", 
                   "Seasonal Max Temperature", "Seasonal Min Temperature")
  units<-c("(mm)", rep("(°C)",3))
  var.num<-c(228, 55, 51, 52)
  years<-c(1993:2015)
  rm(i)
  # ----------------
  # 2 - Calculations and plots
  for (i in 1:length(months)) {
    for (j in 1:length(enss)) { 
      for (k in 1:length(years)) {
        for (l in 1:length(var.num)) {
          # 2.1 - Calculations
          fdir<-paste(datadir[i], years[k], sep="/")
          fname<-paste(var.name[l],"_",enss[j], ".txt", sep="")
          current.file<-read.table(paste(fdir, fname, sep="/"),h=T, sep="\t")
          if (var.name[l]=="PObs") {
            var.mean<-colSums(current.file[,-1])
          } else {
            var.mean<-colMeans(current.file[,-1])
          }
          a1<-slot(subs1, "data")
          a1<-cbind(a1,var.mean)
          slot(subs, "data")<-a1
          classes<-pretty(var.mean,100)
          my.colfun<-colorRampPalette(c("red","yellow","green","deepskyblue","blue"))
          if (j==1) {
            dir.create(outdir[i])
          }
          if( k==1) {
            dir.create(paste(outdir[i], "/ens", enss[j], sep="")) 
          }
          main1<-paste(var.long.name[l], years[k], "Ensemble", enss[j], units[l], sep=" ")
          # ----------------
          # 3.2. Maps (sub-basins)
          png(paste(outdir[i], "/ens", enss[j], "/",
                    paste("map", "_",var.name[l], "_",years[k],".png", sep=""),sep=""),
                    800,800)
          trellis.par.set(par.main.text=list(cex=cexx))
          print(
            spplot(subs, zcol="var.mean", 
                   at=classes, 
                   col.regions=my.colfun(length(classes)-1),
                   main=main1,
                   colorkey=list(labels=list(cex=cexx)),
                   scales=list(draw = TRUE))
          )
          dev.off()
          print(main1)
        }
      }
    }
  }

bias_correct_sf_system5_data_4_niger_hype.R

###---------------------------------------------------------------------------
###Script to bias correct the ECMWF SF System 5 with the WFDEI reanalys data  
###Authors: Bernard Minoungou, Abdou Ali, Mohamed Hamatan                     
###---------------------------------------------------------------------------

#-----------
#1 - General 
  
  #-------------
  # 1.1. Clear work directory and load packages
    rm(list=ls())
    library(raster)
    library(rgdal)
    require(lattice)
    library(qmap)

  #--------------
  # 1.2. Magnification for plots
    cexx<-1.5
    
    
  #--------------
  # 1.3. Define SF system 5 dir
    months<-vector()
    for (i in 1:12) {
      if (i<=9) {
        months<-c(months, paste("0", i, sep=""))
      } 
      else {
        months<-c(months, paste(i))
      }
    }
  
    datadir<-paste("I:/ecmwf_syst5_hype_format/D03_Niger_month", months, sep="")
    
  #-------------- 
  #1.5. Défine corrected data, bias parameters and Q-Q plot dir
    biaisdir<-paste("I:/ecmwf_syst5_hype_format/Biais/D03_Niger_month", 
                    months, sep="")
    plotdir<-paste("I:/analysis_ecmwf_syst5/Biais/plot/D03_Niger_month",
                   months, sep="")
    biasparamdir<-paste("I:/analysis_ecmwf_syst5/Biais/Bias_parameters/D03_Niger_month",
                        months, sep="")
  #--------------
  #1.6. Define SF ensembles index
    ens<-c(0:24)
    enss<-vector()
    for (i in 1:25) {
      if (i<=9) {
        enss<-c(enss, paste("00", i, sep=""))
      } 
      else {
        enss<-c(enss, paste("0", i, sep=""))
      }
    }
  #--------------
  # 1.7. Define variables names, variables attributes and years
    var.name<-c("PObs", "Tobs", "TMAXobs", "TMINobs")
    var.name<-c("PObs", "Tobs", "TMAXobs", "TMINobs")
    var.long.name<-c("Precipitation", "Mean Temperature", 
                 "Max Temperature", "Min Temperature")
    units<-c("[mm]", rep("[°C]",3))
    var.num<-c(228, 55, 51, 52)
    years<-c(1993:2015)
    rm(i)

#-----------
#2 - Get WFDEI Data (Pobs, TObs, TMAXObs, TMINObs) 
  wfdei<-list()
  pcp<-read.table("C:/Niger_2.22/Pobs.txt", sep="\t",h=T)
  tobs<-read.table("C:/Niger_2.22/Tobs.txt", sep="\t",h=T)
  tmaxobs<-read.table("C:/Niger_2.22/TMAXobs.txt", sep="\t",h=T)
  tminobs<-read.table("C:/Niger_2.22/TMINobs.txt", sep="\t",h=T)
  wfdei[[1]]<-pcp
  wfdei[[2]]<-tobs
  wfdei[[3]]<-tmaxobs
  wfdei[[4]]<-tminobs
  
#-----------
#3 - Bias correction and save corrected data, bias parameters and Q-Q Plots
for (i in 1:length(months)) {
  dir.create(biaisdir[i]) # Bias output dir
  dir.create(plotdir[i]) # Bias correction QQ Plot dir by initialization month 
  dir.create(biasparamdir[i]) # Bias parameter sub-dir by initialization month
  for (j in 1:length(var.num)) {
    for (k in 1:length(enss)) {
      value.ens.all.year<-NULL
      #-----------
      #3.1 - loop through years
      for (l in 1:length(years)) {
        if (j==2&k==1) {
          dir.create(paste(biaisdir[i], years[l], sep="/"))
        }
        fdir<-paste(datadir[i], years[l], sep="/")
        fname<-paste(var.name[j],"_",enss[k], ".txt", sep="")
        current.file<-read.table(paste(fdir, fname, sep="/"),h=T, sep="\t")
        value.ens.all.year<-rbind(value.ens.all.year,current.file)
      }
      date.common<-intersect(value.ens.all.year[,1],wfdei[[j]][,1])
      mymatch1<-match(date.common,value.ens.all.year[,1])
      mymatch2<-match(date.common,wfdei[[j]][,1])
      value.ens.intersect<-value.ens.all.year[mymatch1,]
      value.wfdei.intersect<-wfdei[[j]][mymatch2,]
      
      #3.2 - Compute bias parameter for precipitation data
      if (j==1) {
        current.fit<-fitQmap(obs = value.wfdei.intersect[,-1],
                             mod =value.ens.intersect[,-1],
                             method ="QUANT", 
                             type="linear", 
                             wet.day=TRUE)
      } else {
        #3.3 - Compute bias parameter for temperature data
        current.fit<-fitQmap(obs = value.wfdei.intersect[,-1],
                             mod =value.ens.intersect[,-1],
                             method ="QUANT", 
                             type="linear")
      }
      
      
      #3.4 - Apply bias correction to historical data
      value.ens.corr<-round(doQmap(value.ens.intersect[,-1], current.fit),2)
      
      
      ###3.5 Save bias parameter for operationnal correction
      save(current.fit, 
           file=paste(biasparamdir[i],"/",var.name[j], "_ens_", enss[k], ".RData",sep=""))
      
      ###3.6 Plot quantile-quantile map to evaluate the correction
      png(paste(plotdir[i], "/QQplot_", var.name[j], "_ens_", enss[k], ".png", sep=""),  800,800)
      xymin<-min(min(rowMeans(value.wfdei.intersect[,-1])),
                min(rowMeans(value.ens.corr, na.rm=T)),
                min(rowMeans(value.ens.intersect[,-1], na.rm=T)))
      xymax<-max(max(rowMeans(value.wfdei.intersect[,-1])),
                max(rowMeans(value.ens.corr, na.rm=T)),
                max(rowMeans(value.ens.intersect[,-1], na.rm=T)))
      qqplot(rowMeans(value.wfdei.intersect[,-1]), 
             rowMeans(value.ens.corr, na.rm=T),
             xlim=c(xymin,xymax), ylim=c(xymin,xymax),
             col="green",
             xlab=paste("Observed",  var.long.name[j], units[j], sep=" "), 
             ylab=paste("Forecasted",  var.long.name[j], units[j], sep=" "), 
             main=paste("Q-Q Plot", var.long.name[j], "Ensemble", enss[k], sep=" "),
             cex.lab=cexx,cex.main=cexx,cex.axis=cexx
             )
      points(sort(rowMeans(value.wfdei.intersect[,-1], na.rm=T)), 
             sort(rowMeans(value.ens.intersect[,-1], na.rm=T)), 
             col = "red")
      abline(0,1, col="grey", lwd=2)
      legend("topleft", c("Corrected", "Uncorrected"), 
             col=c("green", "red"), pch=c(1,1), cex=1.5)
      dev.off()
      
      ###3.7 Save corrected data year by year
      value.ens.all.year.cor<-round(doQmap(value.ens.all.year[,-1], current.fit),2)
      value.ens.all.year.cor<-cbind(value.ens.all.year[,1],value.ens.all.year.cor)
      colnames(value.ens.all.year.cor)<-c("DATE", sub("X","",colnames(value.ens.all.year)[-1]))
      for (m in 1:length(years)) {
        current.save.file<-value.ens.all.year.cor[which(substr(value.ens.all.year.cor[,1],1,4)==years[m]),]
        write.table(current.save.file,
                    paste(biaisdir[i], years[m], paste(var.name[j], "_", enss[k], ".txt", sep=""), sep="/"),
                    sep="\t", col.names=T, row.names=F, quote=F)
      }
      print(paste(months[i],"_", var.name[j], "_ens_", enss[k], sep=""))
    }
  }
}

extract_sf_system5_in_niger_hype_format.R

# Script to to extract ECMWF Seasonal forecast system 5 into Niger-HYPE V 2.22 format 
# Authors: Bernard Minoungou, Abdou Ali, Mohamed Hamatan

#--------------
# 1 - General
  rm(list=ls())
  # Load librairies
  library(rgdal)
  library(raster)
  library(ncdf4)
  # Get Niger Sub-bassin shapefile and create centroids
  subs<-readOGR(dsn="C:/Niger_2.22/GIS", layer="Current_Niger_onedelta_subbasins_2012-10-17")
  centroid<-getSpPPolygonsLabptSlots(subs)
  centroid<-as.data.frame(centroid)
  colnames(centroid)<-c("x","y")
  centroid<-as.data.frame(centroid)
  coordinates(centroid)<-~x+y
  projection(centroid)<-CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
  
  
#-------------- 
# 2 - Conversion
  #--------------
  #Define SF file dir
  months<-vector()
  for (i in 1:12) {
      if (i<=9) {
      months<-c(months, paste("0", i, sep=""))
    } 
    else {
      months<-c(months, paste(i))
    }
  }
  
  datadir<-paste("I:/ECMWF_system5/D03_Niger_month", months, 
                 "/nobackup/smhid12/sm_jorro/GLORIOUS/Data/PreProcessSF/ECMWFSyst5/D03_Niger/m",
                 months, sep="")
  #--------------
  #Define SF ensembles index
  ens<-c(0:24)
  enss<-vector()
  for (i in 1:25) {
    if (i<=9) {
      enss<-c(enss, paste("00", i, sep=""))
    } 
    else {
      enss<-c(enss, paste("0", i, sep=""))
    }
  }
  #--------------
  #Define variables names, variables attributes and years
  var.name<-c("PObs", "Tobs", "TMAXobs", "TMINobs")
  var.num<-c(228, 55, 51, 52)
  years<-c(1993:2015)
  rm(i)
  #--------------
  #Get Data and save file into Niger-HYPE format
  for (i in 1:length(months)) {
    for (j in 1:length(ens)) {
      for (l in 1:length(years)) {
       for (k in 1:length(var.num)) {
          fname<-paste(var.num[k], "_D03_Niger_ECMWF-sys5_season_03deg_orig_reg_day_r", 
                       ens[j], "_", years[l],"-",months[i], ".nc",sep="")
          rr<-brick(paste(datadir[i], fname, sep="/"))
      
          #--------------
          #Extract and convert cumilated precipitation data into daily precipitation
          if (var.name[k]=="PObs") {
            get_Data.cum<-t(round(extract(rr, centroid)*1000, 5))
            get_Data<-matrix(NA, nrow=nrow(get_Data.cum), ncol=ncol(get_Data.cum))
            get_Data[1,]<-get_Data.cum[1,]
            for (m in 2:nrow(get_Data)) {
              get_Data[m,]<-round((get_Data.cum[m,]-get_Data.cum[m-1,]),2)
            }
            get_Data[which(get_Data<0)]<-0
            colnames(get_Data)<-slot(subs, "data")[, "SUBID"]
            Date<-paste(substr(rownames(get_Data.cum),2,5), substr(rownames(get_Data.cum),7,8),
                        substr(rownames(get_Data.cum),10,11), sep="-")
          } else {
            #--------------
            #Extract Temperature data and convert from Kelvin into °C
            get_Data<-t(round((extract(rr, centroid)-273.15), 2))
            colnames(get_Data)<-slot(subs, "data")[, "SUBID"]
            Date<-paste(substr(rownames(get_Data),2,5), substr(rownames(get_Data),7,8),
                        substr(rownames(get_Data),10,11), sep="-")
          }
          get_Data<-cbind(as.character(Date),get_Data)
          colnames(get_Data)[1]<-"DATE"
          if (k==1) {
            dir.create(paste("I:/ecmwf_syst5_hype_format/", "D03_Niger_month", months[i], sep=""))
          }
          dir.create(paste("I:/ecmwf_syst5_hype_format/", "D03_Niger_month", months[i], "/",years[l],sep=""))
          file.name<-paste(var.name[k], "_", enss[j], ".txt", sep="")
          current_file.dir<-paste("I:/ecmwf_syst5_hype_format/", "D03_Niger_month", months[i], "/",years[l],sep="")
          #--------------
          #Save data into Niger-HYPE format
          write.table(get_Data, paste(current_file.dir, file.name, sep="/"), col.names=T,
                      row.names=F, quote=F, sep="\t")
          process.progress<-paste(months[i], ens[j], years[l], var.name[k], sep="_")
          print(process.progress)
        }
      }
    }
  }

force_hype_with_system5_data.R

### Script to force HYPE with SF system 5 data
### Authors: Bernard Minoungou, Abdou Ali, Mohamed Hamatan

# -----------
# 1 - General
#Remove workspace and set work directory
  rm(list=ls())
  setwd("C:/Niger_2.22")
  
  #--------------
  #Define SF ensembles index
  ens<-c(0:24)
  enss<-vector()
  for (i in 1:25) {
    if (i<=9) {
      enss<-c(enss, paste("00", i, sep=""))
    } 
    else {
      enss<-c(enss, paste("0", i, sep=""))
    }
  }

#--------------
#Define SF file dir and output dir
months<-vector()
for (i in 1:12) {
  if (i<=9) {
    months<-c(months, paste("0", i, sep=""))
  } 
  else {
    months<-c(months, paste(i))
  }
}
datadir<-paste("J:/ecmwf_syst5_hype_format/Biais/D03_Niger_month", months, sep="")
outdir<-paste("J:/ecmwf_syst5_hype_format/sortie_hype_system5/D03_Niger_month",months, sep="")
rm(i)

#--------------
# 2 Simulate automatically NIGER-HYPE over years and ensembles
years<-c(1993:2013)
for (i in 12:length(months)){
  dir.create(outdir[i], recursive=T)
  for (j in 1:length(years)){
    print(years[j])
    dir.create(paste(outdir[i], years[j],sep="/"))
    datee<-as.Date(c(1:215), origin=paste(years[j], months[i], "01", sep="-"))
    list.of.files <- list.files(paste(datadir[i], years[j], sep="/"), all.files=F)
    file.copy(paste(datadir[i], years[j],list.of.files, sep="/"),"C:/Niger_2.22",overwrite = TRUE)
    info<-readLines(paste(getwd(), "info.txt", sep="/"))
    info[6]<-paste("resultdir","\t","\t",  
                   paste("J:\\ecmwf_syst5_hype_format\\sortie_hype_system5\\D03_Niger_month", months[i], "\\", years[j], "\\", sep=""), 
                   "\t!!Windows",
                   sep="")
    info[12]<-paste("bdate","\t","\t",datee[1],sep="")
    info[13]<-paste("cdate","\t","\t",datee[1],sep="")
    info[14]<-paste("edate","\t","\t",datee[length(datee)],sep="")
    writeLines(info, "info.txt")
    for (k in 1:length(enss)) {
      path<-paste("C:/Niger_2.22/HYPE_5_0_0.exe -s", enss[k], sep=" ")
      system(path, wait=T)
      print(path)
    }
  }
}