View source: R/depth_from_unknown.R
depth_from_unknown | R Documentation |
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.
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 )
dat |
A dataframe which defines the FVCOM file date names, hours, and mesh (node) IDs for which layer depths should be calculated (see |
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 |
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 |
corrupt |
A vector of numbers, representing WeStCOMS date names, which define corrupt files (see |
read_fvcom |
A function which is used to load files (see |
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 |
extension |
A string which defines the extension of the FVCOM arrays (see |
cl |
(optional) A cluster objected created by the parallel package (see |
pass2varlist |
A list containing the names of exported objects. This may be required if |
depth_of_specified |
A logical input that defines whether or not to calculate the depths of all layers in |
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 |
assign_layer_method |
A character which defines the method by which layer IDs should be assigned to observed depths, if |
warning_threshold |
A number which defines the number of metres beyond the approximate depth of the seabed below mean sea level (i.e., |
verbose |
A logical input which defines whether or not to display messages to the console detailing function progress. |
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.
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).
Edward Lavender
#### 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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.