# Round 2....


# Uganda

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)
}


uga <- ne_countries(scale = "medium", returnclass = "sf", country="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_updated.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

# issue 2. different recorded densities
# look at how many field sampled per site
hist(table(cassava$field_id), breaks=-.5:25.5, col="grey", xlab="Number of separate fileds in Cassava producing sites", main="")
as.data.frame(table(cassava$Density))


# issue 3. individual plants
barplot(table(plants$tot_ind_plants), xlab="Number of individual plants", ylab="Number of fields")
plot(plants$tot_area_ind_plants, plants$tot_monoculture_area)

# issue 4. cropping
table(cassava$Name)
table(cassava$Name, cassava$field_id)
propTable <- as.matrix(table(cassava$Name, cassava$field_id)) / matrix(colSums(as.matrix(table(cassava$Name, cassava$field_id))), nrow=2, ncol=76, byrow=TRUE)
levelplot(as.matrix(table(cassava$Name, cassava$field_id)), asp=1, xlab="", ylab="Field site", scales=list(y=list(lab=NULL)), main="Number of Fields", col.regions=heat.colors(50))
levelplot(propTable, asp=1, xlab="", ylab="Field site", scales=list(y=list(lab=NULL)), main="Proportion of Fields", col.regions=heat.colors(50))
hist(propTable[1,], col="grey", main="", xlab="Proportion of Intercropping Fields")



cassava$weights <- 1
cassave$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$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)
table(dat_uga$cass_area_presence)
table(dat_uga$tot_cassava_area_w > 0)

ggplot(uga) + geom_sf() + geom_point(data=dat_uga, aes(x=Centre_Lon, y=Centre_Lat)) + theme_bw() + labs(x="Longitude", y="Latitude")






ggplot(dat_uga, aes(x=tot_cassava_area, y=tot_cassava_area_w)) + geom_point()

hist(dat_uga$tot_cassava_area_w, col="grey", breaks=30, xlab="Weighted Cassava Production (per m2)", main="")

# Spatial Cassava production ----
ggplot(dat_uga, aes(x=Centre_Lon, y=Centre_Lat, colour=cass_area_presence)) + geom_point(size=2) + labs(x="Longitude", y="Latitude") + theme_bw() + guides(size="none") + coord_fixed()
ggplot(dat_uga, aes(x=Centre_Lon, y=Centre_Lat, size=tot_cassava_area_w, colour=tot_cassava_area_w)) + geom_point() + labs(x="Longitude", y="Latitude") + theme_bw() + guides(colour="none") + coord_fixed()
ggplot(dat_uga, aes(x=Centre_Lon, y=Centre_Lat,  size=tot_cassava_area_w, colour=log(tot_cassava_area_w))) + geom_point() + labs(x="Longitude", y="Latitude") + theme_bw() + guides(size="none") + coord_fixed()
ggplot(dat_uga, aes(x=Centre_Lon, y=Centre_Lat,  size=tot_intercrop_area_w, colour=log(tot_intercrop_area_w))) + geom_point() + labs(x="Longitude", y="Latitude") + theme_bw() + guides(size="none") + coord_fixed()
ggplot(dat_uga, aes(x=Centre_Lon, y=Centre_Lat,  size=tot_monoculture_area_w, colour=log(tot_monoculture_area_w))) + geom_point() + labs(x="Longitude", y="Latitude") + theme_bw() + guides(size="none") + coord_fixed()
ggplot(dat_uga, aes(x=Centre_Lon, y=Centre_Lat,  size=tot_area_ind_plants, colour=log(tot_area_ind_plants))) + geom_point() + labs(x="Longitude", y="Latitude") + theme_bw() + guides(size="none") + coord_fixed()


# Cassava spatial model...
offset <- min(dat_uga$tot_cassava_area_w[dat_uga$tot_cassava_area_w > 0]) / 2
par(mfrow=c(2, 1))
plot(log(dat_uga$tot_cassava_area_w + offset) ~ dat_uga$Centre_Lon, ylab="Log Cassava production", xlab="Longitude")
plot(log(dat_uga$tot_cassava_area_w + offset) ~ dat_uga$Centre_Lat, ylab="Log Cassava production", xlab="Latitude")

mod1 <- lm(tot_cassava_area_w ~ Centre_Lon + Centre_Lat, data = dat_uga)
par(mfrow=c(2,2))
plot(mod1)
# need to deal with the bimodality... perhaps a hurdle model?
summary(mod1)
# or let's first try some higher order spatial terms...
dat_uga$x <- (dat_uga$Centre_Lon - mean(dat_uga$Centre_Lon)) 
dat_uga$y <- (dat_uga$Centre_Lat - mean(dat_uga$Centre_Lat)) 


mod2 <- lm(log(tot_cassava_area_w + offset)~ x * y + I(x^2) + I(y^2), data = dat_uga)
summary(mod2)
plot(mod2)
mods <- step(mod2)
summary(mods)

# how about fitting a gam
gm <- gam(log(tot_cassava_area_w + offset) ~s(Centre_Lon) + s(Centre_Lat), data=dat_uga)
par(mfrow=c(2,2))
gam.check(gm)
par(mfrow=c(1,2))
plot(gm)

gm2 <- gam(log(tot_cassava_area_w + offset) ~s(Centre_Lon) + s(Centre_Lat) + s(Centre_Lon,Centre_Lat), data=dat_uga)
par(mfrow=c(2,2))
gam.check(gm2)
plot(gm2)
# not really enough data in the middle to do this
# let's be smarter and define a factor splitting up the trajectories

ggplot(dat_uga, aes(x=Centre_Lon, y=Centre_Lat, colour=field_id)) + geom_point(size=2) + labs(x="Longitude", y="Latitude") + theme_bw() + guides(size="none") + coord_fixed()

dat_uga$EastWest <- factor(ifelse(dat_uga$field_id > 8 & dat_uga$field_id < 52, "West", "East"))
dat_uga$NorthSouth <- factor(ifelse(dat_uga$Centre_Lat < 1.75, "South", "North"))
table(dat_uga$EastWest)
table(dat_uga$NorthSouth)
ggplot(dat_uga, aes(x=Centre_Lon, y=Centre_Lat, colour=EastWest)) + geom_point(size=2) + labs(x="Longitude", y="Latitude") + theme_bw() + coord_fixed()
ggplot(dat_uga, aes(x=Centre_Lon, y=Centre_Lat, colour=NorthSouth)) + geom_point(size=2) + labs(x="Longitude", y="Latitude") + theme_bw() + coord_fixed()


gm3 <- gam(log(tot_cassava_area_w + offset) ~s(Centre_Lon, by=NorthSouth) + s(Centre_Lat, by=EastWest), data=dat_uga)
par(mfrow=c(2,2))
gam.check(gm3)
par(mfrow=c(2,2))
plot(gm3)


# what about a variogram?
sp_uga <- dat_uga
coordinates(sp_uga) <- ~ Centre_Lon + Centre_Lat
vg_myData <-variogram(tot_cassava_area_w ~ 1, dat_uga)
N <- length(vg_myData$gamma)
totSS <- sum((vg_myData$gamma - mean(vg_myData$gamma))^2)
par(mfrow=c(1,1))
plot(vg_myData$dist, vg_myData$gamma, pch=19, ylab="Semi-variance", xlab="Distance")
modelFits <- NULL
modelAssess <- data.frame(resSS = NA, r2_adj = NA, r2 = NA, rmse = NA)
ranges <- NULL
m_ind <- 1
for (m in c("Exp", "Cir", "Sph")){
  # Fit variogram Model
  fv <- fit.variogram(vg_myData, vgm(,m))
  # ensure nugget and range are positive
  if (sum(fv$range < 0) == 0){
    if (fv[1,2] < 0){
      fv <- fit.variogram(vg_myData, vgm(, m))
      nparam <- sum(is.na(as.data.frame(vgm(,m))))
      ranges[m_ind] <- fv$range
    } else {
      nparam <- sum(is.na(as.data.frame(vgm(m, range))))
      ranges[m_ind] <- fv$range[2]
    }
  } else {
    fv <- fit.variogram(vg_myData, vgm(, m))
    nparam <- sum(is.na(as.data.frame(vgm(,m))))
    ranges[m_ind] <- fv$range
    
  }
  
  # plot fitted model
  lines_v <- variogramLine(fv, maxdist=max(vg_myData$dist))
  lines(lines_v, col=m_ind)

  
  # Assess variogram model
  fitv <- variogramLine(fv, dist_vector=vg_myData$dist)$gamma
  resSS <- sum((vg_myData$gamma - fitv)^2)
  r2adj <- 1 - ((resSS/(N - nparam)) / (totSS/(N - 1)))
  r2 <- cor(vg_myData$gamma, fitv, method="pearson")^2
  rmse <- sqrt(mean((vg_myData$gamma - fitv)^2))
  
  # store models
  modelAssess[m_ind, ] <- c(resSS, r2adj, r2, rmse)
  modelFits[[m_ind]] <- fv
  
  m_ind <- m_ind + 1
  
}
# hmm doesn't look to fit too well....



# Cassava production by region...
names(dat_uga)
# admin_name , region
table(dat_uga$region)
ggplot(dat_uga, aes(x=region, y=log(tot_cassava_area_w +offset), colour=region)) + geom_boxplot()
a1 <- lm(log(tot_cassava_area_w +offset) ~ region, data=dat_uga)
anova(a1)
par(mfrow=c(2,2))
plot(a1)

table(dat_uga$admin_name)



# Cassava production (modelled vs sampled) -----

g1 <- ggplot(dat_uga, aes(x=cass_area_presence, y=log(Cass_HA))) + geom_boxplot()
g2 <- ggplot(dat_uga, aes(x=cass_area_presence, y=log(Cass_Prod))) + geom_boxplot()
g3 <- ggplot(dat_uga, aes(x=cass_area_presence, y=log(SPAM2010))) + geom_boxplot()
g4 <- ggplot(dat_uga, aes(x=cass_area_presence, y=log(MapSPAM_HA))) + geom_boxplot()
ggarrange(g1,g2,g3,g4)

# t.tests
offset <- min(dat_uga$Cass_HA[dat_uga$Cass_HA > 0], na.rm=TRUE)/2
# do a welch t-test
t.test(log(Cass_HA + offset) ~ cass_area_presence, data=dat_uga)

offset <- min(dat_uga$Cass_Prod[dat_uga$Cass_Prod > 0], na.rm=TRUE)/2
t.test(log(Cass_Prod + offset) ~ cass_area_presence, data=dat_uga)

t.test(log(SPAM2010) ~ cass_area_presence, data=dat_uga)
t.test(log(MapSPAM_HA) ~ cass_area_presence, data=dat_uga)

offset <- min(dat_uga$tot_cassava_area_w[dat_uga$tot_cassava_area_w > 0]) / 2
g1 <- ggplot(dat_uga, aes(x=tot_cassava_area_w, y=Cass_HA)) + geom_point()
g2 <- ggplot(dat_uga, aes(x=tot_cassava_area_w, y=Cass_Prod)) + geom_point()
g3 <- ggplot(dat_uga, aes(x=tot_cassava_area_w, y=SPAM2010)) + geom_point()
g4 <- ggplot(dat_uga, aes(x=tot_cassava_area_w, y=MapSPAM_HA)) + geom_point()
ggarrange(g1,g2,g3,g4, nrow=1, ncol=4)

g1 <- ggplot(dat_uga, aes(x=tot_monoculture_area_w, y=Cass_HA)) + geom_point()
g2 <- ggplot(dat_uga, aes(x=tot_monoculture_area_w, y=Cass_Prod)) + geom_point()
g3 <- ggplot(dat_uga, aes(x=tot_monoculture_area_w, y=SPAM2010)) + geom_point()
g4 <- ggplot(dat_uga, aes(x=tot_monoculture_area_w, y=MapSPAM_HA)) + geom_point()
ggarrange(g1,g2,g3,g4, nrow=1, ncol=4)

g1 <- ggplot(dat_uga, aes(x=tot_intercrop_area_w, y=Cass_HA)) + geom_point()
g2 <- ggplot(dat_uga, aes(x=tot_intercrop_area_w, y=Cass_Prod)) + geom_point()
g3 <- ggplot(dat_uga, aes(x=tot_intercrop_area_w, y=SPAM2010)) + geom_point()
g4 <- ggplot(dat_uga, aes(x=tot_intercrop_area_w, y=MapSPAM_HA)) + geom_point()
ggarrange(g1,g2,g3,g4, nrow=1, ncol=4)

g1 <- ggplot(dat_uga, aes(x=tot_area_ind_plants, y=Cass_HA)) + geom_point()
g2 <- ggplot(dat_uga, aes(x=tot_area_ind_plants, y=Cass_Prod)) + geom_point()
g3 <- ggplot(dat_uga, aes(x=tot_area_ind_plants, y=SPAM2010)) + geom_point()
g4 <- ggplot(dat_uga, aes(x=tot_area_ind_plants, y=MapSPAM_HA)) + geom_point()
ggarrange(g1,g2,g3,g4, nrow=1, ncol=4)

offset <- min(dat_uga$Cass_HA[dat_uga$Cass_HA > 0], na.rm=TRUE)/2
gmod1 <- gam(log(Cass_HA + offset) ~s(Centre_Lon, by=NorthSouth) + s(Centre_Lat, by=EastWest), data=dat_uga)
par(mfrow=c(2,2))
gam.check(gmod1)
par(mfrow=c(2,2))
plot(gmod1)
title(main="Cass HA", outer=TRUE, line=-2)

offset <- min(dat_uga$Cass_Prod[dat_uga$Cass_Prod > 0], na.rm=TRUE)/2
gmod2 <- gam(log(Cass_Prod + offset) ~s(Centre_Lon, by=NorthSouth) + s(Centre_Lat, by=EastWest), data=dat_uga)
par(mfrow=c(2,2))
gam.check(gmod2)
par(mfrow=c(2,2))
plot(gmod2)
title(main="Cass Prod", outer=TRUE, line=-2)

gmod3 <- gam(log(SPAM2010) ~s(Centre_Lon, by=NorthSouth) + s(Centre_Lat, by=EastWest), data=dat_uga)
par(mfrow=c(2,2))
gam.check(gmod3)
par(mfrow=c(2,2))
plot(gmod3)
title(main="SPAM 2010", outer=TRUE, line=-2)

gmod4 <- gam(log(MapSPAM_HA) ~s(Centre_Lon, by=NorthSouth) + s(Centre_Lat, by=EastWest), data=dat_uga)
par(mfrow=c(2,2))
gam.check(gmod4)
par(mfrow=c(2,2))
plot(gmod4)
title(main="MapSPAM_HA", outer=TRUE, line=-2)

plot(log(dat_uga$MapSPAM_HA) ~dat_uga$Centre_Lat)


# Cassava modelled production ----
# HarvArea <- read.csv("X:/Kirsty Hassall/Anna Szyniszewska/Cassava/Outputs/CassavaProduction/CassavaMap_HarvArea_v1.csv")
# Prod <- read.csv("X:/Kirsty Hassall/Anna Szyniszewska/Cassava/Outputs/CassavaProduction/CassavaMap_ Prod_v1.csv")
# 
# names(HarvArea)
# 
# ggplot(HarvArea[HarvArea$Country == "Uganda", ], aes(x=Region, y=HA_PerCapita)) + geom_boxplot()
HarvArea <- read.table("EastTrajectory_HarvArea.txt", header=TRUE, sep=",")
head(HarvArea)
Prod <- read.table("EastTrajectory_Prod.txt", header=TRUE, sep=",")
head(Prod)


hist(HarvArea$grid_code)
summary(HarvArea$grid_code)
offh <- min(HarvArea$grid_code[HarvArea$grid_code > 0])/2
offp <- min(Prod$grid_code[Prod$grid_code > 0])/2

plot(log(HarvArea$grid_code + offh) ~ HarvArea$Lat)
plot(log(Prod$grid_code + offp) ~ Prod$Lat)

hmod <- gam(log(grid_code + offh) ~ s(Lat), data=HarvArea)
par(mfrow=c(2,2))
gam.check(hmod)
par(mfrow=c(1,1))
plot(hmod)

pmod <- gam(log(grid_code + offp) ~ s(Lat), data=Prod)
par(mfrow=c(2,2))
gam.check(pmod)
par(mfrow=c(1,1))
plot(pmod)

# Cassava production and population ----
ggplot(dat_uga, aes(x=cass_area_presence, y=log10(Population))) + geom_boxplot()
t.test(log10(Population) ~ cass_area_presence, data=dat_uga)

offset <- min(dat_uga$tot_cassava_area_w[dat_uga$tot_cassava_area_w > 0]) / 2
ggplot(dat_uga, aes(x=log(tot_cassava_area_w + offset), y=log10(Population))) + geom_point()
ggplot(dat_uga, aes(x=(tot_cassava_area_w ), y=(Population))) + geom_point()

ggplot(dat_uga, aes(y=log(tot_cassava_area_w + offset), x=tot_buildings)) + geom_point()


g1 <- ggplot(dat_uga, aes(x=Cass_HA, y=Population)) + geom_point()
g2 <- ggplot(dat_uga, aes(x=Cass_Prod, y=Population)) + geom_point()
g3 <- ggplot(dat_uga, aes(x=MapSPAM_HA, y=Population)) + geom_point()
g4 <- ggplot(dat_uga, aes(x=SPAM2010, y=Population)) + geom_point()
ggarrange(g1,g2,g3,g4)


ggplot(dat_uga, aes(x=tot_buildings, y=Population)) + geom_point()

dat_uga$tot_cassava_area_w_build <- dat_uga$tot_cassava_area_w * dat_uga$area_m2 / (dat_uga$area_m2 - dat_uga$tot_buildings*50)

ggplot(dat_uga, aes(x=tot_cassava_area_w, y=tot_cassava_area_w_build)) + geom_point()

ggplot(dat_uga, aes(x=(tot_cassava_area_w_build ),y=Population)) + geom_point()