-
Notifications
You must be signed in to change notification settings - Fork 1
/
index_climatology_plot_CSDI.R
285 lines (232 loc) · 11.9 KB
/
index_climatology_plot_CSDI.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
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
# WFP Climate Risk Project
# =----------------------
# Bar graphs
# A. Esquivel - H. Achicanoy - C.Saavedra
# =----------------------
bar_graphs <- function(country, iso, region = 'all'){
path_save <- glue::glue('{root}/7.Results/{country}/results/monthly_g/')
dir.create(path_save, recursive = TRUE)
# Historic indices
h0 <- paste0(root,'/7.Results/',country,'/past/',iso,'_indices_monthly.fst') %>%
tidyft::parse_fst(.) %>%
base::as.data.frame()
h0$model <- 'Historical'
h0$season <- gsub(pattern = 's', replacement = '', h0$season)
h0$period <- '1981-2019'
h0$time <- 'Historical'
ff <- list.files(path = paste0(root,'/7.Results/',country,'/future/'), pattern = '.fst$', full.names = T, recursive = T) %>%
grep(pattern = '_indices_monthly.fst$', x = ., value = T) %>%
grep(pattern = iso, x = ., value = T)
if(length(ff) > 0){
f1 <- ff %>%
grep(pattern = '2021-2040', x = ., value = T) %>%
purrr::map(.f = function(m){
tb <- m %>%
tidyft::parse_fst(path = .) %>%
base::as.data.frame()
tb$model <- m %>% strsplit(x = ., split = '/') %>% purrr::map(10) %>% unlist()
tb$period <- '2021-2040'
tb$season <- gsub(pattern = 's', replacement = '', tb$season)
tb$time <- 'Future'
return(tb)
}) %>%
dplyr::bind_rows()
f2 <- ff %>%
grep(pattern = '2041-2060', x = ., value = T) %>%
purrr::map(.f = function(m){
tb <- m %>%
tidyft::parse_fst(path = .) %>%
base::as.data.frame()
tb$model <- m %>% strsplit(x = ., split = '/') %>% purrr::map(10) %>% unlist()
tb$period <- '2041-2060'
tb$season <- gsub(pattern = 's', replacement = '', tb$season)
tb$time <- 'Future'
return(tb)
}) %>%
dplyr::bind_rows()
}
# Identify intersected pixels
px <- base::intersect(h0$id, f1$id)
h0 <- h0[h0$id %in% px,]
f1 <- f1[f1$id %in% px,]
f2 <- f2[f2$id %in% px,]
all <- dplyr::bind_rows(h0,f1,f2)
# -------------------------------------------------------------------------- #
to_do <- readxl::read_excel(glue::glue('{root}/1.Data/regions_ind.xlsx') )%>%
dplyr::filter(ISO3 == iso) %>%
dplyr::rename('Livehood_z' = 'Livelihood zones', 'NT_X'= "NT-X")
# region <- to_do$Regions[1]
if(region == 'all'){
if('CSDI' %in% names(all)){
var_s <- to_do %>% dplyr::mutate( Regions = 'all', Livehood_z = 'all', Short_Name = 'all') %>%
dplyr::mutate_at(.vars = vars(ATR:CSDI) , .funs = function(x){x <- ifelse(x == '-', 0, x) %>% as.integer()}) %>%
dplyr::group_by(ISO3, Country, Regions, Livehood_z, Short_Name ) %>%
dplyr::summarise_all(. , sum, na.rm = TRUE) %>% ungroup()
} else {
var_s <- to_do %>% dplyr::mutate( Regions = 'all', Livehood_z = 'all', Short_Name = 'all') %>%
dplyr::mutate_at(.vars = vars(ATR:SHI) , .funs = function(x){x <- ifelse(x == '-', 0, x) %>% as.integer()}) %>%
dplyr::group_by(ISO3, Country, Regions, Livehood_z, Short_Name ) %>%
dplyr::summarise_all(. , sum, na.rm = TRUE) %>% ungroup()
}
title = 'Country'
}else{
if('CSDI' %in% names(all)){
var_s <- to_do %>% dplyr::filter(Regions == region) %>%
dplyr::mutate_at(.vars = vars(ATR:CSDI) , .funs = function(x){x <- ifelse(x == '-', 0, x) %>% as.integer()})
} else{
var_s <- to_do %>% dplyr::filter(Regions == region) %>%
dplyr::mutate_at(.vars = vars(ATR:SHI) , .funs = function(x){x <- ifelse(x == '-', 0, x) %>% as.integer()})
}
title = dplyr::filter(to_do, Regions == region)$Short_Name
}
vars <- dplyr::select(var_s, -ISO3, -Country, -Regions, -Livehood_z, -Short_Name) %>%
tidyr::gather(key = 'var',value = 'value') %>%
dplyr::filter(value > 0) %>% dplyr::pull(var)
if(sum(vars == 'NWLD') == 1){vars <- c(vars, "NWLD50", "NWLD90")}
if(sum(vars == 'THI') == 1){
vars <- c(vars[vars != 'THI'], glue::glue('THI_{0:3}'))
vars <- c(vars, 'THI_23')}
if(sum(vars == 'HSI') == 1){
vars <- c(vars[vars != 'HSI'], glue::glue('HSI_{0:3}'))
vars <- c(vars, 'HSI_23')}
if(sum(vars == 'SLGP') == 1){
vars <- c(vars[vars != 'SLGP'], 'SLGP_CV')}
vars <- vars[!(vars %in% c("ISO3","Country","Regions","Livehood_z", "SPI","TAI" ))]
# -------------------------------------------------------------------------- #
allh <- all[all$time == 'Historical',]
allf <- all[all$time == 'Future',]
tsth <- allh %>% as_tibble() %>% ungroup() %>%
dplyr::select(season,period, year, ATR:THI_3) %>%
dplyr::distinct()
tstf <- allf %>% as_tibble() %>% ungroup() %>%
dplyr::select(season,period, year, ATR:THI_3) %>%
dplyr::distinct()
# =-------------------------
all_data <- dplyr::bind_rows(tsth,tstf) %>%
as.tibble() %>%
dplyr::mutate(THI_23 = THI_2 + THI_3, HSI_23 = HSI_2 + HSI_3) %>%
dplyr::group_by(season,period) %>%
dplyr::summarise_all(~mean(. , na.rm = TRUE)) %>%
dplyr::select(-year) %>%
dplyr::mutate_at(.vars = vars[vars %in% c('NDD', 'NT_X', 'NDWS', 'NWLD', 'NWLD50', 'NWLD90','SHI')],
.funs = ~round(. , 0)) %>%
ungroup()
all_data_sd <- dplyr::bind_rows(tsth,tstf) %>%
as.tibble() %>%
dplyr::mutate(THI_23 = THI_2 + THI_3, HSI_23 = HSI_2 + HSI_3) %>%
dplyr::group_by(season,period) %>%
dplyr::summarise_all(~sd(. , na.rm = TRUE)) %>%
dplyr::select(-year) %>%
ungroup()
all_data <- dplyr::full_join(tidyr::pivot_longer(all_data, cols= ATR:HSI_23 ,names_to='Index',values_to='Value') ,
tidyr::pivot_longer(all_data_sd, cols= ATR:HSI_23 ,names_to='Index',values_to='sd'))
p <- all_data %>%
ggplot2::ggplot(aes(x = factor(season, levels = 1:12), y = Value, fill = period)) +
ggplot2::geom_bar(stat = "identity", position = position_dodge()) +
ggplot2::geom_pointrange(aes(ymin=Value-sd, ymax=Value+sd), width=.2, position=position_dodge(.9)) +
ggplot2::facet_wrap(~Index, scales = 'free') +
ggplot2::scale_fill_brewer(palette = "Paired") +
ggplot2::theme_bw() +
ggplot2::xlab('') + ggplot2::ylab('') +
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(paste0(path_save, 'climatology_all.jpeg'), device = 'jpeg', width = 24, height = 12, units = 'in', dpi = 350)
# =-------------------------
complete <- all_data %>% dplyr::filter( Index %in% vars )
erase <- dplyr::distinct(all_data, Index)$Index[dplyr::distinct(all_data, Index)$Index %in% c(glue::glue('HSI_{2:3}'), glue::glue('THI_{2:3}') )]
complete <- complete %>%
dplyr::filter(!(Index %in% c("HSI_2","HSI_3","THI_2","THI_3")) ) # Arreglar la siguiente linea.
filt_graphs <- complete %>%
ggplot2::ggplot(aes(x = factor(season, levels = 1:12), y = Value, fill = period)) +
ggplot2::geom_bar(stat = "identity", position = position_dodge()) +
ggplot2::geom_pointrange(aes(ymin=Value-sd, ymax=Value+sd), width=.2, position=position_dodge(.9)) +
ggplot2::facet_wrap(~Index, scales = 'free') +
ggplot2::scale_fill_brewer(palette = "Paired") +
ggplot2::theme_bw() +
ggplot2::xlab('') + ggplot2::ylab('') +
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(paste0(path_save, 'climatology_subset_I.jpeg'), device = 'jpeg', width = 24, height = 12, units = 'in', dpi = 350)
# =-------------------------------
# By regions
shp <- terra::vect(paste0(root , "/1.Data/shps/", tolower(country), "/",tolower(iso),"_regions/",tolower(iso),"_regions.shp"))
shp$Short_Name <- to_do$Short_Name
ref <- terra::rast(paste0(root,"/1.Data/chirps-v2.0.2020.01.01.tif")) %>%
terra::crop(., terra::ext(shp))
shr <- terra::rasterize(x = shp, y = ref, field= 'region') %>% setNames(c('value'))
allhr <- cbind(terra::extract(x = shr, y = allh[,c('x','y')]),allh)
allfr <- cbind(terra::extract(x = shr, y = allf[,c('x','y')]),allf)
allhr2 <- allhr %>% tidyr::drop_na()
allfr2 <- allfr %>% tidyr::drop_na()
tsth <- allhr2 %>% dplyr::select(value,season,period,ATR:THI_3) %>%
dplyr::distinct() %>%
as.tibble() %>% dplyr::mutate(THI_23 = THI_2 + THI_3, HSI_23 = HSI_2 + HSI_3)
tsth$value <- factor(tsth$value)
levels(tsth$value) <- sort(unique(shp$region))
tstf <- allfr2 %>% dplyr::select(value,season,period,ATR:THI_3) %>%
dplyr::distinct() %>% as.tibble() %>%
dplyr::mutate(THI_23 = THI_2 + THI_3, HSI_23 = HSI_2 + HSI_3)
tstf$value <- factor(tstf$value)
levels(tstf$value) <- sort(unique(shp$region))
# Guardados...
region_mean <- dplyr::bind_rows(tsth,tstf) %>%
dplyr::group_by(value,season,period) %>%
dplyr::summarise_all(~mean(. , na.rm = TRUE)) %>%
dplyr::mutate_at(.vars = vars[vars %in% c('NDD', 'NT_X', 'NDWS', 'NWLD', 'NWLD50', 'NWLD90','SHI')],
.funs = ~round(. , 0)) %>%
dplyr::ungroup()
region_sd <- dplyr::bind_rows(tsth,tstf) %>%
dplyr::group_by(value,season,period) %>%
dplyr::summarise_all(~sd(.)) %>%
dplyr::ungroup()
region_data <- dplyr::full_join(tidyr::pivot_longer(region_mean, cols= ATR:HSI_23 ,names_to='Index',values_to='Value') ,
tidyr::pivot_longer(region_sd, cols= ATR:HSI_23 ,names_to='Index',values_to='sd'))
# =---
region_data %>%
dplyr::group_split(Index) %>%
purrr::map(function(tbl){
mande <- to_do$Short_Name
# Modificar de acuerdo a lo que se encuentre... sino pues...
mande <- str_replace(mande, '-', '\n')
mande <- str_replace(mande, ':', '\n')
mande <- str_replace(mande, '/', '\n')
# mande <- str_replace(mande, '[[:punct:]]', '\n')
names(mande) <- to_do$Regions
mande_label <- labeller(value = mande)
gg <- tbl %>%
ggplot2::ggplot(aes(x = factor(season, levels = 1:12), y = Value, fill = period)) +
ggplot2::geom_bar(stat = "identity", position = position_dodge()) +
ggplot2::geom_pointrange(aes(ymin=Value-sd, ymax=Value+sd), width=.2, position=position_dodge(.9)) +
ggplot2::facet_grid(Index~value, labeller = mande_label ,scales = 'free') +
ggplot2::scale_fill_brewer(palette = "Paired") +
ggplot2::theme_bw() +
ggplot2::xlab('') +
ggplot2::ylab('') +
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),
strip.text.y = element_text(size = 17),
plot.caption = element_text(size = 15, hjust = 0),
legend.position = "bottom")
ggplot2::ggsave(filename = paste0(path_save, '/clim_Ind_',unique(tbl$Index),'.jpeg'), plot = gg, device = 'jpeg', width = 20, height = 7, units = 'in', dpi = 350)
return(print(gg))
})
}