# Code to align data inputs
# create weightings
# calculate cassava production area



library(sf)
library("rnaturalearth")
library("rnaturalearthdata")
#library(rgeos)
library(terra)

# variables that could be tuned for various sensitivity analyses
w_ind <- 1
w_intercrop <- 0.75
w_mono <- 1

CIV_densityweights <- data.frame("dense"=1.5, "dense..?"=1.5, "densea"=1.5,"regular/sparse"=0.75,"sparce"=0.5,"sparse"=0.5,"sparse and uprooted"=0.5,"sparse but"=0.5,"sparse2aqqw1"=0.5,"very dense"=1.75,"very sparce"=0.25,"very sparse"=0.25)

UGA_densityweights <- data.frame("dense"=1.5,"dense at places"=1.5,"high density"=1.5,"low / medium density"=0.75,"low density"=0.5,"medium/low density"=0.75,"regular at times sparse"=0.75,"sparse"=0.25,"sparse estimated"=0.25,"very dense"=1.75,"very low density"=0.25,"very sparse"=0.25)

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", "CIV_densityweights", "UGA_densityweights", "w_ind", "w_intercrop", "w_mono")))])
# 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"] <- CIV_densityweights$dense
  cassava$weights[cassava$Desnity == "dense..?"] <- CIV_densityweights$dense...
  cassava$weights[cassava$Desnity == "densea"] <- CIV_densityweights$densea
  cassava$weights[cassava$Desnity == "regular/sparse"] <- CIV_densityweights$regular.sparse
  cassava$weights[cassava$Desnity == "sparce"] <- CIV_densityweights$sparce
  cassava$weights[cassava$Desnity == "sparse"] <- CIV_densityweights$sparse
  cassava$weights[cassava$Desnity == "sparse and uprooted"] <- CIV_densityweights$sparse.and.uprooted
  cassava$weights[cassava$Desnity == "sparse but"] <- CIV_densityweights$sparse.but
  cassava$weights[cassava$Desnity == "sparse2aqqw1"] <- CIV_densityweights$sparse2aqqw1
  cassava$weights[cassava$Desnity == "very dense"] <- CIV_densityweights$very.dense
  cassava$weights[cassava$Desnity == "very sparce"] <- CIV_densityweights$very.sparce
  cassava$weights[cassava$Desnity == "very sparse"] <- CIV_densityweights$very.sparse
  
  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) / areaSampled
  dat_civ$tot_cassava_area_w <- (dat_civ$tot_area_ind_plants*w_ind + dat_civ$tot_intercrop_area_w*w_intercrop + dat_civ$tot_monoculture_area_w*w_mono) / areaSampled
  
  # create a Cassava presence/absence category
  dat_civ$cass_area_presence <- factor(dat_civ$tot_cassava_area > 0)
  
  fields_civ <- cassava
}

rm(list=ls()[which(!(ls()%in%c("civ","uga","dat_civ","dat_uga", "CIV_densityweights", "UGA_densityweights", "w_ind", "w_intercrop", "w_mono", "fields_civ", "fields_uga")))])
# 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(as.factor(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"] <- UGA_densityweights$dense
  cassava$weights[cassava$Desnity == "dense at places"] <- UGA_densityweights$dense.at.places
  cassava$weights[cassava$Desnity == "high density"] <- UGA_densityweights$high.density
  cassava$weights[cassava$Desnity == "low / medium density"] <- UGA_densityweights$low...medium.density
  cassava$weights[cassava$Desnity == "low density"] <- UGA_densityweights$low.density
  cassava$weights[cassava$Desnity == "medium/low density"] <- UGA_densityweights$medium.low.density
  cassava$weights[cassava$Desnity == "regular at times sparse"] <- UGA_densityweights$regular.at.times.sparse
  cassava$weights[cassava$Desnity == "sparse"] <- UGA_densityweights$sparse
  cassava$weights[cassava$Desnity == "sparse estimated"] <- UGA_densityweights$sparse.estimated
  cassava$weights[cassava$Desnity == "very dense"] <- UGA_densityweights$very.dense
  cassava$weights[cassava$Desnity == "very low density"] <- UGA_densityweights$very.low.density
  cassava$weights[cassava$Desnity == "very sparse"] <- UGA_densityweights$very.sparse
  
  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) / areaSampled
  dat_uga$tot_cassava_area_w <- (dat_uga$tot_area_ind_plants*w_ind + dat_uga$tot_intercrop_area_w*w_intercrop + dat_uga$tot_monoculture_area_w*w_mono) / areaSampled
  
  # create a Cassava presence/absence category
  dat_uga$cass_area_presence <- factor(dat_uga$tot_cassava_area > 0)
  
  fields_uga <- cassava
}

rm(list=ls()[which(!(ls()%in%c("CIV_buff_pointstats","UGA_buff_pointstats","civ","uga","dat_civ","dat_uga", "fields_civ", "fields_uga")))])