###############################################################################################
###############################################################################################
## accesseing via API: AIzaSyB7WZEpto9Jf1Do1SJBF7B8le7yVrEWN_U
###############################################################################################
###############################################################################################

library(httr)
library(jsonlite)
res = GET("https://api.isda-africa.com/v1/layers?key=AIzaSyB7WZEpto9Jf1Do1SJBF7B8le7yVrEWN_U")
layers <- fromJSON(rawToChar(res$content))
layers$property$bedrock_depth

###############################################################################################
## the soil variables
###############################################################################################
## properties for 0-20 and 20-50 cm depth
Property_2depth <- c("aluminium_extractable", "bulk_density", "calcium_extractable","carbon_organic", "carbon_total","clay_content",
                  "cation_exchange_capacity","iron_extractable","magnesium_extractable", "nitrogen_total","phosphorous_extractable",
                  "potassium_extractable", "sand_content", "silt_content", "stone_content", "sulphur_extractable", "texture_class",
                  "zinc_extractable", "ph") 

## properties with no depth
cropCovers <- c("crop_cover_2015", "crop_cover_2016", "crop_cover_2017", "crop_cover_2018","crop_cover_2019")##no depth
landCovers <- c("land_cover_2015","land_cover_2016", "land_cover_2017","land_cover_2018","land_cover_2019")
slopeAngle <- "slope_angle"

## properties with 0 - 200 cm depth
bedRock <- "bedrock_depth"

## properties with 0-50 cm depth
FCC <- "fcc" ##0-50 cm FertilityCapabilityClassification

coordPoints <- data.frame(lat=-11.375, lon=31.525)
coordPoints <- data.frame(lat=-2.425, lon=28.875)

###############################################################################################
## functions to get iSDA data via API: some variables have two depth and some have no depth info or only 21 depth. 
###############################################################################################

pH = GET("https://api.isda-africa.com/v1/soilproperty?key=AIzaSyB7WZEpto9Jf1Do1SJBF7B8le7yVrEWN_U&lat=-11.375&lon=31.525&property=ph&depth=0-20")
ph_latlong = fromJSON(rawToChar(pH$content))
ph_latlong$property$ph$value$value
ph_latlong$property$ph$depth$value
ph_latlong$property$ph$depth$unit
ph_latlong$property$ph$uncertainty[[1]]


## properties for 0-20 and 20-50 cm depth
Property_2depth <- c("aluminium_extractable", "bulk_density", "calcium_extractable","carbon_organic", "carbon_total","clay_content",
                     "cation_exchange_capacity","iron_extractable","magnesium_extractable", "nitrogen_total","phosphorous_extractable",
                     "potassium_extractable", "sand_content", "silt_content", "stone_content", "sulphur_extractable", "texture_class",
                     "zinc_extractable", "ph") 

## properties with no depth
cropCovers <- c("crop_cover_2015", "crop_cover_2016", "crop_cover_2017", "crop_cover_2018","crop_cover_2019")##no depth
landCovers <- c("land_cover_2015","land_cover_2016", "land_cover_2017","land_cover_2018","land_cover_2019")
slopeAngle <- "slope_angle"

## properties with 0 - 200 cm depth
bedRock <- "bedrock_depth"

## properties with 0-50 cm depth
FCC <- "fcc" ##0-50 cm FertilityCapabilityClassification




#' @description a function to get soil data from iSDA 30m res for layers given for 0-20 and 0-50 cm depth
#' @param coordPoints a data frame with lat and lon column
#' @param soilPropAll a list of soil variable with two depth
#' @param depth a vector with "0-20" or "0-50"
#' @return a data frame with lat, lon, Variable, depth, depthUnit, value and lower and upper limits for 50, 68 and
getiSDA_data2Depth <- function(coordPoints, soilPropAll, depth){
  all_SoilProp <- NULL
  for(j in 1:length(soilPropAll)){
    soilProp <- soilPropAll[j]
    print(soilProp)
    soilVarData <- NULL
    for(i in 1:nrow(coordPoints)){
      print(i)
      lat <- coordPoints$lat[i]
      lon <- coordPoints$lon[i]
      url1 <- "https://api.isda-africa.com/v1/soilproperty?key=AIzaSyB7WZEpto9Jf1Do1SJBF7B8le7yVrEWN_U&"
      coord <- paste("lat=", lat, "&lon=", lon, sep="")
      url2 <- paste("&property=",soilProp, "&depth=", depth, sep="")
      apiURL <- paste(url1, coord, url2, sep="")
      soilData <- GET(apiURL)
      pointData = fromJSON(rawToChar(soilData$content))
      names(pointData$property) <- "soilProp"
      if(pointData$property[[1]]$value$value != "No data"){
        if(soilProp == "texture_class"){
          SP <- data.frame(lat=lat, lon=lon, Variable=soilProp, depth=pointData$property[[1]]$depth$value[1],
                           depthUnit = pointData$property[[1]]$depth$unit[1],value=pointData$property[[1]]$value$value[1],
                           CI_50_lower = NA, CI_50_upper = NA,
                           CI_68_lower = NA, CI_68_upper = NA,
                           CI_90_lower = NA, CI_90_upper = NA)
        }else{
          UN <- pointData$property$soilProp$uncertainty[[1]]
          SP <- data.frame(lat=lat, lon=lon, Variable=soilProp, depth=pointData$property[[1]]$depth$value,
                           depthUnit = pointData$property[[1]]$depth$unit,value=pointData$property[[1]]$value$value, 
                           CI_50_lower = UN$lower_bound[1], CI_50_upper = UN$upper_bound[1],
                           CI_68_lower = UN$lower_bound[2], CI_68_upper = UN$upper_bound[2],
                           CI_90_lower = UN$lower_bound[3], CI_90_upper = UN$upper_bound[3])
        }
        soilVarData <- rbind(soilVarData, SP)
      }else{
        soilVarData <- soilVarData
      }
       
    }
    
    all_SoilProp <- rbind(all_SoilProp, soilVarData)
  }
  
  all_SoilProp$depth <- gsub("-", "_", all_SoilProp$depth)

  
  return(all_SoilProp)
}

#' @description a function to get soil data from iSDA 30m res for layers given for one level depth
#' @param coordPoints a data frame with lat and lon column
#' @param soilPropAll a list of soil variable with two depth
#' @param depth a vector with "0-20" or "0-50"
#' @return a data frame with lat, lon, Variable, depth, depthUnit, value and lower and upper limits for 50, 68 and
#' @example coordPoints <- data.frame(lat=-11.375, lon=31.525); getiSDA_data1Depth(coordPoints)
getiSDA_data1Depth <- function(coordPoints){
   soilVarDataB <- NULL
    for(i in 1:nrow(coordPoints)){
      print(i)
      lat <- coordPoints$lat[i]
      lon <- coordPoints$lon[i]
      url1 <- "https://api.isda-africa.com/v1/soilproperty?key=AIzaSyB7WZEpto9Jf1Do1SJBF7B8le7yVrEWN_U&"
      coord <- paste("lat=", lat, "&lon=", lon, sep="")
      soilvars <- c("land_cover_2015","land_cover_2016", "land_cover_2017","land_cover_2018","land_cover_2019", 
                    "crop_cover_2015", "crop_cover_2016", "crop_cover_2017", "crop_cover_2018","crop_cover_2019", 
                    "slope_angle","fcc")
      dd <- NULL
      for (h in 1:length(soilvars)){
        lCovers <- soilvars[h]
        if (lCovers == "fcc"){
          url3 <- paste("&property=",lCovers, "&depth=0-50", sep="")
          apiURL2 <- paste(url1, coord, url3, sep="")
          landCover = GET(apiURL2)
          landCover_latlong = fromJSON(rawToChar(landCover$content))
          names(landCover_latlong$property) <- "lCovers"
          llc <- data.frame(LC = landCover_latlong$property$lCovers$value$value)
          names(llc)[1] <- lCovers
          if(h == 1){
            dd <- llc
          }else{
            dd <- cbind(dd, llc)
          }
        }else{
        
          url3 <- paste("&property=",lCovers, "&depth=0", sep="")
          apiURL2 <- paste(url1, coord, url3, sep="")
          landCover = GET(apiURL2)
          landCover_latlong = fromJSON(rawToChar(landCover$content))
          names(landCover_latlong$property) <- "lCovers"
          if (lCovers == "land_cover_2015"){
            llc <- data.frame(lat=lat, lon=lon, LC = landCover_latlong$property$lCovers$value$value)
            names(llc)[3] <- lCovers
          }else{
            llc <- data.frame(LC = landCover_latlong$property$lCovers$value$value)
            names(llc)[1] <- lCovers
          }
          if(h == 1){
            dd <- llc
          }else{
            dd <- cbind(dd, llc)
          }
        }
      }
      soilVarDataB <- rbind(soilVarDataB, dd)
    }
   colnames(soilVarDataB) <- gsub("-", "_", colnames(soilVarDataB))
   
  return(soilVarDataB)
}



###############################################################################################################
## 
###############################################################################################################
require(tidyverse)
setwd("D:/ACAI_Wrapper/Rwanda")
ACAI_NOT <- read.csv("NOT_Sep2020.csv") 
ACAI_NOT <- unique(select(ACAI_NOT, lat, lon))
ACAI_NOT <- ACAI_NOT[complete.cases(ACAI_NOT),]
str(ACAI_NOT)

MainProp_020cm <- getiSDA_data2Depth(coordPoints = ACAI_NOT, soilPropAll = Property_2depth, depth = "0-20")
MainProp_2050cm <- getiSDA_data2Depth(coordPoints = ACAI_NOT, soilPropAll = Property_2depth, depth = "20-50")
bedRock_200cm <- getiSDA_data2Depth(coordPoints = ACAI_NOT, soilPropAll = "bedrock_depth", depth = "0-200")
landCropCover <- getiSDA_data1Depth(coordPoints=ACAI_NOT)
ACAI_isda <- rbind(MainProp_020cm, MainProp_2050cm,bedRock_200cm)
saveRDS(landCropCover, "ACAI_isda_cover.RDS" )
saveRDS(ACAI_isda, "ACAI_isda.RDS")



  
  
setwd("/home/akilimo/projects/SoilNPK/Coordinates")
RW_NOT <- read.csv("Rwanda_Coordinates.csv") 
unique(RW_NOT$NAME_2)
RW_NOT <- unique(select(RW_NOT, lat, long))
RW_NOT <- RW_NOT[complete.cases(RW_NOT),]
names(RW_NOT) <- c("lat", "lon")
str(RW_NOT)
RW_Main_020cm <- getiSDA_data2Depth(coordPoints = RW_NOT, soilPropAll = Property_2depth, depth = "0-20")
RW_Main_2050cm <- getiSDA_data2Depth(coordPoints = RW_NOT, soilPropAll = Property_2depth, depth = "20-50")
RW_bedRock_200cm <- getiSDA_data2Depth(coordPoints = RW_NOT, soilPropAll = "bedrock_depth", depth = "0-200")
RW_landCropCover <- getiSDA_data1Depth(RW_NOT)
RW_isda <- rbind(RW_Main_020cm, RW_Main_2050cm, RW_bedRock_200cm)
saveRDS(RW_landCropCover, "RW_isda_cover.RDS" )
saveRDS(RW_isda, "RW_isda.RDS")


setwd("/home/akilimo/projects/SoilNPK/Coordinates")
coordPoints <- data.frame(lat=c(-1.575, -1.625, -1.625, -1.625, -2.075, -2.075), 
                          lon=c(29.325, 29.325, 29.375, 29.425, 29.475, 29.525),
                          NAME_1= c("Iburengerazuba", "Iburengerazuba","Iburengerazuba","Iburengerazuba","Iburengerazuba","Iburengerazuba"),
                          NAME_2=c("Rubavu", "Rubavu", "Rubavu", "Rubavu","Karongi","Karongi"))


RW_Main_020cm_missing <- getiSDA_data2Depth(coordPoints = coordPoints, soilPropAll = Property_2depth, depth = "0-20")
RW_Main_2050cm_missing <- getiSDA_data2Depth(coordPoints = coordPoints, soilPropAll = Property_2depth, depth = "20-50")
RW_bedRock_200cm_missing <- getiSDA_data2Depth(coordPoints = coordPoints, soilPropAll = "bedrock_depth", depth = "0-200")
RW_landCropCover_missing <- getiSDA_data1Depth(coordPoints)
RW_isda_missing <- rbind(RW_Main_020cm_missing, RW_Main_2050cm_missing, RW_bedRock_200cm_missing)
saveRDS(RW_landCropCover_missing, "RW_isda_cover_missing.RDS" )
saveRDS(RW_isda_missing, "RW_isda_missing.RDS")



setwd("/home/akilimo/projects/SoilNPK/Coordinates")
BU_NOT <- read.csv("Burundi_Coordinates.csv") 
BU_NOT <- unique(select(BU_NOT, lat, long))
BU_NOT <- BU_NOT[complete.cases(BU_NOT),]
names(BU_NOT) <- c("lat", "lon")
str(BU_NOT)
BU_Main_020cm <- getiSDA_data2Depth(coordPoints = BU_NOT, soilPropAll = Property_2depth, depth = "0-20")
BU_Main_2050cm <- getiSDA_data2Depth(coordPoints = BU_NOT, soilPropAll = Property_2depth, depth = "20-50")
BU_bedRock_200cm <- getiSDA_data2Depth(coordPoints = BU_NOT, soilPropAll = "bedrock_depth", depth = "0-200")
BU_landCropCover <- getiSDA_data1Depth(BU_NOT)
BU_isda <- rbind(BU_Main_020cm, BU_Main_2050cm, BU_bedRock_200cm)
saveRDS(BU_landCropCover, "BU_isda_cover.RDS" )
saveRDS(BU_isda, "BU_isda.RDS")


setwd("/home/akilimo/projects/SoilNPK/Coordinates")
DRC_NOT <- read.csv("DRC_Coordinates.csv") 
DRC_NOT <- unique(select(DRC_NOT, lat, long))
DRC_NOT <- DRC_NOT[complete.cases(DRC_NOT),]
names(DRC_NOT) <- c("lat", "lon")
str(DRC_NOT)
DRC_Main_020cm <- getiSDA_data2Depth(coordPoints = DRC_NOT, soilPropAll = Property_2depth, depth = "0-20")
DRC_Main_2050cm <- getiSDA_data2Depth(coordPoints = DRC_NOT, soilPropAll = Property_2depth, depth = "20-50")
DRC_bedRock_200cm <- getiSDA_data2Depth(coordPoints = DRC_NOT, soilPropAll = "bedrock_depth", depth = "0-200")
DRC_landCropCover <- getiSDA_data1Depth(DRC_NOT)
DRC_isda <- rbind(DRC_Main_020cm, DRC_Main_2050cm, DRC_bedRock_200cm)
saveRDS(DRC_landCropCover, "DRC_isda_cover.RDS" )
saveRDS(DRC_isda, "DRC_isda.RDS")

DRC_landCropCover$cropCover <- ifelse(DRC_landCropCover$crop_cover_2015 > 0, "Yes", 
                                      ifelse(DRC_landCropCover$crop_cover_2016 > 0 , "Yes", 
                                             ifelse(DRC_landCropCover$crop_cover_2017 > 0, "Yes",
                                                    ifelse(DRC_landCropCover$crop_cover_2018 > 0, "Yes",
                                                           ifelse(DRC_landCropCover$crop_cover_2019 > 0,"Yes", "No")))))


DRC_landCropCover_crop <- filter(DRC_landCropCover, cropCover == "Yes")
DRC_landCropCover_crop <- DRC_landCropCover_crop[DRC_landCropCover_crop$cropCover == "Yes"]


par(mfrow=c(2,3))
hist(DRC_landCropCover_crop$crop_cover_2015)
hist(DRC_landCropCover_crop$crop_cover_2016)
hist(DRC_landCropCover_crop$crop_cover_2017)
hist(DRC_landCropCover_crop$crop_cover_2018)
hist(DRC_landCropCover_crop$crop_cover_2019)
summary(DRC_landCropCover_crop)



#### #############################################################################
## use pedo-transfer functions to get soil hydraulic data aggregated for the top 50cm (avw)
### for LINTUL, we need single value of hydraulics data per oixel
## for the random forest data we could consider teh two depth apart ...
##################################################################################
#' get the weighted average and soil hydraulic data
gethydraulics <- function(idata){
  idata$location <- paste(idata$lat, idata$lon, sep="_")
  isda_texClass <- idata[idata$Variable == "texture_class", ]
  isda_texClass_top <- isda_texClass[isda_texClass$depth == "0-20",]
  isda_texClass_bottom <- isda_texClass[isda_texClass$depth == "20-50",]
  isda_texClass_top$texture_class_top <- isda_texClass_top$value
  isda_texClass_bottom$texture_class_bottom <- isda_texClass_bottom$value
  
  isda_texClass_top <- unique(isda_texClass_top[, c("lat","lon","Variable","location","texture_class_top")])
  isda_texClass_bottom <- unique(isda_texClass_bottom[, c("lat","lon","Variable","location","texture_class_bottom")])
  isda_texClass <- unique(merge(isda_texClass_top, isda_texClass_bottom, by=c("lat","lon","Variable","location")))
  
  
  isda_rest <- idata[!idata$Variable %in% c("texture_class", "bedrock_depth"), ]
  
  isda_rest$value <- as.numeric(as.character(isda_rest$value))
  
  
  isda_rest$depth <- revalue(isda_rest$depth, c("0-20" = 20, "20-50" = 30))
  isda_rest$depth <- as.numeric(as.character(isda_rest$depth))
  isda_rest$value2 <- isda_rest$value * isda_rest$depth
  isda_avwt <- ddply(isda_rest, .(lat, lon, Variable, location), summarise, value = (sum(value2)/50))
  
  
  isda_wide <- spread(isda_avwt, "Variable", "value")
  isda_wide <- unique(merge(isda_wide, isda_texClass[,c("location","texture_class_top", "texture_class_bottom")], by = "location"))
  
  isda_wide$PercentSOM <- (isda_wide$carbon_organic * 2)/10
  
  
  ##### WP ######
  isda_wide$wp_wAverage <- (-0.024*isda_wide$sand_content/100) + 0.487*isda_wide$clay_content/100 + 0.006*isda_wide$PercentSOM + 
    0.005*(isda_wide$sand_content/100 *isda_wide$PercentSOM) - 0.013*(isda_wide$clay_content/100 * isda_wide$PercentSOM) + 
    0.068*(isda_wide$sand_content/100 * isda_wide$clay_content/100 ) + 0.031
  isda_wide$wp_wAverage <- (isda_wide$wp_wAverage + (0.14 * isda_wide$wp_wAverage - 0.02))*100
  summary(isda_wide$wp_wAverage)
  
  
  ##### FC ######
  isda_wide$FC_wAverage <- -0.251*isda_wide$sand_content/100 + 0.195*isda_wide$clay_content/100 + 0.011*isda_wide$PercentSOM + 
    0.006*(isda_wide$sand_content/100 * isda_wide$PercentSOM) - 0.027*(isda_wide$clay_content/100 * isda_wide$PercentSOM) + 
    0.452*(isda_wide$sand_content/100 * isda_wide$clay_content/100) + 0.299
  isda_wide$FC_wAverage <- (isda_wide$FC_wAverage + (1.283 * isda_wide$FC_wAverage^2 - 0.374 * isda_wide$FC_wAverage - 0.015))*100
  summary(isda_wide$FC_wAverage)
  
  ##### soil water at saturation######
  isda_wide$sws_wAverage <- 0.278*(isda_wide$sand_content/100)+0.034*(isda_wide$clay_content/100)+0.022*isda_wide$PercentSOM-
    0.018*(isda_wide$sand_content/100*isda_wide$PercentSOM)-
    0.027*(isda_wide$clay_content/100*isda_wide$PercentSOM)-0.584*(isda_wide$sand_content/100*isda_wide$clay_content/100)+0.078
  
  isda_wide$sws_wAverage <- (isda_wide$sws_wAverage +(0.636*isda_wide$sws_wAverage-0.107))*100
  isda_wide$sws_wAverage <- (isda_wide$FC/100+isda_wide$sws_wAverage/100-(0.097*isda_wide$sand_content/100)+0.043)*100
  return(isda_wide)
  
}


setwd("/home/akilimo/projects/SoilNPK/Coordinates")
RW_isda <- readRDS("RW_isda.RDS")
RW_isda_missing <- readRDS("RW_isda_missing.RDS")
RW_isda <- rbind(RW_isda, RW_isda_missing)

RW_isda_avw <- gethydraulics(idata=RW_isda)
head(RW_isda_avw)
saveRDS(RW_isda_avw, "RW_isda_avw.RDS")


BU_isda <- readRDS("BU_isda.RDS")
BU_isda_avw <- gethydraulics(idata=BU_isda)
head(BU_isda_avw)
saveRDS(BU_isda_avw, "BU_isda_avw.RDS")


DRC_isda <- readRDS("DRC_isda.RDS")
DRC_isda_avw <- gethydraulics(idata=DRC_isda)
head(DRC_isda_avw)
saveRDS(DRC_isda_avw, "DRC_isda_avw.RDS")



#### #############################################################################
## use pedo-transfer functions to get soil hydraulic for the two depth apart
#### #############################################################################
##TODO

gethydraulics_bydepth <- function(idata){
  head(idata)
  idata <- droplevels(idata[!idata$Variable == "bedrock_depth", ])
  idata$Variable <- paste(idata$Variable, idata$depth, sep="_")
  isda_wide <- spread(idata[, c("lat","lon","Variable", "value")], "Variable", "value")
  colnames(isda_wide) <- gsub("-", "_", colnames(isda_wide))
  isda_wide$index <- c(1:nrow(isda_wide)) 
 
  isdatex <- isda_wide[, c("texture_class_0_20", "texture_class_20_50", "index")]
  isdachar <- subset(isda_wide, select=-c(texture_class_0_20, texture_class_20_50))
  isdachars <- isdachar %>% mutate_if(is.character,as.numeric)
  isda_wide <- unique(merge(isdachars, isdatex, by="index"))
  str(isda_wide)
  isda_wide$percentSOM_0_20 <- (isda_wide$carbon_organic_0_20 * 2)/10
  isda_wide$percentSOM_20_50 <- (isda_wide$carbon_organic_20_50 * 2)/10
  idata <- isda_wide
  
  ##### WP ######
  idata$wp_0_20A <- (-0.024*idata$sand_content_0_20/100) + 0.487*idata$clay_content_0_20/100 + 0.006*idata$percentSOM_0_20 + 
    0.005*(idata$sand_content_0_20/100 *idata$percentSOM_0_20) - 0.013*(idata$clay_content_0_20/100 * idata$percentSOM_0_20) + 
    0.068*(idata$sand_content_0_20/100 * idata$clay_content_0_20/100 ) + 0.031
  idata$wp_0_20 <- (idata$wp_0_20A + (0.14 * idata$wp_0_20A - 0.02))*100
  
  idata$wp_20_50A <- (-0.024*idata$sand_content_20_50/100) + 0.487*idata$clay_content_20_50/100 + 0.006*idata$percentSOM_20_50 + 
    0.005*(idata$sand_content_20_50/100 *idata$percentSOM_20_50) - 0.013*(idata$clay_content_20_50/100 * idata$percentSOM_20_50) + 
    0.068*(idata$sand_content_20_50/100 * idata$clay_content_20_50/100 ) + 0.031
  idata$wp_20_50 <- (idata$wp_20_50A + (0.14 * idata$wp_20_50A - 0.02))*100
  
  idata <- subset(idata, select=-c(wp_0_20A, wp_20_50A))
  
  head(idata)
  
  ##### FC ######
  idata$FC_0_20_A <- -0.251*idata$sand_content_0_20/100 + 0.195*idata$clay_content_0_20/100 + 0.011*idata$percentSOM_0_20 + 
    0.006*(idata$sand_content_0_20/100 * idata$percentSOM_0_20) - 0.027*(idata$clay_content_0_20/100 * idata$percentSOM_0_20) + 
    0.452*(idata$sand_content_0_20/100 * idata$clay_content_0_20/100) + 0.299
  idata$FC_0_20 <- (idata$FC_0_20_A + (1.283 * idata$FC_0_20_A^2 - 0.374 * idata$FC_0_20_A - 0.015))*100
  
  
  idata$FC_20_50_A <- -0.251*idata$sand_content_20_50/100 + 0.195*idata$clay_content_20_50/100 + 0.011*idata$percentSOM_20_50 + 
    0.006*(idata$sand_content_20_50/100 * idata$percentSOM_20_50) - 0.027*(idata$clay_content_20_50/100 * idata$percentSOM_20_50) + 
    0.452*(idata$sand_content_20_50/100 * idata$clay_content_20_50/100) + 0.299
  idata$FC_20_50 <- (idata$FC_20_50_A + (1.283 * idata$FC_20_50_A^2 - 0.374 * idata$FC_20_50_A - 0.015))*100
  
  
  idata <- subset(idata, select=-c(FC_0_20_A, FC_20_50_A))
  
  
  ##### soil water at ######
  idata$sws_0_20_A <- 0.278*(idata$sand_content_0_20/100)+0.034*(idata$clay_content_0_20/100)+0.022*idata$percentSOM_0_20-
    0.018*(idata$sand_content_0_20/100*idata$percentSOM_0_20)-
    0.027*(idata$clay_content_0_20/100*idata$percentSOM_0_20)-0.584*(idata$sand_content_0_20/100*idata$clay_content_0_20/100)+0.078
  
  idata$sws_0_20_B <- (idata$sws_0_20_A +(0.636*idata$sws_0_20_A-0.107))*100
  idata$sws_0_20 <- (idata$FC_0_20/100+idata$sws_0_20_B/100-(0.097*idata$sand_content_0_20/100)+0.043)*100
  
  
  idata$sws_20_50_A <- 0.278*(idata$sand_content_20_50/100)+0.034*(idata$clay_content_20_50/100)+0.022*idata$percentSOM_20_50-
    0.018*(idata$sand_content_20_50/100*idata$percentSOM_20_50)-
    0.027*(idata$clay_content_20_50/100*idata$percentSOM_20_50)-0.584*(idata$sand_content_20_50/100*idata$clay_content_20_50/100)+0.078
  idata$sws_20_50_B <- (idata$sws_20_50_A +(0.636*idata$sws_20_50_A-0.107))*100
  idata$sws_20_50 <- (idata$FC_20_50/100+idata$sws_20_50_B/100-(0.097*idata$sand_content_20_50/100)+0.043)*100
  
  
  idata <- subset(idata, select=-c(sws_0_20_A, sws_0_20_B, sws_20_50_A, sws_20_50_B))
  str(idata)
  idata$location <- paste(idata$lon, idata$lat, sep="_")
  
  
  return(idata)
  
}


setwd("/home/akilimo/projects/SoilNPK/Coordinates")
RW_isda <- readRDS("RW_isda.RDS")
RW_isda_missing <- readRDS("RW_isda_missing.RDS")


RW_isda_H2O2dapth <- gethydraulics_bydepth(idata=RW_isda)
head(RW_isda_H2O2dapth)
RW_isda_H2O2dapth <- RW_isda_H2O2dapth[complete.cases(RW_isda_H2O2dapth), ]
saveRDS(RW_isda_H2O2dapth, "RW_isda_H2O2dapth.RDS")

RW_isda_H2O2dapth_missing <- gethydraulics_bydepth(idata=RW_isda_missing)
head(RW_isda_H2O2dapth_missing)
RW_isda_H2O2dapth_missing <- RW_isda_H2O2dapth_missing[complete.cases(RW_isda_H2O2dapth_missing), ]
saveRDS(RW_isda_H2O2dapth_missing, "RW_isda_H2O2dapth_missing.RDS")

RW_isda_H2O2dapth <- readRDS("RW_isda_H2O2dapth.RDS")
RW_NOT <- read.csv("Rwanda_Coordinates.csv") 
RW_isda_H2O2dapth$location <- paste(RW_isda_H2O2dapth$lat, RW_isda_H2O2dapth$lon,se="_")

RW_isda_H2O2dapth <- unique(merge(RW_isda_H2O2dapth, RW_NOT[, c("NAME_1", "NAME_2", "location")], by="location"))

head(RW_isda_H2O2dapth)
head(RW_NOT)


setwd("/home/akilimo/projects/SoilNPK/Coordinates")
BU_isda <- readRDS("BU_isda.RDS")
BU_isda_H2O2dapth <- gethydraulics_bydepth(idata=BU_isda)
head(BU_isda_H2O2dapth)
BU_isda_H2O2dapth <- BU_isda_H2O2dapth[complete.cases(BU_isda_H2O2dapth), ]
saveRDS(BU_isda_H2O2dapth, "BU_isda_H2O2dapth.RDS")


setwd("/home/akilimo/projects/SoilNPK/Coordinates")
DRC_isda <- readRDS("DRC_isda.RDS")
DRC_isda_H2O2dapth <- gethydraulics_bydepth(idata=DRC_isda)
head(DRC_isda_H2O2dapth)
DRC_isda_H2O2dapth <- DRC_isda_H2O2dapth[complete.cases(DRC_isda_H2O2dapth), ]
saveRDS(DRC_isda_H2O2dapth, "DRC_isda_H2O2dapth.RDS")


setwd("/home/akilimo/projects/SoilNPK/Coordinates")
ACAI_isda <- readRDS("ACAI_isda.RDS")
ACAI_isda_H2O2dapth <- gethydraulics_bydepth(idata=ACAI_isda)
head(ACAI_isda_H2O2dapth)
ACAI_isda_H2O2dapth <- ACAI_isda_H2O2dapth[complete.cases(ACAI_isda_H2O2dapth), ]
saveRDS(ACAI_isda_H2O2dapth, "ACAI_isda_H2O2dapth.RDS")







#################################################################################################
## get iSDA CILACA trial data 
#################################################################################################
CIALCA_coordinates <- readRDS("CIALCA_coordainte.RDS")
str(CIALCA_coordinates)
names(CIALCA_coordinates) <- c("lat", "lon")

CIALCA_Main_020cm <- getiSDA_data2Depth(coordPoints = CIALCA_coordinates, soilPropAll = Property_2depth, depth = "0-20")
CIALCA_Main_2050cm <- getiSDA_data2Depth(coordPoints = CIALCA_coordinates, soilPropAll = Property_2depth, depth = "20-50")
CIALCA_bedRock_200cm <- getiSDA_data2Depth(coordPoints = CIALCA_coordinates, soilPropAll = "bedrock_depth", depth = "0-200")
CIALCA_landCropCover <- getiSDA_data1Depth(CIALCA_coordinates)
CIALCA_isda <- rbind(CIALCA_Main_020cm, CIALCA_Main_2050cm, CIALCA_bedRock_200cm)
setwd("/home/akilimo/projects/SoilNPK/Coordinates")
saveRDS(CIALCA_landCropCover, "CIALCA_isda_cover.RDS" )
saveRDS(CIALCA_isda, "CIALCA_isda.RDS")

CIALCA_isda <- readRDS("CIALCA_isda.RDS")
CIALCA_landCropCover <- readRDS("CIALCA_isda_cover.RDS")

CIALCA_landCropCover <- CIALCA_landCropCover[, c("lat","lon", "slope_angle", "fcc")]
CIALCA_isda_main <- CIALCA_isda[!CIALCA_isda$Variable == "bedrock_depth", ]
CIALCA_bedRock_200cm <- CIALCA_isda[CIALCA_isda$Variable == "bedrock_depth", ]
bedrockdepth1 <- CIALCA_bedRock_200cm[seq(1,180, by=2), ]
bedrockdepth2 <- CIALCA_bedRock_200cm[seq(2,180, by=2), ]

summary(bedrockdepth1$value)
summary(bedrockdepth2$value)

CIALCA_isda <- rbind(CIALCA_isda_main, bedrockdepth1)
head(CIALCA_isda)
CIALCA_isda$Variable <- paste(CIALCA_isda$Variable, CIALCA_isda$depth, sep="_")
CIALCA_isda_wide <- spread(CIALCA_isda[, c("lat","lon","Variable", "value")], Variable, value, )
head(CIALCA_isda_wide)

CIALCA_isda_wide <- unique(merge(CIALCA_isda_wide, CIALCA_landCropCover, by=c("lat", "lon")))
str(CIALCA_isda_wide)

#################################################################################################
#################################################################################################
## get iSDA data from the bulk download files
#################################################################################################
#################################################################################################

library(raster)
library(sf)
library(rgdal)
require(rgeos)
setwd("E:/iSDASoils/pH")
isda_20_pH <- raster("sol_ph_h2o_md_30m_0..20cm_2001..2017_v0.13_wgs84.tif")
isda_50_pH <- raster("sol_ph_h2o_m_30m_20..50cm_2001..2017_v0.13_wgs84.tif")
data.frame(pH_upper = raster::extract(isda_20_pH, data.frame(x=31.525, y=-11.375)))
data.frame(pH_lower = raster::extract(isda_50_pH, data.frame(x=31.525, y=-11.375)))


setwd("E:/iSDASoils/soiltexture")
isda_20_texClass <- raster("sol_texture.class_m_30m_0..20cm_2001..2017_v0.13_wgs84.tif")
isda_50_texClass <- raster("sol_texture.class_m_30m_20..50cm_2001..2017_v0.13_wgs84.tif")
data.frame(texCalss_upper = raster::extract(isda_20_texClass, data.frame(x=29.625, y=-16)))
data.frame(texCalss_lower = raster::extract(isda_50_texClass, data.frame(x=31.525, y=-11.375)))


isda_20_silt <- raster("sol_silt_tot_psa_m_30m_0..20cm_2001..2017_v0.13_wgs84.tif")
isda_50_silt <- raster("sol_silt_tot_psa_m_30m_20..50cm_2001..2017_v0.13_wgs84.tif")
data.frame(silt_upper = raster::extract(isda_20_silt, data.frame(x=31.525, y=-11.375)))
data.frame(silt_lower = raster::extract(isda_50_silt, data.frame(x=31.525, y=-11.375)))


isda_20_clay <- raster("sol_clay_tot_psa_m_30m_0..20cm_2001..2017_v0.13_wgs84.tif")
isda_50_clay <- raster("sol_clay_tot_psa_m_30m_20..50cm_2001..2017_v0.13_wgs84.tif")
data.frame(clay_upper = raster::extract(isda_20_clay, data.frame(x=31.525, y=-11.375)))
data.frame(clay_lower = raster::extract(isda_50_clay, data.frame(x=31.525, y=-11.375)))


isda_20_sand <- raster("sol_sand_tot_psa_m_30m_0..20cm_2001..2017_v0.13_wgs84.tif")
isda_50_sand <- raster("sol_sand_tot_psa_m_30m_20..50cm_2001..2017_v0.13_wgs84.tif")
data.frame(sand_upper = raster::extract(isda_20_sand, data.frame(x=31.525, y=-11.375)))
data.frame(sand_lower = raster::extract(isda_50_sand, data.frame(x=31.525, y=-11.375)))


setwd("E:/iSDASoils/PMehlich")
isda_20_P <- raster("sol_log.p_mehlich3_m_30m_0..20cm_2001..2017_v0.13_wgs84.tif")
isda_50_P <- raster("sol_log.p_mehlich3_m_30m_20..50cm_2001..2017_v0.13_wgs84.tif")
data.frame(Pmehlich_upper = raster::extract(isda_20_P, data.frame(x=31.525, y=-11.375)))
data.frame(Pmehlich_lower = raster::extract(isda_50_P, data.frame(x=31.525, y=-11.375)))