library(raster)
library(dismo)
library(rgdal)
require(lpSolve)
require(tidyr)
library(deSolve)   
library(plyr)
library(parallel)
library(doParallel)
library(foreach)
library(randomForest)
require(ggplot2)
require(ggrepel)
require(rpart)
require(lme4)
require(MASS)
require(party)
require(tidyverse)


##2.a function to read the rainfall raster data of one month and give pixel values  
#' 
#' @param tifFiles  Daily precipitation from CHIRPs – UCSB (ftp://ftp.chg.ucsb.edu/pub/org/chg/products/CHIRPS-2.0/)


#' @param Nrdays 
#' @param month 
#' @param filename 
#' @param year 
#' @returnType 
#' @return 
#' 
#' @author Meklit
#' @export
getrasterData <- function(tifFiles, Nrdays=31, month, filename, year, p){
  mdata <- stack(tifFiles)
  #tid <- p$trialID
  p <- p[,c("x", "y")]
  month_data <- as.data.frame(raster::extract(mdata, p))
  month_data$long <- p$x
  month_data$lat <- p$y
  colnames(month_data) <- c(as.character(c(1:Nrdays)),"long", "lat") 
  month_data$Year <- year
  month_data$Month <- month
 # month_data$trialID <- tid
  return(month_data)
}


#' read one year rainfall data for the defined region 
#' @param year 
#' @param zone: the country or region name
#' @param coordPoints a two column data frame with c("long", "lat") columns
#' @param outputPath : to be defined by user and it will be where the output will be saved 
#' @return 
#' 
#' @author Meklit
#' @export
#' @example mdata <- NOT_AC_DAP[1, ]
extract_rf <- function(year,  mdata){
  # extract_rf <- function(year, leapYear = FALSE, coordPoints, country, outputPath, outputFormat){
  p <- mdata[, c("lon", "lat")]
  colnames(p) <- c("x", "y")
  leapYears <- seq(1984, 2024, 4)
  notleapYears <- seq(1980, 2024, 1)
  notleapYears <- notleapYears[!notleapYears%in% leapYears]
  if(year %in% notleapYears){
    nrdays <- 28 ## nr of days in Feb
    nrweeks <- 59
  }else if(year %in% leapYears){
    nrdays <- 29
    nrweeks <- 60
  }
  setwd("D:/Daily precipitation/RainfallTIF_1995_2016")
  
  setwd(as.character(year))
  
  suppressWarnings(rm(list_tif))
  
  list_tif <- list.files(getwd(), full.names = T, pattern = ".tif")
  if(nrweeks==59){
    jan_ds <- list_tif[1:31]
    Feb_ds <- list_tif[32:59]
    Mar_ds <- list_tif[60:90]
    Apr_ds <- list_tif[91:120]
    May_ds <- list_tif[121:151]
    Jun_ds <- list_tif[152:181]
    JUl_ds <- list_tif[182:212]
    Aug_ds <- list_tif[213:243]
    Sep_ds <- list_tif[244:273]
    Oct_ds <- list_tif[274:304]
    Nov_ds <- list_tif[305:334]
    Dec_ds <- list_tif[335:365]
    
  }else{
    jan_ds <- list_tif[1:31]
    Feb_ds <- list_tif[32:60]
    Mar_ds <- list_tif[61:91]
    Apr_ds <- list_tif[92:121]
    if(year !=2020){
      May_ds <- list_tif[122:152]
      Jun_ds <- list_tif[153:182]
      JUl_ds <- list_tif[183:213]
      Aug_ds <- list_tif[214:244]
      Sep_ds <- list_tif[245:274]
      Oct_ds <- list_tif[275:305]
      Nov_ds <- list_tif[306:335]
      Dec_ds <- list_tif[336:366]
    }
  }
  
  
  January_df <- getrasterData(tifFiles=jan_ds, Nrdays=31, month=01, filename="January", year=year, p=p)
  February_df <- getrasterData(tifFiles=Feb_ds, Nrdays=nrdays, month=02, filename="February", year=year,p) ## set the nr of days correct here
  March_df <- getrasterData(tifFiles=Mar_ds, Nrdays=31, month=03, filename="March", year=year,p)
  April_df <- getrasterData(tifFiles=Apr_ds, Nrdays=30, month=04, filename="April", year=year,p)
  if(year !=2020){
    May_df <- getrasterData(tifFiles=May_ds, Nrdays=31, month=05, filename="May", year=year,p)
    Jun_df <- getrasterData(tifFiles=Jun_ds, Nrdays=30, month=06, filename="Jun", year=year,p)
    July_df <- getrasterData(tifFiles=JUl_ds, Nrdays=31, month=07, filename="July", year=year,p)
    Aug_df <- getrasterData(tifFiles=Aug_ds, Nrdays=31, month=08, filename="Aug", year=year,p)
    Sep_df <- getrasterData(tifFiles=Sep_ds, Nrdays=30, month=09, filename="Sep", year=year,p)
    Oct_df <- getrasterData(tifFiles=Oct_ds, Nrdays=31, month=10, filename="Oct", year=year,p)
    Nov_df <- getrasterData(tifFiles=Nov_ds, Nrdays=30, month=11, filename="Nov", year=year,p)
    Dec_df <- getrasterData(tifFiles=Dec_ds, Nrdays=31, month=12, filename="Dec", year=year,p)
  }
  
 

  require(tidyr)
  January_long <- gather(January_df, Day, RAIN, as.character(c(1:31)), factor_key=TRUE)
  Feb_long <- gather(February_df, Day, RAIN, as.character(c(1:nrdays)), factor_key=TRUE)
  Mar_long <- gather(March_df, Day, RAIN, as.character(c(1:31)), factor_key=TRUE)
  Apr_long <- gather(April_df, Day, RAIN, as.character(c(1:30)), factor_key=TRUE)
  if(year !=2020){
    May_long <- gather(May_df, Day, RAIN, as.character(c(1:31)), factor_key=TRUE)
    Jun_long <- gather(Jun_df, Day, RAIN, as.character(c(1:30)), factor_key=TRUE)
    July_long <- gather(July_df, Day, RAIN, as.character(c(1:31)), factor_key=TRUE)
    Aug_long <- gather(Aug_df, Day, RAIN, as.character(c(1:31)), factor_key=TRUE)
    Sep_long <- gather(Sep_df, Day, RAIN, as.character(c(1:30)), factor_key=TRUE)
    Oct_long <- gather(Oct_df, Day, RAIN, as.character(c(1:31)), factor_key=TRUE)
    Nov_long <- gather(Nov_df, Day, RAIN, as.character(c(1:30)), factor_key=TRUE)
    Dec_long <- gather(Dec_df, Day, RAIN, as.character(c(1:31)), factor_key=TRUE)
  }

  
  if(nrdays==29){
    January_long$Date  <- as.numeric(January_long$Day)
    Feb_long$Date <- as.numeric(Feb_long$Day) +31 
    Mar_long$Date <- as.numeric(Mar_long$Day) +60
    Apr_long$Date <- as.numeric(Apr_long$Day) + 91
    if(year !=2020){
      May_long$Date <- as.numeric(May_long$Day) + 121
      Jun_long$Date <- as.numeric(Jun_long$Day) + 152
      July_long$Date <- as.numeric(July_long$Day) + 182
      Aug_long$Date <- as.numeric(Aug_long$Day) + 213
      Sep_long$Date <-as.numeric(Sep_long$Day) + 244
      Oct_long$Date <- as.numeric(Oct_long$Day) + 274
      Nov_long$Date <- as.numeric(Nov_long$Day) + 305
      Dec_long$Date <- as.numeric(Dec_long$Day) +335
    }
  }else{
    January_long$Date  <- as.numeric(January_long$Day)
    Feb_long$Date <- as.numeric(Feb_long$Day) +31 
    Mar_long$Date <- as.numeric(Mar_long$Day) +59
    Apr_long$Date <- as.numeric(Apr_long$Day) + 90
    if(year !=2020){
      May_long$Date <- as.numeric(May_long$Day) + 120
      Jun_long$Date <- as.numeric(Jun_long$Day) + 151
      July_long$Date <- as.numeric(July_long$Day) + 181
      Aug_long$Date <- as.numeric(Aug_long$Day) + 212
      Sep_long$Date <- as.numeric(Sep_long$Day) + 243
      Oct_long$Date <- as.numeric(Oct_long$Day) + 273
      Nov_long$Date <- as.numeric(Nov_long$Day) + 304
      Dec_long$Date <- as.numeric(Dec_long$Day) +334
    }
  }
  
   if(year==2020){
    listOfDataFrames <- list(January_long, Feb_long, Mar_long, Apr_long)
    
  }else{
    listOfDataFrames <- list(January_long, Feb_long, Mar_long, Apr_long, May_long, Jun_long, July_long, Aug_long, Sep_long, Oct_long,
                             Nov_long,Dec_long)
  }
 
  df <- do.call("rbind", listOfDataFrames)
  #  setwd(outputPath)
  # 
  # if(outputFormat == "csv"){
  #   fname <- paste(country, "_", year, "_Rainfall.csv", sep="")
  #   write.csv(df, file=fname)
  # }else{
  #   fname <- paste(country, "_", year, "_Rainfall.Rds", sep="")
  #   saveRDS(df, file=fname)
  # }
  return(df)
}


plotRF <- function(m, month){
  rf96P <- sum_96[sum_96$Month==m, ]
  rf99P <- sum_99[sum_99$Month==m, ]
  rf02P <- sum_02[sum_02$Month==m, ]
  
  Ton5 <- rf96P[, c("long", "lat", "MonthlyRain")]
  Ton599 <- rf99P[, c("long", "lat", "MonthlyRain")]
  Ton502 <- rf02P[, c("long", "lat", "MonthlyRain")]
  
  colnames(Ton5) <- c('X', 'Y','Z')
  colnames(Ton599) <- c('X', 'Y','Z')
  colnames(Ton502) <- c('X', 'Y','Z')
  
  Ton5 <- as.matrix(Ton5)
  Ton599 <- as.matrix(Ton599)
  Ton502 <- as.matrix(Ton502)
  
  extentTon5 <- extent(Ton5[,1:2])
  extentTon599 <- extent(Ton599[,1:2])
  extentTon502 <- extent(Ton502[,1:2])
  
  rasOct <- raster(extentTon5, ncol=100, nrow=100, crs='+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0')
  rasOct <- raster(extentTon599, ncol=100, nrow=100, crs='+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0')
  rasOct <- raster(extentTon502, ncol=100, nrow=100, crs='+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0')
  
  
  Nsoil_rast <- rasterize(Ton5[, 1:2], rasOct, Ton5[,3], fun=mean)
  Nsoil_rast99 <- rasterize(Ton599[, 1:2], rasOct, Ton599[,3], fun=mean)
  Nsoil_rast02 <- rasterize(Ton502[, 1:2], rasOct, Ton502[,3], fun=mean)
  
  plot(Nsoil_rast, main=paste(" 1996 monthly Rf", month, sep=" "))
  plot(Nsoil_rast99, main=paste(" 1999 monthly Rf", month, sep=" "))
  plot(Nsoil_rast02, main=paste(" 2002 monthly Rf", month, sep=" "))
  #plot(boundary_N, add=TRUE)
}
# 
# plotRF (m=1, month="Jan")



install.packages("chirps")
install_github("ropensci/chirps", build_vignettes = TRUE)
library("remotes")
library("chirps")

lonlat <- data.frame(lon = c(-55.0281,-54.9857, -55.0714),
                     lat = c(-2.8094, -2.8756, -3.5279))

dates <- c("2017-01-01", "2017-12-31")
dat <- get_chirps(lonlat, dates)
p_ind <- precip_indices(dat, timeseries = TRUE, intervals = 30)






##############################################################################################
## weather data from NASA:
##############################################################################################
## md is a dataframe with Lon and Lat column
#' Title: function to get data from NASA on daily min max temperature, wind speed, day tie radiation, vapour pressure for days as defined
#'
#' @param nrCluster: how many cores of your PC can be assigned to this task, has to be at least 2 
#' @param RangeRows: this function parallelize computation, if you have >1000 locations, better define it by 1000 coordinates in a run  
#' @param startDate: NASA data is requied for dates starting from this date 
#' @param endDate: NASA data is required for dates upto this date 
#' @param coordPoints: a two column data frame with c("long", "lat") columns  
#' @param fNames: if you want to write out the data toyour PC give a fil name for it, it will be saved as .RDS 
#' @param outputPath: to be defined by user and it will be where the output will be saved 
#'
#' @return a data frame with ("long","lat", "YEAR", "MONTH", "DATE", "DateYear", "YYMMDD", "TMAX", "TMIN", "WIND","DTR", "VP","location")
#' @examples: 
#' ##coord <- data.frame(lon=c(38.675, 38.825, 38.875), lat=rep(-4.775,3))
#' ## getNasaData(RangeRows = c(1:3), nrCluster = 8, startDate = "2011-01-01", endDate = "2011-12-31", coordPoints = coord, fNames = "NasaData_2011")
getNasaData <- function(nrCluster, RangeRows=c(1:1000), startDate, endDate, coordPoints, fNames=NULL, outputPath){
  library(nasapower)
  library(parallel)
  if(nrCluster > 1){
    cls = makePSOCKcluster(nrCluster)
    registerDoParallel(cls)
    
    if(length(RangeRows) > nrow(coordPoints) ){
      RangeRows <- c(1:nrow(coordPoints) )
    }
    
    NasaData_year <- foreach(i=RangeRows, .packages = c('plyr','nasapower'), .combine = "rbind") %dopar% {
      startDate <- startDate
      endDate <- endDate
      as.data.frame(nasapower::get_power(community = "AG",
                                         lonlat = c(coordPoints$lon[i], coordPoints$lat[i]),
                                         pars = c("T2M_MAX","T2M_MIN","WS2M","ALLSKY_SFC_SW_DWN"),
                                         dates = c(startDate, endDate),
                                         temporal_average = "DAILY"))
      
    }
    stopCluster(cls)
  }else{
    NasaData_year <- NULL
    for(i in 1:nrow(coordPoints)){
      NasaD <-  as.data.frame(nasapower::get_power(community = "AG",
                                                           lonlat = c(coordPoints$lon[i], coordPoints$lat[i]),
                                                           pars = c("T2M_MAX","T2M_MIN","WS2M","ALLSKY_SFC_SW_DWN"),
                                                           dates = c(startDate, endDate),
                                                           temporal_average = "DAILY"))
      NasaData_year <- rbind(NasaData_year, NasaD)
    }
  }


  colnames(NasaData_year) <- c("lon", "lat", "YEAR", "MONTH","DATE", "DateYear","YYMMDD","TMAX", "TMIN", "WIND", "DTR")
  NasaData_year$VP <- round((6.11*exp(17.4*NasaData_year$TMIN/(NasaData_year$WIND + 239)))/10, digits=2)
  NasaData_year$location <- paste(NasaData_year$lon, NasaData_year$lat, sep="_")
  length(unique(NasaData_year$location))
 
  if(!is.null(fNames)){
    setwd(outputPath)
    saveRDS(NasaData_year, paste(fNames, ".RDS", sep=""))
  }
  return(NasaData_year)
}



#' Title: ro run getNasaData function for all points in a country with 1000 coordiantes at a time
#'
#' @param nrCluster: how many cores of your PC can be assigned to this task, has to be at least 2 
#' @param coordPoints: a two column data frame with c("long", "lat") columns  
#' @param fNames: if you want to write out the data toyour PC give a fil name for it, it will be saved as .RDS 
#' @param outputPath: to be defined by user and it will be where the output will be saved 
#' @param year: the year for which NASA data is needed
#'
#' @return a data frame will be written out to the outputPath with ("long","lat", "YEAR", "MONTH", "DATE", "DateYear", "YYMMDD", "TMAX", "TMIN", "WIND","DTR", "VP","location")
#' @examples: NasaData_parallel(year=1996, coordPoints, 
#' outputPath = "D:/ACAI_Wrapper/cloud_compute/LINTUL/NG_data/NG_solar", nrCluster=8, country = "Nigeria)

NasaData_parallel <- function(year, coordPoints, outputPath, nrCluster, country){
  nrRuns <- data.frame(A=c(seq(1, (nrow(coordPoints)-1000),by=1000)), B=c(seq(1001, nrow(coordPoints), by=1000)))
  ND_Year <- NULL
  for(j in 1:nrow(nrRuns)){
    print(nrRuns[j,])
    lim1 <- nrRuns$A[j]
    if(j == nrow(nrRuns)){
      lim2 <- nrow(coordPoints)
    }else{
      lim2 <- nrRuns$B[j]
    }
    
    NDY <- getNasaData(RangeRows = c(lim1:lim2), nrCluster = nrCluster, 
                       startDate = paste(year, "-01-01", sep=""), endDate = paste(year,"-12-31", sep=""), 
                       coordPoints = coordPoints, fNames = paste("NasaData_", country, "_", year, LETTERS[j],sep=""), 
                       outputPath = outputPath)
    
     ND_Year <- rbind(ND_Year, NDY)
  }
  return(ND_Year)
}


################### get rainfall data, after copying the geotiff files
#' Title get Chirps rainfall data 
#' @param GPS_year  data with lat, lon, plantYear, harvestYear
#'
#' @return rainfall data by lat lon and from planting and harvest year
#' @examples
RF_GPS_year <- function(GPS_year){
  loc_years_rf <- NULL
  for(i in 1:nrow(GPS_year)){
    print(i)
    GPS_year_i <- GPS_year[i, ]
    if(GPS_year_i$plantYear == GPS_year_i$harvestYear){
      year <- GPS_year_i$plantYear
      RWerf <- extract_rf(year, GPS_year_i)
    }else{
      years <- seq(GPS_year_i$plantYear, GPS_year_i$harvestYear, 1)
      RWerf <- NULL
      for(k in unique(years)){
        year <- k
        erfk <- extract_rf(year, mdata=GPS_year_i)
        RWerf <- rbind(RWerf, erfk)
      }
    }
    loc_years_rf <- rbind(loc_years_rf, RWerf)
  }
  return(loc_years_rf)
}




### NOT AC data rainfall
setwd("D:/ACAI_Wrapper/Rwanda")
NOT_AC <- readRDS("NOT_formodel.RDS")
NOT_AC_DAP<- filter(NOT_AC, plantYear %in% c(2016,2017,2018,2019))
NOT_AC_DAP <- unique(NOT_AC_DAP[,c("lat", "lon", "plantYear", "harvestYear")])
str(NOT_AC_DAP)
NOT_AC_RF <- RF_GPS_year(GPS_year=NOT_AC_DAP)
NOT_AC_RF$location <- paste(NOT_AC_RF$lat, NOT_AC_RF$long, sep="_")
length(unique(NOT_AC_RF$location))
head(NOT_AC_RF)
str(NOT_AC_RF)

setwd("/home/akilimo/projects/SoilNPK")
saveRDS(NOT_AC_RF, "NOT_AC_Rainfall.RDS")

## CIALCA cassava
setwd("D:/ACAI_Wrapper/Rwanda/Coordinates")
CIALCA_cassava <- readRDS("CIALCA_isda.RDS")
CIALCA_cassava <- unique(CIALCA_cassava[, c("lat", "lon")])
names(CIALCA_cassava) <- c("lat", "lon")
str(CIALCA_cassava)
CIALCA_cassava$plantYear <- 2018
CIALCA_cassava$harvestYear <- 2020
CIALCA_NOT_RF <- RF_GPS_year(CIALCA_cassava)
setwd("/home/akilimo/projects/SoilNPK")
saveRDS(CIALCA_NOT_RF, "NOT_CIALCA_Rainfall.RDS")





## RW and BU RF
setwd("D:/ACAI_Wrapper/Rwanda/Coordinates")
RWCoor <- read.csv("Rwanda_Coordinates.csv")
BUCoor <- read.csv("Burundi_Coordinates.csv")
RW_BU <- rbind(RWCoor, BUCoor)
RW_BU <- unique(RW_BU[, c("lat", "long")])
names(RW_BU) <- c("lat", "lon")
str(RW_BU) ##1697


DRCCoor <- read.csv("DRC_Coordinates.csv")
DRCCoor <- unique(DRCCoor[, c("lat", "long")])
names(DRCCoor) <- c("lat", "lon")



##2010 to 2020 rainfall data

for(years in 2010:2020){
  RW_BU$plantYear <- years
  RW_BU$harvestYear <- years
  RWBU_RF_years <- RF_GPS_year(RW_BU)
  setwd("/home/akilimo/projects/SoilNPK/dataSource/rainfall")
  saveRDS(RWBU_RF_years, paste("RWBU_RF_", years,".RDS",sep = ""))
}

RWBU_RF_2010 <- readRDS("RWBU_RF_2010.RDS")


RF_GPS_year_par <- function(rf_Coor, years){
  rf_Coor$plantYear <- years
  rf_Coor$harvestYear <- years
  cls = makePSOCKcluster(7)
  registerDoParallel(cls)
  g=0
  print(g)
  rf_Coor_year <- foreach(i=c(1:nrow(rf_Coor)), .packages = c('plyr'), .combine = "rbind") %dopar% {
    g <- g+i
    setwd("D:/ACAI_Wrapper/Rwanda")
    source("RF_NASA_Functions.R")
    extract_rf(year = years, mdata = rf_Coor[i, ])
  }
  stopCluster(cls)
  return(rf_Coor_year)
}


## charles
setwd("D:/ODKdata/Charles")
charlesgps <- read.csv("charlesgps.csv")
CGPS <- unique(charlesgps[, c("Latitude", "Longitude")])
names(CGPS) <- c("lat", "lon")
str(CGPS) 
Charles_2017 <- RF_GPS_year_par(rf_Coor=CGPS, years=2017)
Charles_2018 <- RF_GPS_year_par(rf_Coor=CGPS, years=2018)
Charles_2019 <- RF_GPS_year_par(rf_Coor=CGPS, years=2019)

CharlesRF <- rbind(Charles_2017, Charles_2018, Charles_2019)
head(CharlesRF)
setwd("D:/ODKdata/Charles")
write.csv(CharlesRF, "Rainfall_2018_2019.csv", row.names = FALSE)
##################


CI_RF_2011 <- RF_GPS_year_par(rf_Coor=RW_BU, years=2011)
str(CI_RF_2011)
setwd("D:/ACAI_Wrapper/Rwanda/Rainfall")
saveRDS(CI_RF_2011, "RWBU_RF_2011.RDS")


CI_RF_2012 <- RF_GPS_year_par(rf_Coor=RW_BU, years=2012)
str(CI_RF_2012)
setwd("D:/ACAI_Wrapper/Rwanda/Rainfall")
saveRDS(CI_RF_2012, "RWBU_RF_2012.RDS")

CI_RF_2013 <- RF_GPS_year_par(rf_Coor=RW_BU, years=2013)
str(CI_RF_2013)
setwd("D:/ACAI_Wrapper/Rwanda/Rainfall")
saveRDS(CI_RF_2013, "RWBU_RF_2013.RDS")

CI_RF_2014 <- RF_GPS_year_par(rf_Coor=RW_BU, years=2014)
str(CI_RF_2014)
setwd("D:/ACAI_Wrapper/Rwanda/Rainfall")
saveRDS(CI_RF_2014, "RWBU_RF_2014.RDS")


CI_RF_2015 <- RF_GPS_year_par(rf_Coor=RW_BU, years=2015)
str(CI_RF_2015)
setwd("D:/ACAI_Wrapper/Rwanda/Rainfall")
saveRDS(CI_RF_2015, "RWBU_RF_2015.RDS")

CI_RF_2016 <- RF_GPS_year_par(rf_Coor=RW_BU, years=2016)
str(CI_RF_2016)
setwd("D:/ACAI_Wrapper/Rwanda/Rainfall")
saveRDS(CI_RF_2016, "RWBU_RF_2016.RDS")

CI_RF_2017 <- RF_GPS_year_par(rf_Coor=RW_BU, years=2017)
str(CI_RF_2017)
setwd("D:/ACAI_Wrapper/Rwanda/Rainfall")
saveRDS(CI_RF_2017, "RWBU_RF_2017.RDS")

CI_RF_2018 <- RF_GPS_year_par(rf_Coor=RW_BU, years=2018)
str(CI_RF_2018)
setwd("D:/ACAI_Wrapper/Rwanda/Rainfall")
saveRDS(CI_RF_2018, "RWBU_RF_2018.RDS")

CI_RF_2019 <- RF_GPS_year_par(rf_Coor=RW_BU, years=2019)
str(CI_RF_2019)
setwd("D:/ACAI_Wrapper/Rwanda/Rainfall")
saveRDS(CI_RF_2019, "RWBU_RF_2019.RDS")



#### check years 2010 to 2019 and see which three years would be selected

setwd("D:/ACAI_Wrapper/Rwanda/Coordinates")
RWCoor <- read.csv("Rwanda_Coordinates.csv")
BUCoor <- read.csv("Burundi_Coordinates.csv")
BUCoor$country <- "Bu"
RWCoor$country <- "Rw"
RW_BU <- rbind(RWCoor, BUCoor)
RW_BU$location <- paste(RW_BU$long, RW_BU$lat, sep="_")

setwd("D:/ACAI_Wrapper/Rwanda/Rainfall")

CI_RF_cum <- NULL
for(years in 2010:2019){
  CI_RF <- readRDS(paste("RWBU_RF_", years,".RDS", sep=""))
  CI_RF$location <- paste(CI_RF$long, CI_RF$lat, sep="_")
  CI_RF <- unique(merge(CI_RF, RW_BU[, c("location", "country")], by="location"))
  avRAin <- ddply(CI_RF, .(Date, country), summarise, avRain = mean(RAIN))
  avRAin$year <- years
  avRAinRW <- avRAin[avRAin$country == "Rw", ]
  avRAinBu <- avRAin[avRAin$country == "Bu", ]
  avRAinRW$cumRain <- cumsum(avRAinRW$avRain)
  avRAinBu$cumRain <- cumsum(avRAinBu$avRain)
  
  
  CI_RF_cum <- rbind(CI_RF_cum, avRAinRW, avRAinBu)
}

head(CI_RF_cum)

saveRDS(CI_RF_cum, "RwBu_cum_Rainfall.RDS")
CI_RF_cum <- readRDS("RwBu_cum_Rainfall.RDS")

ggrain <- ggplot(CI_RF_cum, aes(Date, cumRain, col=factor(year)))+
  geom_line(size=2)+
  facet_grid(.~country)+
  theme_bw()

RWtotal <- CI_RF_cum[CI_RF_cum$country == "Rw" & CI_RF_cum$Date == 365, ]
Butotal <- CI_RF_cum[CI_RF_cum$country == "Bu" & CI_RF_cum$Date == 365, ]
RWtotal[order(RWtotal$cumRain),]
Butotal[order(Butotal$cumRain),]

quantile(RWtotal$cumRain, probs=seq(0,1,0.05) ) 
quantile(Butotal$cumRain, probs=seq(0,1,0.05) ) 

rbind(data.frame(country="Rwanda", Q3Year = 2010, Q2Year = 2019, Q1Year = 2016),
      data.frame(country="Burundi", Q3Year = 2015, Q2Year = 2011, Q1Year = 2016))

######################################################################33
## NASA
######################################################################33
##RW: Q1, Q2, Q3 = 2016, 2019, 2010
## Bu: Q1, Q2, Q3 = 2016, 2011, 2015

## RWBU
setwd("/home/akilimo/projects/SoilNPK/Coordinates")
RWCoor <- read.csv("Rwanda_Coordinates.csv")
BUCoor <- read.csv("Burundi_Coordinates.csv")
DRCCoor <- read.csv("DRC_Coordinates.csv")

RW_BU <- rbind(RWCoor, BUCoor)
RWCoor <- unique(RWCoor[, c("long","lat")])
str(RWCoor) ##828

BUCoor <- unique(BUCoor[, c("long","lat")])
str(BUCoor) ##869

DRCCoor <- unique(DRCCoor[, c("long","lat")])
str(DRCCoor) ##2094


setwd("D:/ACAI_Wrapper/Rwanda")
source("RF_NASA_Functions.R")
NasaData_parallel(year=2010, coordPoints=RWCoor, outputPath = "D:/ACAI_Wrapper/Rwanda/solar", nrCluster=6, country = "Rwanda")
NasaData_parallel(year=2019, coordPoints=RWCoor, outputPath = "D:/ACAI_Wrapper/Rwanda/solar", nrCluster=6, country = "Rwanda")
NasaData_parallel(year=2016, coordPoints=RWCoor, outputPath = "D:/ACAI_Wrapper/Rwanda/solar", nrCluster=6, country = "Rwanda")

NasaData_parallel(year=2015, coordPoints=BUCoor, outputPath = "D:/ACAI_Wrapper/Rwanda/solar", nrCluster=6, country = "Burundi")
NasaData_parallel(year=2011, coordPoints=BUCoor, outputPath = "D:/ACAI_Wrapper/Rwanda/solar", nrCluster=6, country = "Burundi")
NasaData_parallel(year=2016, coordPoints=BUCoor, outputPath = "D:/ACAI_Wrapper/Rwanda/solar", nrCluster=6, country = "Burundi")



setwd("D:/ACAI_Wrapper/Rwanda/solar")
NasaData_Rwanda_2010 <- readRDS("NasaData_Rwanda_2010.RDS")
head(NasaData_Rwanda_2010)
length(unique(NasaData_Rwanda_2010$location))


#### next step is running LINTUL