find_cells: Identify the mesh cells enclosing coordinates

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

Identify the mesh cells enclosing coordinates

Description

This function identifies the mesh cells which surround inputted coordinates. Coordinates can be inputted as lat and long vectors, a dataframe (dat) or a SpatialPoints-class object (xysp). The function can return a dataframe or vector of mesh IDs for each pair of coordinates or each unique pair of coordinates inputted.

Usage

find_cells(
  lat,
  long,
  dat = NULL,
  xysp = NULL,
  mesh,
  proj = sp::CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"),
  f = NULL,
  return = 1
)

Arguments

lat	A numeric vector of latitudes (in decimal degrees). Alternatively, dat or xysp can be provided (see below).
long	A numeric vector of longitudes (in decimal degrees). Alternatively, dat or xysp can be provided (see below).
dat	A dataframe with columns 'lat' and 'long' (in decimal degrees) at which mesh IDs will be identified. Alternatively, lat and long or xysp can be provided.
xysp	A SpatialPoints-class object which defines the positions at which to identify mesh IDs. Alternatively, lat and long or dat can be provided (see above).
mesh	A mesh, created by build_mesh, that surrounds nodes or elements.
proj	A projection string of class CRS-class. The default is World Geodetic System 84 (WGS84). This should be the same as for mesh and inputted coordinates.
f	A function to process mesh IDs after identification. For example, function(x) as.numeric(as.character(x)) can be useful.
return	A number (1, 2, 3 or 4) which specifies the format of the return output. return = 1 returns a dataframe with columns 'lat', 'long' and 'mesh_ID', comprising only unique coordinates. (Note that this may contain duplicate mesh IDs because unique coordinates may lie within the same mesh cell.) return = 2 returns only the mesh IDs corresponding to unique coordinates (e.g. to add to an existing dataframe). return = 3 returns a dataframe, as above, but including all inputted coordinate pairs. (This will be the same as outputted with return = 1 if their are no duplicate coordinates.) return = 4 returns only the mesh IDs corresponding to these coordinates.

Details

Ensure that the coordinate system of inputted coordinates, the mesh and the projection are identical (whether directly specified, e.g. if coordinates are specified via a SpatialPoints-class object or otherwise, e.g. if coordinates are supplied via lat and long vectors.

Value

The function returns a vector or a dataframe, depending on the input to return (see above).

Author(s)

Edward Lavender

Examples

#### Define some coordinates for which to identify the mesh cell ID
# In this hypothetical example, we'll choose coordinates for which
# ... we know the mesh ID in advance, so we can verify the function
# ... returns the correct outputs:
set.seed(1)
xy <- sp::coordinates(dat_mesh_around_nodes)
select <- sample(1:nrow(xy), 10)
xy <- xy[select, ]
selected_nodes <- as.numeric(rownames(xy))
lat <- as.numeric(xy[, 2])
long <- as.numeric(xy[, 1])

#### Set the mesh projection appropriately:
proj <- sp::CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
raster::crs(dat_mesh_around_nodes) <- proj

#### Example 1: Find cells from lat/long input:
cells <- find_cells(lat = lat,
                    long = long,
                    mesh = dat_mesh_around_nodes,
                    f = function(x) as.numeric(as.character(x)),
                    return = 1)
cells
identical(cells$mesh_ID, selected_nodes)

#### Example 2: Find cells from dat input:
cells <- find_cells(dat = data.frame(lat = lat, long = long),
                    mesh = dat_mesh_around_nodes,
                    f = function(x) as.numeric(as.character(x)),
                    return = 1)
identical(cells$mesh_ID, selected_nodes)

#### Example 3: Find cells from SpatialPoints input:
cells <- find_cells(xysp = sp::SpatialPoints(xy, proj4string = proj),
                    mesh = dat_mesh_around_nodes,
                    f = function(x) as.numeric(as.character(x)),
                    return = 1)
identical(cells$mesh_ID, selected_nodes)

#### Example 4: Find cells for unique coordinates or all coordinates in a dataframe:
# Imagine we have duplicate coordinates:
xy <- rbind(xy, xy)
rownames(xy) <- NULL
# Loop over each return type (1:4) and compare outputs:
cells_ls <-
  lapply(1:4, function(i){
    cells <- find_cells(xysp = sp::SpatialPoints(xy, proj4string = proj),
                        long = long,
                        mesh = dat_mesh_around_nodes,
                        f = function(x) as.numeric(as.character(x)),
                        return = i)
    return(cells)
  })
utils::str(cells_ls)

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