Nothing
#' Check inputs for temporalBehaviour function
#'
#' @noRd
#' @param sts StormsList object
#' @param points data.frame
#' @param product character
#' @param windThreshold numeric
#' @param method character
#' @param asymmetry character
#' @param tempRes numeric
#' @param verbose logical
#' @return NULL
checkInputsTemporalBehaviour <- function(sts, points, product, windThreshold, method, asymmetry,
empiricalRMW, tempRes, verbose) {
# Checking sts input
stopifnot("no data found" = !missing(sts))
# Checking points input
stopifnot("no data found" = !missing(points))
stopifnot("points must be data.frame" = identical(class(points), "data.frame"))
stopifnot(
"colnames of points must be \"x\" (Eastern degree), \"y\" (Northern degree)" = colnames(points) == c("x", "y")
)
stopifnot("Invalid points coordinates" = points$x > -180 & points$x <= 360 &
points$y >= -90 & points$y <= 90)
# Checking product input
stopifnot("Invalid product" = product %in% c("TS", "PDI", "Exposure"))
stopifnot("Only one product must be chosen" = length(product) == 1)
# Checking windThreshold input
if ("Exposure" %in% product) {
stopifnot("windThreshold must be numeric" = identical(class(windThreshold), "numeric"))
stopifnot("invalid value(s) for windThreshold input (must be > 0)" = windThreshold > 0)
}
# Checking method input
stopifnot("Invalid method input" = method %in% c("Willoughby", "Holland", "Boose"))
stopifnot("Only one method must be chosen" = length(method) == 1)
# Checking asymmetry input
stopifnot("Invalid asymmetry input" = asymmetry %in% c("None", "Chen", "Miyazaki"))
stopifnot("Only one asymmetry must be chosen" = length(asymmetry) == 1)
# Checking empiricalRMW input
stopifnot("empiricalRMW must be logical" = identical(class(empiricalRMW), "logical"))
# Checking tempRes input
stopifnot("tempRes must be numeric" = identical(class(tempRes), "numeric"))
stopifnot("tempRes must be length 1" = length(tempRes) == 1)
stopifnot("invalid tempRes: must be either 60, 30 or 15" = tempRes %in% c(60, 30, 15))
# Checking verbose input
stopifnot("verbose must be numeric" = identical(class(verbose), "numeric"))
stopifnot("verbose must length 1" = length(verbose) == 1)
stopifnot("verbose must be either 0, 1 or 2" = verbose %in% c(0, 1, 2))
}
#' rasterizePDI counterpart function for non raster data
#'
#' @noRd
#' @param wind numeric vector. Wind speed values
#' @param tempRes numeric. Time resolution, used for the numerical integration
#' over the whole track
#'
#' @return numeric. PDI computed using the wind speed values in wind
computePDI <- function(wind, tempRes) {
# Computing surface drag coefficient
rho <- 1
cd <- 0.002
# Raising to power 3
pdi <- wind**3
# Applying both rho and surface drag coefficient
pdi <- pdi * rho * cd
# Integrating over the whole track
pdi <- sum(pdi, na.rm = TRUE) * tempRes
return(round(pdi, 3))
}
#' rasterizeExposure counterpart function for non raster data
#'
#' @noRd
#' @param wind numeric vector. Wind speed values
#' @param tempRes numeric. Time resolution, used for the numerical integration
#' over the whole track
#' @param threshold numeric vector. Wind threshold
#'
#' @return numeric vector of length 5 (for each category). Exposure computed
#' using the wind speed values in wind
computeExposure <- function(wind, tempRes, threshold) {
exposure <- c()
for (t in threshold) {
ind <- which(wind >= t)
expo <- rep(0, length(wind))
expo[ind] <- 1
exposure <- c(exposure, sum(expo, na.rm = TRUE) * tempRes)
}
return(exposure)
}
#' rasterizeProduct counterpart function for non raster data
#'
#' @noRd
#' @param product character. Product input from temporalBehaviour
#' @param wind numeric vector. Wind speed values
#' @param direction numeric vector. Wind direction
#' @param tempRes numeric. Time resolution, used for the numerical integration
#' over the whole track
#' @param result numeric array. Similar as finalStack, i.e where to add the
#' computed product
#' @param threshold numeric vector. Wind threshold
#'
#' @return numeric array of dimension:
#' \itemize{
#' \item number of observations: number of points. If product == "MSW"
#' \item 1 : number of points. If product == "PDI"
#' \item 5 : number of points. If product == "Exposure"
#' }
computeProduct <- function(product, wind, direction, tempRes, result, threshold) {
if (product == "TS") {
prod <- cbind(wind, direction)
} else if (product == "PDI") {
prod <- computePDI(wind, tempRes)
} else if (product == "Exposure") {
prod <- computeExposure(wind, tempRes, threshold)
}
return(cbind(result, prod))
}
#' Arrange result before the end of temporalBehaviour
#'
#' @noRd
#' @param finalResult list of data.frame. Where to add the computed product
#' @param result result output from computeProduct function
#' @param product character. Product input from temporalBehaviour
#' @param points points input from temporalBehaviour
#' @param isoT numeric vector. Iso Times of observations
#' @param indices numeric vector. Indices of observations
#' @param st Storm object.
#' @param threshold numeric vector. Wind threshold
#'
#' @return finalResult
finalizeResult <- function(finalResult, result, product, points, isoT, indices, st, threshold) {
if (product == "TS") {
l <- list()
for (i in seq_len(dim(points)[1])) {
j <- 2 * (i - 1) + 1
df <- data.frame(result[, j], result[, j:j + 1], indices = indices, isoTimes = isoT)
colnames(df) <- c("speed", "direction", "indices", "isoTimes")
l <- append(l, list(df))
}
names(l) <- rownames(points)
} else if (product == "PDI") {
l <- data.frame(result, row.names = "PDI")
colnames(l) <- rownames(points)
} else {
l <- data.frame(result, row.names = paste("Min threshold:", threshold, "m/s"))
colnames(l) <- rownames(points)
}
dfn <- list(l)
names(dfn) <- st@name
finalResult <- append(finalResult, dfn)
return(finalResult)
}
############################
# temporalBehaviour function#
############################
#' Computing wind behaviour time series and summary statistics at given point locations
#'
#' The `temporalBehaviour()` function allows computing wind speed and direction
#' for a given location or set of locations along the lifespan of a tropical cyclone.
#' It also allows to compute three associated summary statistics.
#'
#' @param sts `StormsList` object.
#' @param points data.frame. Consisting of two columns names as "x" (for the longitude) and
#' "y" (for the latitude), providing the coordinates in decimal degrees of the point locations. Row names
#' can also be provided to named the locations.
#' @param product character. Desired output statistics:
#' \itemize{
#' \item `"TS"`, for time series of wind speeds and directions (default setting),
#' \item `"PDI"`, for power dissipation index, or
#' \item `"Exposure"`, for the duration of exposure to defined wind thresholds.
#' }
#' @param windThreshold numeric vector. Minimal wind threshold(s) (in \eqn{m.s^{-1}}) used to
#' compute the duration of exposure when `product="Exposure"`. Default value is to set NULL, in this
#' case, the windthresholds are the one used in the scale defined in the stromsList.
#' @param method character. Model used to compute wind speed and direction.
#' Three different models are implemented:
#' \itemize{
#' \item `"Willoughby"`, for the symmetrical model developed by Willoughby et al. (2006) (default setting),
#' \item `"Holland"`, for the symmetrical model developed by Holland (1980), or
#' \item `"Boose"`, for the asymmetrical model developed by Boose et al. (2004).
#' }
#' @param asymmetry character. If `method="Holland"` or `method="Willoughby"`,
#' this argument specifies the method used to add asymmetry. Can be:
#' \itemize{
#' \item `"Chen"`, for the model developed by Chen (1994) (default setting),
#' \item `"Miyazaki"`, for the model developed by Miyazaki et al. (1962), or
#' \item `"None"`, for no asymmetry.
#' }
#' @param empiricalRMW logical. Whether (TRUE) or not (FALSE) to compute
#' the radius of maximum wind (`rmw`) empirically using the model developed by
#' Willoughby et al. (2006). If `empiricalRMW==FALSE` (default setting) then the
#' `rmw` provided in the `StormsList` is used.
#' @param tempRes numeric. Temporal resolution (min). Can be `60` (default setting),
#' `30` or `15`.
#' @param verbose numeric. Information displayed. Can be:
#' \itemize{
#' \item `2`, information about the processes and outputs are displayed (default setting),
#' \item `1`, information about the processes are displayed, or
#' \item `0`, no information displayed.
#' }
#' @returns For each storm and each point location, the `temporalBehaviour()` function returns
#' a data.frame. The data frames are organised into named lists. Depending on the `product` argument
#' different data.frame are returned:
#' \itemize{
#' \item if `product == "TS"`, the function returns a data.frame with
#' one row for each observation (or interpolated observation) and
#' four columns for wind speed (in \eqn{m.s^{-1}}), wind direction (in degree),
#' the observation number, and the ISO time of observations,
#' \item if `product == "PDI"`, the function returns a data.frame with one row
#' for each point location and one column for the PDI,
#' \item if `product == "Exposure"`, the function returns a data.frame with one
#' row for the duration of exposure to winds above each wind speed threshold defined
#' by the `windThreshold` argument and one column for each point location.
#' }
#'
#' @details Storm track data sets, such as those extracted from IBRTrACKS (Knapp et
#' al., 2010), usually provide observation at a 3- or 6-hours temporal
#' resolution. In the temporalBehaviour() function, linear interpolations are
#' used to reach the temporal resolution specified in the `tempRes` argument
#' (default value = 60 min).
#'
#' The Holland (1980) model, widely used in the literature, is based on the
#' 'gradient wind balance in mature tropical cyclones. The wind speed distribution
#' is computed from the circular air pressure field, which can be derived from
#' the central and environmental pressure and the radius of maximum winds.
#'
#' \eqn{v_r = \sqrt{\frac{b}{\rho} \times \left(\frac{R_m}{r}\right)^b \times (p_{oci} - p_c)
#' \times e^{-\left(\frac{R_m}{r}\right)^b} + \left(\frac{r \times f}{2}\right)^2} -
#' \left(\frac{r \times f}{2}\right)}
#'
#' with,
#'
#' \eqn{b = \frac{\rho \times e \times v_m^2}{p_{oci} - p_c}}
#'
#' \eqn{f = 2 \times 7.29 \times 10^{-5} \sin(\phi)}
#'
#' where, \eqn{v_r} is the tangential wind speed (in \eqn{m.s^{-1}}),
#' \eqn{b} is the shape parameter,
#' \eqn{\rho} is the air density set to \eqn{1.15 kg.m^{-3}},
#' \eqn{e} being the base of natural logarithms (~2.718282),
#' \eqn{v_m} the maximum sustained wind speed (in \eqn{m.s^{-1}}),
#' \eqn{p_{oci}} is the pressure at outermost closed isobar of the storm (in \eqn{Pa}),
#' \eqn{p_c} is the pressure at the centre of the storm (in \eqn{Pa}),
#' \eqn{r} is the distance to the eye of the storm (in \eqn{km}),
#' \eqn{R_m} is the radius of maximum sustained wind speed (in \eqn{km}),
#' \eqn{f} is the Coriolis force (in \eqn{N.kg^{-1}}, and
#' \eqn{\phi} being the latitude).
#'
#' The Willoughby et al. (2006) model is an empirical model fitted to aircraft
#' observations. The model considers two regions: inside the eye and at external
#' radii, for which the wind formulations use different exponents to better match
#' observations. In this model, the wind speed increases as a power function of the
#' radius inside the eye and decays exponentially outside the eye after a smooth
#' polynomial transition across the eyewall.
#'
#' \eqn{\left\{\begin{aligned}
#' v_r &= v_m \times \left(\frac{r}{R_m}\right)^{n} \quad if \quad r < R_m \\
#' v_r &= v_m \times \left((1-A) \times e^{-\frac{|r-R_m|}{X1}} +
#' A \times e^{-\frac{|r-R_m|}{X2}}\right) \quad if \quad r \geq R_m \\
#' \end{aligned}
#' \right.
#' }
#'
#' with,
#'
#' \eqn{n = 2.1340 + 0.0077 \times v_m - 0.4522 \times \ln(R_m) - 0.0038 \times |\phi|}
#'
#' \eqn{X1 = 287.6 - 1.942 \times v_m + 7.799 \times \ln(R_m) + 1.819 \times |\phi|}
#'
#' \eqn{A = 0.5913 + 0.0029 \times v_m - 0.1361 \times \ln(R_m) - 0.0042 \times |\phi|} and \eqn{A\ge0}
#'
#' where, \eqn{v_r} is the tangential wind speed (in \eqn{m.s^{-1}}),
#' \eqn{v_m} is the maximum sustained wind speed (in \eqn{m.s^{-1}}),
#' \eqn{r} is the distance to the eye of the storm (in \eqn{km}),
#' \eqn{R_m} is the radius of maximum sustained wind speed (in \eqn{km}),
#' \eqn{\phi} is the latitude of the centre of the storm,
#' \eqn{X2 = 25}.
#'
#' Asymmetry can be added to Holland (1980) and Willoughby et al. (2006) wind fields as follows,
#'
#' \eqn{\vec{V} = \vec{V_c} + C \times \vec{V_t}}
#'
#' where, \eqn{\vec{V}} is the combined asymmetric wind field,
#' \eqn{\vec{V_c}} is symmetric wind field,
#' \eqn{\vec{V_t}} is the translation speed of the storm, and
#' \eqn{C} is function of \eqn{r}, the distance to the eye of the storm (in \eqn{km}).
#'
#' Two formulations of C proposed by Miyazaki et al. (1962) and Chen (1994) are implemented.
#'
#' Miyazaki et al. (1962)
#' \eqn{C = e^{(-\frac{r}{500} \times \pi)}}
#'
#' Chen (1994)
#' \eqn{C = \frac{3 \times R_m^{\frac{3}{2}} \times
#' r^{\frac{3}{2}}}{R_m^3 + r^3 +R_m^{\frac{3}{2}} \times r^{\frac{3}{2}}}}
#'
#' where, \eqn{R_m} is the radius of maximum sustained wind speed (in \eqn{km})
#'
#' The Boose et al. (2004) model, or “HURRECON” model, is a modification of the
#' Holland (1980) model. In addition to adding
#' asymmetry, this model treats of water and land differently, using different
#' surface friction coefficient for each.
#'
#' \eqn{v_r = F\left(v_m - S \times (1 - \sin(T)) \times \frac{v_h}{2} \right)
#' \times \sqrt{\left(\frac{R_m}{r}\right)^b \times e^{1 - \left(\frac{R_m}{r}\right)^b}}}
#'
#' with,
#'
#' \eqn{b = \frac{\rho \times e \times v_m^2}{p_{oci} - p_c}}
#'
#' where, \eqn{v_r} is the tangential wind speed (in \eqn{m.s^{-1}}),
#' \eqn{F} is a scaling parameter for friction (\eqn{1.0} in water, \eqn{0.8} in land),
#' \eqn{v_m} is the maximum sustained wind speed (in \eqn{m.s^{-1}}),
#' \eqn{S} is a scaling parameter for asymmetry (usually set to \eqn{1}),
#' \eqn{T} is the oriented angle (clockwise/counter clockwise in Northern/Southern Hemisphere) between
#' the forward trajectory of the storm and a radial line from the eye of the storm to point $r$
#' \eqn{v_h} is the storm velocity (in \eqn{m.s^{-1}}),
#' \eqn{R_m} is the radius of maximum sustained wind speed (in \eqn{km}),
#' \eqn{r} is the distance to the eye of the storm (in \eqn{km}),
#' \eqn{b} is the shape parameter,
#' \eqn{\rho = 1.15} is the air density (in \eqn{kg.m^{-3}}),
#' \eqn{p_{oci}} is the pressure at outermost closed isobar of the storm (in \eqn{Pa}), and
#' \eqn{p_c} is the pressure at the centre of the storm (\eqn{pressure} in \eqn{Pa}).
#'
#' @references
#' Boose, E. R., Serrano, M. I., & Foster, D. R. (2004). Landscape and regional impacts of hurricanes in Puerto Rico.
#' Ecological Monographs, 74(2), Article 2. https://doi.org/10.1890/02-4057
#'
#' Chen, K.-M. (1994). A computation method for typhoon wind field. Tropic Oceanology, 13(2), 41–48.
#'
#' Holland, G. J. (1980). An Analytic Model of the Wind and Pressure
#' Profiles in Hurricanes. Monthly Weather Review, 108(8), 1212–1218.
#' https://doi.org/10.1175/1520-0493(1980)108<1212:AAMOTW>2.0.CO;2
#'
#' Knapp, K. R., Kruk, M. C., Levinson, D. H., Diamond, H. J., & Neumann, C. J. (2010). The
#' International Best Track Archive for Climate Stewardship (IBTrACS).
#' Bulletin of the American Meteorological Society, 91(3), Article 3. https://doi.org/10.1175/2009bams2755.1
#'
#' Miyazaki, M., Ueno, T., & Unoki, S. (1962). The theoretical investigations of typhoon surges
#' along the Japanese coast (II).
#' Oceanographical Magazine, 13(2), 103–117.
#'
#' Willoughby, H. E., Darling, R. W. R., & Rahn, M. E. (2006). Parametric Representation of the
#' Primary Hurricane Vortex. Part II: A New Family of Sectionally Continuous Profiles.
#' Monthly Weather Review, 134(4), 1102–1120. https://doi.org/10.1175/MWR3106.1
#'
#' @examples
#' \donttest{
#' # Creating a stormsDataset
#' sds <- defStormsDataset()
#'
#' # Geting storm track data for tropical cyclone Pam (2015) near Vanuatu
#' pam <- defStormsList(sds = sds, loi = "Vanuatu", names = "PAM")
#'
#' pts <- data.frame(x = c(168.5, 168), y = c(-17.9, -16.3))
#' row.names(pts) <- c("point_1", "point_2")
#'
#' # Computing time series of wind speed and direction for Pam
#' # over points 1 and 2 defined above
#' ts.pam <- temporalBehaviour(pam, points = pts)
#'
#' # Computing PDI for Pam over points 1 and 2 defined above
#' pdi.pam <- temporalBehaviour(pam, points = pts, product = "PDI")
#'
#' # Computing the duration of exposure to wind speeds above the thresholds
#' # used by the Saffir-Simpson hurricane wind scale for Pam
#' # over points 1 and 2 defined above
#' exp.pam <- temporalBehaviour(pam, points = pts, product = "Exposure")
#' }
#'
#' @export
temporalBehaviour <- function(sts,
points,
product = "TS",
windThreshold = NULL,
method = "Willoughby",
asymmetry = "Chen",
empiricalRMW = FALSE,
tempRes = 60,
verbose = 1) {
if (is.null(windThreshold)) {
windThreshold <- sts@scale
}
checkInputsTemporalBehaviour(sts, points, product, windThreshold, method, asymmetry, empiricalRMW, tempRes, verbose)
if (verbose > 0) {
cat("=== temporalBehaviour processing ... ===\n\n")
cat("Initializing data ...")
}
if (method == "Boose") {
# Map for intersection
world <- rworldmap::getMap(resolution = "high")
world <- sf::st_as_sf(world)
world <- sf::st_transform(world, crs = wgs84)
indCountries <- which(sf::st_intersects(sts@spatialLoiBuffer, world$geometry, sparse = FALSE) == TRUE)
asymmetry <- "None"
} else {
indCountries <- NULL
}
points$x[points$x < 0] <- points$x[points$x < 0] + 360
buffer <- 2.5
# Initializing final result
finalResult <- list()
if (verbose > 0) {
cat(" Done\n\n")
cat("Computation settings:\n")
cat(" (*) Temporal resolution: Every", tempRes, "min\n")
cat(" (*) Method used:", method, "\n")
cat(" (*) Product(s) to compute:", product, "\n")
cat(" (*) Asymmetry used:", asymmetry, "\n")
if (empiricalRMW) {
cat(" (*) rmw computed according to the empirical formula (See Details section)")
}
cat(" (*) Points: lon-lat\n")
for (i in seq_len(dim(points)[1])) {
cat(" ", points$x[i], " ", points$y[i], "\n")
}
cat("\n")
cat("\nStorm(s):\n")
cat(" (", getNbStorms(sts), ") ", getNames(sts), "\n\n")
}
for (st in sts@data) {
# Handling indices inside loi.buffer or not
ind <- getIndices(st, 2, "none")
# Getting data associated with storm st
dataTC <- getDataInterpolate(st, ind, tempRes, empiricalRMW, method)
# Computing distances from the eye of storm for every observations x, and
# every points y
distEye <- terra::distance(
x = cbind(dataTC$lon, dataTC$lat),
y = cbind(points$x, points$y),
lonlat = TRUE
)
res <- c()
# For each point
for (i in seq_len(dim(points)[1])) {
# Coordinates
pt <- points[i, ]
# Computing distance between eye of storm and point P
x <- pt$x - dataTC$lon
y <- pt$y - dataTC$lat
# Computing wind speed/direction
dist2p <- distEye[, i]
output <- computeWindProfile(
dataTC, i, method, asymmetry, x, y, pt, dist2p,
buffer, sts@spatialLoiBuffer,
world, indCountries
)
vr <- output$wind
dir <- output$direction
# Computing product
res <- computeProduct(product, vr, dir, tempRes, res, windThreshold)
}
finalResult <- finalizeResult(
finalResult, res, product, points,
dataTC$isoTimes, dataTC$indices, st,
windThreshold
)
}
if (verbose > 0) {
cat("\n=== DONE ===\n\n")
if (verbose > 1) {
cat("Output:\n")
cat("DataFrame with", length(finalResult), "storm:\n")
cat("index - name of the storm - Observation Point - Point indices\n")
stormNames <- names(finalResult)
for (i in seq_along(stormNames)) {
stormName <- stormNames[i]
pointNames <- names(finalResult[[stormName]])
for (j in seq_along(pointNames)) {
pointName <- pointNames[j]
cat(
" ", i, " ",
stormName, " ",
pointName, " ",
finalResult[[stormName]][[pointName]]$indices, "\n"
)
}
}
cat("\n")
}
}
return(finalResult)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.