update threat pollution layers

R/fnc_diffusive_model.R

@@ -0,0 +1,68 @@
+#' Diffusive model
+#'
+#' Function to evaluate the influence area of a stressor using a passive diffusive model
+#'
+#' @param dat sf object with points characterizing stressor intensity
+#' @param field character, name of field containing stresseor intensity
+#' @param threshold numeric, minimum threshold in percent of the global maximum at which to stop the diffusive model
+#' @param globalmaximum global maximum value of the stressor in the study area
+#' @param decay percent by which the value of the stressor is reduced when diffusing to the adjacent cells
+#' @param increment numeric, distance between successive steps in the model. Units are in the units of the spatial object provided.
+#'
+#' @keywords cumulative footprint
+#'
+#' @export
+#'
+#' @details
+#'
+
+diffusive <- function(dat, field, threshold, globalmaximum, decay, increment) {
+  # -----
+  val <- res <- list()
+  for (i in 1:nrow(dat)) {
+    val[[i]] <- seq(
+      from = as.numeric(dat[i, field, drop = TRUE]),
+      to = as.numeric(globalmaximum * (threshold / 100)),
+      by = -(decay / 100)
+    )
+    nVal <- length(val[[i]])
+    temp <- dplyr::select(dat[i, ], geometry)
+
+
+    rings <- circles <- list()
+    if (nVal > 1) {
+      # Buffers and concentric circles for passive diffusion
+      circles[[1]] <- rings[[1]] <- sf::st_buffer(temp, increment)
+      for (j in 2:nVal) {
+        circles[[j]] <- sf::st_buffer(circles[[(j - 1)]], increment)
+        rings[[j]] <- sf::st_difference(circles[[j]], circles[[(j - 1)]])
+      }
+
+      # Add intensity values and rasterize
+      res[[i]] <- dplyr::bind_rows(rings) |>
+        dplyr::mutate(intensity = val[[i]])
+    } else if (nVal == 1) {
+      res[[i]] <- sf::st_buffer(temp, increment) |>
+        dplyr::mutate(intensity = val[[i]])
+    }
+  }
+
+  # Combine and rasterize
+  r <- as(
+    stars::read_stars("project-data/grid/grid.tif"),
+    "Raster"
+  )
+  res <- res |>
+    dplyr::bind_rows() |>
+    sf::st_transform(sf::st_crs(r)) |>
+    fasterize::fasterize(
+      r,
+      field = "intensity",
+      fun = "sum"
+    ) |>
+    stars::st_as_stars() |>
+    mask_aoi()
+
+  # Return
+  res
+}

R/threat_pollution.R

@@ -81,15 +81,17 @@ threat_pollution <- function() {
   neec <- neec[iid, ]
 
   # Remove downloaded data
-  fs::dir_delete(tmp)
+  # fs::dir_delete(tmp)
   # ----------------------------------------------------------------
   # To explore:
   # Impact to Waterbody field.  Consider if we should treat “actual” vs. “potential” spills differently.  Data on “potential” spills only started being recorded in April 2021 onward (no differentiation before that).  Some “Potential” spills may report very large amounts of substance that was not actually spilled (e.g.  Incident Number NL-20200611-03045-20 was a disabled vessel with 400000 L of fuel on board, but nothing spilled).  These “Potential” spills still represent risk even if no true pollution occurred.
 
   # =~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~= #
   # NASP data
   # =~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~= #
-  nasp <- pipedat::importdat("376f0891", "format")[[1]]
+  nasp <- pipedat::importdat("376f0891", "format")[[1]] |>
+    dplyr::filter(!is.na(volume_l)) |>
+    dplyr::rename(geometry = geom)
 
   # =~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~= #
   # ISTOP data
@@ -98,12 +100,119 @@ threat_pollution <- function() {
   istop <- pipedat::importdat("48ea8a05", "format")[[1]] |>
     dplyr::filter(!is.na(date)) |>
     dplyr::distinct() |>
-    dplyr::filter(category %in% c("1A", "1B", "B"))
+    dplyr::filter(category %in% c("1A", "1B", "B")) |>
+    dplyr::rename(geometry = geom) |>
+    dplyr::mutate(
+      area = sf::st_area(geometry),
+      area = units::set_units(area, km^2),
+      area = as.numeric(area)
+    )
 
   # =~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~= #
-  # Threat layer TODO
+  # Threat layer - diffusive model
   # =~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~= #
-  neec
-  nasp
-  istop
+  # Using diffusive model function to generate threat layer for each dataset
+  # https://github.com/EffetsCumulatifsNavigation/ceanav/blob/main/R/fnc_diffusive_model.R
+  # Parameters
+  threshold <- .05
+  decay <- 2
+  increment <- 100
+
+  # NEEC
+  neec_tmp <- neec |>
+    dplyr::mutate(
+      quantity = log(quantity + 1),
+      quantity = quantity / max(quantity)
+    ) |>
+    sf::st_transform(3348)
+  neec_threat <- diffusive(
+    dat = neec_tmp,
+    field = "quantity",
+    threshold = threshold,
+    globalmaximum = min(neec_tmp$quantity),
+    decay = decay,
+    increment = increment
+  )
+
+  # NASP
+  nasp_tmp <- nasp |>
+    dplyr::mutate(
+      volume_l = log(volume_l + 1),
+      volume_l = volume_l / max(volume_l)
+    ) |>
+    sf::st_transform(3348)
+  nasp_threat <- diffusive(
+    dat = nasp_tmp,
+    field = "volume_l",
+    threshold = threshold,
+    globalmaximum = min(nasp_tmp$volume_l),
+    decay = decay,
+    increment = increment
+  )
+
+  # ISTOP
+  istop_tmp <- istop |>
+    dplyr::mutate(
+      area = log(area + 1),
+      area = area / max(area)
+    ) |>
+    sf::st_transform(3348)
+  istop_threat <- diffusive(
+    dat = istop_tmp,
+    field = "area",
+    threshold = threshold,
+    globalmaximum = min(istop_tmp$area),
+    decay = decay,
+    increment = increment
+  )
+
+  # Export
+  out <- here::here("output", "threats")
+  pipedat::chk_create(out)
+  pipedat::masterwrite(neec_threat, here::here(out, "pollution_neec"))
+  pipedat::masterwrite(nasp_threat, here::here(out, "pollution_nasp"))
+  pipedat::masterwrite(istop_threat, here::here(out, "pollution_istop"))
+
+
+  # # ----
+  # KERNEL EXPLORATION
+  # make_kernel <- function(dat, sig, wts) {
+  #   # Window
+  #   global_parameters()
+  #   bbox <- unlist(param$bbox$base)[c("xmin", "xmax", "ymin", "ymax")] |>
+  #     pipedat::bbox_poly(param$crs) |>
+  #     sf::st_transform(3348) |>
+  #     sf::st_bbox()
+  #   wdw <- spatstat.geom::owin(xrange = c(bbox["xmin"], bbox["xmax"]), yrange = c(bbox["ymin"], bbox["ymax"]))
+
+  #   # Points
+  #   xy <- sf::st_transform(dat, 3348) |>
+  #     sf::st_coordinates()
+  #   pts <- spatstat.geom::ppp(
+  #     x = xy[, "X"],
+  #     y = xy[, "Y"],
+  #     window = wdw,
+  #     marks = dat[, wts, drop = TRUE]
+  #   )
+
+  #   # Kernel
+  #   kern <- spatstat.explore::density.ppp(
+  #     x = pts,
+  #     sigma = sig,
+  #     # weights = log(pts$marks + 1),
+  #     weights = pts$marks,
+  #     edge = TRUE,
+  #     kernel = "gaussian"
+  #   )
+
+  #   # Spatial objet and mask to study area
+  #   kern |>
+  #     terra::rast() |>
+  #     stars::st_as_stars() |>
+  #     sf::st_set_crs(3348) #|>
+  #   # pipedat::masteringrid()
+  # }
+
+  # x <- make_kernel(neec, 10000, "quantity")
+  # plot(x)
 }

R/zzz.R

@@ -24,9 +24,9 @@ rep_hyperlien <- function(texte, url) {
   nl <- length(texte)
   hyperlien <- character(nl)
 
-  for(i in 1:nl) {
-    if(!is.null(url[i])) {
-      hyperlien[i] <- paste0("[",texte[i],"](",url[i],")")
+  for (i in 1:nl) {
+    if (!is.null(url[i])) {
+      hyperlien[i] <- paste0("[", texte[i], "](", url[i], ")")
     } else {
       hyperlien[i] <- texte[i]
     }
@@ -49,3 +49,16 @@ clean <- function() {
 chk_create <- function(path) {
   if (!file.exists(path)) dir.create(path, recursive = TRUE)
 }
+
+# ------------------------------------------------------------------------------
+# Function to mask on area of interest
+mask_aoi <- function(dat) {
+  # Area of interest
+  aoi <- sf::st_read("project-data/aoi/aoi.gpkg", quiet = TRUE)
+  tmp <- dat
+  tmp[[1]][] <- NA
+  aoi$val_ras <- 1
+  aoi_mask <- stars::st_rasterize(aoi["val_ras"], template = tmp)
+  stars::st_dimensions(aoi_mask) <- stars::st_dimensions(dat)
+  dat * aoi_mask
+}