-
Notifications
You must be signed in to change notification settings - Fork 4
/
plotBehaviour.R
330 lines (278 loc) · 11.1 KB
/
plotBehaviour.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
#' Check inputs for plotBehaviour function
#'
#' @noRd
#' @param sts StormsList object
#' @param rasterProduct Spatraster
#' @param xlim numeric vector
#' @param ylim numeric vector
#' @param labels logical
#' @param by numeric
#' @param pos numeric
#' @param colorPalette character vector
#' @param main character
#' @param legends character
#' @param dynamicPlot logical
#' @return NULL, just stops the function if an error is found
checkInputsPlotBehaviour <- function(sts, rasterProduct, xlim, ylim, labels, by, pos, colorPalette, main, legends, dynamicPlot) {
# Checking sts input
stopifnot("no data to plot" = !missing(sts))
# Checking rasterProduct
stopifnot("no data to plot" = !missing(rasterProduct))
stopifnot("Raster stack are not allowed. Please subset your desired layer" = terra::nlyr(rasterProduct) == 1)
# Checking xlim input
if (!is.null(xlim)) {
stopifnot("xlim must be numeric" = identical(class(xlim), "numeric"))
stopifnot("xlim must length 2" = length(xlim) == 2)
stopifnot("xlim must have valid longitude coordinates" = xlim >= 0 &
xlim <= 360)
}
# Checking ylim input
if (!is.null(ylim)) {
stopifnot("ylim must be numeric" = identical(class(ylim), "numeric"))
stopifnot("ylim must length 2" = length(ylim) == 2)
stopifnot("ylim must have valid latitude coordinates" = ylim >= -90 &
ylim <= 90)
}
# Checking labels input
stopifnot("labels must be logical" = identical(class(labels), "logical"))
# Checking by input
stopifnot("by must be numeric" = identical(class(by), "numeric"))
stopifnot("by must length 1" = length(by) == 1)
stopifnot("by must be as integer" = round(by) == by)
# Checking pos input
stopifnot("pos must be numeric" = identical(class(pos), "numeric"))
stopifnot("pos must length 1" = length(pos) == 1)
stopifnot("pos must be between either 1, 2, 3 or 4" = pos %in% c(1, 2, 3, 4))
# Checking colorPalette input
if (!is.null(colorPalette)) {
stopifnot("colorPalette must be character" = identical(class(colorPalette), "character"))
}
# Checking main input
if (!is.null(main)) {
stopifnot("main must be character" = identical(class(main), "character"))
stopifnot("main must be length 1" = length(main) == 1)
}
# Checking legend
stopifnot("legends must be character" = identical(class(legends), "character"))
stopifnot("legends must be length 1" = length(legends) == 1)
stopifnot(
"legends must be either topright, topleft, bottomleft, bottomright, or none" =
legends %in% c("topright", "topleft", "bottomleft", "bottomright", "none")
)
#Checking mode input
stopifnot("dynamicPlot must be logical" = identical(class(dynamicPlot), "logical"))
stopifnot("dynamicPlot must length 1" = length(dynamicPlot) == 1)
}
#' Plotting spatial wind behaviour
#'
#' The `plotBehaviour()` function allows plotting spatial statistics generated using
#' the `spatialBehaviour()` function and stored in `SpatRaster` objects.
#'
#' @param sts `StormsList` object.
#' @param rasterProduct layer name in a `SpatRaster` object. The names of the layers follow
#' the following terminology:
#' \itemize{
#' \item for "MSW" or "PDI", the name of the storm in capital letters and the name of the
#' statistic separated by underscores (e.g., "PAM_MSW", "PAM_PDI"),
#' \item for duration of exposure, the name of the storm in capital letters, "Exposure",
#' and the threshold value separated by underscores (e.g., "PAM_Exposure_18", "PAM_Exposure_33", ...).
#' \item for wind profiles, the name of the storm in capital letters, "Speed" or "Direction",
#' and the indices of the observation separated by underscores (e.g., "PAM_Speed_41", "PAM_Direction_41",...).
#' }
#' @param colorPalette character vector. The color palette used to plot the raster layer.
#' If `colorPalette=NULL` (default setting), the default color palette is used.
#' @param main character. Title of the plot. If `main=NULL` (default setting),
#' a default title is generated based on the name of the layer.
#' @param xlim numeric vector. The x limits of the plot.
#' @param ylim numeric vector. The y limits of the plot.
#' @param labels logical. Whether (TRUE) or not (FALSE, default setting) to add labels with the name
#' of the storm and the indices and ISO times of the observation.
#' @param by numeric. If `labels=TRUE`, defines the frequency at which labels are plotted.
#' Default value is set to `8` which corresponds to a 24h (or 48h) time interval between the labelled observations
#' when observations are made every 3 (or 6) hours.
#' @param pos numeric. If `labels=TRUE`, defines the position of the labels, `1` (above the observation),
#' `2` (on the left), `3` (below, default setting), and `4` (on the right).
#' @param legends character. Indicates where to plot the legend, `"topright"`(default setting), `"topleft"`,
#' `"bottomleft"`, `"bottomright"`, or `"none"` (legend not plotted).
#' @param dynamicPlot logical. Whether (FALSE, default setting) or (TRUE) to plot the
#' data dynamicaly using leaflet library
#' @returns A plot of the storm track data with the raster layer.
#'
#' @examples
#' \donttest{
#' # Creating a stormsDataset
#' sds <- defStormsDataset()
#'
#' # Getting storm track data for tropical cyclone Pam (2015)
#' pam <- defStormsList(sds = sds, loi = "Vanuatu", names = "PAM")
#'
#' # Plotting maximum sustained wind speed for Pam (2015) near Vanuatu
#' pam.msw <- spatialBehaviour(pam, verbose = 0)
#' plotBehaviour(pam, pam.msw)
#'
#' # dynamicPlot mode
#' plotBehaviour(pam, pam.msw, dynamicPlot = TRUE)
#'
#' # Plotting 2D wind speed profile for Pam (2015) near Vanuatu
#' pam.prof <- spatialBehaviour(pam, product = "Profiles", verbose = 0)
#' plotBehaviour(pam, pam.prof$PAM_Speed_37, labels = TRUE, pos = 4)
#'
#
#' }
#' @export
plotBehaviour <- function(sts,
rasterProduct,
colorPalette = NULL,
main = NULL,
xlim = NULL,
ylim = NULL,
labels = FALSE,
by = 8,
pos = 3,
legends = "topright",
dynamicPlot = FALSE) {
checkInputsPlotBehaviour(
sts, rasterProduct, xlim, ylim, labels, by, pos, colorPalette,
main, legends, dynamicPlot
)
name <- strsplit(names(rasterProduct), split = "_", fixed = TRUE)[[1]][1]
product <- strsplit(names(rasterProduct), split = "_", fixed = TRUE)[[1]][2]
if (!(name %in% getNames(sts))) {
stop("Imcompatibility between rasterProduct and sts (name not found in sts)")
}
# Handling spatial extent
xmin <- terra::ext(rasterProduct)$xmin
xmax <- terra::ext(rasterProduct)$xmax
ymin <- terra::ext(rasterProduct)$ymin
ymax <- terra::ext(rasterProduct)$ymax
if (!is.null(xlim)) {
xlim <- xlim[order(xlim)]
xmin <- xlim[1]
xmax <- xlim[2]
}
if (!is.null(ylim)) {
ylim <- ylim[order(ylim)]
ymin <- ylim[1]
ymax <- ylim[2]
}
# Choose color, range and legends depending on input raster
if (product == "MSW") {
col <- sts@scalePalette
range <- c(17, 95)
leg <- ifelse(dynamicPlot, "MSW (m.s <sup>-1</sup>)", expression(paste("MSW (m.s"^"-1", ")")))
} else if (product == "PDI") {
col <- pdiPalette
range <- c(0, max(terra::values(rasterProduct), na.rm = TRUE))
leg <- ifelse(dynamicPlot, "PDI (J.m<sup>2</sup>)",expression(paste("PDI (J.m"^"2", ")")))
} else if (product == "Exposure") {
col <- exposurePalette
range <- c(0, max(terra::values(rasterProduct), na.rm = TRUE))
leg <- ifelse(dynamicPlot, "Duration of exposure (h)",expression(paste("Duration of exposure (h)")))
} else if (product == "Speed") {
col <- sts@scalePalette
range <- c(17, 95)
leg <- ifelse(dynamicPlot, "Radial wind speed (m.s <sup>-1</sup>)",expression(paste("Radial wind speed (m.s"^"-1", ")")))
} else if (product == "Direction") {
col <- exposurePalette
range <- c(0, 360)
leg <- ifelse(dynamicPlot, "Wind direction (degree)",expression(paste("Wind direction (degree)")))
}
if (!is.null(colorPalette)) {
col <- colorPalette
}
if (!is.null(main)) {
leg <- main
}
if (!dynamicPlot) {
# Plotting track
plotStorms(
sts = sts, names = name, xlim = c(xmin, xmax), ylim = c(ymin, ymax),
legends = legends
)
# Adding title
graphics::title(leg)
terra::plot(rasterProduct,
col = col,
type = "continuous",
xlim = c(xmin - 1, xmax + 1), # we extend W & E by 1°. Needs to be in agreement with plotStorm
ylim = c(ymin - 1, ymax + 1), # we extend S & N by 1°. Needs to be in agreement with plotStorm
alpha = 0.7,
axes = FALSE,
range = range,
legend = TRUE,
add = TRUE
)
# Adding track again (to emphazise)
plotTrack(sts@data[[name]], sts@scale, sts@scalePalette)
# Adding labels
if (labels && product != "Profiles" && product != "WindDirection") {
plotLabels(sts@data[[name]], by, pos)
}
if (labels && (product == "Profiles" || product == "WindDirection")) {
ind <- as.numeric(strsplit(names(rasterProduct), split = "_", fixed = TRUE)[[1]][3])
if (round(ind) == ind) {
# It is a real observation
graphics::text(sts@data[[name]]@obs.all$lon[ind],
sts@data[[name]]@obs.all$lat[ind],
labels = paste0(
name, "\n", sts@data[[name]]@obs.all$iso.time[ind],
"\n(", ind, ")"
),
pos = pos,
cex = 0.6
)
} else {
# It is an interpolated observation
indf <- floor(ind)
indc <- ceiling(ind)
pos2 <- switch(pos,
"1" = 3,
"2" = 4,
"3" = 1,
"4" = 2
)
graphics::text(sts@data[[name]]@obs.all$lon[indf],
sts@data[[name]]@obs.all$lat[indf],
labels = paste0(
name, "\n", sts@data[[name]]@obs.all$iso.time[indf],
"\n(", indf, ")"
),
pos = pos,
cex = 0.6
)
graphics::text(sts@data[[name]]@obs.all$lon[indc],
sts@data[[name]]@obs.all$lat[indc],
labels = paste0(
name, "\n", sts@data[[name]]@obs.all$iso.time[indc],
"\n(", indc, ")"
),
pos = pos2,
cex = 0.6
)
}
}
}else{
# dynamicPlot plot
map <- plotStorms(sts = sts,
names = name,
xlim = c(xmin, xmax),
ylim = c(ymin, ymax),
legends = legends, dynamicPlot = TRUE)
pal <- leaflet::colorNumeric(col,
terra::values(rasterProduct),
na.color = "transparent")
map <- leaflet::addRasterImage(map,
rasterProduct,
colors = pal,
opacity = 0.8
)
#Adding legends
map <- leaflet::addLegend(map,
legends,
pal = pal,
values = terra::values(rasterProduct),
title = leg,
opacity = 0.8)
map
}
}