In BigelowLab/fvcom: fvcom

knitr::opts_chunk$set(echo = TRUE)

Convenient access to FVCOM datasets from R.

Background

FVCOM produces ocean circulation models on an irregular mesh. Here is the user manual.

As an example, NECOFS leverages this mesh to deliver models for the Gulf of Maine, Mass coastal waters and a small number of estuaries. NECOFS provides an THREDDS catalog of its OPeNDAP resources.

Requirements

• R 4.1+
• rlang
• ncdf4
• tibble
• dplyr
• sf
• stars
• locate # Not available on CRAN

Installation

remotes::install_github("BigelowLab/locate")`
remotes::install_github("BigelowLab/fvcom")

Gulf of Maine centric products

There are a number of models with a Gulf of Maine focus including forecasts and hindcasts. For each of these OPeNDAP resources we have provided a simple interface class.

• GOM GOMPhysics()

• Boston BostonPhysics()

• Scituate ScituatePhysics()

• Mass Bay MassBayPhysics()

• Saco Bay SacoBayPhysics()

• Hampton HamptonPhysics()

• Global GlobalPhysics()

Accessing FVCOM data

suppressPackageStartupMessages({
  library(sf) 
  library(fvcom) 
  library(ncdf4) 
  library(dplyr)
})

Open the FVCOM resource as you would any NetCDF file.

uri_base <- "http://www.smast.umassd.edu:8080/thredds/dodsC/models/fvcom/NECOFS/Archive/Seaplan_33_Hindcast_v1"
uri <- file.path(uri_base, "gom5_201812.nc")
x <- nc_open(uri)
x

Data structure

Nodes and Elements

The mesh is defined by non-intersecting triangular elements bounded by three nodes. Nodes are shown below as the solid dots, three nodes define the boundary of element. Scalar values, like temperature, salinty and height are defined at nodes. Vector values, such as velocity are defined at the element centroids. Within the NetCDF object you can make complete identification by examining which variables use the node dimensions versus those that use the nele dimension.

Node and Element locations

Use the functions fvcom_nodes and fvcom_elems to extract location information as either xy or lonlat pairs. Note we get both xy and lonlat here for nodes. The same can be had for elements.

dplyr::left_join(fvcom::fvcom_nodes(x, what = 'lonlat'),                   
                 fvcom::fvcom_nodes(x, what = 'xy'), by = "node") 

Timestamps

Retrieving the timestamps provides you with a choice for assumptions you want to make about the hourly model output. Timestamps stored internally do not land on each hour.

fvcom_time(x, internal = TRUE) |> head()

We can get a proxy for these times, but prettily settled on each hour. The choice is yours.

fvcom_time(x, internal = FALSE) |> head()

Variables

Variables (the oceanographic ones) can be extract by node or element. It is possible to select a subset, but the operation on the whole dataset is quick and it is just as easy to subset after you have the table in hand.

v <- get_node_var(x, var = 'zeta')
v

Mesh

Some computation is required (not a lot) to produce the mesh which is comprised of non-intersecting polygons (triangles). Meshes can be derived from the lists of nodes or the list of elements. We have chosen to use elements by default as they are simple to construct from the adjancency lists provided in the NetCDF resource.

The each element in the mesh is defined by the three nodes identified by index, “p1”, “p2” and “p3”. The geometry is defined by either “lonlat” coordinates or projected mercator coordinates, “xy”.

mesh <- get_mesh_geometry(x, what = 'lonlat') 
mesh

plot(sf::st_geometry(mesh), axes = TRUE)

Mesh with variables

We can assign variable values to the polygons by reusing the mesh table. These variables are associated with either node or elements. If node-referenced variables are requested the mean of three neighboring nodes (which define an element) is computed. Where element-referenced variables are requested no averaging is done - the values are simply assigned to the mesh element.

mesh <- get_mesh(x, vars = c("zeta", "u", "v"), mesh = mesh)
plot(mesh[c("u", "v")], lty = 'blank', main = c("u", "v"))

You can request variables at various dimensions such as times and sigma levels/layers - the default is the first of each dimension. While you can request one of these in ‘real’ values (such as POSIXct time), you can also provide a 1-based index into that dimension. The example below requests the same variables as above but at the 24th time interval. See the functions get_node_var and get_elem_var for details.

mesh <- get_mesh(x, vars = c("zeta", "u", "v"), mesh = mesh, time = 24)
plot(mesh[c("u", "v")], lty = 'blank', main = c("u", "v"))

Rasterize

The mesh can be interpolated on to a regular grid (“rasterize”).

template = default_template(mesh)
uv <- sapply(c("u", "v"), 
    function(f) {
        fvcom::rasterize(mesh[f], template = template) 
    }, simplify = FALSE)
par(mfrow = c(1,2))
plot(uv[['u']], key.pos = NULL, reset = FALSE)
plot(uv[['v']], key.pos = NULL, reset = FALSE)

ncdf4::nc_close(x)

BigelowLab/fvcom documentation built on Nov. 8, 2024, 2:24 p.m.