Read the FVCOM mesh itself, and draw it

accessFVCOM() returns values at points - nodes for scalars, element
centroids for velocities. That is what matchData() needs and it is not the
grid: an FVCOM mesh is a triangulation, and which nodes form each triangle
lives in a connectivity array a point fetch never reads. Plotting what came
back therefore showed a scatter of dots where the model has cells, with no
way to draw a boundary, shade a cell, or see where the mesh refines - which
is most of why anyone looks at an unstructured model.

fvcom_mesh() reads that array and returns the triangles as POLYGONs with a
DEPTH per cell, carrying no time dimension because it is the grid rather
than anything measured on it. what = "nodes" and "elements" return the two
point sets alone.

plot_mesh() draws it, bare or shaded, and does the join when given a fetch.
That join is the part worth automating: a fetch is subset to the bounding
box, so its row order says nothing about the mesh's own element numbering,
and joining by position would mislabel every cell without complaining.

Two things only came out by trying it. Node-centred values sit at triangle
*corners*, shared between the cells meeting there - and st_contains
excludes boundaries, so a node join returned nothing whatsoever. It tests
intersection instead, which catches corners and centroids alike. And
plot.sf wants a palette function it can call with the number of breaks it
chose, not a vector of colours; a vector fails with "must have one more
break than colour".

A triangle is kept when any vertex falls inside the box, so the mesh covers
what was asked for rather than stopping short, and the edge is ragged by
design - which is also why an element join leaves the outermost triangles
NA rather than borrowing a neighbour's value. GOM7 remeshes between months,
so date says which month's mesh to read.

Also describes the five access functions up front in the README, which had
been listed in a table but never actually explained, and adds the mesh
plotting to it.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

Author: chross22

NAMESPACE

@@ -31,6 +31,7 @@ export(forecast_variables)
 export(fvcom_archive)
 export(fvcom_archives)
 export(fvcom_dictionary)
+export(fvcom_mesh)
 export(fvcom_variables)
 export(grid_resolution)
 export(hycom_archives)
@@ -42,6 +43,7 @@ export(matchData)
 export(plot_coverage)
 export(plot_env)
 export(plot_matched)
+export(plot_mesh)
 export(plot_series)
 export(product_url)
 export(refresh_climate_index)

NEWS.md

@@ -2,6 +2,31 @@
 
 ## New features
 
+* **The FVCOM mesh itself, through `fvcom_mesh()`.** `accessFVCOM()` returns
+  values at points — nodes for scalars, element centroids for velocities — which
+  is what matching needs but is not the grid. Drawing it shows dots where the
+  model has triangles.
+
+  `fvcom_mesh()` reads the connectivity array a point fetch never touches and
+  returns the triangles as `POLYGON`s, with `DEPTH` per cell. It carries no time
+  dimension: it is the grid. `what = "nodes"` and `"elements"` return the two
+  point sets on their own.
+
+  **`plot_mesh()`** draws it — bare, or shaded by a covariate. Pass a fetch as
+  `values` and the spatial join is done for you, which matters more than it
+  sounds: a fetch is subset to the bounding box, so its row order says nothing
+  about the mesh's element numbering, and joining by position would silently
+  mislabel every cell.
+
+  Node values sit at triangle *corners* and element values at centroids, so the
+  join tests intersection rather than containment — a corner lies on the
+  boundary between cells, and nothing contains it. Joining by containment
+  returned no node values at all.
+
+  A triangle is kept when any vertex is inside the box, so the mesh covers what
+  was asked for rather than stopping short of it, and the edge is ragged by
+  design. `GOM7` remeshes between months, so `date` says which month's mesh.
+
 * **EML metadata, through `write_eml()`.** Writes an Ecological Metadata
   Language document for a matched table — the standard EDI, LTER and DataONE
   expect alongside a deposited dataset.

R/fvcom.R

@@ -911,3 +911,174 @@ fvcom_dictionary <- function(mesh = c("all", "node", "element")) {
   class(dictionary) <- c("datamatch_dictionary", "data.frame")
   dictionary
 }
+
+#' Read an FVCOM mesh on its own, with no data on it
+#'
+#' [accessFVCOM()] returns values at points — nodes for scalars, element
+#' centroids for velocities. That is what [matchData()] needs, but it is not the
+#' mesh: it is a scatter of the mesh's points, and drawing it shows dots where
+#' the grid is triangles. This returns the triangles themselves, so the grid can
+#' actually be plotted.
+#'
+#' @section Why the points are not enough:
+#' An FVCOM mesh is a triangulation, and which nodes form each triangle lives in
+#' a connectivity array (`nv`) that a point fetch never reads. Without it there
+#' is no way to draw a cell boundary, shade a cell by its value, or show how the
+#' resolution changes across a shelf — all of which are the usual reasons for
+#' looking at an unstructured model in the first place.
+#'
+#' The result carries no time dimension and no covariates. It is the grid.
+#'
+#' @section What comes back:
+#' \itemize{
+#'   \item `"polygons"` (the default) — one `POLYGON` per element, with the
+#'     element index and `DEPTH`, the mean of its three nodes' bathymetry. This
+#'     is what to plot.
+#'   \item `"nodes"` — the mesh nodes as `POINT`s, with `DEPTH`. Scalars such as
+#'     `SST` and `BOTS` live here.
+#'   \item `"elements"` — the element centroids as `POINT`s. Velocities and
+#'     stresses live here.
+#' }
+#'
+#' @section Bounding box:
+#' A triangle is kept when **any** of its three vertices falls inside the box, so
+#' the returned mesh covers the box rather than stopping short of it. The edge is
+#' therefore ragged, and a triangle may extend a little beyond what was asked
+#' for. Clipping to the box exactly would cut triangles into shapes the model
+#' does not have.
+#'
+#' @section Meshes that move:
+#' `GOM3` has one mesh for its whole record. `GOM7` does not — it is operational
+#' output and remeshes between files, carrying 198,594 nodes in January 2025 and
+#' 207,081 from March. For an archive like that the mesh belongs to a particular
+#' month, so `date` says which; it defaults to the start of the record.
+#'
+#' @param archive which archive's mesh: a name from [fvcom_archives()], or a spec
+#'   from [fvcom_archive()]
+#' @param bounding_box <list> optional named list with `xmin`, `xmax`, `ymin`,
+#'   `ymax`, or an `sf`/`sfc` object. `NULL` returns the whole mesh, which for
+#'   GOM3 is 90,415 triangles.
+#' @param what <char> `"polygons"`, `"nodes"`, or `"elements"`
+#' @param date the month whose mesh to read, for an archive that remeshes.
+#'   Ignored for a single-mesh archive.
+#' @return an `sf` object: `POLYGON` for `"polygons"`, `POINT` otherwise, in
+#'   EPSG:4326
+#' @examples
+#' \dontrun{
+#' bb <- list(xmin = -70, xmax = -66, ymin = 41, ymax = 44)
+#'
+#' mesh <- fvcom_mesh(bounding_box = bb)
+#' plot(sf::st_geometry(mesh))                 # the grid itself
+#' plot(mesh["DEPTH"], border = NA)            # shaded by bathymetry
+#'
+#' # Shade the triangles by a fetched value. Join spatially rather than by
+#' # position: accessFVCOM() returns only the points inside the box, so its row
+#' # order does not correspond to the mesh's element numbering.
+#' currents <- accessFVCOM(vars = "UBAR", years = 2010, months = 6,
+#'                         bounding_box = bb)
+#' shaded <- sf::st_join(mesh, currents["UBAR"], join = sf::st_contains)
+#' plot(shaded["UBAR"], border = NA)
+#' }
+#' @seealso [accessFVCOM()] for values on the mesh, [fvcom_archives()]
+#' @export
+fvcom_mesh <- function(archive = "GOM3", bounding_box = NULL,
+                       what = c("polygons", "nodes", "elements"), date = NULL) {
+  what <- match.arg(what)
+
+  if (is.list(archive)) {
+    spec <- archive
+    if (is.null(spec$url)) {
+      stop("An archive given as a list must carry a `url`. Build one with ",
+           "fvcom_archive().", call. = FALSE)
+    }
+  } else {
+    archives <- fvcom_archives()
+    if (!archive %in% names(archives)) {
+      stop("Unknown archive '", archive, "'. Built in: ",
+           paste(names(archives), collapse = ", "),
+           "\nTo read any other, describe it with fvcom_archive(url).",
+           call. = FALSE)
+    }
+    spec <- archives[[archive]]
+  }
+
+  # A per-month archive's url is a template, and its mesh can differ between
+  # months, so some month has to be named. See the Meshes that move section.
+  if (identical(spec$layout, "per_month")) {
+    when <- if (is.null(date)) spec$start else parse_dates(date)[1]
+    spec$url <- sprintf(spec$url, as.integer(format(when, "%Y")),
+                        as.integer(format(when, "%Y")),
+                        as.integer(format(when, "%m")))
+  }
+
+  handle <- fvcom_open(spec)
+  on.exit(ncdf4::nc_close(handle), add = TRUE)
+
+  if (identical(what, "elements")) {
+    coords <- fvcom_coordinates(handle, "element")
+    coords$element <- seq_len(nrow(coords))
+    keep <- if (is.null(bounding_box)) {
+      seq_len(nrow(coords))
+    } else {
+      fvcom_in_box(coords, bounding_box)
+    }
+    return(sf::st_as_sf(coords[keep, , drop = FALSE], coords = c("x", "y"),
+                        crs = sf::st_crs(4326)))
+  }
+
+  nodes <- fvcom_coordinates(handle, "node")
+  nodes$node <- seq_len(nrow(nodes))
+  depth <- if (!is.null(handle$var[["h"]])) {
+    as.numeric(ncdf4::ncvar_get(handle, "h"))
+  } else {
+    rep(NA_real_, nrow(nodes))
+  }
+  nodes$DEPTH <- depth
+
+  if (identical(what, "nodes")) {
+    keep <- if (is.null(bounding_box)) {
+      seq_len(nrow(nodes))
+    } else {
+      fvcom_in_box(nodes, bounding_box)
+    }
+    return(sf::st_as_sf(nodes[keep, , drop = FALSE], coords = c("x", "y"),
+                        crs = sf::st_crs(4326)))
+  }
+
+  if (is.null(handle$var[["nv"]])) {
+    stop("This archive has no `nv` connectivity array, so its triangles cannot ",
+         "be built.\nWithout it the mesh is only a set of points; use ",
+         "what = \"nodes\".", call. = FALSE)
+  }
+  # [nele, 3], one-based node indices - FVCOM is written in Fortran.
+  nv <- ncdf4::ncvar_get(handle, "nv")
+
+  # A triangle is kept when any vertex is inside, so the mesh covers the box
+  # rather than stopping short of it.
+  keep <- if (is.null(bounding_box)) {
+    seq_len(nrow(nv))
+  } else {
+    inside <- rep(FALSE, nrow(nodes))
+    inside[fvcom_in_box(nodes, bounding_box)] <- TRUE
+    which(inside[nv[, 1]] | inside[nv[, 2]] | inside[nv[, 3]])
+  }
+  if (length(keep) == 0) {
+    stop("No mesh triangles fall inside the bounding box.", call. = FALSE)
+  }
+
+  x <- nodes$x
+  y <- nodes$y
+  triangles <- lapply(keep, function(i) {
+    vertices <- nv[i, ]
+    # Closed ring: the first vertex repeated at the end, as sf requires.
+    ring <- cbind(x[c(vertices, vertices[1])], y[c(vertices, vertices[1])])
+    sf::st_polygon(list(ring))
+  })
+
+  out <- sf::st_sf(
+    element = keep,
+    DEPTH = rowMeans(matrix(depth[nv[keep, ]], ncol = 3)),
+    geometry = sf::st_sfc(triangles, crs = sf::st_crs(4326)))
+  attr(out, "datamatch_mesh") <- spec$mesh %||% "unstructured"
+  out
+}

R/plot.R

@@ -370,3 +370,125 @@ pretty_steps <- function(n) {
 numeric_for_plot <- function(x) {
   if (is.numeric(x)) x else as.integer(as.factor(x))
 }
+
+#' Map covariates on an unstructured mesh, or anywhere else
+#'
+#' [plot_env()] rasterises, which is right for a regular grid and wrong for an
+#' FVCOM mesh: a triangulation has no rows and columns to rasterise onto, so the
+#' result is either blocky or interpolated over cells the model does not have.
+#' This draws the geometry it is given.
+#'
+#' @section What it is for:
+#' Two jobs. With a mesh from [fvcom_mesh()] and nothing else, it draws the grid
+#' itself — which is how you see where a model resolves a shelf finely and where
+#' it does not. With `var` naming a column, it shades each cell by that value.
+#'
+#' @section Shading a mesh with fetched values:
+#' [accessFVCOM()] returns points and [fvcom_mesh()] returns triangles, so the
+#' two are joined spatially rather than by position — a fetch is subset to the
+#' bounding box, so its row order says nothing about the mesh's own numbering.
+#' Pass `values` and that join is done for you:
+#'
+#' ```
+#' mesh <- fvcom_mesh(bounding_box = bb)
+#' sst  <- accessFVCOM(vars = "SST", years = 2010, months = 6, bounding_box = bb)
+#' plot_mesh(mesh, "SST", values = sst)
+#' ```
+#'
+#' Element-centred values (`UO`, `UBAR`, `TAUX`) sit at the centroid, one per
+#' triangle. Node-centred ones (`SST`, `BOTS`) sit at the *corners*, shared
+#' between the triangles meeting there, so each triangle takes the mean of its
+#' three — which is why the join tests intersection rather than containment: a
+#' corner is on the boundary, and nothing contains it.
+#'
+#' @param x an `sf` object with polygon or point geometry — typically from
+#'   [fvcom_mesh()], but any `sf` object works
+#' @param var which column to shade by; `NULL` draws the bare geometry
+#' @param values optional `sf` object holding the values, from [accessFVCOM()]
+#'   or any access function. Joined to `x` spatially.
+#' @param time which time step of `values` to use, as an index or a named vector
+#'   of `YEAR`/`MONTH`/`DAY`. Only needed when `values` carries more than one.
+#' @param palette a [grDevices::hcl.colors()] palette name
+#' @param border colour for cell edges. `NA` hides them, which is what you want
+#'   on a fine mesh where the edges would otherwise be all you see.
+#' @param main plot title
+#' @param ... passed to [plot()]
+#' @return the `sf` object that was drawn, invisibly, with the joined column if
+#'   `values` was given
+#' @examples
+#' \dontrun{
+#' bb <- list(xmin = -70, xmax = -66, ymin = 41, ymax = 44)
+#' mesh <- fvcom_mesh(bounding_box = bb)
+#'
+#' plot_mesh(mesh)                      # the grid itself
+#' plot_mesh(mesh, "DEPTH")             # shaded by bathymetry
+#'
+#' sst <- accessFVCOM(vars = "SST", years = 2010, months = 6, bounding_box = bb)
+#' plot_mesh(mesh, "SST", values = sst)
+#' }
+#' @seealso [fvcom_mesh()] for the grid, [plot_env()] for regular grids
+#' @export
+plot_mesh <- function(x, var = NULL, values = NULL, time = 1,
+                      palette = "viridis", border = NA, main = NULL, ...) {
+  if (!inherits(x, "sf")) {
+    stop("`x` must be an sf object, such as fvcom_mesh() returns.",
+         call. = FALSE)
+  }
+
+  label <- NULL
+  if (!is.null(values)) {
+    if (is.null(var)) {
+      stop("`values` was given but `var` was not, so there is nothing to take ",
+           "from it.", call. = FALSE)
+    }
+    if (!var %in% names(values)) {
+      stop("`values` has no column '", var, "'.\nIt has: ",
+           paste(covariate_columns(values), collapse = ", "), call. = FALSE)
+    }
+
+    # One time step: a mesh cell can hold one value, so which step must be
+    # settled before the join rather than averaged across silently.
+    steps <- time_steps(values)
+    if (nrow(steps$table) > 1) {
+      step <- select_time_step(values, time)
+      values <- values[step$rows, ]
+      label <- step$label
+    }
+
+    # st_join returns one row per (cell, point) pair, so a cell containing
+    # several points appears several times. The row is marked first so those can
+    # be collapsed back onto the cell they came from.
+    marked <- x
+    marked$.cell <- seq_len(nrow(x))
+    # st_intersects rather than st_contains: node-centred values sit at triangle
+    # *corners*, on the boundary between cells, and st_contains excludes
+    # boundaries - so a node join returned nothing at all. Intersects catches
+    # both, since an element centroid strictly inside a triangle intersects it
+    # too.
+    joined <- sf::st_drop_geometry(
+      sf::st_join(marked, values[var], join = sf::st_intersects))
+
+    # A node-centred variable puts several points in one triangle; an
+    # element-centred one puts exactly its centroid in it. Averaging handles
+    # both, and a triangle containing no point stays NA.
+    averaged <- tapply(joined[[var]], joined$.cell,
+                       function(v) if (all(is.na(v))) NA_real_ else mean(v, na.rm = TRUE))
+    x[[var]] <- as.numeric(averaged[as.character(seq_len(nrow(x)))])
+  }
+
+  if (is.null(var)) {
+    plot(sf::st_geometry(x), border = if (is.na(border)) "grey40" else border,
+         main = main %||% "FVCOM mesh", ...)
+    return(invisible(x))
+  }
+
+  check_columns(x, var)
+  # plot.sf wants a palette *function* it can call with the number of breaks it
+  # chose, not a fixed vector of colours - passing a vector gives "must have one
+  # more break than colour".
+  plot(x[var], pal = function(n) grDevices::hcl.colors(n, palette),
+       border = border,
+       main = main %||% paste0(var, if (!is.null(label)) paste0("  ", label)),
+       ...)
+  invisible(x)
+}

R/provenance.R

@@ -19,7 +19,7 @@
 #' most: which model the value came from at all.
 #'
 #' @param x an object from [accessEnvDat()], [accessFVCOM()], [accessHYCOM()],
-#'   [accessCCMP()]
+#'   [accessCCMP()] or [accessERDDAP()]
 #' @return <char> the source tag, or `NA` if the object carries none — which is
 #'   the case for anything built by hand or produced before this was recorded
 #' @examples