### 1.1 ==== Create a Mesh Based on Gradient Analysis ====
#' @title Creates a mesh based on gradient analysis
#'
#' @description A function to generate a mesh, placing points preferentially
#' where there are larger changes in the values of a set of environmental
#' covariates
#'
#' @param covars A \code{\link[terra]{SpatRaster}} object that contains the
#' covariate information
#' @param numCovarPoints An integer scalar that provides the number of extra
#' covariate-related triangulation points that will be produced and added to the
#' list of parameters appended to the \code{loc} argument before running the
#' \code{\link[fmesher]{fm_mesh_2d_inla}} mesh creation function.
#' @param ... Parameters to be passed to the
#' \code{\link[fmesher]{fm_mesh_2d_inla}} mesh creation function
#' @param funcApply A function to apply to the \code{covars} object before the
#' priority analysis is performed. A value of \code{NULL} means that no
#' function will be applied to the covariate information. The default option is
#' that the slope of the covariates will be calculated through the application
#' of the \code{\link[terra]{terrain}} function.
#' @param filename An optional character scalar giving the location to store the
#' calculated intensity raster object. If this value is NULL (the default) then
#' the raster is not stored in a file.
#' @param writeOptions A list of named options to be passed to the
#' \code{\link[terra]{writeRaster}} function when producing raster outputs.
#'
#' @return An \code{inla.mesh} object
#'
#' @author Joseph D. Chipperfield, \email{joechip90@@googlemail.com}
#' @seealso \code{\link[fmesher]{fm_mesh_2d_inla}},
#' \code{\link[terra]{SpatRaster}}, \code{\link[terra]{terrain}},
#' \code{\link[terra]{writeRaster}}
#' @export
<- (covars = , numCovarPoints = 100, , funcApply = (curRaster, extraOptions = writeOptions) {
(terra::terrain, ((x = curRaster, v = "slope", filename = (fileext = ".tif")), (extraOptions)))
}, filename = , writeOptions = ()) {
### 1.1.1 ---- Sanity check the inputs ----
ellipsisArgs <- ((()))
# Sanity check the filename input
inFileName <- filename
(!(inFileName)) {
inFileName <- ((inFileName), error = (err) {
("invalid argument given for the intensity raster file name: ", err)
})
((inFileName) <= 0) {
inFileName <-
} ((inFileName) > 1) {
("argument for the intensity raster file name has a length greater than one: only the first element will be used")
inFileName <- inFileName[1]
}
((inFileName) || inFileName == "") {
inFileName <-
}
}
# Function to import arguments that can be a matrix or a set of sf points
importPointArg <- (inArg, crsToUse) {
outOb <-
(!(inArg)) {
((inArg) || (inArg)) {
# If the input is a matrix then convert it to an sf object using the supplied
# coordinate reference system
outOb <- (inArg)
((outOb) < 2) {
("invalid entry for coordinate matrix: fewer than two columns")
} ((outOb) > 2) {
("coordinate matrix input has more than two columns: only the first two will be used")
outOb <- outOb, 1:2]
}
(outOb) <- ("x", "y")
outOb <- sf::st_as_sf(outOb, coords = ("x", "y"), crs = crsToUse)
} {
# Otherwise convert the object to a sf object
outOb <- (sf::st_as_sf(inArg), error = (err) {
("invalid entry for coordinate specification: ", err)
})
((sf::st_crs(outOb))) {
# Use the input CRS if none is provided
sf::st_crs(outOb) <- crsToUse
} (!(crsToUse) && sf::st_crs(outOb) != crsToUse) {
("coordinate matrix input does not share a cordinate reference system with the other parameters: performing conversion between reference systems")
outOb <- sf::st_transform(outOb, crsToUse)
}
}
}
outOb
}
# Retrieve the covariate raster object
covarsIn <- covars
crsToUse <-
(!(covars)) {
(!(covars, "SpatRaster") && !(covars, "SpatRasterDataset")) {
# Coerce to a raster object if it isn't already one
covarsIn <- (terra::rast(covarsIn), error = (err) {
("invalid entry for the input covariates: ", err)
})
}
crsToUse <- sf::st_crs(covarsIn)
}
# Retrieve the number of covariate-related triangulation points
numPointsIn <- ((numCovarPoints), error = (err) {
("invalid entry for the number of covariate-related triangulation points: ", err)
})
((numPointsIn) <= 0) {
numPointsIn <- 100
} ((numPointsIn) > 1) {
("more than one value entered for the number of covariate-related triangulation points: only the first value will be used")
numPointsIn <- numPointsIn[1]
}
((numPointsIn)) {
numPointsIn <- 100
} (numPointsIn < 0) {
("invalid entry for the number of covariate-related triangulation points: argument must be a positive integer")
}
# Retrieve the initial triangulation location points provided by the user (if any)
locIn <-
(!((ellipsisArgs)) && ("loc" == (ellipsisArgs))) {
locIn <- importPointArg(ellipsisArgs"loc"]], crsToUse)
crsToUse <- sf::st_crs(locIn)
}
# Retrieve the cutoff argument
cutoffIn <- (fmesher::fm_mesh_2d_inla)$cutoff
(!((ellipsisArgs)) && ("cutoff" == (ellipsisArgs))) {
cutoffIn <- ((ellipsisArgs"cutoff"]]), error = (err) {
("error processing the cutoff argument: ", err)
})
}
((cutoffIn) <= 0) {
cutoffIn <- (fmesher::fm_mesh_2d_inla)$cutoff
} ((cutoffIn) == 1) {
cutoffIn <- (cutoffIn, 2)
} ((cutoffIn) > 2) {
("more than two values entered for the cut-off distance: only the first two values will be used")
cutoffIn <- cutoffIn[1:2]
}
(((cutoffIn))) {
cutoffIn[is.na(cutoffIn)] <- (fmesher::fm_mesh_2d_inla)$cutoff
}
((cutoffIn <= 0)) {
("invalid entry for the cut-off distance: argument must be positive")
}
(!((ellipsisArgs)) && ("cutoff" == (ellipsisArgs))) {
# If the cutoff argument is present in the ellipsis arguments then set the
# cutoff argument to only be the second argument in the vector
ellipsisArgs"cutoff"]] <- cutoffIn[2]
}
### 1.1.2 ---- Apply the processing function to the covariates ----
intensityRast <-
(!(covarsIn)) {
# Sanity check the covariate processing function
inFuncApply <- funcApply
(!(funcApply)) {
inFuncAppply <- ((funcApply), error = (err) {
("invalid entry for the covaraite processing function: ", err)
})
}
# Function to apply the processing function to each layer of a SpatRaster object
applyFuncApply <- (inRast, retrCoords, inFuncApply) {
# Apply the function to each layer in the object
tempOut <- (X = inRast, FUN = (curRast, retrCoords, inFuncApply) {
# Apply the relevant processing function
tempRaster <- curRast
(!(inFuncApply)) {
tempRaster <- inFuncApply(curRast)
}
# Retrieve the elements of the processed raster at the requested coordinates
terra::extract(tempRaster, retrCoords), 1]
}, retrCoords = retrCoords, inFuncApply = inFuncApply)
# Sum the output matrix
(X = tempOut, MARGIN = 1, FUN = )
}
# Retrieve the coordinates of the cells of the finest-scale raster in the covariates
covarSum <-
((covarsIn, "SpatRasterDataset")) {
# If the input is a spatial raster dataset then go through each raster brick
# and apply the function to each layer within
retrCoords <- (X = covarsIn, FUN = terra::crds)
retrCoords <- retrCoords[which.max((X = retrCoords, FUN = ))]]
covarSum <- (X = (X = covarsIn, FUN = applyFuncApply, retrCoords = retrCoords, inFuncApply = inFuncApply), MARGIN = 1, FUN = )
} {
# Apply the function to each layer of the raster
retrCoords <- terra::crds(covarsIn)
covarSum <- applyFuncApply(covarsIn, retrCoords = retrCoords, inFuncApply = inFuncApply)
}
# Make spatial objects around the calculated intensity values
intensityFrame <- ((retrCoords), (intensity = covarSum))
intensityPoints <- sf::st_as_sf(intensityFrame, coords = ("x", "y"), crs = crsToUse)
intensityRast <- terra::rast(intensityFrame, = "xyz", crs = crsToUse$wkt)
(!(inFileName)) {
(terra::writeRaster, ((
x = intensityRast,
filename = inFileName
), (writeOptions)))
}
(!((ellipsisArgs)) && ("boundary" != (ellipsisArgs))) {
ellipsisArgs <- (ellipsisArgs, (
# If the user hasn't set a "boundary" argument then use the boundary of the non-zero
# elements of the covariate raster as the boundary
boundary = fmesher::fm_as_segm(sf::st_as_sf(terra::as.polygons(intensityRast > -)))
))
}
### 1.1.3 ---- Thin the points based on intensity ----
# First remove those points that have NA intensities
intensityPoints <- intensityPoints!(intensityPoints$intensity), ]
(!(locIn)) {
# If there is already a set of input locations then augment the intensity points frame to include those
sf::st_geometry(intensityPoints) <- (locIn, "sf_column")
nonGeomColnames <- (locIn)[names(locIn) != (locIn, "sf_column")]
((nonGeomColnames) > 0) {
intensityPoints <- (intensityPoints, ((, = (intensityPoints), = (nonGeomColnames), = (, nonGeomColnames))))
}
# If there is already a set of input locations then augment the data frame containing those entries
locIn <- (locIn, (intensity = (, (locIn))))
numPointsIn <- numPointsIn + (locIn)
distMat <- sf::st_distance(locIn, intensityPoints)
(cutoffIn) <- (distMat)
# units::units(cutoffIn) <- units::units(distMat)
# Remove those entries in the intensity points that are closer than the cutoff distance to the input points
intensityPoints <- intensityPoints[apply(X = distMat > cutoffIn[1], MARGIN = 2, FUN = ), ]
}
((intensityPoints) > 0) {
# Create a distance matrix between each of the putative points
distMat <- sf::st_distance(intensityPoints)
(cutoffIn) <- (distMat)
# units::units(cutoffIn) <- units::units(distMat)
((numPointsIn)) {
numPointsIn <-
}
}
(((locIn) || (locIn) < numPointsIn) && (intensityPoints) > 0) {
# Find the index of the cell value that has th highest intensity
highestIndex <- (intensityPoints$intensity)
# Add the point corresponding to the highest intensity and pass it to the locations
locIn <- (locIn, intensityPoints[highestIndex, ])
# Remove those points
isNotNear <- distMat[highestIndex, ] > cutoffIn[1]
intensityPoints <- intensityPoints[isNotNear, ]
distMat <- distMat[isNotNear, isNotNear]
}
}
### 1.1.4 ---- Process the generated points ----
# Final rejigging of outputs to conform with possible user-supplied coordinate reference systems
(((ellipsisArgs)) || ((ellipsisArgs) != "crs")) {
# There is no CRS provided by the user so use the one applied derived from the
# spatial input
ellipsisArgs <- (ellipsisArgs, (crs = crsToUse))
} {
# Otherwise overwrite the CRS with the user-specified one
crsToUse <- (sf::st_crs(ellipsisArgs"crs"]]), error = (err) {
("invalid entry for the coordinate reference system: ", err)
})
(!(intensityRast)) {
# Overwrite the CRS of the intensity raster if needed
intensityRast <- terra::project(intensityRast, crsToUse$wkt)
}
(!(locIn)) {
((sf::st_crs(locIn))) {
# Input points don't have a coordinate reference system - use the user-specified one
sf::st_crs(locIn) <- crsToUse
} {
# Input point already have a coordinate reference system - transform to the user-specified one
locIn <- sf::st_transform(locIn, crs = crsToUse)
}
}
(("boundary" == (ellipsisArgs))) {
((sf::st_crs(ellipsisArgs"boundary"]]))) {
# Utility function to overwrite CRS on existing segm objects
overwriteCRS <- (curSegm, crsToUse) {
tempOut <- fmesher::fm_as_segm(curSegm)
tempOut$crs <- fmesher::fm_crs(crsToUse)
tempOut
}
# Boundary doesn't have a coordinate reference system - use the user-specified one
((ellipsisArgs"boundary"]], "fm_segm_list")) {
ellipsisArgs"boundary"]] <- fmesher::fm_as_segm_list((X = (ellipsisArgs"boundary"]]), FUN = overwriteCRS, crsToUse = crsToUse))
} {
ellipsisArgs"boundary"]] <- overwriteCRS(ellipsisArgs"boundary"]], crsToUse = crsToUse)
}
} {
# Utility function to transform CRS on boundary objects
transformCRS <- (curSegm, crsToUse) {
tempOut <- fmesher::fm_as_segm(curSegm)
fmesher::fm_transform(tempOut, crs = fmesher::fm_crs(crsToUse))
}
# Boundary already has a coordinate reference system - transform the objects
((ellipsisArgs"boundary"]], "fm_segm_list")) {
ellipsisArgs"boundary"]] <- fmesher::fm_as_segm_list((X = (ellipsisArgs"boundary"]]), FUN = transformCRS, crsToUse = crsToUse))
} {
ellipsisArgs"boundary"]] <- transformCRS(ellipsisArgs"boundary"]], crsToUse = crsToUse)
}
}
}
}
# Retrieve the boundary if it has been set
outBoundary <-
(!((ellipsisArgs)) && ((ellipsisArgs) == "boundary")) {
outBoundary <- (fmesher::fm_as_sfc(ellipsisArgs"boundary"]]))
}
# Add the generated locations to the ellipsis arguments
ellipsisArgs$loc <- sf::st_as_sfc(locIn)
# Run the mesh making algorithm with the updated parameters
meshOut <- (fmesher::fm_mesh_2d_inla, ellipsisArgs)
meshClass <- (meshOut)
# Add the intensity and initial point placement information and define
# the class of the object to be of type "fm_mesh_2d_intensity"
(!(intensityRast)) {
intensityRast <- terra::wrap(intensityRast)
}
meshOut <- (meshOut, (intensity = intensityRast, initialLocs = locIn, boundary = outBoundary))
(meshOut) <- ("fm_mesh_2d_intensity", meshClass)
meshOut
}
#' @title Plot a Mesh Object
#'
#' @description Function to plot a \code{fm_mesh_2d_intensity} object. This
#' object type is what is returned from the \code{\link{priorityMesh}} function.
#'
#' @param object A \code{fm_mesh_2d_intensity} object created from the
#' \code{\link{priorityMesh}} function.
#' @param ... A series of named parameters prefixed with \code{"intensity."},
#' \code{"mesh."}, \code{"initialLocs."}, \code{"boundary."} to pass parameters
#' to the \code{\link[ggplot2]{geom_sf}} or
#' \code{\link[tidyterra]{geom_spatraster}} geometry functions that determine
#' the plotting behaviour of the intensity surface, mesh, initial triangulation
#' nodes, and boundary polygon respectively. Is is also possible to specify
#' other arguments that control the look of elements related to the aesthetic
#' properties of the plot, see 'Details' below.
#'
#' @details
#' In addition to the arguments outlined above, the user may also set the
#' following optional arguments that control the look of elements related to the
#' aesthetic properties of the plot:
#' \describe{
#' \item{\code{show.intensity}}{A logical that if \code{TRUE} (the default)
#' displays the intensity surface.}
#' \item{\code{show.initialLocs}}{A logical that if \code{TRUE} (the default)
#' displays the initial triangulation points.}
#' \item{\code{show.boundary}}{A logical that if \code{TRUE} (the default)
#' displays the boundary.}
#' \item{\code{low.intensity.col}}{A character scalar containing the colour to
#' use for the cells with the lowest intensity values in the intensity surface.
#' }
#' \item{\code{high.intensity.col}}{A character scalar containing the colour to
#' use for the cells with the highest intensity values in the intensity
#' surface.}
#' \item{\code{na.intensity.col}}{A character scalar containing the colour to
#' use for the cells with an NA value in the intensity surface.}
#' \item{\code{provided.initialLocs.col}}{A character scalar containing the
#' colour to use to plot the points that were provided by the user as initial
#' triangulation nodes.}
#' \item{\code{generated.initialLocs.col}}{A character scalar containing the
#' colour to use to plot the points that were generated by the
#' \code{\link{priorityMesh}} function.}
#' }
#'
#' @return A \code{ggplot} object
#'
#' @author Joseph D. Chipperfield, \email{joechip90@@googlemail.com}
#' @seealso \code{\link[fmesher]{fm_mesh_2d_inla}} \code{\link{priorityMesh}}
#' \code{\link[ggplot2]{geom_sf}} \code{\link[tidyterra]{geom_spatraster}}
#' @export
<- (object, ) {
# Default values for the visualisation of the aesthetic components
low.intensity.col <- retrieveOtherArgs("low.intensity.col", , defaultVal = ::(255, 255, 240, maxColorValue = 255))
high.intensity.col <- retrieveOtherArgs("high.intensity.col", , defaultVal = ::(110, 139, 61, maxColorValue = 255))
na.intensity.col <- retrieveOtherArgs("na.intensity.col", , defaultVal = )
provided.initialLocs.col <- retrieveOtherArgs("provided.initialLocs.col", , defaultVal = ::(125, 158, 192, maxColorValue = 255))
generated.initialLocs.col <- retrieveOtherArgs("generated.initialLocs.col", , defaultVal = ::(198, 113, 113, maxColorValue = 255))
show.intensity <- retrieveOtherArgs("show.intensity", , defaultVal = )
show.initialLocs <- retrieveOtherArgs("show.initialLocs", , defaultVal = )
show.boundary <- retrieveOtherArgs("show.boundary", , defaultVal = )
# Create a set of default elements for the geometries
defaultArgs <- (
intensity.data = terra::unwrap(object$intensity),
intensity.show.legend = ,
mesh.fill = ,
mesh.colour = ::(170, 170, 170, maxColorValue = 255),
initialLocs.mapping = ggplot2::aes(colour = (((("intensity")), "Provided", "Generated"), = ("Provided", "Generated"))),
initialLocs.data = object$initialLocs,
initialLocs.show.legend = ,
boundary.fill = ,
boundary.colour = ::(218, 112, 214, maxColorValue = 255),
boundary.linewidth = 1,
boundary.data = object$boundary
)
# Create a temporary sf object from the mesh polygons
tempMesh <- fmesher::fm_as_sfc(object)
# Initialise an output plot object
outPlot <- ggplot2::ggplot(tempMesh)
(!(object$intensity) && (show.intensity)) {
# If the intensity surface exists then plot that object
outPlot <- outPlot + (tidyterra::geom_spatraster, retrievePrefixArgs("intensity", , defaultArgs = defaultArgs)) +
ggplot2::scale_fill_gradient(low = low.intensity.col, high = high.intensity.col, na.value = na.intensity.col)
}
# Add the mesh polygons
outPlot <- outPlot + (ggplot2::geom_sf, retrievePrefixArgs("mesh", , defaultArgs = defaultArgs))
(!(object$initialLocs) && (show.initialLocs)) {
# If initial locations have been provided then plot those points
outPlot <- outPlot + (ggplot2::geom_sf, retrievePrefixArgs("initialLocs", , defaultArgs = defaultArgs)) +
ggplot2::scale_color_manual(values = (provided.initialLocs.col, generated.initialLocs.col))
}
(!(object$boundary) && (show.boundary)) {
# If the boundary exists then add that line to the plot
outPlot <- outPlot + (ggplot2::geom_sf, retrievePrefixArgs("boundary", , defaultArgs = defaultArgs))
}
# Add the information to colour the aesthetics in the plot
outPlot <- outPlot + ggplot2::theme_classic()
outPlot
}
