# Figures for paper

library(ggplot2)
library(ggpubr)
library(ggcorrplot)
library(mgcv)
library(tidyr)
library(dplyr)

setwd('D:/validating_cassava_spatial_disaggregation_model_in_SSAfrica')


# get aligned data
source('D:/validating_cassava_spatial_disaggregation_model_in_SSAfrica/1_alignData.R')

res_names <- c("tot_cassava_area", "tot_monoculture_area", "tot_intercrop_area", "tot_area_ind_plants")

# ****************************************************
# Figure 4 - Cassava distribution characteristics ----
g1 <- ggplot(dat_civ[dat_civ$tot_cassava_area > 0, ], aes(x=tot_monoculture_area)) + 
  geom_histogram(bins=10) + 
  theme_bw() + 
  labs(x=expression(atop("Total monoculture area (m"^2*")", "per location")), y="Count") + 
  ggtitle("C\U00F4te d'Ivoire") + 
  coord_cartesian(xlim=c(0,30000),ylim=c(0,35))

g2 <- ggplot(dat_civ[dat_civ$tot_cassava_area > 0, ], aes(x=tot_intercrop_area)) + 
  geom_histogram(bins=10) + theme_bw() + 
  labs(x=expression(atop("Total intercropped area (m"^2*")", "per location")), y="Count") + 
  ggtitle(" ")+ 
  coord_cartesian(xlim=c(0,30000),ylim=c(0,35))

civ_corr <- cor(dat_civ[, c(res_names)], use="pairwise", method="spearman")
colnames(civ_corr) <- c("Total", "Monoculture", "Intercrop", "Individual plants")
row.names(civ_corr) <- c("Total", "Monoculture", "Intercrop", "Individual plants")

g3 <- ggcorrplot(civ_corr, type="upper")


g4 <- ggplot(dat_uga[dat_uga$tot_cassava_area > 0, ], aes(x=tot_monoculture_area)) + 
  geom_histogram(bins=10) + theme_bw() + labs(x=expression(atop("Total monoculture area (m"^2*")", "per location")), y="Count") + ggtitle("Uganda") + 
  coord_cartesian(xlim=c(0,30000),ylim=c(0,35))

g5 <- ggplot(dat_uga[dat_uga$tot_cassava_area > 0, ], aes(x=tot_intercrop_area)) + 
  geom_histogram(bins=10) + theme_bw() + labs(x=expression(atop("Total intercropped area (m"^2*")", "per location")), y="Count") + ggtitle(" ") + 
  coord_cartesian(xlim=c(0,30000),ylim=c(0,35))

uga_corr <- cor(dat_uga[, c(res_names)], use="pairwise", method="spearman")
colnames(uga_corr) <- c("Total", "Monoculture", "Intercrop", "Individual plants")
row.names(uga_corr) <- c("Total", "Monoculture", "Intercrop", "Individual plants")

g6 <- ggcorrplot(uga_corr, type="upper")

ggsave("fig4.tiff", ggarrange(g1, g2, g3, g4, g5, g6), scaling=1, width=10)

tiff("fig4.tif", width=800)
#tiff("fig4.tif", width=720*5, height=480*5, res=72*5)
print(
  ggarrange(g1, g2, g3, g4, g5, g6)
)
dev.off()



# ****************************************************
# Figure 5 - Survey vs model summaries ----

dat_civ$cass_area_presence <- as.character(dat_civ$cass_area_presence)
dat_civ$cass_area_presence[dat_civ$cass_area_presence == "TRUE"] <- "Present"
dat_civ$cass_area_presence[dat_civ$cass_area_presence == "FALSE"] <- "Absent"
dat_civ$cass_area_presence <- as.factor(dat_civ$cass_area_presence)

c_civ <- min(dat_civ$Cass_Prod[dat_civ$Cass_Prod > 0], na.rm=TRUE) / 2
cha_civ <- min(dat_civ$Cass_Prod[dat_civ$Cass_HA > 0], na.rm=TRUE) / 2
c_uga <- min(dat_civ$Cass_Prod[dat_uga$Cass_Prod > 0], na.rm=TRUE) / 2
cha_uga <- min(dat_uga$Cass_HA[dat_uga$Cass_HA > 0], na.rm=TRUE) / 2

dat_uga$cass_area_presence <- as.character(dat_uga$cass_area_presence)
dat_uga$cass_area_presence[dat_uga$cass_area_presence == "TRUE"] <- "Present"
dat_uga$cass_area_presence[dat_uga$cass_area_presence == "FALSE"] <- "Absent"
dat_uga$cass_area_presence <- as.factor(dat_uga$cass_area_presence)

g1 <- ggplot(dat_civ, aes(y=log(Cass_Prod + c_civ), x=cass_area_presence)) + geom_boxplot() + labs(y=expression(atop("CassavaMap predicted", "production (log t/km"^2*")")), x="Cassava production") + geom_point() + ggtitle("C\U00F4te d'Ivoire")

g2 <- ggplot(dat_civ, aes(y=log(Cass_HA + cha_civ), x=cass_area_presence)) + geom_boxplot() + labs(y=expression(atop("CassavaMap predicted","harvest area (log ha/km"^2*")")), x="Cassava production") + geom_point() + ggtitle(" ")

g3 <- ggplot(dat_civ, aes(y=log(SPAM2010), x=cass_area_presence)) + geom_boxplot() + labs(y=expression(atop("MapSPAM predicted","production (log t/km"^2*")")), x="Cassava production") + geom_point() + ggtitle(" ")

g4 <- ggplot(dat_uga, aes(y=log(Cass_Prod + c_uga), x=cass_area_presence)) + geom_boxplot() + labs(y=expression(atop("CassavaMap predicted", "production (log t/km"^2*")")), x="Cassava production") + geom_point() + ggtitle("Uganda")

g5 <- ggplot(dat_uga, aes(y=log(Cass_HA + cha_uga), x=cass_area_presence)) + geom_boxplot() + labs(y=expression(atop("CassavaMap predicted","harvest area (log ha/km"^2*")")), x="Cassava production") + geom_point() + ggtitle(" ")

g6 <- ggplot(dat_uga, aes(y=log(SPAM2010), x=cass_area_presence)) + geom_boxplot() + labs(y=expression(atop("MapSPAM predicted","production (log t/km"^2*")")), x="Cassava production") + geom_point() + ggtitle(" ")

tiff("fig5.tif", width=800)
#tiff("fig5.tif",width=720*5, height=480*5, res=72*5)
print(
  ggarrange(g1, g2, g3, g4, g5, g6, nrow=2, ncol=3)
)
dev.off()


# ******************************************************
# Figure 6 - Spatial trends in Ivory Coast --------------

# ggplot(civ) + geom_sf() + geom_point(data=dat_civ, aes(x=Centre_Lon, y=Centre_Lat)) + theme_bw() + labs(x="Longitude", y="Latitude")

civ_admin <- read_sf("D://validating_cassava_spatial_disaggregation_model_in_SSAfrica//CIV_GADM_Adm1_Adm2//CIV_adm2.shp")


# ggplot(civ_admin) + geom_sf() + geom_point(data=dat_civ, aes(x=Centre_Lon, y=Centre_Lat)) + theme_bw() + labs(x="Longitude", y="Latitude") 

# ggplot(civ_admin) + geom_sf() + geom_point(data=dat_civ, aes(x=Centre_Lon, y=Centre_Lat)) + theme_bw() + labs(x="Longitude", y="Latitude") + geom_point(data = dat_civ, aes(x=Centre_Lon, y=Centre_Lat, size=tot_cassava_area_w, colour=tot_cassava_area_w)) +  labs(x="Longitude", y="Latitude") + theme_bw() + guides(size="none", colour=guide_legend(title="Total Cassava Production\n (weighted per m2)"))


g1 <- ggplot(civ_admin) + geom_sf() + geom_point(data=dat_civ, aes(x=Centre_Lon, y=Centre_Lat)) +
  theme_bw() + labs(x="Longitude", y="Latitude") +
  geom_point(data = dat_civ, aes(x=Centre_Lon, y=Centre_Lat, size=tot_cassava_area_w, colour=tot_cassava_area_w)) +
  labs(x="Longitude", y="Latitude") +
  theme_bw() +
  guides(size="none", colour=guide_legend(title=expression(atop("Weighted cassava", "production (log t/km"^2*")")))) +
  xlim(-6, -3) + ylim(5, 8)

offset <- min(dat_civ$tot_cassava_area_w[dat_civ$tot_cassava_area_w > 0]) / 2
gm <- gam(log(tot_cassava_area_w + offset) ~s(Centre_Lon) + s(Centre_Lat), data=dat_civ)

newdata1 <- data.frame(Centre_Lon = seq(from=min(dat_civ$Centre_Lon), to=max(dat_civ$Centre_Lon), length=100), Centre_Lat = rep(mean(dat_civ$Centre_Lat), 100))

newdata2 <- data.frame(Centre_Lat = seq(from=min(dat_civ$Centre_Lat), to=max(dat_civ$Centre_Lat), length=100), Centre_Lon = rep(mean(dat_civ$Centre_Lon), 100))

pred1 <- predict(gm, newdata1, se=TRUE)
newdata1$log_cassava <- pred1$fit
newdata1$lower <- pred1$fit - 1.96*pred1$se.fit
newdata1$upper <- pred1$fit + 1.96*pred1$se.fit

pred2 <- predict(gm, newdata2, se=TRUE)
newdata2$log_cassava <- pred2$fit
newdata2$lower <- pred2$fit - 1.96*pred2$se.fit
newdata2$upper <- pred2$fit + 1.96*pred2$se.fit

g2 <- ggplot(newdata1, aes(x=Centre_Lon, y=log_cassava)) +
  geom_ribbon(data = newdata1, alpha = 0.3,
              aes(ymin = lower, ymax = upper, fill = "confidence interval")) +
  geom_line(data = newdata1, aes(color = "GAM")) +
  scale_fill_manual(values = "lightblue", name = NULL) +
  scale_color_manual(values = "darkblue", name = NULL) +
  theme_bw() + labs(x= "Longitude", y=expression(atop("Cassava production (log t/km"^2*")"))) +
  guides(fill="none", colour="none")

g3 <- ggplot(newdata1, aes(x=Centre_Lat, y=log_cassava)) +
  geom_ribbon(data = newdata2, alpha = 0.3,
              aes(ymin = lower, ymax = upper, fill = "confidence interval")) +
  geom_line(data = newdata2, aes(color = "GAM")) +
  scale_fill_manual(values = "lightblue", name = NULL) +
  scale_color_manual(values = "darkblue", name = NULL) +
  theme_bw() +
  labs(x= "Latitude", y=expression(atop("Cassava production (log t/km"^2*")"))) +
  guides(fill="none", colour="none")


civ_prod_10000 <- read.csv("D:/Cassava_files/Kirsty/Stats/raster_to_points/civ_prod_10000.txt")

offset <- min(civ_prod_10000$GRID_CODE[civ_prod_10000$GRID_CODE > 0], na.rm=TRUE) / 2
gm_prod <- gam(log(GRID_CODE + offset) ~ s(Longitude) + s(Latitude), data=civ_prod_10000)

g4 <-  ggplot(civ_prod_10000, aes(x=Longitude, y=Latitude, colour=GRID_CODE)) +
  geom_point() +
  theme_bw() +
  guides(colour=guide_legend(title=expression(atop("CassavaMap predicted", "production (log t/km"^2*")")))) +
  coord_fixed()

newdata1 <- data.frame(Longitude = seq(from=min(civ_prod_10000$Longitude), to=max(civ_prod_10000$Longitude), length=100), Latitude = rep(mean(civ_prod_10000$Latitude), 100))

newdata2 <- data.frame(Latitude = seq(from=min(civ_prod_10000$Latitude), to=max(civ_prod_10000$Latitude), length=100), Longitude = rep(mean(civ_prod_10000$Longitude), 100))

pred1 <- predict(gm_prod, newdata1, se=TRUE)
newdata1$log_cassava <- pred1$fit
newdata1$lower <- pred1$fit - 1.96*pred1$se.fit
newdata1$upper <- pred1$fit + 1.96*pred1$se.fit

pred2 <- predict(gm_prod, newdata2, se=TRUE)
newdata2$log_cassava <- pred2$fit
newdata2$lower <- pred2$fit - 1.96*pred2$se.fit
newdata2$upper <- pred2$fit + 1.96*pred2$se.fit

g5 <- ggplot(newdata1, aes(x=Longitude, y=log_cassava)) +
  geom_ribbon(data = newdata1, alpha = 0.3,
              aes(ymin = lower, ymax = upper, fill = "confidence interval")) +
  geom_line(data = newdata1, aes(color = "GAM")) +
  scale_fill_manual(values = "lightblue", name = NULL) +
  scale_color_manual(values = "darkblue", name = NULL) +
  theme_bw() +
  labs(x= "Longitude", y=expression(atop("CassavaMap production (log t/km"^2*")"))) +
  guides(fill="none", colour="none")

g6 <- ggplot(newdata1, aes(x=Latitude, y=log_cassava)) +
  geom_ribbon(data = newdata2, alpha = 0.3,
              aes(ymin = lower, ymax = upper, fill = "confidence interval")) +
  geom_line(data = newdata2, aes(color = "GAM")) +
  scale_fill_manual(values = "lightblue", name = NULL) +
  scale_color_manual(values = "darkblue", name = NULL) +
  theme_bw() +
  labs(x= "Latitude", y=expression(atop("CassavaMap production (log t/km"^2*")"))) +
  guides(fill="none", colour="none")

#tiff("fig6.tif",width=880*5, height=480*5, res=72*5)
tiff("fig6.tif", width=800)
print(
  ggarrange(g1, g2, g3, g4, g5, g6, nrow=2, ncol=3, widths=c(1.5, 1, 1))
)
dev.off()

# ******************************************************
# Figure 7 - Spatial trends in Uganda --------------


uga_admin <- read_sf("D://validating_cassava_spatial_disaggregation_model_in_SSAfrica//Uganda_2014-Regions//Uganda_regions_2014.shp")

# gs1 <- ggplot(uga_admin) + geom_sf() + geom_point(data=dat_uga, aes(x=Centre_Lon, y=Centre_Lat)) + theme_bw() + labs(x="Longitude", y="Latitude") + geom_point(data = dat_uga, aes(x=Centre_Lon, y=Centre_Lat, size=tot_cassava_area_w, colour=tot_cassava_area_w)) +  labs(x="Longitude", y="Latitude") + theme_bw() + guides(size="none", colour=guide_legend(title="Total Cassava ")) + xlim(30, 35) + ylim(0,4)

gs2 <- ggplot(uga_admin) +
  geom_sf() +
  geom_point(data=dat_uga, aes(x=Centre_Lon, y=Centre_Lat)) +
  theme_bw() + labs(x="Longitude", y="Latitude") +
  geom_point(data = dat_uga, aes(x=Centre_Lon, y=Centre_Lat, size=tot_monoculture_area, colour=tot_monoculture_area)) +
  labs(x="Longitude", y="Latitude") + theme_bw() + guides(size="none", colour=guide_legend(title="Area in monoculture")) +
  xlim(30, 35) +
  ylim(0,4)

gs3 <- ggplot(uga_admin) +
  geom_sf() +
  geom_point(data=dat_uga, aes(x=Centre_Lon, y=Centre_Lat)) +
  theme_bw() + labs(x="Longitude", y="Latitude") +
  geom_point(data = dat_uga, aes(x=Centre_Lon, y=Centre_Lat, size=tot_intercrop_area, colour=tot_intercrop_area)) +
  labs(x="Longitude", y="Latitude") +
  theme_bw() +
  guides(size="none", colour=guide_legend(title="Area in intercropped")) +
  xlim(30, 35) +
  ylim(0,4)

ggarrange(gs2, gs3, nrow=1)


g1 <- ggplot(uga_admin) +
  geom_sf() +
  geom_point(data=dat_uga, aes(x=Centre_Lon, y=Centre_Lat)) +
  theme_bw() +
  labs(x="Longitude", y="Latitude") +
  geom_point(data = dat_uga, aes(x=Centre_Lon, y=Centre_Lat, size=tot_cassava_area_w, colour=tot_cassava_area_w)) +
  labs(x="Longitude", y="Latitude") +
  theme_bw() +
  guides(size="none", colour=guide_legend(title=expression(atop("Weighted cassava", "production (log t/km"^2*")")))) +
  xlim(30, 35) +
  ylim(0,4)

offset <- min(dat_uga$tot_cassava_area_w[dat_uga$tot_cassava_area_w > 0]) / 2
gm <- gam(log(tot_cassava_area_w + offset) ~s(Centre_Lon) + s(Centre_Lat), data=dat_uga)

newdata1 <- data.frame(Centre_Lon = seq(from=min(dat_uga$Centre_Lon), to=max(dat_uga$Centre_Lon), length=100), Centre_Lat = rep(mean(dat_uga$Centre_Lat), 100))

newdata2 <- data.frame(Centre_Lat = seq(from=min(dat_uga$Centre_Lat), to=max(dat_uga$Centre_Lat), length=100), Centre_Lon = rep(mean(dat_uga$Centre_Lon), 100))

pred1 <- predict(gm, newdata1, se=TRUE)
newdata1$log_cassava <- pred1$fit
newdata1$lower <- pred1$fit - 1.96*pred1$se.fit
newdata1$upper <- pred1$fit + 1.96*pred1$se.fit

pred2 <- predict(gm, newdata2, se=TRUE)
newdata2$log_cassava <- pred2$fit
newdata2$lower <- pred2$fit - 1.96*pred2$se.fit
newdata2$upper <- pred2$fit + 1.96*pred2$se.fit

g2 <- ggplot(newdata1, aes(x=Centre_Lon, y=log_cassava)) +
  geom_ribbon(data = newdata1, alpha = 0.3,
              aes(ymin = lower, ymax = upper, fill = "confidence interval")) +
  geom_line(data = newdata1, aes(color = "GAM")) +
  scale_fill_manual(values = "lightblue", name = NULL) +
  scale_color_manual(values = "darkblue", name = NULL) +
  theme_bw() + labs(x= "Longitude", y=expression(atop("Cassava production (log t/km"^2*")"))) +
  guides(fill="none", colour="none")

g3 <- ggplot(newdata1, aes(x=Centre_Lat, y=log_cassava)) +
  geom_ribbon(data = newdata2, alpha = 0.3,
              aes(ymin = lower, ymax = upper, fill = "confidence interval")) +
  geom_line(data = newdata2, aes(color = "GAM")) +
  scale_fill_manual(values = "lightblue", name = NULL) +
  scale_color_manual(values = "darkblue", name = NULL) +
  theme_bw() + labs(x= "Latitude", y=expression(atop("Cassava production (log t/km"^2*")"))) +
  guides(fill="none", colour="none")

uga_prod_10000 <- read.csv("D:/Cassava_files/Kirsty/Stats/raster_to_points/uga_prod_10000.txt")

offset <- min(uga_prod_10000$GRID_CODE[uga_prod_10000$GRID_CODE > 0], na.rm=TRUE) / 2
gm_prod <- gam(log(GRID_CODE + offset) ~ s(Longitude) + s(Latitude), data=uga_prod_10000)

g4 <-  ggplot(uga_prod_10000, aes(x=Longitude, y=Latitude, colour=GRID_CODE)) +
  geom_point() +
  theme_bw() +
  guides(colour=guide_legend(title=expression(atop("CassavaMap predicted", "production (log t/km"^2*")")))) +
  coord_fixed()

newdata1 <- data.frame(Longitude = seq(from=min(uga_prod_10000$Longitude), to=max(uga_prod_10000$Longitude), length=100), Latitude = rep(mean(uga_prod_10000$Latitude), 100))

newdata2 <- data.frame(Latitude = seq(from=min(uga_prod_10000$Latitude), to=max(uga_prod_10000$Latitude), length=100), Longitude = rep(mean(uga_prod_10000$Longitude), 100))

pred1 <- predict(gm_prod, newdata1, se=TRUE)
newdata1$log_cassava <- pred1$fit
newdata1$lower <- pred1$fit - 1.96*pred1$se.fit
newdata1$upper <- pred1$fit + 1.96*pred1$se.fit

pred2 <- predict(gm_prod, newdata2, se=TRUE)
newdata2$log_cassava <- pred2$fit
newdata2$lower <- pred2$fit - 1.96*pred2$se.fit
newdata2$upper <- pred2$fit + 1.96*pred2$se.fit

g5 <- ggplot(newdata1, aes(x=Longitude, y=log_cassava)) +
  geom_ribbon(data = newdata1, alpha = 0.3,
              aes(ymin = lower, ymax = upper, fill = "confidence interval")) +
  geom_line(data = newdata1, aes(color = "GAM")) +
  scale_fill_manual(values = "lightblue", name = NULL) +
  scale_color_manual(values = "darkblue", name = NULL) +
  theme_bw() + labs(x= "Longitude", y=expression(atop("CassavaMap production (log t/km"^2*")"))) +
  guides(fill="none", colour="none")

g6 <- ggplot(newdata1, aes(x=Latitude, y=log_cassava)) +
  geom_ribbon(data = newdata2, alpha = 0.3,
              aes(ymin = lower, ymax = upper, fill = "confidence interval")) +
  geom_line(data = newdata2, aes(color = "GAM")) +
  scale_fill_manual(values = "lightblue", name = NULL) +
  scale_color_manual(values = "darkblue", name = NULL) +
  theme_bw() + labs(x= "Latitude", y=expression(atop("CassavaMap production (log t/km"^2*")"))) + guides(fill="none", colour="none")

tiff("fig7.tif", width=800)
print(
  ggarrange(g1, g2, g3, g4, g5, g6, nrow=2, ncol=3, widths=c(1.4,1,1))
)
dev.off()


# ******************************************************
# Figure 8 - Regression model outputs --------------





# given main effect is gam, let's only plot the gam models
load("3_modelVsSurvey_output.rdata")

# some summary figures of the models fitted
getSummaries <- function(r1){
  r1$cass_type <- NA
  r1$cass_type[grep("Harv", r1$cass_model)] <- "Harv"
  r1$cass_type[grep("HA", r1$cass_model)] <- "Harv"
  r1$cass_type[grep("Prod", r1$cass_model)] <- "Prod"
  table(r1$cass_type, useNA="ifany")
  r1$cass_summary <- sapply(r1$cass_model, FUN=function(x)strsplit(as.character(x), "_", fixed=TRUE)[[1]][1])
  r1$cass_summary[r1$cass_summary == "Cass"] <- "mean" # mean of a single point = point estimate
  table(r1$cass_summary, useNA="ifany")
  r1$cass_dist <- NA
  r1$cass_dist[grep("5", r1$cass_model)] <- "5km"
  r1$cass_dist[grep("2", r1$cass_model)] <- "2km"
  r1$cass_dist[grep("10", r1$cass_model)] <- "10km"
  r1$cass_dist[is.na(r1$cass_dist)] <- "0km"
  table(r1$cass_dist, useNA="ifany")
  # table(r1$cass_dist,r1$cass_model, useNA="ifany")
  
  r1$population_type <- NA
  r1$population_type[grep("not included", r1$population)] <- "not included"
  r1$population_type[grep("Population", r1$population)] <- "population"
  r1$population_type[grep("lspop", r1$population)] <- "population"
  table(r1$population_type, useNA="ifany")
  r1$population_summary <- sapply(r1$population, FUN=function(x)strsplit(as.character(x), "_", fixed=TRUE)[[1]][1])
  r1$population_summary[r1$population_summary %in% c("Population")] <- "mean" # mean of a single point = point estimate
  table(r1$population_summary, useNA="ifany")
  r1$population_dist <- NA
  r1$population_dist[grep("5", r1$population)] <- "5km"
  r1$population_dist[grep("2", r1$population)] <- "2km"
  r1$population_dist[grep("10", r1$population)] <- "10km"
  r1$population_dist[is.na(r1$population_dist)] <- "0km"
  r1$population_dist[r1$population == "not included"] <- "not included"
  table(r1$population_dist, useNA="ifany")
  
  r1$settlement_type <- NA
  r1$settlement_type[grep("not included", r1$settlement)] <- "not included"
  r1$settlement_type[grep("settle", r1$settlement)] <- "settlement"
  table(r1$settlement_type, useNA="ifany")
  r1$settlement_summary <- sapply(r1$settlement, FUN=function(x)strsplit(as.character(x), "_", fixed=TRUE)[[1]][1])
  table(r1$settlement_summary, useNA="ifany")
  r1$settlement_dist <- NA
  r1$settlement_dist[grep("5", r1$settlement)] <- "5km"
  r1$settlement_dist[grep("2", r1$settlement)] <- "2km"
  r1$settlement_dist[grep("10", r1$settlement)] <- "10km"
  r1$settlement_dist[is.na(r1$settlement_dist)] <- "0km"
  r1$settlement_dist[r1$settlement == "not included"] <- "not included"
  table(r1$settlement_dist, useNA="ifany")
  
  r1$modeltype <- factor(r1$modeltype)
  r1$cass_type <- factor(r1$cass_type)
  r1$cass_summary <- factor(r1$cass_summary)
  r1$cass_dist <- factor(r1$cass_dist)
  r1$population_type <- factor(r1$population_type)
  r1$population_summary <- factor(r1$population_summary)
  r1$population_dist <- factor(r1$population_dist)
  r1$settlement_type <- factor(r1$settlement_type)
  r1$settlement_summary <- factor(r1$settlement_summary)
  r1$settlement_dist <- factor(r1$settlement_dist)
  
  return(r1)
}

r1_CIV <- getSummaries(r1_CIV)
r2_CIV <- getSummaries(r2_CIV)
r3_CIV <- getSummaries(r3_CIV)
r4_CIV <- getSummaries(r4_CIV)
r5_CIV <- getSummaries(r5_CIV)
r6_CIV <- getSummaries(r6_CIV)
r7_CIV <- getSummaries(r7_CIV)

r1_CIV_pres <- getSummaries(r1_CIV_pres)
r2_CIV_pres <- getSummaries(r2_CIV_pres)
r3_CIV_pres <- getSummaries(r3_CIV_pres)
r4_CIV_pres <- getSummaries(r4_CIV_pres)
r5_CIV_pres <- getSummaries(r5_CIV_pres)
r6_CIV_pres <- getSummaries(r6_CIV_pres)
r7_CIV_pres <- getSummaries(r7_CIV_pres)

r1_UGA <- getSummaries(r1_UGA)
r2_UGA <- getSummaries(r2_UGA)
r3_UGA <- getSummaries(r3_UGA)
r4_UGA <- getSummaries(r4_UGA)
r5_UGA <- getSummaries(r5_UGA)
r6_UGA <- getSummaries(r6_UGA)
r7_UGA <- getSummaries(r7_UGA)

r1_UGA_pres <- getSummaries(r1_UGA_pres)
r2_UGA_pres <- getSummaries(r2_UGA_pres)
r3_UGA_pres <- getSummaries(r3_UGA_pres)
r4_UGA_pres <- getSummaries(r4_UGA_pres)
r5_UGA_pres <- getSummaries(r5_UGA_pres)
r6_UGA_pres <- getSummaries(r6_UGA_pres)
r7_UGA_pres <- getSummaries(r7_UGA_pres)


r1_CIV$cass_dist <- factor(r1_CIV$cass_dist, levels=c("not included","0km", "2km", "5km", "10km"))
r1_UGA$cass_dist <- factor(r1_UGA$cass_dist, levels=c("not included","0km", "2km", "5km", "10km"))

r1_CIV$settlement_dist <- factor(r1_CIV$settlement_dist, levels=c("not included","0km", "2km", "5km", "10km"))
r1_UGA$settlement_dist <- factor(r1_UGA$settlement_dist, levels=c("not included","0km", "2km", "5km", "10km"))
r1_CIV$population_dist <- factor(r1_CIV$population_dist, levels=c("not included","0km", "2km", "5km", "10km"))
r1_UGA$population_dist <- factor(r1_UGA$population_dist, levels=c("not included","0km", "2km", "5km", "10km"))

g1 <- ggplot(r1_CIV %>% filter(modeltype == "gam"), aes(x=cass_type, y=AIC, colour=cass_dist)) + geom_boxplot() + ggtitle("C\U00F4te d'Ivoire") + guides(colour=guide_legend(title="Distance")) +  scale_colour_discrete(drop = FALSE)

# ggplot(r1_CIV %>% filter(modeltype == "gam"), aes(x=cass_summary, y=AIC, colour=cass_summary)) + geom_boxplot() + ggtitle("Ivory Coast")

# ggplot(r1_CIV %>% filter(modeltype == "gam"), aes(x=population_summary, y=AIC, colour=population_summary)) + geom_boxplot() + ggtitle("Ivory Coast")

g2 <- ggplot(r1_CIV %>% filter(modeltype == "gam"), aes(x=settlement_type, y=AIC, colour=settlement_dist)) + geom_boxplot() + ggtitle(" ") + guides(colour=guide_legend(title="Distance")) +   scale_colour_discrete(drop = FALSE)

g3 <- ggplot(r1_CIV %>% filter(modeltype == "gam"), aes(x=population_type, y=AIC, colour=population_dist)) + geom_boxplot() + ggtitle(" ") + guides(colour=guide_legend(title="Distance")) +   scale_colour_discrete(drop = FALSE)

g4 <- ggplot(r1_UGA %>% filter(modeltype == "gam"), aes(x=cass_type, y=AIC, colour=cass_dist)) + geom_boxplot() + ggtitle("Uganda") + guides(colour=guide_legend(title="Distance")) +   scale_colour_discrete(drop = FALSE)

# ggplot(r1_UGA %>% filter(modeltype == "gam"), aes(x=cass_summary, y=AIC, colour=cass_summary)) + geom_boxplot() + ggtitle("Uganda")

# ggplot(r1_UGA %>% filter(modeltype == "gam"), aes(x=population_summary, y=AIC, colour=population_summary)) + geom_boxplot() + ggtitle("Uganda")

g5 <- ggplot(r1_UGA %>% filter(modeltype == "gam"), aes(x=settlement_type, y=AIC, colour=settlement_dist)) + geom_boxplot() + ggtitle(" ") + guides(colour=guide_legend(title="Distance")) +   scale_colour_discrete(drop = FALSE)

g6 <- ggplot(r1_UGA %>% filter(modeltype == "gam"), aes(x=population_type, y=AIC, colour=population_dist)) + geom_boxplot() + ggtitle(" ") + guides(colour=guide_legend(title="Distance")) +   scale_colour_discrete(drop = FALSE)

tiff("fig8.tif", width=720)
print(
  ggarrange(g1, g2, g3, g4, g5, g6, nrow=2, ncol=3, common.legend=TRUE, legend="right")
)
dev.off()

# Figure 9 - Fitted model -----------------------

source('D:/validating_cassava_spatial_disaggregation_model_in_SSAfrica/2_alignRasterOutputs.R')

dat_civ$Centre_Lon <- round(dat_civ$Centre_Lon, 7)
dat_civ$Centre_Lat <- round(dat_civ$Centre_Lat, 7)
CIV_buff_pointstats$Centre_Lon <- round(CIV_buff_pointstats$Centre_Lon, 7)
CIV_buff_pointstats$Centre_Lat <- round(CIV_buff_pointstats$Centre_Lat, 7)
CIVdata <- merge(dat_civ, CIV_buff_pointstats, by=c("Centre_Lon", "Centre_Lat"), all.x=TRUE, all.y=TRUE)

dat_uga$Centre_Lon <- round(dat_uga$Centre_Lon, 7)
dat_uga$Centre_Lat <- round(dat_uga$Centre_Lat, 7)
UGA_buff_pointstats$Centre_Lon <- round(UGA_buff_pointstats$Centre_Lon, 7)
UGA_buff_pointstats$Centre_Lat <- round(UGA_buff_pointstats$Centre_Lat, 7)
UGAdata <- merge(dat_uga, UGA_buff_pointstats, by=c("Centre_Lon", "Centre_Lat"), all.x=TRUE, all.y=TRUE)


mod_CIV <- gam(tot_cassava_area ~ s(mean_Prod.10000,min_lspop.2000)+s(mean_settle.2000), data=CIVdata)

mod_UGA <- gam(tot_cassava_area ~ s(Cass_HA,median_lspop.5000)+s(mean_settle.2000), data=UGAdata)

par(mfrow=c(2,2), mar=c(4, 8, 4, 1))
plot(mod_CIV, select=1, scheme=2, ylim=c(0, 30), xlim=c(0, 200), xlab="CassavaMap production (mean at 10km)", ylab="Population density\n(minimum at 2km)", main="C\U00F4te d'Ivoire")
plot(mod_CIV, select=2, main="",  xlim=c(0, 0.15), xlab="Settlement density (2km)", ylab="Centred spline")

plot(mod_UGA, select=1, scheme=2, ylim=c(0, 250), xlim=c(0, 40), ylab="Population Density\n(Median at 5km)", xlab="CassavaMap Harvest area (point location)", main="Uganda")
plot(mod_UGA, select=2, main="",  xlim=c(0, 0.2), xlab="Settlement density (2km)", ylab="Centred spline")