113 incorrect units in pdi calculation and figures

[thomasarsouze]

@thomaspibanez : can you double-check that :

• the bugs described are correct
• the fixes for PDI and Exposure are correct

[thomaspibanez]

@thomasarsouze

For PDI calculation, indeed, it looks like PDI values were too high (in our paper with average annual values computed by Christophe we obtained PDI values up to 1000: https://besjournals.onlinelibrary.wiley.com/doi/10.1111/1365-2745.13039?utm_medium=article&utm_source=researchgate.net). Yes it should be J.m-2 and not J.m

For the duration of exposure, I don't know which wind speed was used in the example provided but yes it looks like 250 hours (>10 days) is overly too long. I checked in the code (below) and yes it is min (tempRes = 15, 30, or 60 min). So we need to change the unit and also specify in the documentation that the output is in min.

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)
}

[thomasarsouze]

@thomaspibanez :

• for PDI, I'm not sure to get it, if we correct and use "s" instead of "min", as it should be, this will increase PDI by one order of magnitude, right ? cf. last line of the code snippet:

 # Surface sea-level air density in kg.m-3
  rho <- 1
  # Surface drag coefficient
  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 and converting minutes to seconds
  pdi <- sum(pdi, na.rm = TRUE) * tempRes * 60

• for exposure time, in the documentation it was stated that the exposure time was given in hours, but the result was in minutes. I've chosen to correct the code and put everything in "h". You suggest that we should keep "min". I don't care, but you're the one more used to community standards, so should we keep "h" and correct the calculation, or keep "min" and don't change the code but change the documentation ?

[thomasarsouze]

@thomaspibanez : can you confirm the points mentioned above so that we can merge the fix ?

[thomaspibanez]

@thomasarsouze : I've discussed with C. Menkes, we would prefer exposure time in hours, so we have to divide by 60 in the calculation

Regarding the PDI Christophe is ok with the formula currently implemented in StromR and tempRes is in min so we have to multiply by 60 to have it in seconds.

That said I agree that there is an issue with PDI. If we take the example we show in the article for spatialBehaviour for PAM in Vanuatu we have exposure time of about 200 min for wind if 58 m.s if we apply the formula we have 20060(58**3)*0.002 = 4, 682, 688 which is much higher than the 250.000 (and this a conservative estimate)

[thomaspibanez]

@thomasarsouze

I was looking at the code to find the issue with the PDI. I saw that some part of the code are duplicated is it normal?

#' Compute PDI raster
#'
#' @nord
#' @param finalStack list of SpatRaster. Where to add the computed MSW raster
#' @param tempRes numeric. Time resolution, used for the numerical integration
#' over the whole track
#' @param stack SpatRaster stack. All the PDI rasters used to compute MSW
#' @param name character. Name of the storm. Used to give the correct layer name
#' in finalStack
#' @param spaceRes character. spaceRes input from spatialBehaviour
#' @param threshold numeric vector. Wind threshold
#'
#'
#' @return list of SpatRaster
rasterizePDI <- function(finalStack, stack, tempRes, spaceRes, name, threshold) {
nbg <- switch(spaceRes,
"30sec" = 59,
"2.5min" = 11,
"5min" = 5,
"10min" = 3
)

Integrating over the whole track

prod <- sum(stack, na.rm = TRUE) * tempRes

Applying focal function to smooth results

prod <- terra::focal(prod, w = matrix(1, nbg, nbg), mean, na.rm = TRUE, pad = TRUE)
names(prod) <- paste0(name, "_PDI")

return(c(finalStack, prod))
}

#' Compute PDI raster
#'
#' @nord
#' @param finalStack list of SpatRaster. Where to add the computed MSW raster
#' @param tempRes numeric. Time resolution, used for the numerical integration
#' over the whole track
#' @param stack SpatRaster stack. All the PDI rasters used to compute MSW
#' @param name character. Name of the storm. Used to give the correct layer name
#' in finalStack
#' @param spaceRes character. spaceRes input from spatialBehaviour
#' @param threshold numeric vector. Wind threshold
#'
#'
#' @return list of SpatRaster
rasterizeExp <- function(finalStack, stack, tempRes, spaceRes, name, threshold) {
nbg <- switch(spaceRes,
"30sec" = 59,
"2.5min" = 11,
"5min" = 5,
"10min" = 3
)

for (l in seq_along(threshold)) {
ind <- seq(l, terra::nlyr(stack), length(threshold))
# Integrating over the whole track
prod <- sum(terra::subset(stack, ind), na.rm = TRUE) * tempRes
# Applying focal function to smooth results
prod <- terra::focal(prod, w = matrix(1, nbg, nbg), mean, na.rm = TRUE)
names(prod) <- paste0(name, "Exposure", threshold[l])
finalStack <- c(finalStack, prod)
}

return(finalStack)
}