Add barrier.R

R/barrier.R

@@ -0,0 +1,155 @@
+#' Transform a mesh object into a mesh with correlation barriers
+#'
+#' @param spde_obj Output from [make_mesh()].
+#' @param barrier_sf An sf object with polygons defining the barriers. For
+#'   example, a coastline dataset for ocean data. **Note that this object must
+#'   have the same projection as the data used to generate the x and y columns
+#'   in `spde_obj`.**
+#' @param range_fraction The fraction of the spatial range that barrier
+#'   triangles have.
+#' @param proj_scaling If `spde_obj` was created with scaling of the coordinates
+#'   after the projection (e.g., dividing UTMs by 1000 so the spatial range is
+#'   on a reasonable scale) the x and y values in `spde_obj` are multiplied by
+#'   this scaling factor before applying the projection from `barrier_sf`.
+#' @param plot Logical.
+#'
+#' @return A list similar to [make_mesh()] but with `spde_barrier` and a
+#' couple other helper list elements added.
+#'
+#' If `plot = TRUE`, then a basic plot will be created as a side effect. Each
+#' grey dot represents the center of a "normal" mesh triangle. Each red cross
+#' represents the center of a "barrier" mesh triangle.
+#' @export
+#' @references
+#' Bakka, H., Vanhatalo, J., Illian, J., Simpson, D., and Rue, H. 2019.
+#' Non-stationary Gaussian models with physical barriers.
+#' <https://arxiv.org/abs/1608.03787>
+#'
+#' <https://sites.google.com/a/r-inla.org/www/barrier-model>
+#'
+#' <https://haakonbakkagit.github.io/btopic107.html>
+#' @examples
+#' if (require("sf", quietly = TRUE) &&
+#'   require("ggplot2", quietly = TRUE) &&
+#'   require("dplyr", quietly = TRUE) &&
+#'   require("rnaturalearth", quietly = TRUE)) {
+#'
+#' # For applied situations on finer scales, you may with to use scale = "large".
+#' # For that, first: remotes::install_github("ropensci/rnaturalearthhires")
+#' map_data <- rnaturalearth::ne_countries(
+#'   scale = "medium",
+#'   returnclass = "sf", country = "canada")
+#'
+#' # Crop the polygon for plotting and efficiency:
+#' sf::st_bbox(map_data)
+#' bc_coast <- suppressWarnings(suppressMessages(
+#'   sf::st_crop(map_data,
+#'     c(xmin = -134, ymin = 46, xmax = -120, ymax = 57))))
+#'
+#' crs_utm9 <- 3156 # Pick a projection, here UTM9
+#'
+#' sf::st_crs(bc_coast) <- 4326 # 'WGS84'
+#' bc_coast <- sf::st_transform(bc_coast, crs_utm9)
+#'
+#' # Project our survey data coordinates:
+#' survey <- pcod %>% dplyr::select(lon, lat, density) %>%
+#'   sf::st_as_sf(crs = 4326, coords = c("lon", "lat")) %>%
+#'   sf::st_transform(crs_utm9)
+#'
+#' # Plot our coast and survey data:
+#' ggplot(bc_coast) +
+#'   geom_sf() +
+#'   geom_sf(data = survey, size = 0.5)
+#'
+#' # Note that a barrier mesh won't do much here for this
+#' # example data set, but we nonetheless use it as an example.
+#'
+#' # Prepare for making the mesh
+#' # First, we will extract the coordinates:
+#' surv_utm_coords <- sf::st_coordinates(survey)
+#'
+#' # Then we will scale coordinates to km so the range parameter
+#' # is on a reasonable scale for estimation:
+#' pcod$X1000 <- surv_utm_coords[,1] / 1000
+#' pcod$Y1000 <- surv_utm_coords[,2] / 1000
+#'
+# # Construct our mesh:
+#' mesh <- make_mesh(pcod, xy_cols = c("X1000", "Y1000"),
+#'   n_knots = 200, type = "kmeans")
+#' plot(mesh)
+#'
+#' # Add on the barrier mesh component:
+#' bspde <- add_barrier_mesh(
+#'   mesh, bc_coast, range_fraction = 0.1,
+#'   proj_scaling = 1000, plot = TRUE
+#' )
+#'
+#' # In the above, the grey dots are the centre of triangles that are in the
+#' # ocean. The red crosses are centres of triangles that are over land. The
+#' # spatial range will be assumed to be 0.1 (`range_fraction`) over land
+#' # compared to over water.
+#'
+#' # We can make a more advanced plot if we want:
+#' mesh_df_water <- bspde$mesh_sf[bspde$normal_triangles, ]
+#' mesh_df_land <- bspde$mesh_sf[bspde$barrier_triangles, ]
+#'
+#' ggplot(bc_coast) +
+#'   geom_sf() +
+#'   geom_sf(data = mesh_df_water, size = 1, colour = "blue") +
+#'   geom_sf(data = mesh_df_land, size = 1, colour = "green")
+#'
+#' # Now, when we fit our model with the new mesh, it will automatically
+#' # include a barrier structure in the spatial correlation:
+#' fit <- sdmTMB(density ~ s(depth, k = 3), data = pcod, mesh = bspde,
+#'   family = tweedie(link = "log"))
+#' fit
+#' }
+
+add_barrier_mesh <- function(spde_obj, barrier_sf, range_fraction = 0.2,
+  proj_scaling = 1, plot = FALSE) {
+
+  if (!requireNamespace("sf", quietly = TRUE)) {
+    cli_abort("The sf package must be installed to use this function.")
+  }
+
+  assert_that(
+    is.numeric(range_fraction), range_fraction <= 1,
+    range_fraction > 0, length(range_fraction) == 1
+  )
+  assert_that(is.numeric(proj_scaling), length(proj_scaling) == 1)
+  assert_that(inherits(barrier_sf, "sf"))
+  assert_that(inherits(spde_obj, "sdmTMBmesh"))
+  assert_that(is.logical(plot))
+
+  mesh <- spde_obj$mesh
+  tl <- length(mesh$graph$tv[, 1]) # the number of triangles in the mesh
+  pos_tri <- matrix(0, tl, 2)
+  for (i in seq_len(tl)) {
+    temp <- mesh$loc[mesh$graph$tv[i, ], ]
+    pos_tri[i, ] <- colMeans(temp)[c(1, 2)]
+  }
+  mesh_sf <- as.data.frame(pos_tri)
+  mesh_sf$X <- mesh_sf$V1 * proj_scaling
+  mesh_sf$Y <- mesh_sf$V2 * proj_scaling
+  mesh_sf <- sf::st_as_sf(mesh_sf, crs = sf::st_crs(barrier_sf), coords = c("X", "Y"))
+
+  intersected <- sf::st_intersects(mesh_sf, barrier_sf)
+  water.triangles <- which(lengths(intersected) == 0)
+  land.triangles <- which(lengths(intersected) > 0)
+
+  if (plot) {
+    plot(pos_tri[water.triangles, ],
+      col = "grey40", asp = 1,
+      xlab = spde_obj$xy_cols[1], ylab = spde_obj$xy_cols[2]
+    )
+    points(pos_tri[land.triangles, ], col = "red", pch = 4)
+  }
+  # barrier_spde <- INLA::inla.barrier.fem(mesh, barrier.triangles = land.triangles)
+  barrier_spde <- INLAspacetime::mesh2fem.barrier(mesh, barrier.triangles = land.triangles)
+  spde_obj$spde_barrier <- barrier_spde
+  spde_obj$barrier_scaling <- c(1, range_fraction)
+  spde_obj$mesh_sf <- mesh_sf
+  spde_obj$barrier_triangles <- land.triangles
+  spde_obj$normal_triangles <- water.triangles
+  spde_obj
+}