Newer
Older
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
---
output:
pdf_document: default
html_document: default
---
Insurance Risk Indices
========================
Insurance Risk Indices are an ensemble of indices relevant for insurance industry. These indices are based on Expert Team on Climate Change Detection Indices (ETCCDI). There are currently 5 available indices to compute for extreme: heat (tx90p), cold (tn10p), wind (wx), drought (ccd) and flooding (rx5day). The individual indices can be combined into a single index with or without weighting for each component. This combined index is roughly analogous to the Actuaris Climate Risk Index.
### 1- Load dependencies
This example requires the following system libraries:
- libssl-dev
- libnecdf-dev
- cdo
The **ClimProjDiag R package** should be loaded by running the following lines in R once it's integrated into CRAN mirror.
```r
library(ClimProjDiag)
```
All the other R packages involved can be installed directly from CRAN and loaded as follows:
```r
library(s2dverification)
library(abind)
library(multiApply)
library(ggplot2)
library(parallel)
```
### 2- Synthetic data
Daily maximum and minimum temperature, wind speed and precipitation are necessary to compute the different indices in both, the reference period (1971 - 2000) and the future projection (2006 - 2100). The defined region will be in the northern hemisphere between -40 - 20 ºE and 25 - 60 ºN.
Maximum temperature is generated considering the annual cycle:
```r
lat <- seq(25, 60, 5)
lon <- seq(-35, 20 ,5)
tmax_historical <- NULL
grid1 <- 293 - 10 * cos(2 * pi / 365 * (1 : 10958)) + rnorm(10958)
gridlon <- NULL
for (i in 1 : 12) {
gridlon <- cbind(gridlon,
grid1 + rnorm(10958, sd = 5) * cos(2 * pi / 365 * (1 : 10958)))
}
for (j in 1 : 8) {
gridnew <- apply(gridlon, 2, function(x) {x - rnorm(10958, mean = j * 0.5, sd = 3)})
tmax_historical <- abind(tmax_historical, gridnew, along = 3)
}
tmax_historical <- InsertDim(InsertDim(tmax_historical, posdim = 1, lendim = 1),
posdim = 1, lendim = 1)
names(dim(tmax_historical)) <- c("model", "var", "time", "lon", "lat")
time <- seq(ISOdate(1971, 1, 1), ISOdate(2000, 12, 31), "day")
metadata <- list(time = list(standard_name = 'time', long_name = 'time',
calendar = 'proleptic_gregorian',
units = 'days since 1970-01-01 00:00:00', prec = 'double',
dim = list(list(name = 'time', unlim = FALSE))))
attr(time, "variables") <- metadata
attr(tmax_historical, 'Variables')$dat1$time <- time
```
A similar procedure is considered to build the synthetic data for the future projections. However, a little trend is added.
```r
tmax_projection <- NULL
grid1 <- 298 - 10 * cos(2 * pi / 365 * (1 : 34698)) + rnorm(34698) +
(1 : 34698) * rnorm(1, mean = 4) / 34698
gridlon <- NULL
for (i in 1 : 12) {
gridlon <- cbind(gridlon, grid1 + rnorm(34698, sd = 5) *
cos(2 * pi / 365 * (1 : 34698)))
}
for (j in 1 : 8) {
gridnew <- apply(gridlon, 2, function(x) {x -
rnorm(34698, mean = j * 0.5, sd = 3)})
tmax_projection <- abind(tmax_projection, gridnew, along = 3)
}
tmax_projection <- InsertDim(InsertDim(tmax_projection, posdim = 1, lendim = 1), posdim = 1, lendim = 1)
names(dim(tmax_projection)) <- c("model", "var", "time", "lon", "lat")
time <- seq(ISOdate(2006, 1, 1), ISOdate(2100, 12, 31), "day")
metadata <- list(time = list(standard_name = 'time', long_name = 'time',
calendar = 'proleptic_gregorian',
units = 'days since 1970-01-01 00:00:00', prec = 'double',
dim = list(list(name = 'time', unlim = FALSE))))
attr(time, "variables") <- metadata
attr(tmax_projection, 'Variables')$dat1$time <- time
```
To build synthetic precipitation data, a lognormal distribution is considered:
```r
ppt_historical <- rlnorm(10958 * 12 * 8)
dim(ppt_historical) <- c(model = 1, var = 1, time = 10958, lon = 12, lat = 8)
time <- seq(ISOdate(1971, 1, 1), ISOdate(2000, 12, 31), "day")
metadata <- list(time = list(standard_name = 'time', long_name = 'time',
calendar = 'proleptic_gregorian',
units = 'days since 1970-01-01 00:00:00', prec = 'double',
dim = list(list(name = 'time', unlim = FALSE))))
attr(time, "variables") <- metadata
attr(ppt_historical, 'Variables')$dat1$time <- time
ppt_projection <- rlnorm(34698 * 12 * 8)
dim(ppt_projection) <- c(model = 1, var = 1, time = 34698, lon = 12, lat = 8)
time <- seq(ISOdate(2006, 1, 1), ISOdate(2100, 12, 31), "day")
metadata <- list(time = list(standard_name = 'time', long_name = 'time',
calendar = 'proleptic_gregorian',
units = 'days since 1970-01-01 00:00:00', prec = 'double',
dim = list(list(name = 'time', unlim = FALSE))))
attr(time, "variables") <- metadata
attr(ppt_projection, 'Variables')$dat1$time <- time
```
### 3- Computing the Extreme Heat Index
The Extreme Heat Index (*t90p*) is defined as the percentage of days when the maximum temperature exceeds the 90th percentile.
In order to evaluate the future projections, it is necessary to compute the index during a reference historical period.
The next steps should be followed:
To remove seasonality effects, the anomaly is computed for each day and gridpoint by applying the `DailyAno` function. The name of the first dimensions is defined as 'time' dimension.
```r
anomaly_data <- apply(tmax_historical, c(1,2,4,5), DailyAno, dates = attributes(tmax_historical)$Variables$dat1$time)
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
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
names(dim(anomaly_data))[1] <- "time"
```
This data can be detrended by applying `Trend` function from **s2dverification package**. In order to remove the trend from the `tmax_historical`, the correction is calculated by subtracting the `detrended_data` to the `anomaly_data`.
```r
detrended_data <- Trend(anomaly_data,
posTR = which(names(dim(anomaly_data)) == "time"))
diff <- anomaly_data - detrended_data$detrended
diff <- aperm(diff, c(2,3,1,4,5))
detrended_data <- tmax_historical - diff
```
For each gridpoint and day of the year (from the 1st of January to the 31st of December), the maximum temperature on the position at the 90 percent of the series will be calculated as the threshold.
```r
quantile <- 0.9
thresholds <- Threshold(detrended_data, qtiles = quantile,
ncores = detectCores() -1)
```
```r
> dim(thresholds)
jdays model var lon lat
366 1 1 12 8
```
By indicating the metric and introducing the threshold, `Climdex()` function will return the extreme heat index during the reference period.
```r
base_index <- Climdex(data = detrended_data, metric = 't90p',
threshold = thresholds, ncores = detectCores() - 1)
```
The output of ´Climdex´ function will be a ´list()´ object. Index values are saved in the ´base_index$result´ label.
```r
> str(base_index)
List of 2
$ result: num [1:30, 1, 1, 1:12, 1:8] 11.23 8.74 10.41 11.78 10.14 ...
$ years : num [1:30] 1971 1972 1973 1974 1975 ...
> dim(base_index$result)
year model var lon lat
30 1 1 12 8
```
Now, the standard deviation is computed in order to standardize the index. Notice that, by definition, the mean of the percentage of the number of days exceeding the 90th percentile is 10. Only standard deviation is computed.
```r
base_sd <- Apply(list(base_index$result), target_dims = list(c(1)),
AtomicFun = "sd")$output1
```
The index can be computed by considering the threshold obtain for the reference period.
```r
projection_index <- Climdex(data = tmax_projection, metric = 't90p',
threshold = thresholds, ncores = detectCores() - 1)
```
Its normalized with mean 10 and the standard deviation of the reference period.
```r
base_mean <- 10
base_sd <- InsertDim(base_sd, 1, dim(projection_index$result)[1])
HeatExtremeIndex <- (projection_index$result - base_mean) / base_sd
```
A spatial representation of the mean index values is obtained and save in PNG format in the working directory with the name: "SpatialExtremeHeatIndex.png". The matrix `masc` is build and shown as dots in the plot indicating wich pixels are considered land.
```r
masc <- rep(0, 8 * 12)
masc[c(5 : 12, 18 : 24, 31 : 34, 43, 44, 47, 56 : 60, 67 : 72, 79,
82 : 84, 93 : 96)] <- 1
dim(masc) <- c(12, 8)
PlotEquiMap(Mean1Dim(HeatExtremeIndex,
which(names(dim(HeatExtremeIndex)) == "year")),
lon = lon, lat = lat, filled.continents = FALSE,
toptitle = "Extreme Heat Index", dots = masc,
fileout = "SpatialExtremeHeatIndex.png")
```
![Spatial distribution of Extreme Heat Index](SpatialExtremeHeatIndex.png)
The inland average of the Extreme Heat Index can be computed to plot its time evolution using `WeigthedMean` function. `Smoothing()` returns the smoothed time series for a 3 year moving window which can be modified using `runmeanlen` parameter.
```r
temporal <- WeightedMean(HeatExtremeIndex, lon = lon, lat = lat, mask = drop(masc))
temporal_3ysmooth <- Smoothing(temporal, runmeanlen = 3, numdimt = 1)
```
The next code should be run to plot and save the original average and the 3 year smoothed data.
```r
png("Temporal_Inland_ExtremeHeatIndex.png", width = 8, height = 5, units = 'in',
res = 100, type = "cairo")
plot(2006 : 2100, temporal, type = "l", lty = 5, lwd = 2, bty = 'n',
xlab = "Time (years)", ylab = "Extreme Heat Index",
main = "Inland average Extreme Heat Index")
lines(2006 : 2100, temporal_3ysmooth, col = "darkgreen", lwd = 2)
legend('bottomright', c('Anual', '3 years smooth'), col = c(1, 'darkgreen'),
lty = c(5, 1), lwd = 2, bty = 'n')
dev.off()
```
![Temporal Inland Extreme Heat Index](Temporal_Inland_ExtremeHeatIndex.png)
### 4- Extreme Drought Index
The Extreme Drought Index (*cdd*), which measures the maximum length of a dry spell, is defined as the maximum number of consecutive days with the daily precipitation amount lower than 1 mm.
To compute the Extreme Drought Index during the reference period and its standar deviation and mean:
*Note: Precipitation data is not detrended. Furthermore, this index doesn't require to compute a threshold as `Climdex` function integrates the threshold of precipitation amount lower than 1 mm internally. However, this case requires the calculation of the mean.*
```r
base_index <- Climdex(data = ppt_historical, metric = 'cdd',
ncores = detectCores() - 1)
base_mean <- Apply(list(base_index$result), target_dims = list(c(1)),
AtomicFun = "mean")$output1
base_sd <- Apply(list(base_index$result), target_dims = list(c(1)),
AtomicFun = "sd")$output1
```
The object `base_index` contains the output of the `Climdex` function as two list with the next dimensions:
```r
> str(base_index)
List of 2
$ result: num [1:30, 1, 1, 1:12, 1:8] 6 11 8 8 8 12 9 10 6 8 ...
$ years : num [1:30] 1971 1972 1973 1974 1975 ...
```
The Extreme Drought Index is computed and standardized:
```r
projection_index <- Climdex(data = ppt_projection, metric = 'cdd',
ncores = detectCores() - 1)
base_mean <- InsertDim(base_mean, 1, dim(projection_index$result)[1])
base_sd <- InsertDim(base_sd, 1, dim(projection_index$result)[1])
DroughtExtremeIndex <- (projection_index$result - base_mean) / base_sd
```
Spatial representation of the Extreme Drought Index:
```r
PlotEquiMap(Mean1Dim(DroughtExtremeIndex,
which(names(dim(DroughtExtremeIndex)) == "year")),
lon = lon, lat = lat, filled.continents = FALSE,
toptitle = "Drought Index", brks = seq(-1, 1, 0.01),
fileout = "SpatialDroughtIndex.png")
```
![Spatial distribution of Extreme Drought Index](SpatialDroughtIndex.png)
Evolution of inland average of the Extreme Drought Index:
```r
temporal <- WeightedMean(DroughtExtremeIndex, lon = lon, lat = lat,
mask = drop(masc))
temporal_5ysmooth <- Smoothing(temporal, runmeanlen = 5, numdimt = 1)
png("Temporal_Inland_ExtremeDroughtIndex.png", width = 8, height = 5, units= 'in',
res = 100, type = "cairo")
plot(2006: 2100, temporal, type = "l", lty = 5, lwd = 2, bty = 'n',
xlab = "Time (years)", ylab = "Extreme Drought Index",
main = "Inland average Extreme Drought Index")
lines(2006 : 2100, temporal_5ysmooth, col = "darkgreen",lwd = 2)
legend('bottomleft', c('Anual', '3 years smooth'), col= c(1, 'darkgreen'),
lty = c(5, 1), lwd = 2, bty = 'n')
dev.off()
```
![Temporal evolution of Inland Extreme Drought Index](Temporal_Inland_ExtremeDroughtIndex.png)
### 5- Extreme Flooding Index
The Extreme Flooding Index (*rx5day*) is defined as the maximum precipitation amount in 5 consecutive days.
The Extreme Flooding Index during the reference period and its standard deviation and mean can be calculated by executing:
```r
base_index <- Climdex(data = ppt_historical, metric = 'rx5day',
ncores = detectCores() - 1)
base_mean <- Apply(list(base_index$result), target_dims = list(c(1)),
AtomicFun = "mean")$output1
base_sd <- Apply(list(base_index$result), target_dims = list(c(1)),
AtomicFun = "sd")$output1
```
The Extreme Flooding Index is computed and standardized:
```r
projection_index <- Climdex(data = ppt_projection, metric = 'rx5day',
ncores = detectCores() - 1)
base_mean <- InsertDim(base_mean, 1, dim(projection_index$result)[1])
base_sd <- InsertDim(base_sd, 1, dim(projection_index$result)[1])
FloodingExtremeIndex <- (projection_index$result - base_mean) / base_sd
```
Spatial representation of the Extreme Flooding Index:
```r
PlotEquiMap(Mean1Dim(FloodingExtremeIndex,
which(names(dim(FloodingExtremeIndex)) == "year")), lon = lon,
lat = lat, filled.continents = FALSE,
toptitle = "Extreme Flooding Index",
brks = seq(-1, 1, 0.1), fileout = "SpatialFloodingIndex.png")
```
![Spatial distribution of Extreme Flooding Index](SpatialFloodingIndex.png)
- Evolution of inland average of the Extreme Flooding Index:
```r
temporal <- WeightedMean(FloodingExtremeIndex, lon = lon, lat = lat,
mask = drop(masc))
temporal_3ysmooth <- Smoothing(temporal, runmeanlen = 3, numdimt = 1)
png("Temporal_Inland_ExtremeFloodingIndex.png", width = 8, height = 5,
units= 'in', res = 100, type = "cairo")
plot(2006 : 2100, temporal, type = "l", lty = 5, lwd = 2, bty = 'n',
xlab = "Time (years)", ylab = "Extreme Flooding Index",
main = "Inland average Extreme Flooding Index")
lines(2006 : 2100, temporal_3ysmooth, col = "darkgreen",lwd = 2)
legend('bottomleft', c('Anual', '3 years smooth'), col= c(1, 'darkgreen'),
lty = c(5, 1), lwd = 2, bty = 'n')
dev.off()
```
![Temporal evolution of Inland Extreme Flooding Index](Temporal_Inland_ExtremeFloodingIndex.png)
### 6- Combining Indices
The individual indices can be combined into a single index with or without weighting for each component. This combined index is roughly analogous to the Actuaris Climate Risk Index (see http://actuariesclimateindex.org/home/). Extreme Indices should be saved in the same `list` object.
```r
indices <- list()
indices[[1]] <- HeatExtremeIndex
indices[[2]] <- DroughtExtremeIndex
indices[[3]] <- FloodingExtremeIndex
```
If the `weights` parameter is defined as `NULL`, all indices will be equally weighted if `operation` parameter is set as `mean` (by default). To define other `weights` a vector of length equal to the number of considered indices (5 in this example) and with total sum equal to 1.
```r
aci <- CombineIndices(indices = indices, weights = NULL)
```
A spatial visulitzation can be performs by executing:
```r
PlotEquiMap(Mean1Dim(aci, which(names(dim(aci)) == "year")), lon = lon,
lat = lat, filled.continents = FALSE, toptitle = "Indices Combination",
fileout = "CombinedIndices.png")
```
![Spatial distribution of Combined Indices](CombinedIndices.png)
*Note: This vignette shows the computation of three indices, however, five different indices can be computed with `Climdex` function. To consider other combination settings run `?CombinedIndices`.*