-
Notifications
You must be signed in to change notification settings - Fork 1
/
time_series_plot.R
146 lines (137 loc) · 7.89 KB
/
time_series_plot.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
# -------------------------------------------------- #
# Climate Risk Profiles -- Time series plot
# A. Esquivel & H. Achicanoy
# Alliance Bioversity-CIAT, 2021
# -------------------------------------------------- #
# Input parameters:
# iso: ISO 3 code in capital letters
# country: country name first character in capital letter
# Output:
# .jpeg graphs per index and season
time_series_plot <- function(country = 'Haiti', iso = 'HTI', seasons){
# Load packages
if(!require(pacman)){install.packages('pacman'); library(pacman)} else {suppressMessages(library(pacman))}
suppressMessages(pacman::p_load(tidyverse,fst))
OSys <<- Sys.info()[1]
root <<- switch(OSys,
'Linux' = '/dapadfs/workspace_cluster_14/WFP_ClimateRiskPr',
'Windows' = '//CATALOGUE/Workspace14/WFP_ClimateRiskPr')
outdir <- paste0(root,'/7.Results/',country,'/results/time_series')
if(!dir.exists(outdir)){ dir.create(outdir, F, T) }
# Load historical time series
pst <- paste0(root,"/7.Results/",country,"/past/",iso,"_indices.fst") %>%
tidyft::parse_fst() %>%
base::as.data.frame()
# Si piden series de tiempo de gSeason, revisar aqui ----
pst$gSeason <- pst$SLGP <- pst$LGP <- NULL
pst$model <- 'Historical'
pst <- pst %>% dplyr::mutate(THI_23 = THI_2 + THI_3, HSI_23 = HSI_2 + HSI_3)
pst_smm <- pst %>% dplyr::group_by(season, year, model) %>% dplyr::summarise_all(median, na.rm = T)
pst_max <- pst %>% dplyr::group_by(season, year, model) %>% dplyr::summarise(NWLD_max = max(NWLD, na.rm = T), NWLD50_max = max(NWLD50, na.rm = T), NWLD90_max = max(NWLD90, na.rm = T))
pst_smm <- dplyr::left_join(x = pst_smm, y = pst_max, by = c('season', 'year', 'model')); rm(pst_max)
pst_smm$id <- pst_smm$x <- pst_smm$y <- NULL
fut_fls <- paste0(root,'/7.Results/',country,'/future') %>%
list.files(., pattern = paste0(iso,'_indices.fst'), full.names = T, recursive = T)
if(length(fut_fls) > 0){
models <- fut_fls %>% strsplit(x = ., split = '/') %>% purrr::map(9) %>% unlist()
fut <- 1:length(fut_fls) %>% purrr::map(.f = function(i){
df <- fut_fls[i] %>%
tidyft::parse_fst() %>%
base::as.data.frame()
df$model <- models[i]; return(df)
}) %>% dplyr::bind_rows()
fut$gSeason <- fut$SLGP <- fut$LGP <- NULL
fut <- fut %>% dplyr::mutate(THI_23 = THI_2 + THI_3, HSI_23 = HSI_2 + HSI_3)
fut_smm <- fut %>% dplyr::group_by(season, year, model) %>% dplyr::summarise_all(median, na.rm = T)
fut_max <- fut %>% dplyr::group_by(season, year, model) %>% dplyr::summarise(NWLD_max = max(NWLD, na.rm = T), NWLD50_max = max(NWLD50, na.rm = T), NWLD90_max = max(NWLD90, na.rm = T))
fut_smm <- dplyr::left_join(x = fut_smm, y = fut_max, by = c('season', 'year', 'model')); rm(fut_max)
fut_smm$id <- fut_smm$x <- fut_smm$y <- NULL
}
ifelse(exists('fut_smm'), tbl <- dplyr::bind_rows(pst_smm,fut_smm), tbl <- pst_smm)
tbl <- tbl %>% tidyr::drop_na()
ss <- sort(unique(tbl$season))
for(s in ss){
dnm <- length(seasons[names(seasons) == s][[1]])
tbl$NDD[tbl$season == s] <- tbl$NDD[tbl$season == s]/dnm
tbl$NT_X[tbl$season == s] <- tbl$NT_X[tbl$season == s]/dnm
tbl$NDWS[tbl$season == s] <- tbl$NDWS[tbl$season == s]/dnm
tbl$NWLD[tbl$season == s] <- tbl$NWLD[tbl$season == s]/dnm
tbl$NWLD50[tbl$season == s] <- tbl$NWLD50[tbl$season == s]/dnm
tbl$NWLD90[tbl$season == s] <- tbl$NWLD90[tbl$season == s]/dnm
tbl$NWLD_max[tbl$season == s] <- tbl$NWLD_max[tbl$season == s]/dnm
tbl$NWLD50_max[tbl$season == s] <- tbl$NWLD50_max[tbl$season == s]/dnm
tbl$NWLD90_max[tbl$season == s] <- tbl$NWLD90_max[tbl$season == s]/dnm
}; rm(ss, s)
# Obtain number of seasons
seasons <- as.character(unique(tbl$season))
for(i in 1:length(seasons)){
dir.create(path = paste0(outdir,'/all_s',i), F, T)
tbl_lng <- tbl %>%
tidyr::pivot_longer(cols = c(TAI:NWLD90_max), names_to = 'Indices', values_to = 'Value') %>%
dplyr::ungroup() %>%
dplyr::group_by(Indices, add = T) %>%
dplyr::group_split()
tbl_lng %>%
purrr::map(.f = function(df){
df <- unique(df)
df <- df[df$season == seasons[i],]
df <- tidyr::drop_na(df)
vr <- df$Indices %>% unique()
if(vr == 'AMT'){ ylb <- expression('Temperature ('*~degree*C*')') }
if(vr == 'ATR'){ ylb <- 'mm/season' }
if(vr == 'TAI'){ ylb <- 'Aridity' }
if(vr %in% c(paste0('HSI_',0:3),'HSI_23',paste0('THI_',0:3),'THI_23')){ ylb <- 'Probability' }
if(vr == 'IRR'){ ylb <- '' }
if(vr %in% c('NDD','NT_X','NDWS','NWLD','NWLD50','NWLD90','NWLD_max','NWLD50_max','NWLD90_max')){ df <- df %>% tidyr::drop_na(); ylb <- 'days/month' }
if(vr == 'P5D'){ ylb <- 'mm/5 days' } # df <- df[df$Value < 125,]
if(vr == 'P95'){ ylb <- 'mm/day' } # df <- df[df$Value < 50,]
if(vr == 'SHI'){ ylb <- 'days/season' }
if(vr == 'CSDI'){ ylb <- 'days' }
if(length(models) > 1){
gg <- df %>%
dplyr::filter(model == 'Historical') %>%
ggplot2::ggplot(aes(x = year, y = Value, group = model)) +
# ggplot2::geom_line(alpha = .05, colour = 'gray') +
ggplot2::theme_bw() +
ggplot2::stat_summary(fun = median, geom = "line", lwd = 1.2, colour = 'black', aes(group=1)) +
ggplot2::stat_summary(data = df[df$model!='Historical',], fun = median, geom = "line", lwd = 1.2, colour = 'black', aes(group=1)) +
ggplot2::stat_summary(data = df[df$model!='Historical',] %>% dplyr::group_by(model), fun = median, geom = "line", lwd = 0.5, aes(colour = model)) +
ggplot2::scale_color_brewer(palette = 'Set1') +
ggplot2::ggtitle(paste0(vr,': ',seasons[i])) +
ggplot2::xlab('Year') +
ggplot2::ylab(ylb) +
ggplot2::geom_vline(xintercept = 2020, size = 1, linetype = "dashed", color = "red") +
ggplot2::theme(axis.text = element_text(size = 17),
axis.title = element_text(size = 20),
legend.text = element_text(size = 15),
legend.title = element_blank(),
plot.title = element_text(size = 25),
plot.subtitle = element_text(size = 17),
strip.text.x = element_text(size = 17),
plot.caption = element_text(size = 15, hjust = 0),
legend.position = "bottom")
ggplot2::ggsave(filename = paste0(outdir,'/all_s',i,'/',vr,'.png'), plot = gg, device = 'jpeg', width = 10, height = 8, units = 'in', dpi = 350)
} else {
gg <- df %>%
ggplot2::ggplot(aes(x = year, y = Value, group = model)) +
ggplot2::geom_line(alpha = .05, colour = 'gray') +
ggplot2::theme_bw() +
ggplot2::stat_summary(fun = median, geom = "line", lwd = 1.2, colour = 'black', aes(group=1)) +
ggplot2::ggtitle(paste0(vr,': ',seasons[i])) +
ggplot2::xlab('Year') +
ggplot2::ylab(ylb) +
ggplot2::geom_vline(xintercept = 2020, size = 1, linetype = "dashed", color = "blue") +
ggplot2::theme(axis.text = element_text(size = 17),
axis.title = element_text(size = 20),
legend.text = element_text(size = 17),
legend.title = element_blank(),
plot.title = element_text(size = 25),
plot.subtitle = element_text(size = 17),
strip.text.x = element_text(size = 17),
plot.caption = element_text(size = 15, hjust = 0),
legend.position = "bottom")
ggplot2::ggsave(filename = paste0(outdir,'/all_s',i,'/',vr,'.png'), plot = gg, device = 'jpeg', width = 10, height = 8, units = 'in', dpi = 350)
}
})
}
}