R/barrier.R

In tulpaMesh: Constrained Delaunay Triangulation Meshes for Spatial 'SPDE' Models

#' Identify Barrier Triangles
#'
#' Determines which mesh triangles fall inside barrier regions (e.g.,
#' coastlines, rivers, lakes). A triangle is marked as a barrier if its
#' centroid falls inside any barrier polygon.
#'
#' @param mesh A `tulpa_mesh` object.
#' @param barriers An sf/sfc object defining barrier regions, or a list
#'   of N x 2 coordinate matrices (each defining a closed polygon).
#'
#' @return A logical vector of length `n_triangles`. `TRUE` indicates
#'   the triangle is inside a barrier region.
#'
#' @export
barrier_triangles <- function(mesh, barriers) {
  if (!inherits(mesh, "tulpa_mesh")) {
    stop("mesh must be a tulpa_mesh object")
  }

  # Compute triangle centroids
  v <- mesh$vertices
  tri <- mesh$triangles
  centroids <- cbind(
    (v[tri[, 1], 1] + v[tri[, 2], 1] + v[tri[, 3], 1]) / 3,
    (v[tri[, 1], 2] + v[tri[, 2], 2] + v[tri[, 3], 2]) / 3
  )

  if (inherits(barriers, c("sf", "sfc", "sfg"))) {
    if (!requireNamespace("sf", quietly = TRUE)) {
      stop("sf package required for sf barrier objects")
    }
    barrier_geom <- sf::st_geometry(barriers)
    pts_sf <- sf::st_as_sf(
      data.frame(x = centroids[, 1], y = centroids[, 2]),
      coords = c("x", "y"),
      crs = sf::st_crs(barrier_geom)
    )
    inside <- as.logical(sf::st_intersects(pts_sf, sf::st_union(barrier_geom),
                                            sparse = FALSE))
  } else {
    # List of coordinate matrices — use ray casting
    if (is.matrix(barriers)) {
      barriers <- list(barriers)
    }
    inside <- rep(FALSE, nrow(centroids))
    for (poly in barriers) {
      inside <- inside | points_in_polygon(centroids[, 1], centroids[, 2], poly)
    }
  }

  inside
}