depth_from_unknown: Compute the depth of Sigma layers at specified location(s)...

In edwardlavender/fvcom.tbx: An R Toolbox for the Unstructured Grid Finite Volume Community Ocean Model (FVCOM)

Compute the depth of Sigma layers at specified location(s) and time(s)

Description

In FVCOM, three-dimensional hydrodynamic conditions are predicted for each node (or element) at 11 Sigma layers (or between these layers) for every hour of every day. The depth of each layer depends on the depth of the seabed below mean sea level (in a given location), the tidal elevation (in a given location and at a given time) and a constant which adjusts this depth for each layer, according to the equation depth_{t, l, n} = s_{l} \times (h_n + el_{n, t}) where s is a constant multiplier for layer l, h is the depth of node n below mean sea level and el is the tidal elevation predicted at that node at hour t. If values for h_n and el_{n, t} are already known, depth_from_known calculates the depths of the layers accordingly. Otherwise, depth_from_unknown first extracts necessary parameters for these variables using extract and then uses these to calculate the depths of all layers specified in siglev at specified times and locations. To implement (see depth_from_unknown), the user must supply a dataframe which contains the times (dates and hours) for which depths should be calculated, as well as a vector of constants (siglev) which are used to calculate the depths of corresponding layers. If requested, the function can use the computed depths of each layer to assign user-supplied depths Sigma layer IDs using nearest neighbour interpolation or fractional layer IDs using linear interpolation.

Usage

depth_from_unknown(
  dat,
  h,
  siglev,
  match_hour = data.frame(hour = 0:23, index = 1:24),
  match_mesh = NULL,
  corrupt = NULL,
  read_fvcom = function(con) R.matlab::readMat(con)$data,
  dir2load,
  extension = ".mat",
  cl = NULL,
  pass2varlist = NULL,
  depth_of_specified = FALSE,
  assign_layer = FALSE,
  assign_layer_method = "nearest",
  warning_threshold = 0,
  verbose = TRUE
)

Arguments

dat	A dataframe which defines the FVCOM file date names, hours, and mesh (node) IDs for which layer depths should be calculated (see extract). The dataframe may contain column named 'depth', containing (absolute) depths (m), if you want to assign layers IDs to depth observations (see Description). The dataframe may also contain a column named 'layer' if you want to compute the depths of a particular layer at each row in dat, rather than all layers in siglev (see below). The dataframe should be arranged by 'date_name'.
h	A dataframe which, for each mesh ID, defines the (absolute) depth of the seabed in that cell below mean sea level. Columns should be named 'ID' and 'h' respectively.
siglev	A dataframe which, for each Sigma layer, defines the (absolute value of the) siglev constant. Columns should be named 'layer' and 'siglev' respectively.
match_hour	A dataframe with two integer columns named 'hour' and 'index' which defines the index in FVCOM arrays (i.e. the row) which corresponds to each hour (see extract).
match_mesh	(optional) A dataframe with two columns named 'mesh' and 'index' which defines the index in FVCOM arrays (columns or sheets for 2- and 3-dimensional arrays respectively) which corresponds to each mesh cell (see extract).
corrupt	A vector of numbers, representing WeStCOMS date names, which define corrupt files (see extract).
read_fvcom	A function which is used to load files (see extract).
dir2load	A string which defines the directory from which to load the FVCOM arrays that contain the tidal elevation predictions required to calculate depth (see extract).
extension	A string which defines the extension of the FVCOM arrays (see extract).
cl	(optional) A cluster objected created by the parallel package (see extract).
pass2varlist	A list containing the names of exported objects. This may be required if cl is supplied. This is passed to the varlist argument of clusterExport. Exported objects must be located in the global environment (see extract).
depth_of_specified	A logical input that defines whether or not to calculate the depths of all layers in siglev for each row in dat (depth_of_specified = FALSE), or just the depths of specified layers, specified via dat$layer (depth_of_specified = TRUE).
assign_layer	A logical input that defines whether or not layer IDs should be assigned to observed depths, based on the computed depths of Sigma layers. This method requires a column named 'depth' in dat.
assign_layer_method	A character which defines the method by which layer IDs should be assigned to observed depths, if assign_layer = TRUE. Implemented options are "nearest" or "fractional". If assign_layer_method = "nearest", then each depth observation is assigned the layer ID of the nearest layer. If assign_layer_method = "fractional", a fractional layer ID is computed based on the depths of the two layers which surround the observed depth.
warning_threshold	A number which defines the number of metres beyond the approximate depth of the seabed below mean sea level (i.e., h$h, see above) after which the function will warn the user if there are any deeper depth observations. This may be due to inaccuracies in recorded depths or locations. However, note that in areas of complex bathymetry, there may well be areas in each mesh cell that are much deeper than the depth of the WeStCOMS layer at the seabed. This is why a warning_threshold is provided.
verbose	A logical input which defines whether or not to display messages to the console detailing function progress.

Details

The function currently only supports nearest neighbour interpolation, extracting tidal predictions for specified integer hours. To calculate layer depths for non integer hours, a custom approach is necessary.

Value

The function returns a dataframe, as inputted but with the following columns added: 'index_hour', the row in the FVCOM files which corresponds to the inputted hour; 'index_mesh', the column in the FVCOM files which corresponds to the inputted mesh ID; 'h', the depth (m) of the seabed in that cell below mean sea level; and 'el', the tidal elevation on the inputted hour and for the inputted mesh cell. If depth_of_specified = TRUE and dat$depth is supplied, then the dataframe also contains 'siglev', the siglev value for the specified layers and 'depth_layer', the depth of each inputted layer. Otherwise, the depth of every layer specified in siglev are provided in columns named 'l1',..., 'ln' where 'n' is the deepest layer in the inputted siglev dataframe. If assign_layer = TRUE, a column, 'layer_ID', is also returned; this is the ID of the layer assigned to each inputted depth observation. This column has one or more attributes: 'method' is "nearest" or "fractional" and, if assign_layer_method = "fractional", then the column also has the "computational_details" attribute which returns the information required to compute fractional layer IDs, such as the IDs of the surrounding layers, their depths and the weights applied to each of these to compute the fractional layer ID (see Examples).

Author(s)

Edward Lavender

Examples

#### Define dat
# Imagine we want to compute the depths of layers at the following timestamps
timestamp <- as.POSIXct(c("2016-03-01 00:00:00",
                          "2016-03-01 01:00:00",
                          "2016-03-02 00:00:00",
                          "2016-03-02 01:00:00"), tz = "UTC")
# For all of these nodes:
nodes <- as.integer(as.character(dat_mesh_around_nodes$ID))
dat <- expand.grid(timestamp = timestamp, mesh_ID = nodes)
# dat should be arranged by timestamp
dat <- dplyr::arrange(dat, timestamp)
dat$date_name <- date_name(dat$timestamp, define = "date_name")
dat$hour <- lubridate::hour(dat$timestamp)
head(dat)

#### Define other function inputs
h <- dat_h; colnames(h) <- c("ID", "h")
match_mesh <- data.frame(mesh = dat_nodexy$node_id, index = 1:length(dat_nodexy$node_id))
dir2load <- system.file("WeStCOMS_files/tidal_elevation",
                        package = "fvcom.tbx", mustWork = TRUE)
dir2load <- paste0(dir2load, "/")

#### Example (1) Compute layer depths for all layers using default options
dol1 <- depth_from_unknown(dat = dat,
                            h = h,
                            siglev = dat_siglev,
                            match_hour = data.frame(hour = 0:23, index = 1:24),
                            match_mesh = match_mesh,
                            dir2load = dir2load,
                            extension = ".mat",
                            cl = NULL,
                            pass2varlist = NULL
)
# Examine outputs
head(dol1, 3)
tail(dol1, 3)
# The depths of h and layer 10 (the seabed) should be similar:
dol1$h_from_h <- h$h[match(dat$mesh_ID, h$ID)]
head(dol1[, c("mesh_ID", "h", "h_from_h", "l10")])

#### Example (2) Incorporate parallel loading of tidal elevation files
dol2 <- depth_from_unknown(dat = dat,
                            h = h,
                            siglev = dat_siglev,
                            match_hour = data.frame(hour = 0:23, index = 1:24),
                            match_mesh = match_mesh,
                            dir2load = dir2load,
                            extension = ".mat",
                            cl = parallel::makeCluster(2L),
                            pass2varlist = NULL
)
head(dol2, 3)

#### Example (3) Assign observed depths to the nearest layer
# Imagine we have observations at the following depths.
# (Here, the use of match is simply to constrain hypothetical depths
# ... to be below the maximum depth of the seabed)
dat$depth <- h$h[match(dat$mesh_ID, h$ID)] - 5
dol3 <- depth_from_unknown(dat = dat,
                            h = h,
                            siglev = dat_siglev,
                            match_hour = data.frame(hour = 0:23, index = 1:24),
                            match_mesh = match_mesh,
                            dir2load = dir2load,
                            extension = ".mat",
                            assign_layer = TRUE,
                            assign_layer_method = "nearest"
)
head(dol3)
table(dol3$layer_ID)
# The layer_ID column contains the "method" attribute = "nearest"
attributes(dol3$layer_ID)

#### Example (4) Assign observed depths fractional layer IDs
# ... based on the depths of surrounding layers
dol4 <- depth_from_unknown(dat = dat,
                            h = h,
                            siglev = dat_siglev,
                            match_hour = data.frame(hour = 0:23, index = 1:24),
                            match_mesh = match_mesh,
                            dir2load = dir2load,
                            extension = ".mat",
                            assign_layer = TRUE,
                            assign_layer_method = "fractional"
)
head(dol4)
range(dol4$layer_ID)
# The layer_ID column contains the "method" attribute and the "computation_details" attribute
# ... The latter is a dataframe containing information needed to calculate fractional layer IDs.
str(attributes(dol4$layer_ID))
# For the computational_details attribute, columns are as follows:
# layer_nearest: the ID of the nearest layer
# layer_2ndnearest: the ID of the other surrounding layer
# depth_nearest: the depth of the nearest layer
# depth_2ndnearest: the depth of the other surrounding layer
# diff1: the absolute difference in depth between the observation and the first layer
# diff2: the absolute difference in depth between the observation and the second layer
# w1: the weight applied to the depth of the first layer
# w2 the weight applied to the depth of the second layer
# diff_layer: the absolute difference in depth between the two layers

#### Example (5) When layer IDs are being assigned to depths, potentially erroneous depth
# ... observations (i.e., observations deeper than the seabed) produce a warning:
## Not run: 
dat$depth[1:10] <- dat$depth[1:10] + 30
dol5 <- depth_from_unknown(dat = dat,
                            h = h,
                            siglev = dat_siglev,
                            match_hour = data.frame(hour = 0:23, index = 1:24),
                            match_mesh = match_mesh,
                            dir2load = dir2load,
                            extension = ".mat",
                            assign_layer = TRUE,
                            assign_layer_method = "fractional",
                            warning_threshold = 0
)
head(dol5)

## End(Not run)

#### Example 6) Compute the depths of specified layers with depth_of_specified = TRUE,
# ... rather than all layers in siglev.
set.seed(1)
dat$layer <- sample(1:10, nrow(dat), replace = TRUE)
dol6 <- depth_from_unknown(dat = dat,
                            h = h,
                            siglev = dat_siglev,
                            match_hour = data.frame(hour = 0:23, index = 1:24),
                            match_mesh = match_mesh,
                            dir2load = dir2load,
                            depth_of_specified = TRUE
)
utils::head(dol6, 3)

#### Example (7) Compute the depths of specified layers and assign layer_IDs to 'observed' depths
dat$depth <- 1
set.seed(1)
dat$layer <- sample(1:10, nrow(dat), replace = TRUE)
dol7 <- depth_from_unknown(dat = dat,
                            h = h,
                            siglev = dat_siglev,
                            match_hour = data.frame(hour = 0:23, index = 1:24),
                            match_mesh = match_mesh,
                            dir2load = dir2load,
                            depth_of_specified = TRUE,
                            assign_layer = TRUE,
                            assign_layer_method = "fractional"
)
utils::head(dol7, 3)
# dat$depth_layer is the depth of dat$layer
# dat$layer_ID is the layer ID assigned to the 'observed' depth (dat$depth)

edwardlavender/fvcom.tbx documentation built on Nov. 26, 2022, 10:28 p.m.