library(dplyr)
library(ggplot2)
library(tidyr)
library(lbscripts)
library(sprawl)
library(data.table)
library(tibble)
library(sf)
library(raster)
load("/home/lb/Google_Drive/IREA/Conference&paper/PhenoRice_SRV/Datasets/statistics_modis/extract_Phenorice_final_maskirrig_newvalvalley.RData")
load("/home/lb/Google_Drive/IREA/Conference&paper/PhenoRice_SRV/Datasets/statistics_modis/extract_Phenorice_final_maskirrig_newvalvalley_november_1000_valleybig.RData")
load("/home/lb/Google_Drive/IREA/Conference&paper/PhenoRice_SRV/Datasets/statistics_modis/extract_Phenorice_final_maskirrig_newvalvalley_november_1000_valleybig_alldata.RData")
# load("/home/lb/Google_Drive/IREA/Conference&paper/PhenoRice_SRV/Datasets/statistics_modis/extract_Phenorice_srv_nov_17.RData")
gd_folder <- "/home/lb/Google_Drive/IREA/Conference&paper/PhenoRice_SRV/"
validation_delta <- read_vect(
file.path(gd_folder, "Datasets/Pheno_Validation/validation_delta.shp")) %>%
st_transform((proj4string(r))) %>%
dplyr::select(Name)
validation_valley <- read_vect(
file.path(gd_folder, "Datasets/Pheno_Validation/validation_valley.shp")) %>%
st_transform((proj4string(r))) %>%
dplyr::select(id)
inpheno <- read.csv("/home/lb/Google_Drive/IREA/Conference&paper/PhenoRice_SRV/Datasets/Pheno_Validation/Pheno_Validation_data_arrondissements.csv", stringsAsFactors = FALSE)
inpheno <- as_tibble(inpheno) %>%
dplyr::mutate(sowdate = as.Date(sowdate, format = "%d/%m/%y")) %>%
dplyr::mutate(harvdate = as.Date(harvdate, format = "%d/%m/%y")) %>%
dplyr::filter(season == "Hivern") %>%
dplyr::mutate(year = as.numeric(year)) %>%
dplyr::mutate(sowdoy = datetodoy(sowdate, abort = FALSE), harvdoy = datetodoy(harvdate, abort = FALSE))
inpheno$harvdoy[which(inpheno$harvdoy < 200)] <- inpheno$harvdoy[which(inpheno$harvdoy < 200)] + 365
inpheno <- inpheno %>%
dplyr::mutate(sowddate_tmp = doytodate(sowdoy, 2007), harvdate_tmp = doytodate(harvdoy, 2007)) %>%
dplyr::mutate(village = factor(village))
valid_data <- inpheno %>%
group_by(season, year, site) %>%
summarize(sowdate = mean(sowddate_tmp, na.rm = T),
harvdate = mean(harvdate_tmp, na.rm = T),
sow_sd = sd(sowddate_tmp, na.rm = T),
harv_sd = sd(harvdate_tmp, na.rm = T),
obs = n())
valid_data$site = factor(valid_data$site, labels = c("Boundoum (delta)", "Podor (middle valley)"))
valid_data <- valid_data %>%
dplyr::filter(season == "Hivern")
mod_data_delta <- data_delta %>%
dplyr::filter(band_name %in% c("SOS_1", "SOS_2", "EOS_1", "EOS_2")) %>%
dplyr::select(c(3, 5, 6, 8, 10, 11)) %>%
separate(band_name, sep = "_", into = c("Metric", "Season")) %>%
mutate(Season = factor(Season , labels = c("Dry", "Wet"))) %>%
mutate(site = "Boundoum (delta)")
mod_data_valley <- data_valley %>%
dplyr::filter(band_name %in% c("SOS_1", "SOS_2", "EOS_1", "EOS_2")) %>%
dplyr::select(c(3, 5, 6, 8, 10, 11)) %>%
separate(band_name, sep = "_", into = c("Metric", "Season")) %>%
mutate(Season = factor(Season , labels = c("Dry", "Wet"))) %>%
mutate(site = "Podor (middle valley)")
mod_data <- rbind(mod_data_delta, mod_data_valley)
names(mod_data)[4:5] = c("avg_mod", "sd_mod")
# graph for SOS
sos_valid <- valid_data %>%
ungroup() %>%
dplyr::select(year, sowdate, sow_sd, site) %>%
dplyr::rename(avg = sowdate, sd = sow_sd) %>%
dplyr::mutate(avg = doytodate(datetodoy(avg), 2007))
mod_sos <- mod_data %>%
dplyr::filter(Season == "Wet" & Metric == "SOS") %>%
dplyr::select(year, avg_mod, sd_mod, site) %>%
dplyr::mutate(avg_mod = doytodate(avg_mod, 2007))
sos_valid_data <- left_join(sos_valid, mod_sos) %>%
dplyr::mutate(diff = as.numeric(avg - avg_mod)) %>%
dplyr::group_by( site)
lb_linregr(sos_valid_data, "avg", "avg_mod", "site")
# diff <- sos_valid_data$avg - sos_valid_data$avg_mod
sos_valid_data <- sos_valid_data %>%
filter(year > 2002 & year < 2010) %>%
mutate(year = as.numeric(year))
pp <- ggplot(sos_valid_data, aes(x = year - 0.2, y = avg)) +
geom_point(color = "darkblue") +
geom_errorbar(aes(x = year - 0.2, ymin = avg - sd, ymax = avg + sd), color = "darkblue", width = 0.3) +
geom_point(aes(x = year + 0.2, y = avg_mod), color = "darkred") +
geom_errorbar(aes(x = year + 0.2, ymin = avg_mod - sd_mod, ymax = avg_mod + sd_mod), color = "darkred",
, width = 0.3) +
theme_light() +
scale_y_date(limits = c(as.Date("2007-06-23") , as.Date("2007-10-07")), date_labels = "%d-%b") +
scale_x_continuous("Year", limits = c(2002.6, 2009.4),
breaks = 2003:2009) +
scale_color_manual("", values = c("darkblue", "darkred")) +
ggtitle("Observed vs PhenoRice Sowing Dates") +
theme(plot.title = element_text(family = "AvantGarde")) +
facet_wrap(~site) +
ylab("Sowing Date") + theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
ggplot2::theme(panel.grid.minor = element_line(color = "transparent"))
# graph for EOS
eos_valid <- valid_data %>%
ungroup() %>%
dplyr::select(year, harvdate, harv_sd, site) %>%
dplyr::rename(avg = harvdate, sd = harv_sd) %>%
dplyr::mutate(avg = doytodate(datetodoy(avg), 2007)) %>%
mutate(type = "Observed")
mod_eos <- mod_data %>%
dplyr::filter(Season == "Wet" & Metric == "EOS") %>%
dplyr::select(year, avg_mod, sd_mod, site) %>%
dplyr::mutate(avg_mod = doytodate(avg_mod, 2007)) %>%
# mutate(type = "PhenoRice") %>%
mutate(obs = n())
# eos_valid_data <- rbind(eos_valid, mod_eos) %>%
# filter(year > 2002 & year < 2010)
eos_valid_data <- left_join(eos_valid, mod_eos) %>%
dplyr::mutate(diff = as.numeric(avg - avg_mod)) %>%
dplyr::group_by( site)
regr_eos <- lb_linregr(eos_valid_data, "avg", "avg_mod", "site") %>%
dplyr::mutate(Metric = "Harvesting Date") %>%
dplyr::select(site, Metric, slope, intercept, r2, p, ME, RMSE) %>%
dplyr::rename(Site = site, Intercept = intercept, Slope = slope)
regr_sos <- lb_linregr(sos_valid_data, "avg", "avg_mod", "site")%>%
dplyr::mutate(Metric = "Sowing Date") %>%
dplyr::select(site, Metric, slope, intercept, r2, p, ME, RMSE) %>%
dplyr::rename(Site = site, Intercept = intercept, Slope = slope)
regr_table <- rbind(regr_sos, regr_eos) %>%
xtable::xtable(digits = 3) %>%
print(include.rownames = FALSE)
%>%
print(include.rownames = FALSE)
pp2 <- ggplot(eos_valid_data, aes(x = year - 0.2, y = avg)) +
geom_point(color = "darkblue") +
geom_errorbar(aes(x = year - 0.2, ymin = avg - sd, ymax = avg + sd), color = "darkblue", width = 0.3) +
geom_point(aes(x = year + 0.2, y = avg_mod), color = "darkred") +
geom_errorbar(aes(x = year + 0.2, ymin = avg_mod - sd_mod, ymax = avg_mod + sd_mod), color = "darkred",
, width = 0.3) +
theme_light() +
scale_y_date("Harvesting Date", limits = c(as.Date("2007-09-23") , as.Date("2008-02-01")), date_labels = "%d-%b") +
scale_x_continuous("Year", limits = c(2002.6, 2009.4),
breaks = 2003:2009) +
scale_color_manual("", values = c("darkblue", "darkred")) +
ggtitle("Observed vs PhenoRice Harvesting Dates") +
theme(plot.title = element_text(family = "AvantGarde")) +
facet_wrap(~site) +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
ggplot2::theme(panel.grid.minor = element_line(color = "transparent"))
pp <- ggplot(eos_valid_data, aes(x = factor(year), y = avg)) +
geom_point(aes(color = type), position = position_dodge(width = 0.6)) +
geom_errorbar(aes(color = type, ymin = avg - sd, ymax = avg + sd),
position = position_dodge(width = 0.6)) +
theme_light() +
scale_y_date(limits = c(as.Date("2007-09-30") , as.Date("2008-01-30")), date_labels = "%d-%b") +
scale_color_manual("", values = c("darkblue", "darkred")) +
ggtitle("Observed vs PhenoRice Harvesting Dates") +
theme(plot.title = element_text(family = "AvantGarde")) +
facet_wrap(~site) +
xlab("Year") + ylab("Harvesting Date")
ggplot(sos_valid_data, aes(x = avg, y = avg_mod, color = site, label = year)) +
geom_point() + theme_bw() + geom_text(nudge_y = 2) +
geom_smooth(method = "lm") +
geom_hline(yintercept = 0, colour = "grey70", size = 0.5)+
geom_vline(xintercept = 0, colour = "grey70", size = 0.5) +
scale_colour_manual(values = c("darkred", "darkgreen")) +
ylab("PhenoRice average sos") + xlab("SEAD average sos") +
facet_wrap(~site, scales = "free_y") + geom_abline(slope = 1) +
coord_cartesian(xlim = as.Date(c("2007-07-01", "2007-09-15")),
ylim = as.Date(c("2007-07-01", "2007-09-15"))) +
theme(legend.position = "bottom")
ggplot(eos_valid_data, aes(x = avg, y = avg_mod, color = site, label = year)) +
geom_point() + theme_bw() + geom_text(nudge_y = 2) +
geom_smooth(method = "lm") +
geom_hline(yintercept = 0, colour = "grey70", size = 0.5)+
geom_vline(xintercept = 0, colour = "grey70", size = 0.5) +
scale_colour_manual(values = c("darkred", "darkgreen")) +
ylab("PhenoRice average eos") + xlab("SEAD average eos") +
facet_wrap(~site, scales = "free_y") + geom_abline(slope = 1) +
coord_cartesian(xlim = as.Date(c("2007-10-15", "2007-12-31")),
ylim = as.Date(c("2007-10-15", "2007-12-31"))) +
theme(legend.position = "bottom")
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.