-
Notifications
You must be signed in to change notification settings - Fork 105
/
04_map_figures.R
224 lines (180 loc) · 10.9 KB
/
04_map_figures.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
# 04_map_figures.R
#
# Copyright (C) 2021 Santander Meteorology Group (http://meteo.unican.es)
#
# This work is licensed under a Creative Commons Attribution 4.0 International
# License (CC BY 4.0 - http://creativecommons.org/licenses/by/4.0)
#' @title Generate map figures from ensemble NcMLs
#' @description Generate map figures from the NcMLs created
#' using 03_ensemble_building.R, for Atlas Product Reproducibility.
#' @author M. Iturbide
# This script builds on the climate4R framework
# https://github.com/SantanderMetGroup/climate4R
# Climate4R package for data loading
library(loadeR)
# Climate4R package for data visualization
# <https://doi.org/10.1016/j.envsoft.2017.09.008>
library(visualizeR)
# climate 4R package for geoprocessing
library(geoprocessoR)
# Other utilities for spatial data handling and geoprocessing:
library(sp)
library(rgdal)
# Dev utilities used for source_url
library(devtools)
# Load function for latitudinal chunking
source_url("https://github.com/SantanderMetGroup/climate4R/blob/devel/R/climate4R.chunk.R?raw=TRUE")
# Load functions for computing uncertainty (signal, signal.ens and agreement)
source("../datasets-interactive-atlas/hatching-functions/hatching-functions.R")
# USER PARAMETER SETTINGS ------------------------------------------------------
# Select number of chunks
# Note: chunking sequentially splits the task into manageable data chunks to avoid memory problems
# Chunking operates by spliting the data into a predefined number latitudinal slices (n=2 in this example).
# Further details: https://github.com/SantanderMetGroup/climate4R/tree/master/R
n.chunks <- 100
# Index, scenario, season and reference and future period(s) of interest, e.g.:
AtlasIndex <- "meanpr" # index (monthly mean daily precipitation)
project <- "CMIP6"
scenario <- "ssp585" # scenario
season <- c(12, 1, 2) # (entire year: season = 1:12; boreal winter (DJF): season = c(12, 1, 2); boreal summer (JJA): season = 6:8, and so on...)
signal.period <- 1850:1900
base.period <- 1986:2005
future.period <- 2041:2060
# Path of the shapefile of World coastlines and referece regions, e.g.
## downloable from https://github.com/IPCC-WG1/Atlas/tree/devel/notebooks/auxiliary-material
coast.dir <- "auxiliary-material/WORLD_coastline.shp"
## downloable from https://github.com/IPCC-WG1/Atlas/tree/devel/reference-regions
regions.dir <- "../reference-regions/IPCC-WGI-reference-regions-v4.geojson"
# Source directory where the NcMLs (those created by script3_ensemble_building.R) are located
source.dir <- ""
# Output directory to save the .rda object of the computed deltas and to export the pdf of the figure
out.dir <- ""
# Color key graphical parameter:
# n = min value
# m = max value
# s = cut value frequency
# ct = Brewer color code (see 'RColorBrewer::display.brewer.all()')
if (AtlasIndex != "meanpr") { # case of precipitation
m <- 8
n <- 0
s <- 0.5
ct <- "Reds"
revc <- FALSE
} else {
m <- 50
n <- -50
s <- 5
ct <- "BrBG"
revc <- FALSE
}
## LOAD DATA --------------------------------------------------------------
dataset.hist <- list.files(source.dir, pattern = paste0("historical_", AtlasIndex, ".ncml"), full.names = TRUE)
dataset.ssp <- list.files(source.dir, pattern = paste0(scenario,"_", AtlasIndex, ".ncml"), full.names = TRUE)
# Retain common members among scenarios
hist.members <- dataInventory(dataset.hist)[[AtlasIndex]][["Dimensions"]][["member"]][["Values"]]
ssp.members <- dataInventory(dataset.ssp)[[AtlasIndex]][["Dimensions"]][["member"]][["Values"]]
hist.m.ind <- which(hist.members %in% ssp.members)
ssp.m.ind <- which(ssp.members %in% hist.members)
membernames <- ssp.members[ssp.m.ind]
## Data loading and aggregation
# For convenience, the auxiliary wrapper 'aux.fun' is next defined, performing the following actions in one single step:
# 1. Performs data subsetting along the member dimension to extract the common members defined in the previous step
# 2. Computes the climatology, this is the temporal mean for the whole target period.
aux.fun <- function(grid, members) {
climatology(
subsetGrid(grid, members = members), clim.fun = list(FUN = mean, na.rm = TRUE)
)
}
## A bit more complex auxiliary wrapper is aux.fun.signal, used for the historical data used for computing the signal.
# 1. Performs data subsetting along the member dimension to extract the common members defined in the previous step
# 2. Agregate the data annualy.
aux.fun.signal <- function(grid, members) {
aggregateGrid(
subsetGrid(grid, members = members), aggr.y = list(FUN = mean, na.rm = TRUE)
)
}
## Historical experiment data are loaded sequentially according to the chunking definition and applying the wrapper functions defined above.
hist.s <- climate4R.chunk(n.chunks = n.chunks,
C4R.FUN.args = list(FUN = "aux.fun.signal",
grid = list(dataset = dataset.hist, var = AtlasIndex),
members = hist.m.ind),
loadGridData.args = list(years = signal.period, season = season))
hist <- climate4R.chunk(n.chunks = n.chunks,
C4R.FUN.args = list(FUN = "aux.fun",
grid = list(dataset = dataset.hist, var = AtlasIndex),
members = hist.m.ind),
loadGridData.args = list(years = base.period, season = season))
# redim is a helper function to ensure array structure consistency among the different models
hist.s <- redim(hist.s, drop = TRUE); hist.s <- redim(hist.s)
hist <- redim(hist, drop = TRUE); hist <- redim(hist)
## RCP/SSP experiment future data are loaded sequentially according to the chunking definition and applying the wrapper function defined above.
ssp <- climate4R.chunk(n.chunks = n.chunks,
C4R.FUN.args = list(FUN = "aux.fun",
grid = list(dataset = dataset.ssp, var = AtlasIndex),
members = ssp.m.ind),
loadGridData.args = list(years = future.period, season = season))
# redim is a helper function to ensure array structure consistency among the different models
ssp <- redim(ssp, drop = TRUE); ssp <- redim(ssp)
## CALCULATE DELTAS ------------------------------------------------------------
# Deltas are the arithmetic difference between future and historical time slices
# NOTE: Relative deltas are calculated in the case of precipitation (i.e.: future / historical, in %)
delta <- gridArithmetics(climatology(ssp), climatology(hist), operator = "-")
rel.delta <- gridArithmetics(aggregateGrid(delta, aggr.mem = list(FUN = mean, na.rm = T)),
aggregateGrid(climatology(hist), aggr.mem = list(FUN = mean, na.rm = T)),
100,
operator = c("/", "*"))
# The default projection of the data is WGS84, but the Robinson projection is usually preferred, the nex lines
# reprojects the grids to Robinson
attr(delta$xyCoords, "projection") <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
attr(rel.delta$xyCoords, "projection") <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
delta.rob <- warpGrid(delta, new.CRS = "+proj=robin +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs")
rel.delta.rob <- warpGrid(rel.delta, new.CRS = "+proj=robin +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs")
# The intermediate object can be optionally stored as a R data object:
# save(delta, file = paste0(out.dir, "delta_", AtlasIndex, "_", scenario, "_", paste(season, collapse = "-"),".rda"))
# load(paste0(out.dir, "delta_", AtlasIndex, "_", scenario, "_", paste(season, collapse = "-"),".rda"), verbose = TRUE)
# CALCULATE UNCERTAINTY AND PRODUCE HATCHING ---------------------------------------------------
# Next uncertainty is calculated (check "../datasets-interactive-atlas/hatching-functions/hatching-functions.R"
# for further details on the uncertainty calculation). The outputs are binary C4R grids (0 = uncertain, 1 = certain)
# run e.g. spatialPlot(uncer2) to check the uncertainty areas.
# (1) signal
signal.grid <- signal(hist.s, delta)
uncer1 <- aggregateGrid(signal.grid, aggr.mem = list(FUN = signal.ens, th = 66))
# (2) agreement
uncer2 <- aggregateGrid(delta, aggr.mem = list(FUN = agreement, th = 80))
# The uncertainty is illustrated in the final figure with hatched areas. To create the hatches
# first the grids are reprojected to Robinson
attr(uncer1$xyCoords, "projection") <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
attr(uncer2$xyCoords, "projection") <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
uncer1.rob <- warpGrid(uncer1, new.CRS = "+proj=robin +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs")
uncer2.rob <- warpGrid(uncer2, new.CRS = "+proj=robin +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs")
# Hatches are created. The outputs are a list containing a SpatialLines object (form package sp)
# together with other graphical parameters. This list is later passed to the plotting function.
uncer1.hatch <- map.hatching(clim = climatology(uncer1.rob), threshold = "0.5", angle = "45",
condition = "LT", density = 4, lwd = 0.6,
upscaling.aggr.fun = list(FUN = mean))
uncer2.hatch <- map.hatching(clim = climatology(uncer2.rob), threshold = "0.5", angle = "-45",
condition = "LT", density = 4, lwd = 0.6,
upscaling.aggr.fun = list(FUN = mean))
# LOAD AND PROJECT ADDITIONAL MAP LINES ---------------------------------------
coast <- readOGR(coast.dir)
coast.rob <- spTransform(coast, CRSobj = CRS("+proj=robin +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"))
regions <- readOGR(regions.dir)
regions.rob <- spTransform(regions, CRSobj = CRS("+proj=robin +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"))
# PLOT MAP --------------------------------------------------------------------
# Calculate the ensemble mean (or use the relative delta)
out <- rel.delta.rob
#out <- aggregateGrid(delta.rob, aggr.mem = list(FUN = mean))
pl <- spatialPlot(out,
color.theme = ct,
rev.colors = revc,
at = seq(n, m, s),
set.max = m, set.min = n,
main = list(paste0(AtlasIndex, " mean delta change"), cex = 0.8),
xlab = list(paste0("Period: ", paste(range(future.period), collapse = "-"), ", Season: ", paste(month.abb[season], collapse = "-")), cex = 0.8),
sp.layout = list(uncer1.hatch, uncer2.hatch, list(coast.rob, col = "purple4", first = FALSE), list(regions.rob, col = "black", first = FALSE)),
par.settings = list(axis.line = list(col = 'transparent')))
pl
# Export the figure as PDF file
pdf(paste0(out.dir, "/Delta_", AtlasIndex, "_", scenario, "_", paste(range(future.period), collapse = "-"), "_", paste(month.abb[season], collapse = "-"),".pdf"), width = 10, height = 10)
pl
dev.off()