-
Notifications
You must be signed in to change notification settings - Fork 4
/
plotStorms.R
434 lines (359 loc) · 12.8 KB
/
plotStorms.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
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
#' Get the right color associated with a wind observation
#'
#' This function returns the Saffir Simpson Hurricane Scale color associated
#' with a maximum sustained wind speed
#'
#' @noRd
#' @param msw numeric. Maximum Sustained Wind (m/s)
#' @param scale list of values that defines the scale intensity of the storm, e.g. `sshs`
#' @param palette list of colors that defines the scale intensity of the storm
#' @return color associated with the observation
getColors <- function(msw, scale, palette) {
if (is.na(msw)) {
color <- NA
} else {
i <- findInterval(msw, scale)
color <- palette[i + 1]
}
return(color)
}
#' Plot track of a storm
#'
#' This function plots the track of a storm on a map that should be previously
#' plotted
#'
#' @noRd
#' @param st Storm object
#' @param scale numeric vector. Thresholds used for the scale
#' @param scale named character vector. colors used for the scale
#'
#' @return A plot with the track of the storm
plotTrack <- function(st, scale, scalePalette) {
cexL <- 1
cexP <- 0.6
lon <- st@obs.all$lon
lat <- st@obs.all$lat
msw <- st@obs.all$msw
colors <- unlist(lapply(msw, getColors, scale = scale, palette = scalePalette))
graphics::lines(
lon,
lat,
col = "black",
lty = 2,
lwd = 1,
cex = cexL
)
graphics::points(
lon,
lat,
col = colors,
pch = 19,
lwd = 1,
cex = cexP
)
}
#' Add labels on plot
#'
#' Add names labels, indices of observations and ISO Times for a given storm on
#' a map that should be previsouly plotted
#'
#' @noRd
#' @param st Storm object
#' @param by numeric. Increment of the sequence for the labels to plot
#' @param pos numeric. Must be between 1 and 4. Corresponds to the position of
#' labels according to the observation: 1 (up), 2 (left), 3 (down), 4 (right)
#'
#' @return A plot with added labels
plotLabels <- function(st, by, pos) {
cex <- 0.5
ind <- round(seq(1, getNbObs(st), by))
for (i in ind) {
lon <- st@obs.all$lon[i]
lat <- st@obs.all$lat[i]
graphics::text(lon,
lat,
labels = paste0(st@name, " (", i, ")\n", st@obs.all$iso.time[i]),
pos = pos,
cex = cex
)
}
}
#' Check inputs for plotStorms function
#'
#' @noRd
#' @param sts StormsList object
#' @param names character vector
#' @param category numeric vector
#' @param labels logical
#' @param by numeric
#' @param pos numeric
#' @param legends logical
#' @param loi logical
#' @param xlim numeric vector
#' @param ylim numeric vector
#' @param dynamicPlot logical
#' @return NULL, just stops the function if an error is found
checkInputsPlotStorms <- function(sts, names, category, labels, by,
pos, legends, loi, xlim, ylim, dynamicPlot) {
# Checking sts input
stopifnot("no data to plot" = !missing(sts))
# Checking names input
if (!is.null(names)) {
stopifnot("names must be characters" = identical(class(names), "character"))
stopifnot("Invalid storm names (storm not found)" = names %in% getNames(sts))
}
# Checking category input
if (!is.null(category)) {
stopifnot("category must be numeric(s)" = identical(class(category), "numeric"))
stopifnot("Invalid category input" = category %in% seq(0, length(sts@scale)))
}
# 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 logical inputs
stopifnot("loi must be logical" = identical(class(loi), "logical"))
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 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 storm track data
#'
#' This `plotStorms()` function allows plotting storm track data stored in a `StormsList` object.
#'
#' @param sts `StormsList` object
#' @param names character vector. Name(s) of the storm(s) in capital letters.
#' If `names = NULL` (default setting), all storms are plotted.
#' @param category numeric vector. Category of storms to be plotted among the level in
#' the windscale provided in `sts` input. If `category=NULL` (default setting),
#' all storms are plotted.
#' @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 loi logical. Whether (TRUE, default setting) or not (FALSE) to plot the
#' extent of the buffered location of interest on the map.
#' @param dynamicPlot logical. Whether (FALSE, default setting) or (TRUE) to plot the
#' data dynamicaly using leaflet library
#'
#' @return A plot of the storm track data.
#'
#' @examples
#' #' #Creating a stormsDataset
#' \donttest{
#' dev.off()
#' sds <- defStormsDataset()
#'
#' # Getting storm track data for tropical cyclone Pam (2015)
#' pam <- defStormsList(sds = sds, loi = "Vanuatu", names = "PAM")
#'
#' # Plotting Pam over Vanuatu with labels every 24h
#' plotStorms(sts =pam, labels = TRUE)
#'
#' # Plotting Pam over Vanuatu with labels every 6h on the right side of the observations
#' plotStorms(pam, labels = TRUE, by = 2, pos = 4)
#'
#' # dynamicPlot mode
# # plotStorms(pam, dynamicPlot=TRUE)
#'
#' }
#' @export
plotStorms <- function(sts,
names = NULL,
category = NULL,
xlim = NULL,
ylim = NULL,
labels = FALSE,
by = 8,
pos = 3,
legends = "topright",
loi = TRUE,
dynamicPlot = FALSE) {
checkInputsPlotStorms(sts, names, category, labels, by, pos, legends,loi, xlim, ylim, dynamicPlot)
# Handling spatial extent
ext <- terra::ext(
sf::st_bbox(sts@spatialLoiBuffer)$xmin,
sf::st_bbox(sts@spatialLoiBuffer)$xmax,
sf::st_bbox(sts@spatialLoiBuffer)$ymin,
sf::st_bbox(sts@spatialLoiBuffer)$ymax
)
if (!is.null(xlim)) {
xlim <- xlim[order(xlim)]
ext$xmin <- xlim[1]
ext$xmax <- xlim[2]
}
if (!is.null(ylim)) {
ylim <- ylim[order(ylim)]
ext$ymin <- ylim[1]
ext$ymax <- ylim[2]
}
# Handling categories and names
if (!is.null(category) && is.null(names)) {
if (length(category) == 2) {
category <- category[order(category)]
catInf <- category[1]
catSup <- category[2]
ind <- which(getScale(sts) >= catInf & getScale(sts) <= catSup)
} else {
# length category == 1
ind <- which(getScale(sts) == category)
}
stsAux <- unlist(sts@data)[ind]
} else {
if (!is.null(names)) {
if (!is.null(category)) {
warning("category input ignored\n")
}
ind <- which(getNames(sts) %in% names)
stsAux <- unlist(sts@data)[ind]
} else {
stsAux <- sts@data
}
}
if (!dynamicPlot) {
# Plotting base map
world <- rworldmap::getMap(resolution = "high")
maps::map(world,
fill = TRUE,
col = groundColor,
bg = oceanColor,
wrap = c(0, 360),
xlim = c(ext$xmin - 1, ext$xmax + 1), # we extend W & E by 1°
ylim = c(ext$ymin - 1, ext$ymax + 1)
) # we extend S & N by 1°
maps::map.axes(cex.axis = 1)
# Plotting loi
if (loi) {
plot(sts@spatialLoiBuffer, lwd = 1, border = "blue", add = TRUE)
}
if (loi && as.character(sf::st_geometry_type(sts@spatialLoi)) == "POINT") {
plot(sts@spatialLoi, lwd = 2, col = "blue", pch = 4, add = TRUE)
}
# Plotting track(s) and labels
lapply(stsAux, plotTrack, sts@scale, sts@scalePalette)
if (labels) {
lapply(stsAux, plotLabels, by, pos)
}
# Adding legends
if (legends != "none") {
leg <- names(sts@scalePalette)
col <- sts@scalePalette
lty <- rep(0, length(col))
pch <- rep(19, length(col))
lwd <- rep(1, length(col))
graphics::legend(legends,
title = "Scale",
legend = leg,
col = col,
lty = lty,
pch = pch,
lwd = lwd,
cex = 0.6,
bty = "n"
)
}
}else {
#Init map
map <- leaflet::fitBounds(
leaflet::addProviderTiles(
leaflet::leaflet(options =
leaflet::leafletOptions(worldCopyJump = F,
minZoom = 2,
maxZoom = 12),
width = 650,
height = 700),
leaflet::providers$Esri.NatGeoWorldMap,
group = "Satellite",
options = leaflet::providerTileOptions(errorTileUrl = "Tile not found")),
lng1 = as.numeric(ext$xmin),
lng2 = as.numeric(ext$xmax),
lat1 = as.numeric(ext$ymin),
lat2 = as.numeric(ext$ymax)
)
#Plotting loi
if(loi)
map <- leaflet::addPolylines(map,
data = sts@spatialLoiBuffer,
fillColor = "transparent",
label = "Buffer limit")
if(loi & as.character(sf::st_geometry_type(sts@spatialLoi)) == "POINT")
map <- leaflet::addCircleMarkers(map,
data = sts@spatialLoi,
fillColor = "transparent",
label = "LOI")
#Plotting track(s) and labels
for (st in stsAux) {
data <- st@obs.all
data$type <- unlist(lapply(data$msw,getColors, sts@scale, sts@scalePalette))
labels <- paste0("<b>",st@name, "</b><br>",
"<i>observation ",row.names(data), "</i><br>",
data$iso.time,
"<ul>",
"<li>", "scale: ", data$scale, "</li>",
"<li>", "msw: ", data$msw, "m/s</li>",
"<li>", "rmw: ", data$rmw, "km</li>",
"<li>", "pressure: ", data$pressure, "pa</li>",
"<li>", "poci: ", data$poci, "pa</li>",
"</ul>")
#Plot track
map <- leaflet::addPolylines(map,
data = data,
lng = ~lon,
lat = ~lat,
color = "black",
weight = 2)
map <- leaflet::addCircleMarkers(map,
data = data,
lng = ~lon,
lat = ~lat,
radius = 5,
color = ~type,
stroke = FALSE,
fillOpacity = 1,
popup = labels
)
}
#Adding legends
map <- leaflet::addLegend(map,
legends,
colors = sts@scalePalette,
labels = names(sts@scalePalette),
title = "Scale",
opacity = 1)
map
return(map)
}
}