-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathRSCRIPTBATCH_OVERSEE_indivual_diversity_uncertainties_mappe.R
401 lines (367 loc) · 28.8 KB
/
RSCRIPTBATCH_OVERSEE_indivual_diversity_uncertainties_mappe.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
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
# --------------------------------------------------------------------------------------------------------------------------------
library("tidyverse")
library("RColorBrewer")
library("reshape2")
library("viridis")
library("scales")
library("maps")
library("geosphere")
library("rgdal")
library("raster")
library("betapart")
world2 <- map_data("world2")
# --------------------------------------------------------------------------------------------------------------------------------
### Set the working directories, vectors etc.
WD <- getwd()
rcp <- "rcp85"
# Vectors
SDMs <- c('GAM','GLM','RF','ANN')
# Vector of months
months <- c("jan","feb","mar","apr","may","jun","jul","aug","sep","oct","nov","dec")
# Vector of pools
pools <- c("p1","p2","p3","p4")
# Vector of earth system models
ESMs <- c("CESM-BEC","CNRM-PISCES","GFDL-TOPAZ","IPSL-PISCES","MRI-NEMURO")
# --------------------------------------------------------------------------------------------------------------------------------
### 1°) Load individual projections in perc change in SR for phyto and zoo
setwd("/net/kryo/work/fabioben/OVERSEE/data/tables_composition_ensemble_rcp85/Individual_projections")
# format: "table_ann_changes_",esm,"_",sdm,"_",p,".Rdata"
files <- dir()[grep("table_ann_changes_",dir())]#; files
# Separate the ones for richness (no "beta.div") from those with beta.div
files2 <- files[grep("beta.div", files)]#; files2
files1 <- files[!(files %in% files2)]; files1
require("parallel")
res <- mclapply(files1, function(f) {
# Useless message
message(paste("Loading results from ",f, sep = ""))
t <- get(load(f))
# Get attributes from filename
terms <- do.call(cbind, strsplit(as.character(f),"_"))
t$ESM <- terms[4,1]
t$SDM <- terms[5,1]
p <- terms[6,1]
t$pool <- str_replace(as.character(p),".Rdata","")
# Return
return(t)
}, mc.cores = 25
) # eo mclapply
# Rbind
table <- dplyr::bind_rows(res)
rm(res); gc()
### Computing ensembles (all, SDM/ESM/pool)
ens <- data.frame(table %>%
group_by(id) %>%
summarize(x = unique(x), y = unique(y),
rich_tot = mean(perc_rich_tot, na.rm = T), sd_rich_tot = sd(perc_rich_tot, na.rm = T),
rich_phyto = mean(perc_rich_phyto, na.rm = T), sd_rich_phyto = sd(perc_rich_phyto, na.rm = T),
rich_zoo = mean(perc_rich_zoo, na.rm = T), sd_rich_zoo = sd(perc_rich_zoo, na.rm = T)
)
) # ddf
# summary(ens)
### 2°) First, map stdev (Figs. A and B of the overall panel)
mapA <- ggplot() + geom_raster(aes(x = x, y = y, fill = sd_rich_phyto), data = ens[!is.na(ens$sd_rich_phyto),]) +
geom_contour(colour = "grey25", binwidth = 33, size = 0.25, aes(x = x, y = y, z = sd_rich_phyto),
data = ens[!is.na(ens$sd_rich_phyto),]) +
scale_fill_distiller(name = "Standard\ndeviation", palette = "YlOrRd", direction = 1) +
geom_polygon(aes(x = long, y = lat, group = group), data = world2, fill = "grey70", colour = "black", size = 0.3) +
coord_quickmap() + scale_x_continuous(name = "", breaks = c(0,60,120,180,240,300,360),
labels = c("0°W","60°W","120°W","180°W","-120°W","-60°W","0°W"), expand = c(0,0)) +
scale_y_continuous(name = "", breaks = c(-90,-60,-30,0,30,60,90),
labels = c("-90°N","-60°N","-30°N","0°N","30°N","60°N","90°N"), expand = c(0,0)) +
theme(panel.background = element_rect(fill = "white"),legend.key = element_rect(fill = "grey50"),
panel.grid.major = element_line(colour = "grey70",linetype = "dashed") )
#
mapB <- ggplot() + geom_raster(aes(x = x, y = y, fill = sd_rich_zoo), data = ens[!is.na(ens$sd_rich_zoo),]) +
geom_contour(colour = "grey25", binwidth = 33, size = 0.25, aes(x = x, y = y, z = sd_rich_zoo),
data = ens[!is.na(ens$sd_rich_zoo),]) +
scale_fill_distiller(name = "Standard\ndeviation", palette = "YlOrRd", direction = 1) +
geom_polygon(aes(x = long, y = lat, group = group), data = world2, fill = "grey70", colour = "black", size = 0.3) +
coord_quickmap() + scale_x_continuous(name = "", breaks = c(0,60,120,180,240,300,360),
labels = c("0°W","60°W","120°W","180°W","-120°W","-60°W","0°W"), expand = c(0,0)) +
scale_y_continuous(name = "", breaks = c(-90,-60,-30,0,30,60,90),
labels = c("-90°N","-60°N","-30°N","0°N","30°N","60°N","90°N"), expand = c(0,0)) +
theme(panel.background = element_rect(fill = "white"),legend.key = element_rect(fill = "grey50"),
panel.grid.major = element_line(colour = "grey70",linetype = "dashed") )
require("ggpubr")
ggarrange(mapA, mapB, ncol = 2, nrow = 1, labels = c("A1","A2"))
### 3°) The 8 maps (LETTERS[3:10]) for the each SDM type
# Per SDM (for facet)
ens.SDM <- data.frame(table %>%
group_by(id,SDM) %>%
summarize(x = unique(x), y = unique(y),
rich_tot = mean(perc_rich_tot, na.rm = T),
rich_phyto = mean(perc_rich_phyto, na.rm = T),
rich_zoo = mean(perc_rich_zoo, na.rm = T)
)
) # ddf
summary(ens.SDM)
min <- floor(min(ens.SDM[,c("rich_phyto","rich_zoo")], na.rm = T) )
min
mapC <- ggplot() + geom_raster(aes(x = x, y = y, fill = rich_phyto),
data = ens.SDM[which(ens.SDM$rich_phyto < 55 & ens.SDM$SDM == "GLM"),]) +
geom_raster(aes(x = x, y = y), data = ens.SDM[which(ens.SDM$rich_phyto >= 55 & ens.SDM$SDM == "GLM"),], fill = "#b2182b") +
geom_contour(colour = "grey60", binwidth = 25, size = 0.25, aes(x = x, y = y, z = rich_phyto),
data = ens.SDM[which(ens.SDM$rich_phyto < 55 & ens.SDM$SDM == "GLM"),] ) +
scale_fill_gradient2(name = "Richness\ndifference", low = "#3288bd", high = "#d53e4f", mid = "white", limits = c(min,55)) +
geom_polygon(aes(x = long, y = lat, group = group), data = world2, fill = "grey70", colour = "black", size = 0.3) +
coord_quickmap() + scale_x_continuous(name = "", breaks = c(0,60,120,180,240,300,360),
labels = c("0°W","60°W","120°W","180°W","-120°W","-60°W","0°W"), expand = c(0,0)) +
scale_y_continuous(name = "", breaks = c(-90,-60,-30,0,30,60,90),
labels = c("-90°N","-60°N","-30°N","0°N","30°N","60°N","90°N"), expand = c(0,0)) +
theme(panel.background = element_rect(fill = "white"),legend.key = element_rect(fill = "grey50"),
panel.grid.major = element_line(colour = "grey70",linetype = "dashed") )
mapD <- ggplot() + geom_raster(aes(x = x, y = y, fill = rich_zoo),
data = ens.SDM[which(ens.SDM$rich_zoo < 55 & ens.SDM$SDM == "GLM"),]) +
geom_raster(aes(x = x, y = y), data = ens.SDM[which(ens.SDM$rich_zoo >= 55 & ens.SDM$SDM == "GLM"),], fill = "#b2182b") +
geom_contour(colour = "grey60", binwidth = 25, size = 0.25, aes(x = x, y = y, z = rich_zoo),
data = ens.SDM[which(ens.SDM$rich_zoo < 55 & ens.SDM$SDM == "GLM"),] ) +
scale_fill_gradient2(name = "Richness\ndifference", low = "#3288bd", high = "#d53e4f", mid = "white", limits = c(min,55)) +
geom_polygon(aes(x = long, y = lat, group = group), data = world2, fill = "grey70", colour = "black", size = 0.3) +
coord_quickmap() + scale_x_continuous(name = "", breaks = c(0,60,120,180,240,300,360),
labels = c("0°W","60°W","120°W","180°W","-120°W","-60°W","0°W"), expand = c(0,0)) +
scale_y_continuous(name = "", breaks = c(-90,-60,-30,0,30,60,90),
labels = c("-90°N","-60°N","-30°N","0°N","30°N","60°N","90°N"), expand = c(0,0)) +
theme(panel.background = element_rect(fill = "white"),legend.key = element_rect(fill = "grey50"),
panel.grid.major = element_line(colour = "grey70",linetype = "dashed") )
#
mapE <- ggplot() + geom_raster(aes(x = x, y = y, fill = rich_phyto),
data = ens.SDM[which(ens.SDM$rich_phyto < 55 & ens.SDM$SDM == "GAM"),]) +
geom_raster(aes(x = x, y = y), data = ens.SDM[which(ens.SDM$rich_phyto >= 55 & ens.SDM$SDM == "GAM"),], fill = "#b2182b") +
geom_contour(colour = "grey60", binwidth = 25, size = 0.25, aes(x = x, y = y, z = rich_phyto),
data = ens.SDM[which(ens.SDM$rich_phyto < 55 & ens.SDM$SDM == "GAM"),] ) +
scale_fill_gradient2(name = "Richness\ndifference", low = "#3288bd", high = "#d53e4f", mid = "white", limits = c(min,55)) +
geom_polygon(aes(x = long, y = lat, group = group), data = world2, fill = "grey70", colour = "black", size = 0.3) +
coord_quickmap() + scale_x_continuous(name = "", breaks = c(0,60,120,180,240,300,360),
labels = c("0°W","60°W","120°W","180°W","-120°W","-60°W","0°W"), expand = c(0,0)) +
scale_y_continuous(name = "", breaks = c(-90,-60,-30,0,30,60,90),
labels = c("-90°N","-60°N","-30°N","0°N","30°N","60°N","90°N"), expand = c(0,0)) +
theme(panel.background = element_rect(fill = "white"),legend.key = element_rect(fill = "grey50"),
panel.grid.major = element_line(colour = "grey70",linetype = "dashed") )
mapF <- ggplot() + geom_raster(aes(x = x, y = y, fill = rich_zoo),
data = ens.SDM[which(ens.SDM$rich_zoo < 55 & ens.SDM$SDM == "GAM"),]) +
geom_raster(aes(x = x, y = y), data = ens.SDM[which(ens.SDM$rich_zoo >= 55 & ens.SDM$SDM == "GAM"),], fill = "#b2182b") +
geom_contour(colour = "grey60", binwidth = 25, size = 0.25, aes(x = x, y = y, z = rich_zoo),
data = ens.SDM[which(ens.SDM$rich_zoo < 55 & ens.SDM$SDM == "GAM"),] ) +
scale_fill_gradient2(name = "Richness\ndifference", low = "#3288bd", high = "#d53e4f", mid = "white", limits = c(min,55)) +
geom_polygon(aes(x = long, y = lat, group = group), data = world2, fill = "grey70", colour = "black", size = 0.3) +
coord_quickmap() + scale_x_continuous(name = "", breaks = c(0,60,120,180,240,300,360),
labels = c("0°W","60°W","120°W","180°W","-120°W","-60°W","0°W"), expand = c(0,0)) +
scale_y_continuous(name = "", breaks = c(-90,-60,-30,0,30,60,90),
labels = c("-90°N","-60°N","-30°N","0°N","30°N","60°N","90°N"), expand = c(0,0)) +
theme(panel.background = element_rect(fill = "white"),legend.key = element_rect(fill = "grey50"),
panel.grid.major = element_line(colour = "grey70",linetype = "dashed") )
#
mapG <- ggplot() + geom_raster(aes(x = x, y = y, fill = rich_phyto),
data = ens.SDM[which(ens.SDM$rich_phyto < 55 & ens.SDM$SDM == "ANN"),]) +
geom_raster(aes(x = x, y = y), data = ens.SDM[which(ens.SDM$rich_phyto >= 55 & ens.SDM$SDM == "ANN"),], fill = "#b2182b") +
geom_contour(colour = "grey60", binwidth = 25, size = 0.25, aes(x = x, y = y, z = rich_phyto),
data = ens.SDM[which(ens.SDM$rich_phyto < 55 & ens.SDM$SDM == "ANN"),] ) +
scale_fill_gradient2(name = "Richness\ndifference", low = "#3288bd", high = "#d53e4f", mid = "white", limits = c(min,55)) +
geom_polygon(aes(x = long, y = lat, group = group), data = world2, fill = "grey70", colour = "black", size = 0.3) +
coord_quickmap() + scale_x_continuous(name = "", breaks = c(0,60,120,180,240,300,360),
labels = c("0°W","60°W","120°W","180°W","-120°W","-60°W","0°W"), expand = c(0,0)) +
scale_y_continuous(name = "", breaks = c(-90,-60,-30,0,30,60,90),
labels = c("-90°N","-60°N","-30°N","0°N","30°N","60°N","90°N"), expand = c(0,0)) +
theme(panel.background = element_rect(fill = "white"),legend.key = element_rect(fill = "grey50"),
panel.grid.major = element_line(colour = "grey70",linetype = "dashed") )
mapH <- ggplot() + geom_raster(aes(x = x, y = y, fill = rich_zoo),
data = ens.SDM[which(ens.SDM$rich_zoo < 55 & ens.SDM$SDM == "ANN"),]) +
geom_raster(aes(x = x, y = y), data = ens.SDM[which(ens.SDM$rich_zoo >= 55 & ens.SDM$SDM == "ANN"),], fill = "#b2182b") +
geom_contour(colour = "grey60", binwidth = 25, size = 0.25, aes(x = x, y = y, z = rich_zoo),
data = ens.SDM[which(ens.SDM$rich_zoo < 55 & ens.SDM$SDM == "ANN"),] ) +
scale_fill_gradient2(name = "Richness\ndifference", low = "#3288bd", high = "#d53e4f", mid = "white", limits = c(min,55)) +
geom_polygon(aes(x = long, y = lat, group = group), data = world2, fill = "grey70", colour = "black", size = 0.3) +
coord_quickmap() + scale_x_continuous(name = "", breaks = c(0,60,120,180,240,300,360),
labels = c("0°W","60°W","120°W","180°W","-120°W","-60°W","0°W"), expand = c(0,0)) +
scale_y_continuous(name = "", breaks = c(-90,-60,-30,0,30,60,90),
labels = c("-90°N","-60°N","-30°N","0°N","30°N","60°N","90°N"), expand = c(0,0)) +
theme(panel.background = element_rect(fill = "white"),legend.key = element_rect(fill = "grey50"),
panel.grid.major = element_line(colour = "grey70",linetype = "dashed") )
#
mapI <- ggplot() + geom_raster(aes(x = x, y = y, fill = rich_phyto),
data = ens.SDM[which(ens.SDM$rich_phyto < 55 & ens.SDM$SDM == "RF"),]) +
geom_raster(aes(x = x, y = y), data = ens.SDM[which(ens.SDM$rich_phyto >= 55 & ens.SDM$SDM == "RF"),], fill = "#b2182b") +
geom_contour(colour = "grey60", binwidth = 25, size = 0.25, aes(x = x, y = y, z = rich_phyto),
data = ens.SDM[which(ens.SDM$rich_phyto < 55 & ens.SDM$SDM == "RF"),] ) +
scale_fill_gradient2(name = "Richness\ndifference", low = "#3288bd", high = "#d53e4f", mid = "white", limits = c(min,55)) +
geom_polygon(aes(x = long, y = lat, group = group), data = world2, fill = "grey70", colour = "black", size = 0.3) +
coord_quickmap() + scale_x_continuous(name = "", breaks = c(0,60,120,180,240,300,360),
labels = c("0°W","60°W","120°W","180°W","-120°W","-60°W","0°W"), expand = c(0,0)) +
scale_y_continuous(name = "", breaks = c(-90,-60,-30,0,30,60,90),
labels = c("-90°N","-60°N","-30°N","0°N","30°N","60°N","90°N"), expand = c(0,0)) +
theme(panel.background = element_rect(fill = "white"),legend.key = element_rect(fill = "grey50"),
panel.grid.major = element_line(colour = "grey70",linetype = "dashed") )
mapJ <- ggplot() + geom_raster(aes(x = x, y = y, fill = rich_zoo),
data = ens.SDM[which(ens.SDM$rich_zoo < 55 & ens.SDM$SDM == "RF"),]) +
geom_raster(aes(x = x, y = y), data = ens.SDM[which(ens.SDM$rich_zoo >= 55 & ens.SDM$SDM == "RF"),], fill = "#b2182b") +
geom_contour(colour = "grey60", binwidth = 25, size = 0.25, aes(x = x, y = y, z = rich_zoo),
data = ens.SDM[which(ens.SDM$rich_zoo < 55 & ens.SDM$SDM == "RF"),] ) +
scale_fill_gradient2(name = "Richness\ndifference", low = "#3288bd", high = "#d53e4f", mid = "white", limits = c(min,55)) +
geom_polygon(aes(x = long, y = lat, group = group), data = world2, fill = "grey70", colour = "black", size = 0.3) +
coord_quickmap() + scale_x_continuous(name = "", breaks = c(0,60,120,180,240,300,360),
labels = c("0°W","60°W","120°W","180°W","-120°W","-60°W","0°W"), expand = c(0,0)) +
scale_y_continuous(name = "", breaks = c(-90,-60,-30,0,30,60,90),
labels = c("-90°N","-60°N","-30°N","0°N","30°N","60°N","90°N"), expand = c(0,0)) +
theme(panel.background = element_rect(fill = "white"),legend.key = element_rect(fill = "grey50"),
panel.grid.major = element_line(colour = "grey70",linetype = "dashed") )
### Test panel
labels <- paste(expand.grid(c(1,2), LETTERS[1:10])$Var2, expand.grid(c(1,2), LETTERS[1:10])$Var1, sep = "")
panel1 <- ggarrange(mapA, mapB, mapC, mapD, mapE, mapF, mapG, mapH, mapI, mapJ, ncol = 2, nrow = 5, labels = labels)
# Cool
### Same, with ESM
ens.ESM <- data.frame(table %>%
group_by(id,ESM) %>%
summarize(x = unique(x), y = unique(y),
rich_tot = mean(perc_rich_tot, na.rm = T),
rich_phyto = mean(perc_rich_phyto, na.rm = T),
rich_zoo = mean(perc_rich_zoo, na.rm = T)
)
) # ddf
head(ens.ESM)
min <- floor(min(ens.ESM[,c("rich_phyto","rich_zoo")], na.rm = T) )
min
unique(ens.ESM$ESM)
mapC <- ggplot() + geom_raster(aes(x = x, y = y, fill = rich_phyto),
data = ens.ESM[which(ens.ESM$rich_phyto < 55 & ens.ESM$ESM == "CESM-BEC"),]) +
geom_raster(aes(x = x, y = y), data = ens.ESM[which(ens.ESM$rich_phyto >= 55 & ens.ESM$ESM == "CESM-BEC"),], fill = "#b2182b") +
geom_contour(colour = "grey60", binwidth = 25, size = 0.25, aes(x = x, y = y, z = rich_phyto),
data = ens.ESM[which(ens.ESM$rich_phyto < 55 & ens.ESM$ESM == "CESM-BEC"),] ) +
scale_fill_gradient2(name = "Richness\ndifference", low = "#3288bd", high = "#d53e4f", mid = "white", limits = c(min,55)) +
geom_polygon(aes(x = long, y = lat, group = group), data = world2, fill = "grey70", colour = "black", size = 0.3) +
coord_quickmap() + scale_x_continuous(name = "", breaks = c(0,60,120,180,240,300,360),
labels = c("0°W","60°W","120°W","180°W","-120°W","-60°W","0°W"), expand = c(0,0)) +
scale_y_continuous(name = "", breaks = c(-90,-60,-30,0,30,60,90),
labels = c("-90°N","-60°N","-30°N","0°N","30°N","60°N","90°N"), expand = c(0,0)) +
theme(panel.background = element_rect(fill = "white"),legend.key = element_rect(fill = "grey50"),
panel.grid.major = element_line(colour = "grey70",linetype = "dashed") )
mapD <- ggplot() + geom_raster(aes(x = x, y = y, fill = rich_zoo),
data = ens.ESM[which(ens.ESM$rich_zoo < 55 & ens.ESM$ESM == "CESM-BEC"),]) +
geom_raster(aes(x = x, y = y), data = ens.ESM[which(ens.ESM$rich_zoo >= 55 & ens.ESM$ESM == "CESM-BEC"),], fill = "#b2182b") +
geom_contour(colour = "grey60", binwidth = 25, size = 0.25, aes(x = x, y = y, z = rich_zoo),
data = ens.ESM[which(ens.ESM$rich_zoo < 55 & ens.ESM$ESM == "CESM-BEC"),] ) +
scale_fill_gradient2(name = "Richness\ndifference", low = "#3288bd", high = "#d53e4f", mid = "white", limits = c(min,55)) +
geom_polygon(aes(x = long, y = lat, group = group), data = world2, fill = "grey70", colour = "black", size = 0.3) +
coord_quickmap() + scale_x_continuous(name = "", breaks = c(0,60,120,180,240,300,360),
labels = c("0°W","60°W","120°W","180°W","-120°W","-60°W","0°W"), expand = c(0,0)) +
scale_y_continuous(name = "", breaks = c(-90,-60,-30,0,30,60,90),
labels = c("-90°N","-60°N","-30°N","0°N","30°N","60°N","90°N"), expand = c(0,0)) +
theme(panel.background = element_rect(fill = "white"),legend.key = element_rect(fill = "grey50"),
panel.grid.major = element_line(colour = "grey70",linetype = "dashed") )
#
mapE <- ggplot() + geom_raster(aes(x = x, y = y, fill = rich_phyto),
data = ens.ESM[which(ens.ESM$rich_phyto < 55 & ens.ESM$ESM == "CNRM-PISCES"),]) +
geom_raster(aes(x = x, y = y), data = ens.ESM[which(ens.ESM$rich_phyto >= 55 & ens.ESM$ESM == "CNRM-PISCES"),], fill = "#b2182b") +
geom_contour(colour = "grey60", binwidth = 25, size = 0.25, aes(x = x, y = y, z = rich_phyto),
data = ens.ESM[which(ens.ESM$rich_phyto < 55 & ens.ESM$ESM == "CNRM-PISCES"),] ) +
scale_fill_gradient2(name = "Richness\ndifference", low = "#3288bd", high = "#d53e4f", mid = "white", limits = c(min,55)) +
geom_polygon(aes(x = long, y = lat, group = group), data = world2, fill = "grey70", colour = "black", size = 0.3) +
coord_quickmap() + scale_x_continuous(name = "", breaks = c(0,60,120,180,240,300,360),
labels = c("0°W","60°W","120°W","180°W","-120°W","-60°W","0°W"), expand = c(0,0)) +
scale_y_continuous(name = "", breaks = c(-90,-60,-30,0,30,60,90),
labels = c("-90°N","-60°N","-30°N","0°N","30°N","60°N","90°N"), expand = c(0,0)) +
theme(panel.background = element_rect(fill = "white"),legend.key = element_rect(fill = "grey50"),
panel.grid.major = element_line(colour = "grey70",linetype = "dashed") )
mapF <- ggplot() + geom_raster(aes(x = x, y = y, fill = rich_zoo),
data = ens.ESM[which(ens.ESM$rich_zoo < 55 & ens.ESM$ESM == "CNRM-PISCES"),]) +
geom_raster(aes(x = x, y = y), data = ens.ESM[which(ens.ESM$rich_zoo >= 55 & ens.ESM$ESM == "CNRM-PISCES"),], fill = "#b2182b") +
geom_contour(colour = "grey60", binwidth = 25, size = 0.25, aes(x = x, y = y, z = rich_zoo),
data = ens.ESM[which(ens.ESM$rich_zoo < 55 & ens.ESM$ESM == "CNRM-PISCES"),] ) +
scale_fill_gradient2(name = "Richness\ndifference", low = "#3288bd", high = "#d53e4f", mid = "white", limits = c(min,55)) +
geom_polygon(aes(x = long, y = lat, group = group), data = world2, fill = "grey70", colour = "black", size = 0.3) +
coord_quickmap() + scale_x_continuous(name = "", breaks = c(0,60,120,180,240,300,360),
labels = c("0°W","60°W","120°W","180°W","-120°W","-60°W","0°W"), expand = c(0,0)) +
scale_y_continuous(name = "", breaks = c(-90,-60,-30,0,30,60,90),
labels = c("-90°N","-60°N","-30°N","0°N","30°N","60°N","90°N"), expand = c(0,0)) +
theme(panel.background = element_rect(fill = "white"),legend.key = element_rect(fill = "grey50"),
panel.grid.major = element_line(colour = "grey70",linetype = "dashed") )
#
mapG <- ggplot() + geom_raster(aes(x = x, y = y, fill = rich_phyto),
data = ens.ESM[which(ens.ESM$rich_phyto < 55 & ens.ESM$ESM == "GFDL-TOPAZ"),]) +
geom_raster(aes(x = x, y = y), data = ens.ESM[which(ens.ESM$rich_phyto >= 55 & ens.ESM$ESM == "GFDL-TOPAZ"),], fill = "#b2182b") +
geom_contour(colour = "grey60", binwidth = 25, size = 0.25, aes(x = x, y = y, z = rich_phyto),
data = ens.ESM[which(ens.ESM$rich_phyto < 55 & ens.ESM$ESM == "GFDL-TOPAZ"),] ) +
scale_fill_gradient2(name = "Richness\ndifference", low = "#3288bd", high = "#d53e4f", mid = "white", limits = c(min,55)) +
geom_polygon(aes(x = long, y = lat, group = group), data = world2, fill = "grey70", colour = "black", size = 0.3) +
coord_quickmap() + scale_x_continuous(name = "", breaks = c(0,60,120,180,240,300,360),
labels = c("0°W","60°W","120°W","180°W","-120°W","-60°W","0°W"), expand = c(0,0)) +
scale_y_continuous(name = "", breaks = c(-90,-60,-30,0,30,60,90),
labels = c("-90°N","-60°N","-30°N","0°N","30°N","60°N","90°N"), expand = c(0,0)) +
theme(panel.background = element_rect(fill = "white"),legend.key = element_rect(fill = "grey50"),
panel.grid.major = element_line(colour = "grey70",linetype = "dashed") )
mapH <- ggplot() + geom_raster(aes(x = x, y = y, fill = rich_zoo),
data = ens.ESM[which(ens.ESM$rich_zoo < 55 & ens.ESM$ESM == "GFDL-TOPAZ"),]) +
geom_raster(aes(x = x, y = y), data = ens.ESM[which(ens.ESM$rich_zoo >= 55 & ens.ESM$ESM == "GFDL-TOPAZ"),], fill = "#b2182b") +
geom_contour(colour = "grey60", binwidth = 25, size = 0.25, aes(x = x, y = y, z = rich_zoo),
data = ens.ESM[which(ens.ESM$rich_zoo < 55 & ens.ESM$ESM == "GFDL-TOPAZ"),] ) +
scale_fill_gradient2(name = "Richness\ndifference", low = "#3288bd", high = "#d53e4f", mid = "white", limits = c(min,55)) +
geom_polygon(aes(x = long, y = lat, group = group), data = world2, fill = "grey70", colour = "black", size = 0.3) +
coord_quickmap() + scale_x_continuous(name = "", breaks = c(0,60,120,180,240,300,360),
labels = c("0°W","60°W","120°W","180°W","-120°W","-60°W","0°W"), expand = c(0,0)) +
scale_y_continuous(name = "", breaks = c(-90,-60,-30,0,30,60,90),
labels = c("-90°N","-60°N","-30°N","0°N","30°N","60°N","90°N"), expand = c(0,0)) +
theme(panel.background = element_rect(fill = "white"),legend.key = element_rect(fill = "grey50"),
panel.grid.major = element_line(colour = "grey70",linetype = "dashed") )
#
mapI <- ggplot() + geom_raster(aes(x = x, y = y, fill = rich_phyto),
data = ens.ESM[which(ens.ESM$rich_phyto < 55 & ens.ESM$ESM == "IPSL-PISCES"),]) +
geom_raster(aes(x = x, y = y), data = ens.ESM[which(ens.ESM$rich_phyto >= 55 & ens.ESM$ESM == "IPSL-PISCES"),], fill = "#b2182b") +
geom_contour(colour = "grey60", binwidth = 25, size = 0.25, aes(x = x, y = y, z = rich_phyto),
data = ens.ESM[which(ens.ESM$rich_phyto < 55 & ens.ESM$ESM == "IPSL-PISCES"),] ) +
scale_fill_gradient2(name = "Richness\ndifference", low = "#3288bd", high = "#d53e4f", mid = "white", limits = c(min,55)) +
geom_polygon(aes(x = long, y = lat, group = group), data = world2, fill = "grey70", colour = "black", size = 0.3) +
coord_quickmap() + scale_x_continuous(name = "", breaks = c(0,60,120,180,240,300,360),
labels = c("0°W","60°W","120°W","180°W","-120°W","-60°W","0°W"), expand = c(0,0)) +
scale_y_continuous(name = "", breaks = c(-90,-60,-30,0,30,60,90),
labels = c("-90°N","-60°N","-30°N","0°N","30°N","60°N","90°N"), expand = c(0,0)) +
theme(panel.background = element_rect(fill = "white"),legend.key = element_rect(fill = "grey50"),
panel.grid.major = element_line(colour = "grey70",linetype = "dashed") )
mapJ <- ggplot() + geom_raster(aes(x = x, y = y, fill = rich_zoo),
data = ens.ESM[which(ens.ESM$rich_zoo < 55 & ens.ESM$ESM == "IPSL-PISCES"),]) +
geom_raster(aes(x = x, y = y), data = ens.ESM[which(ens.ESM$rich_zoo >= 55 & ens.ESM$ESM == "IPSL-PISCES"),], fill = "#b2182b") +
geom_contour(colour = "grey60", binwidth = 25, size = 0.25, aes(x = x, y = y, z = rich_zoo),
data = ens.ESM[which(ens.ESM$rich_zoo < 55 & ens.ESM$ESM == "IPSL-PISCES"),] ) +
scale_fill_gradient2(name = "Richness\ndifference", low = "#3288bd", high = "#d53e4f", mid = "white", limits = c(min,55)) +
geom_polygon(aes(x = long, y = lat, group = group), data = world2, fill = "grey70", colour = "black", size = 0.3) +
coord_quickmap() + scale_x_continuous(name = "", breaks = c(0,60,120,180,240,300,360),
labels = c("0°W","60°W","120°W","180°W","-120°W","-60°W","0°W"), expand = c(0,0)) +
scale_y_continuous(name = "", breaks = c(-90,-60,-30,0,30,60,90),
labels = c("-90°N","-60°N","-30°N","0°N","30°N","60°N","90°N"), expand = c(0,0)) +
theme(panel.background = element_rect(fill = "white"),legend.key = element_rect(fill = "grey50"),
panel.grid.major = element_line(colour = "grey70",linetype = "dashed") )
#
mapK <- ggplot() + geom_raster(aes(x = x, y = y, fill = rich_phyto),
data = ens.ESM[which(ens.ESM$rich_phyto < 55 & ens.ESM$ESM == "MRI-NEMURO"),]) +
geom_raster(aes(x = x, y = y), data = ens.ESM[which(ens.ESM$rich_phyto >= 55 & ens.ESM$ESM == "MRI-NEMURO"),], fill = "#b2182b") +
geom_contour(colour = "grey60", binwidth = 25, size = 0.25, aes(x = x, y = y, z = rich_phyto),
data = ens.ESM[which(ens.ESM$rich_phyto < 55 & ens.ESM$ESM == "MRI-NEMURO"),] ) +
scale_fill_gradient2(name = "Richness\ndifference", low = "#3288bd", high = "#d53e4f", mid = "white", limits = c(min,55)) +
geom_polygon(aes(x = long, y = lat, group = group), data = world2, fill = "grey70", colour = "black", size = 0.3) +
coord_quickmap() + scale_x_continuous(name = "", breaks = c(0,60,120,180,240,300,360),
labels = c("0°W","60°W","120°W","180°W","-120°W","-60°W","0°W"), expand = c(0,0)) +
scale_y_continuous(name = "", breaks = c(-90,-60,-30,0,30,60,90),
labels = c("-90°N","-60°N","-30°N","0°N","30°N","60°N","90°N"), expand = c(0,0)) +
theme(panel.background = element_rect(fill = "white"),legend.key = element_rect(fill = "grey50"),
panel.grid.major = element_line(colour = "grey70",linetype = "dashed") )
mapL <- ggplot() + geom_raster(aes(x = x, y = y, fill = rich_zoo),
data = ens.ESM[which(ens.ESM$rich_zoo < 55 & ens.ESM$ESM == "MRI-NEMURO"),]) +
geom_raster(aes(x = x, y = y), data = ens.ESM[which(ens.ESM$rich_zoo >= 55 & ens.ESM$ESM == "MRI-NEMURO"),], fill = "#b2182b") +
geom_contour(colour = "grey60", binwidth = 25, size = 0.25, aes(x = x, y = y, z = rich_zoo),
data = ens.ESM[which(ens.ESM$rich_zoo < 55 & ens.ESM$ESM == "MRI-NEMURO"),] ) +
scale_fill_gradient2(name = "Richness\ndifference", low = "#3288bd", high = "#d53e4f", mid = "white", limits = c(min,55)) +
geom_polygon(aes(x = long, y = lat, group = group), data = world2, fill = "grey70", colour = "black", size = 0.3) +
coord_quickmap() + scale_x_continuous(name = "", breaks = c(0,60,120,180,240,300,360),
labels = c("0°W","60°W","120°W","180°W","-120°W","-60°W","0°W"), expand = c(0,0)) +
scale_y_continuous(name = "", breaks = c(-90,-60,-30,0,30,60,90),
labels = c("-90°N","-60°N","-30°N","0°N","30°N","60°N","90°N"), expand = c(0,0)) +
theme(panel.background = element_rect(fill = "white"),legend.key = element_rect(fill = "grey50"),
panel.grid.major = element_line(colour = "grey70",linetype = "dashed") )
### Arrange in second panel
labels <- paste(expand.grid(c(1,2), LETTERS[6:10])$Var2, expand.grid(c(1,2), LETTERS[6:10])$Var1, sep = "")
panel2 <- ggarrange(mapC, mapD, mapE, mapF, mapG, mapH, mapI, mapJ, mapK, mapL, ncol = 2, nrow = 5, labels = labels)
### Save both panels
setwd("/net/kryo/work/fabioben/OVERSEE/data/")
ggsave(plot = panel1, filename = "panel_maps_uncertainty_A-E.jpg", dpi = 300, width = 10, height = 10)
ggsave(plot = panel2, filename = "panel_maps_uncertainty_F-J.jpg", dpi = 300, width = 10, height = 10)
# --------------------------------------------------------------------------------------------------------------------------------
# --------------------------------------------------------------------------------------------------------------------------------
# --------------------------------------------------------------------------------------------------------------------------------