# Investigating whether aggregated predictions of cassava production better correlate
# using buffer zones from Had

library(ggplot2)
library(sf)
library("rnaturalearth")
library("rnaturalearthdata")
library(rgeos)
library(ggpubr)
library(lattice)
library(gstat)
library(sp)
library(mgcv)


panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...)
{
  usr <- par("usr"); on.exit(par(usr))
  par(usr = c(0, 1, 0, 1))
  r <- (cor(x, y, use="pairwise"))
  txt <- format(c(r, 0.123456789), digits = digits)[1]
  txt <- paste0(prefix, txt)
  if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
  text(0.5, 0.5, txt, cex = cex.cor)
}



civ <- ne_countries(scale = "medium", returnclass = "sf", country="Ivory Coast")
uga <- ne_countries(scale = "medium", returnclass = "sf", country="Uganda")


# input and merge all cassava sample data to define total cassava production

rm(list=ls()[which(!(ls()%in%c("civ","uga","dat_civ","dat_uga","panel.cor")))])
# CIV
{
  
  
  # area of individual plants in each study area
  indplants <- read.csv("KLH_Civ_individual_plants_poly.csv")
  head(indplants)
  table(indplants$field_id)
  
  # aggregate number of plants per field
  plants <- data.frame(aggregate(area_m2 ~ field_id, FUN=sum, data=indplants), table(indplants$field_id))
  names(plants) <- c("field_id", "tot_area_ind_plants", "x", "tot_ind_plants")
  plants <- plants[, c("field_id", "tot_area_ind_plants", "tot_ind_plants")]
  
  # number of buildings / inhabitants in each study area
  inhabit <- read.csv("KLH_Civ_inhabited_pts.csv")
  head(inhabit)
  table(inhabit$Subcategor)
  table(inhabit$FloorNo)
  table(inhabit$field_id)
  
  # aggregate number of buildings per field
  temp <- data.frame(table(inhabit$field_id, inhabit$FloorNo))
  names(temp) <- c("field_id", "FloorNo", "Freq")
  temp2 <- reshape(temp, direction="wide", v.names="Freq", timevar="FloorNo", idvar="field_id")
  temp2$tot <- temp2$Freq.1 + temp2$Freq.2
  buildings <- temp2
  names(buildings)<- c("field_id", "1story_buildings", "2story_buildings", "tot_buildings")
  buildings$field_id <- as.numeric(buildings$field_id)
  
  
  # cassava production models. also includes area of each study area
  models <- read.csv("Civ_study_area_centre_sample_aug20.csv")
  head(models)
  
  # cassava cultivation
  cassava <- read.csv("KLH_Civ_cassava_fields_poly.csv")
  head(cassava)
  table(cassava$Name)
  table(cassava$Name, cassava$Subcategor)
  table(cassava$Density)
  table(cassava$Density, cassava$field_id)
  table(cassava$Density, cassava$Subcategor)
  
  # aggregate to area in monoculture and area in intercropping per field, total cassava area, what is intercropped, number of fields - bit dodgy due to intersecting polygons...
  
  c1 <- aggregate(area_m2 ~ field_id, data=cassava, FUN=sum)
  names(c1) <- c("field_id", "tot_cassava_area")
  c2 <- aggregate(area_m2 ~ field_id, data=cassava[cassava$Name == "Cassava Monoculture", ], FUN=sum)
  names(c2) <- c("field_id", "tot_monoculture_area")
  c3 <- aggregate(area_m2 ~ field_id, data=cassava[cassava$Name == "Cassava", ], FUN=sum)
  names(c3) <- c("field_id", "tot_intercrop_area")
  c4 <- data.frame(table(cassava$field_id))
  names(c4) <- c("field_id", "no_fields")
  c5 <- aggregate(cbind(Banana, Pineapple, Palm.oil, Maize, Papaya, Avocado, Sugar.cane, Cacao, Ivory.Palm, Guava, Eggplant, Cashew, Coffee, Peanut, Okra, Pepper, Mango, Coconut, Rubber, Taro, Lemon, Chilli, Palm.Tree, Soursop, Popo) ~ field_id, data=cassava, FUN=sum, na.rm=TRUE, na.action=na.pass)
  
  fields <- merge(c1, c2, by="field_id", all.x=TRUE, all.y=TRUE)
  fields <- merge(fields, c3, by="field_id", all.x=TRUE, all.y=TRUE)
  fields <- merge(fields, c4, by="field_id", all.x=TRUE, all.y=TRUE)
  fields <- merge(fields, c5, by="field_id", all.x=TRUE, all.y=TRUE)
  
  
  # so, can we refine the estimate of cassava production from the study sites...
  
  # issue 1. area sampled is different
  # scale area by size sampled
  areaSampled <- models$area_m2
  
  
  cassava$weights <- 1
  cassava$weights[cassava$Density == "dense"] <- 1.5
  cassava$weights[cassava$Desnity == "dense..?"] <- 1.5
  cassava$weights[cassava$Desnity == "densea"] <- 1.5
  cassava$weights[cassava$Desnity == "regular/sparse"] <- 0.75
  cassava$weights[cassava$Desnity == "sparce"] <- 0.5
  cassava$weights[cassava$Desnity == "sparse"] <- 0.5
  cassava$weights[cassava$Desnity == "sparse and uprooted"] <- 0.5
  cassava$weights[cassava$Desnity == "sparse but"] <- 0.5
  cassava$weights[cassava$Desnity == "sparse2aqqw1"] <- 0.5
  cassava$weights[cassava$Desnity == "very dense"] <- 1.75
  cassava$weights[cassava$Desnity == "very sparce"] <- 0.25
  cassava$weights[cassava$Desnity == "very sparse"] <- 0.25
  
  cassava$area_w <- cassava$area_m2 * cassava$weights
  
  c1 <- aggregate(area_w ~ field_id, data=cassava, FUN=sum)
  names(c1) <- c("field_id", "tot_cassava_area_w")
  c2 <- aggregate(area_w ~ field_id, data=cassava[cassava$Name == "Cassava Monoculture", ], FUN=sum)
  names(c2) <- c("field_id", "tot_monoculture_area_w")
  c3 <- aggregate(area_w ~ field_id, data=cassava[cassava$Name == "Cassava", ], FUN=sum)
  names(c3) <- c("field_id", "tot_intercrop_area_w")
  
  fields_w <- merge(c1, c2, by="field_id", all.x=TRUE, all.y=TRUE)
  fields_w <- merge(fields_w, c3, by="field_id", all.x=TRUE, all.y=TRUE)
  
  
  # check field ids
  plants$field_id[which(!plants$field_id %in% models$field_id)]
  buildings$field_id[which(!buildings$field_id %in% models$field_id )]
  fields$field_id[which(!fields$field_id %in% models$field_id )]
  fields_w$field_id[which(!fields_w$field_id %in% models$field_id )]
  
  
  dim(models)
  dim(plants)
  dat_civ <- merge(models, plants, by="field_id", all.x=TRUE, all.y=TRUE)
  dim(dat_civ)
  dim(buildings)
  dat_civ <- merge(dat_civ, buildings, by="field_id", all.x=TRUE, all.y=TRUE)
  dim(dat_civ)
  dat_civ <- merge(dat_civ, fields, by="field_id", all.x=TRUE, all.y=TRUE)
  dim(dat_civ)
  dat_civ <- merge(dat_civ, fields_w, by="field_id", all.x=TRUE, all.y=TRUE)
  dim(dat_civ)
  
  
  
  # where no Cassava production was recorded, replace NAs by 0
  names(dat_civ)
  
  which(!models$field_id %in% plants$field_id)
  which(is.na(dat_civ$tot_area_ind_plants))
  dat_civ[which(!models$field_id %in% plants$field_id), "tot_area_ind_plants"] <- 0
  dat_civ[which(!models$field_id %in% plants$field_id), "tot_ind_plants"] <- 0
  which(!models$field_id %in% buildings$field_id)
  which(is.na(dat_civ$`1story_buildings`))
  which(is.na(dat_civ$`2story_buildings`))
  which(is.na(dat_civ$tot_buildings))
  dat_civ[which(!models$field_id %in% buildings$field_id), "1story_buildings"] <- 0
  dat_civ[which(!models$field_id %in% buildings$field_id), "2story_buildings"] <- 0
  dat_civ[which(!models$field_id %in% buildings$field_id), "tot_buildings"] <- 0
  which(is.na(dat_civ$tot_cassava_area))
  which(is.na(dat_civ$tot_intercrop_area))
  which(is.na(dat_civ$tot_monoculture_area))
  which(!models$field_id %in% fields$field_id)
  dat_civ[which(!models$field_id %in% fields$field_id), c("tot_cassava_area", "tot_monoculture_area", "tot_intercrop_area", "no_fields")] <- 0
  dat_civ[which(is.na(dat_civ$tot_intercrop_area)), "tot_intercrop_area"] <- 0
  dat_civ[which(is.na(dat_civ$tot_monoculture_area)), "tot_monoculture_area"] <- 0
  dat_civ[which(!models$field_id %in% fields$field_id), c("Banana", "Pineapple", "Palm.oil","Maize", "Papaya","Avocado", "Scivr.cane", "Cacao", "Ivory.Palm", "Guava","Eggplant","Cashew","Coffee", "Peanut","Okra","Pepper" ,"Mango", "Coconut", "Rubber","Taro","Lemon","Chilli","Palm.Tree","Soursop","Popo")] <- 0
  dat_civ[which(!models$field_id %in% fields_w$field_id), c("tot_cassava_area_w", "tot_monoculture_area_w", "tot_intercrop_area_w")] <- 0
  dat_civ[which(is.na(dat_civ$tot_intercrop_area_w)), "tot_intercrop_area_w"] <- 0
  dat_civ[which(is.na(dat_civ$tot_monoculture_area_w)), "tot_monoculture_area_w"] <- 0
  
  
  dat_civ$tot_cassava_area <- (dat_civ$tot_area_ind_plants + dat_civ$tot_intercrop_area + dat_civ$tot_monoculture_area) / dat_civ$area_m2
  dat_civ$tot_cassava_area_w <- (dat_civ$tot_area_ind_plants + dat_civ$tot_intercrop_area_w*0.75 + dat_civ$tot_monoculture_area_w) / dat_civ$area_m2
  
  # create a Cassava presence/absence category
  dat_civ$cass_area_presence <- factor(dat_civ$tot_cassava_area > 0)
}

rm(list=ls()[which(!(ls()%in%c("civ","uga","dat_civ","dat_uga","panel.cor")))])
# Uganda
{
  # area of individual plants in each study area
  indplants <- read.csv("KLH_Uga_individual_plants_poly.csv")
  head(indplants)
  table(indplants$field_id)
  
  # aggregate number of plants per field
  plants <- data.frame(aggregate(area_m2 ~ field_id, FUN=sum, data=indplants), table(indplants$field_id))
  names(plants) <- c("field_id", "tot_area_ind_plants", "x", "tot_ind_plants")
  plants <- plants[, c("field_id", "tot_area_ind_plants", "tot_ind_plants")]
  plants
  
  # number of buildings / inhabitants in each study area
  inhabit <- read.csv("KLH_Uga_inhabited_pts.csv")
  head(inhabit)
  table(inhabit$Subcategor)
  table(inhabit$FloorNo)
  table(inhabit$field_id)
  
  # aggregate number of buildings per field
  temp <- data.frame(table(inhabit$field_id, as.numeric(inhabit$Subcategor)))
  names(temp) <- c("field_id", "FloorNo", "Freq")
  temp2 <- reshape(temp, direction="wide", v.names="Freq", timevar="FloorNo", idvar="field_id")
  temp2$tot <- temp2$Freq.1 + temp2$Freq.2
  buildings <- temp2
  names(buildings)<- c("field_id", "1story_buildings", "2story_buildings", "3story_buildings","tot_buildings")
  buildings$field_id <- as.numeric(buildings$field_id)
  
  
  # cassava production models. also includes area of each study area
  models <- read.csv("Uga_study_area_centre_sample_aug20.csv")
  head(models)
  models$Cass_HA[models$Cass_HA == -9999] <- NA
  models$Cass_Prod[models$Cass_Prod == -9999] <- NA
  
  # cassava cultivation
  cassava <- read.csv("KLH_Uga_cassava_fields_poly.csv")
  head(cassava)
  table(cassava$Name)
  table(cassava$Name, cassava$Subcategor)
  table(cassava$Density)
  table(cassava$Density, cassava$field_id)
  table(cassava$Density, cassava$Subcategor)
  
  # aggregate to area in monoculture and area in intercropping per field, total cassava area, what is intercropped, number of fields - bit dodgy due to intersecting polygons...
  
  c1 <- aggregate(area_m2 ~ field_id, data=cassava, FUN=sum)
  names(c1) <- c("field_id", "tot_cassava_area")
  c2 <- aggregate(area_m2 ~ field_id, data=cassava[cassava$Name == "Cassava Monoculture", ], FUN=sum)
  names(c2) <- c("field_id", "tot_monoculture_area")
  c3 <- aggregate(area_m2 ~ field_id, data=cassava[cassava$Name == "Cassava", ], FUN=sum)
  names(c3) <- c("field_id", "tot_intercrop_area")
  c4 <- data.frame(table(cassava$field_id))
  names(c4) <- c("field_id", "no_fields")
  c5 <- aggregate(cbind(Banana,Coffee,Maize,Sugarcane,Sweet.Potato,Eucalyptus,Avocado,Taro.Cocoyam,Sorghum,Papaya,Pineapple,Tabacco,	Orange,Jackfruit,Beans,Pumpkin,Pigeonpea,Sesame,Guava,Mango,Cacia,Sambia,Eggplant,Citrus,Fig,Elephant.grass,Hibiscus,Cotton,Millet,Pine.tree,Yam,Graville,Tiak.tree) ~ field_id, data=cassava, FUN=sum, na.rm=TRUE, na.action=na.pass)
  
  fields <- merge(c1, c2, by="field_id", all.x=TRUE, all.y=TRUE)
  fields <- merge(fields, c3, by="field_id", all.x=TRUE, all.y=TRUE)
  fields <- merge(fields, c4, by="field_id", all.x=TRUE, all.y=TRUE)
  fields <- merge(fields, c5, by="field_id", all.x=TRUE, all.y=TRUE)
  
  
  # so, can we refine the estimate of cassava production from the study sites...
  
  # issue 1. area sampled is different
  # scale area by size sampled
  areaSampled <- models$area_m2
  
  
  cassava$weights <- 1
  cassava$weights[cassava$Density == "dense"] <- 1.5
  cassava$weights[cassava$Desnity == "dense at places"] <- 1.5
  cassava$weights[cassava$Desnity == "high density"] <- 1.5
  cassava$weights[cassava$Desnity == "low / medium density"] <- 0.75
  cassava$weights[cassava$Desnity == "low density"] <- 0.5
  cassava$weights[cassava$Desnity == "medium/low density"] <- 0.75
  cassava$weights[cassava$Desnity == "regular at times sparse"] <- 0.75
  cassava$weights[cassava$Desnity == "sparse"] <- 0.25
  cassava$weights[cassava$Desnity == "sparse estimated"] <- 0.25
  cassava$weights[cassava$Desnity == "very dense"] <- 1.75
  cassava$weights[cassava$Desnity == "very low density"] <- 0.25
  cassava$weights[cassava$Desnity == "very sparse"] <- 0.25
  
  cassava$area_w <- cassava$area_m2 * cassava$weights
  
  c1 <- aggregate(area_w ~ field_id, data=cassava, FUN=sum)
  names(c1) <- c("field_id", "tot_cassava_area_w")
  c2 <- aggregate(area_w ~ field_id, data=cassava[cassava$Name == "Cassava Monoculture", ], FUN=sum)
  names(c2) <- c("field_id", "tot_monoculture_area_w")
  c3 <- aggregate(area_w ~ field_id, data=cassava[cassava$Name == "Cassava", ], FUN=sum)
  names(c3) <- c("field_id", "tot_intercrop_area_w")
  
  fields_w <- merge(c1, c2, by="field_id", all.x=TRUE, all.y=TRUE)
  fields_w <- merge(fields_w, c3, by="field_id", all.x=TRUE, all.y=TRUE)
  
  
  # check field ids
  plants$field_id[which(!plants$field_id %in% models$field_id)]
  buildings$field_id[which(!buildings$field_id %in% models$field_id )]
  fields$field_id[which(!fields$field_id %in% models$field_id )]
  fields_w$field_id[which(!fields_w$field_id %in% models$field_id )]
  
  
  dim(models)
  dim(plants)
  dat_uga <- merge(models, plants, by="field_id", all.x=TRUE, all.y=TRUE)
  dim(dat_uga)
  dim(buildings)
  dat_uga <- merge(dat_uga, buildings, by="field_id", all.x=TRUE, all.y=TRUE)
  dim(dat_uga)
  dat_uga <- merge(dat_uga, fields, by="field_id", all.x=TRUE, all.y=TRUE)
  dim(dat_uga)
  dat_uga <- merge(dat_uga, fields_w, by="field_id", all.x=TRUE, all.y=TRUE)
  dim(dat_uga)
  
  
  # where no Cassava production was recorded, replace NAs by 0
  names(dat_uga)
  
  which(!models$field_id %in% plants$field_id)
  which(is.na(dat_uga$tot_area_ind_plants))
  dat_uga[which(!models$field_id %in% plants$field_id), "tot_area_ind_plants"] <- 0
  dat_uga[which(!models$field_id %in% plants$field_id), "tot_ind_plants"] <- 0
  which(!models$field_id %in% buildings$field_id)
  which(is.na(dat_uga$`1story_buildings`))
  which(is.na(dat_uga$`2story_buildings`))
  which(is.na(dat_uga$`3story_buildings`))
  which(is.na(dat_uga$tot_buildings))
  dat_uga[which(!models$field_id %in% buildings$field_id), "1story_buildings"] <- 0
  dat_uga[which(!models$field_id %in% buildings$field_id), "2story_buildings"] <- 0
  dat_uga[which(!models$field_id %in% buildings$field_id), "3story_buildings"] <- 0
  dat_uga[which(!models$field_id %in% buildings$field_id), "tot_buildings"] <- 0
  which(is.na(dat_uga$tot_cassava_area))
  which(is.na(dat_uga$tot_intercrop_area))
  which(is.na(dat_uga$tot_monoculture_area))
  which(!models$field_id %in% fields$field_id)
  dat_uga[which(!models$field_id %in% fields$field_id), c("tot_cassava_area", "tot_monoculture_area", "tot_intercrop_area", "no_fields")] <- 0
  dat_uga[which(is.na(dat_uga$tot_intercrop_area)), "tot_intercrop_area"] <- 0
  dat_uga[which(is.na(dat_uga$tot_monoculture_area)), "tot_monoculture_area"] <- 0
  dat_uga[which(!models$field_id %in% fields$field_id), c("Banana","Coffee","Maize","Sugarcane","Sweet.Potato","Eucalyptus","Avocado","Taro.Cocoyam","Sorghum","Papaya","Pineapple","Tabacco",	"Orange","Jackfruit","Beans","Pumpkin","Pigeonpea","Sesame","Guava","Mango","Cacia","Sambia","Eggplant","Citrus","Fig","Elephant.grass","Hibiscus","Cotton","Millet","Pine.tree","Yam","Graville","Tiak.tree")] <- 0
  dat_uga[which(!models$field_id %in% fields_w$field_id), c("tot_cassava_area_w", "tot_monoculture_area_w", "tot_intercrop_area_w")] <- 0
  dat_uga[which(is.na(dat_uga$tot_intercrop_area_w)), "tot_intercrop_area_w"] <- 0
  dat_uga[which(is.na(dat_uga$tot_monoculture_area_w)), "tot_monoculture_area_w"] <- 0
  
  dat_uga$tot_cassava_area <- (dat_uga$tot_area_ind_plants + dat_uga$tot_intercrop_area + dat_uga$tot_monoculture_area) / dat_uga$area_m2
  dat_uga$tot_cassava_area_w <- (dat_uga$tot_area_ind_plants + dat_uga$tot_intercrop_area_w*0.75 + dat_uga$tot_monoculture_area_w) / dat_uga$area_m2
  
  # create a Cassava presence/absence category
  dat_uga$cass_area_presence <- factor(dat_uga$tot_cassava_area > 0)
}

rm(list=ls()[which(!(ls()%in%c("civ","uga","dat_civ","dat_uga","panel.cor")))])

# load in Had's buffer statistics
uga_zonalstats_2km <- read.csv("ZonalSt_UGA_2km.csv")
uga_zonalstats_5km <- read.csv("ZonalSt_UGA_5km.csv")
uga_zonalstats_10km <- read.csv("ZonalSt_UGA_10km.csv")
head(uga_zonalstats_10km)

# merge in
dat_uga2km <- merge(dat_uga, uga_zonalstats_2km, by.x="X", by.y="OBJECTID")
# are survey points representative of model predictions
ggplot(dat_uga2km, aes(x=X, y=tot_cassava_area_w*100)) + geom_point() + geom_point(aes(y=Cass_HA), colour="red")
ggplot(dat_uga2km, aes(x=X, y=tot_cassava_area*100)) + geom_point() + geom_point(aes(y=Cass_HA), colour="red")

ggplot(dat_uga2km, aes(x=X, y=tot_cassava_area*100)) + geom_point() + geom_point(aes(y=HarvV1_MIN), colour="red") + geom_point(aes(y=HarvV1_MAX), colour="red")