| Title: | Modeling Species Distributions in Three Dimensions |
|---|---|
| Description: | Facilitates modeling species' ecological niches and geographic distributions based on occurrences and environments that have a vertical as well as horizontal component, and projecting models into three-dimensional geographic space. Working in three dimensions is useful in an aquatic context when the organisms one wishes to model can be found across a wide range of depths in the water column. The package also contains functions to automatically generate marine training model training regions using machine learning, and interpolate and smooth patchily sampled environmental rasters using thin plate splines. Davis Rabosky AR, Cox CL, Rabosky DL, Title PO, Holmes IA, Feldman A, McGuire JA (2016) <doi:10.1038/ncomms11484>. Nychka D, Furrer R, Paige J, Sain S (2021) <doi:10.5065/D6W957CT>. Pateiro-Lopez B, Rodriguez-Casal A (2022) <https://CRAN.R-project.org/package=alphahull>. |
| Authors: | Hannah L. Owens [aut, cre, cph] (ORCID: <https://orcid.org/0000-0003-0071-1745>), Emmaline Sheahan [aut] (ORCID: <https://orcid.org/0000-0003-3358-9758>), Carsten Rahbek [aut] (ORCID: <https://orcid.org/0000-0003-4585-0300>) |
| Maintainer: | Hannah L. Owens <[email protected]> |
| License: | GPL-3 |
| Version: | 0.2.4.9000 |
| Built: | 2026-05-24 15:06:04 UTC |
| Source: | https://github.com/hannahlowens/volumodel |
Samples deepest depth values from a
SpatVector data frame and generates a SpatRaster.
bottomRaster(rawPointData)bottomRaster(rawPointData)
rawPointData |
A |
rawPointData is a SpatVector object that
contains measurements of a single environmental variable (e.g.
salinity, temperature, etc.) with x, y, and z coordinates. The
measurements in the data.frame should be organized so that each
column is a depth slice, increasing in depth from left to right. The
function was designed around the oceanographic data shapefiles supplied
by the World Ocean Atlas
(https://www.ncei.noaa.gov/access/world-ocean-atlas-2018/).
The function selects the "deepest" (rightmost) measurement at each
x, y coordinate pair that contains data. These measurements are then
rasterized at the resolution and extent of the x,y coordinates, under
the assumption that x and y intervals are equal and represent the center
of a cell.
A SpatRaster designed to approximate sea bottom
measurements for modeling species' distributions and/or niches.
library(terra) # Create point grid coords <- data.frame(x = rep(seq(1:5), times = 5), y = unlist(lapply(1:5, FUN = function(x) { rep(x, times = 5)}))) # Create data and add NAs to simulate uneven bottom depths dd <- data.frame(SURFACE = 1:25, d5M = 6:30, d10M = 11:35, d25M = 16:40) dd$d25M[c(1:5, 18:25)] <- NA dd$d10M[c(3:5, 21:23)] <- NA dd$d5M[c(4, 22)] <- NA dd[,c("x","y")] <- coords # Create SpatialPointsDataFrame sp <- vect(dd, geom = c("x", "y")) # Here's the function result <- bottomRaster(rawPointData = sp) plot(result)library(terra) # Create point grid coords <- data.frame(x = rep(seq(1:5), times = 5), y = unlist(lapply(1:5, FUN = function(x) { rep(x, times = 5)}))) # Create data and add NAs to simulate uneven bottom depths dd <- data.frame(SURFACE = 1:25, d5M = 6:30, d10M = 11:35, d25M = 16:40) dd$d25M[c(1:5, 18:25)] <- NA dd$d10M[c(3:5, 21:23)] <- NA dd$d5M[c(4, 22)] <- NA dd[,c("x","y")] <- coords # Create SpatialPointsDataFrame sp <- vect(dd, geom = c("x", "y")) # Here's the function result <- bottomRaster(rawPointData = sp) plot(result)
Creates a SpatRaster template from a
SpatVector point object in which the raster cells
are centered on the vector points.
centerPointRasterTemplate(rawPointData)centerPointRasterTemplate(rawPointData)
rawPointData |
A |
rawPointData is a SpatVector object that
contains x and y coordinates.
An empty SpatRaster designed to serve as a template for
rasterizing SpatVector objects.
library(terra) # Create point grid coords <- data.frame(x = rep(seq(1:5), times = 5), y = unlist(lapply(1:5, FUN = function(x) { rep(x, times = 5)}))) # Create data and add NAs to simulate uneven bottom depths dd <- data.frame(SURFACE = 1:25, d5M = 6:30, d10M = 11:35, d25M = 16:40) dd$d25M[c(1:5, 18:25)] <- NA dd$d10M[c(3:5, 21:23)] <- NA dd$d5M[c(4, 22)] <- NA dd[,c("x","y")] <- coords # Create SpatialPointsDataFrame sp <- vect(dd, geom = c("x", "y")) # Here's the function template <- centerPointRasterTemplate(rawPointData = sp) class(template)library(terra) # Create point grid coords <- data.frame(x = rep(seq(1:5), times = 5), y = unlist(lapply(1:5, FUN = function(x) { rep(x, times = 5)}))) # Create data and add NAs to simulate uneven bottom depths dd <- data.frame(SURFACE = 1:25, d5M = 6:30, d10M = 11:35, d25M = 16:40) dd$d25M[c(1:5, 18:25)] <- NA dd$d10M[c(3:5, 21:23)] <- NA dd$d5M[c(4, 22)] <- NA dd[,c("x","y")] <- coords # Create SpatialPointsDataFrame sp <- vect(dd, geom = c("x", "y")) # Here's the function template <- centerPointRasterTemplate(rawPointData = sp) class(template)
function to remove or flag coordinates on land, coordinates which are deeper than the known bathymetry at a given xy location, and/or coordinates outside of a specified depth range when cleaning coordinates for species in marine environments
cleanDepth( occs, bathy, land_poly, depth_range, flag = FALSE, bottom_correct = FALSE )cleanDepth( occs, bathy, land_poly, depth_range, flag = FALSE, bottom_correct = FALSE )
occs |
dataframe containing columns named "longitude", "latitude", and "depth," where depth values are positive |
bathy |
a spatRaster of bathymetry in the same units as the occurrence depth, where values are negative |
land_poly |
optional polygon of continents in same projection as occurrence data |
depth_range |
optional range of depth values in the form of c(upper, lower) in which all points should fall between. depth value should be positive, as in the case of the occs |
flag |
default is F, if T, points are not removed, but flagged if they intersect, and a column is added to the output indicating what points have been flagged |
bottom_correct |
default is F, but if T, points are not removed, and a column is added to the output displaying the bathymetry value at the flagged point for corrective purposes |
assumes depths are listed as positive in the occs and depth_range, while the bathymetry layer lists them as negative when below the sea surface
an object of class dataframe containing cleaned or flagged occurrences
library(terra) library(dplyr) # Create sample bathymetry r1 <- rast(ncol=10, nrow=10) values(r1) <- c(-100:0) # Create test occurrences set.seed(0) longitude <- sample(ext(r1)[1]:ext(r1)[2], size = 10, replace = FALSE) set.seed(0) latitude <- sample(ext(r1)[3]:ext(r1)[4], size = 10, replace = FALSE) set.seed(0) depth <- sample(1:100, size = 10, replace = TRUE) occs <- data.frame(longitude, latitude, depth) #Heres the function result <- cleanDepth(bathy = r1, occs = occs)library(terra) library(dplyr) # Create sample bathymetry r1 <- rast(ncol=10, nrow=10) values(r1) <- c(-100:0) # Create test occurrences set.seed(0) longitude <- sample(ext(r1)[1]:ext(r1)[2], size = 10, replace = FALSE) set.seed(0) latitude <- sample(ext(r1)[3]:ext(r1)[4], size = 10, replace = FALSE) set.seed(0) depth <- sample(1:100, size = 10, replace = TRUE) occs <- data.frame(longitude, latitude, depth) #Heres the function result <- cleanDepth(bathy = r1, occs = occs)
assigns coordinates to depth slices given a user provided raster template. if any points are deeper than the lowest depth slice, they are assigned to the lowest depth slice.
depthMatch(occs, rasterTemplate)depthMatch(occs, rasterTemplate)
occs |
dataframe of occurrence records with colnames "longitude," "latitude," and "depth." Depth values should be positive. |
rasterTemplate |
a spatRaster stack where the layer names correspond to the depth value, as a positive number |
rasterTemplate names and values in the depth column of occs should be positive. Points lower than the deepest depth slice are assigned to the deepest depth slice
object of class dataframe, the same occs fed into the function but with the depth values adjusted to match the spatRaster stack depth slices
library(terra) # Create test raster brick r1 <- rast(ncol = 100, nrow = 100) r2 <- rast(ncol = 100, nrow = 100) r3 <- rast(ncol = 100, nrow = 100) rbrick <- c(r1, r2, r3) names(rbrick) <- c(1, 40, 90) # Create test occurrences set.seed(0) longitude <- sample(ext(r1)[1]:ext(r1)[2], size = 10, replace = FALSE) set.seed(0) latitude <- sample(ext(r1)[3]:ext(r1)[4], size = 10, replace = FALSE) set.seed(0) depth <- sample(1:100, size = 10, replace = TRUE) occs <- data.frame(longitude, latitude, depth) # Here's the function result <- depthMatch(occs = occs, rasterTemplate = rbrick)library(terra) # Create test raster brick r1 <- rast(ncol = 100, nrow = 100) r2 <- rast(ncol = 100, nrow = 100) r3 <- rast(ncol = 100, nrow = 100) rbrick <- c(r1, r2, r3) names(rbrick) <- c(1, 40, 90) # Create test occurrences set.seed(0) longitude <- sample(ext(r1)[1]:ext(r1)[2], size = 10, replace = FALSE) set.seed(0) latitude <- sample(ext(r1)[3]:ext(r1)[4], size = 10, replace = FALSE) set.seed(0) depth <- sample(1:100, size = 10, replace = TRUE) occs <- data.frame(longitude, latitude, depth) # Here's the function result <- depthMatch(occs = occs, rasterTemplate = rbrick)
Takes list of rasters of species distributions (interpreted as 1 = presence, 0 = absence), which do not have to have the same extents, and stack them to create an estimate of species richness that matches the extent and resolution of a template.
diversityStack(rasterList, template)diversityStack(rasterList, template)
rasterList |
A |
template |
A |
A SpatRaster
library(terra) rast1 <- rast(ncol=10, nrow=10) values(rast1) <- rep(0:1, 50) rast2 <- rast(ncol=10, nrow=10) values(rast2) <- c(rep(0, 50), rep(1,50)) rastList <- list(rast1, rast2) result <- diversityStack(rasterList = rastList, template = rast2) result plot(result)library(terra) rast1 <- rast(ncol=10, nrow=10) values(rast1) <- rep(0:1, 50) rast2 <- rast(ncol=10, nrow=10) values(rast2) <- c(rep(0, 50), rep(1,50)) rastList <- list(rast1, rast2) result <- diversityStack(rasterList = rastList, template = rast2) result plot(result)
Reduces number of occurrences to resolution of input raster
downsample(occs, rasterTemplate, verbose = TRUE)downsample(occs, rasterTemplate, verbose = TRUE)
occs |
A |
rasterTemplate |
A |
verbose |
|
A data.frame with two columns named "longitude"
and "latitude" or with names that were used when coercing
input data into this format.
library(terra) # Create sample raster r <- rast(ncol=10, nrow=10) values(r) <- 1:100 # Create test occurrences set.seed(0) longitude <- sample(ext(r)[1]:ext(r)[2], size = 10, replace = FALSE) set.seed(0) latitude <- sample(ext(r)[3]:ext(r)[4], size = 10, replace = FALSE) occurrences <- as.data.frame(cbind(longitude,latitude)) # Here's the function result <- downsample(occs = occurrences, rasterTemplate = r)library(terra) # Create sample raster r <- rast(ncol=10, nrow=10) values(r) <- 1:100 # Create test occurrences set.seed(0) longitude <- sample(ext(r)[1]:ext(r)[2], size = 10, replace = FALSE) set.seed(0) latitude <- sample(ext(r)[3]:ext(r)[4], size = 10, replace = FALSE) occurrences <- as.data.frame(cbind(longitude,latitude)) # Here's the function result <- downsample(occs = occurrences, rasterTemplate = r)
A user will likely already have environmental data organized into a list of 'spatRaster' stacks where each element is a different environmental variable and each layer corresponds to a depth slice for background sampling and extracting. However, the maxent_3D function requires the data be input as a list of 'spatRaster' stacks where each element corresponds to a depth slice and each layer is an environmental variable for the purpose of projecting the model back into geographic space, and this function has been developed for ease of conversion.
env_stack_transform(envs_all, envs_names)env_stack_transform(envs_all, envs_names)
envs_all |
a list of 'spatRaster' stacks where each element is a different environmental variable and each layer is a depth slice. |
envs_names |
A 'character' vector of the names of the environmental variables to be applied to the layers of the resulting list of 'spatRaster' stacks |
The extents of each environmental variable layer should match on a depth slice by depth slice basis.
A list of 'spatRaster' stacks where each list element is a depth slice and each layer is an environmental variable.
library(terra) # creating a list of spatRaster stacks where each element is an environmental variable # and each layer is a depth env1_d1 <- rast(ncol = 50, nrow = 50) values(env1_d1) <- sample(c(1:100), size = 2500, replace = TRUE) env2_d1 <- rast(ncol = 50, nrow = 50) values(env2_d1) <- sample(c(1:100), size = 2500, replace = TRUE) env1_d2 <- rast(ncol = 50, nrow = 50) values(env1_d2) <- sample(c(1:100), size = 2500, replace = TRUE) env2_d2 <- rast(ncol = 50, nrow = 50) values(env2_d2) <- sample(c(1:100), size = 2500, replace = TRUE) env1 <- c(env1_d1, env1_d2) env2 <- c(env2_d1, env2_d2) envs <- list(env1, env2) envnames <- c("env1", "env2") # Here's the function result <- env_stack_transform(envs_all = envs, envs_names = envnames)library(terra) # creating a list of spatRaster stacks where each element is an environmental variable # and each layer is a depth env1_d1 <- rast(ncol = 50, nrow = 50) values(env1_d1) <- sample(c(1:100), size = 2500, replace = TRUE) env2_d1 <- rast(ncol = 50, nrow = 50) values(env2_d1) <- sample(c(1:100), size = 2500, replace = TRUE) env1_d2 <- rast(ncol = 50, nrow = 50) values(env1_d2) <- sample(c(1:100), size = 2500, replace = TRUE) env2_d2 <- rast(ncol = 50, nrow = 50) values(env2_d2) <- sample(c(1:100), size = 2500, replace = TRUE) env1 <- c(env1_d1, env1_d2) env2 <- c(env2_d1, env2_d2) envs <- list(env1, env2) envnames <- c("env1", "env2") # Here's the function result <- env_stack_transform(envs_all = envs, envs_names = envnames)
Uses thin plate spline regression from
fields package to interpolate missing two-dimensional
raster values.
interpolateRaster(inputRaster, fast = FALSE, ...)interpolateRaster(inputRaster, fast = FALSE, ...)
inputRaster |
An object of class |
fast |
A logical operator. Setting to |
... |
For any additional arguments passed to |
Missing data values from original raster
are replaced with interpolated values. User has the
option of choosing fastTps to speed calculation,
but be advised that this is only an approximation
of a true thin plate spline.
An object of class raster
library(terra) library(fields) # Create sample raster r <- rast(ncol=50, nrow=50) values(r) <- 1:2500 # Introduce a "hole" values(r)[c(117:127, 167:177, 500:550)] <- NA plot(r) # Patch hole with interpolateRaster interpolatedRaster <- interpolateRaster(r) plot(interpolatedRaster) fastInterp <- interpolateRaster(r, fast = TRUE, aRange = 3.0) plot(fastInterp)library(terra) library(fields) # Create sample raster r <- rast(ncol=50, nrow=50) values(r) <- 1:2500 # Introduce a "hole" values(r)[c(117:127, 167:177, 500:550)] <- NA plot(r) # Patch hole with interpolateRaster interpolatedRaster <- interpolateRaster(r) plot(interpolatedRaster) fastInterp <- interpolateRaster(r, fast = TRUE, aRange = 3.0) plot(fastInterp)
Automatically generates background shapefiles for sampling pseudoabsences and/or background points for niche modeling or species distribution modeling. Delineating background sampling regions can be one of the trickiest parts of generating a good model. Automatically generated background shapefiles should be inspected carefully prior to model use.
Useful references, among others:
Barve N, Barve V, Jiménez-Valverde A, Lira-Noriega A, Maher SP, Peterson AT, Soberón J, Villalobos F. 2011. The crucial role of the accessible area in ecological niche modeling and species distribution modeling. Ecological modelling 222:1810-9.
Merow, C, Smith MJ, Silander JA. 2013. A practical guide to MaxEnt for modeling species’ distributions: what it does, and why inputs and settings matter." Ecography 36: 1058-69.
Murphy SJ. 2021. Sampling units derived from geopolitical boundaries bias biodiversity analyses. Global Ecology and Biogeography 30: 1876-88.
marineBackground(occs, clipToOcean = TRUE, verbose = TRUE, ...)marineBackground(occs, clipToOcean = TRUE, verbose = TRUE, ...)
occs |
A |
clipToOcean |
|
verbose |
|
... |
Additional optional arguments to pass to
|
The meat of this function is a special-case wrapper
around getDynamicAlphaHull() from the rangeBuilder package.
The function documented here is especially useful in cases where
one wants to automatically generate training regions that overlap
the international date line. Regions that exceed the line are cut
and pasted into the appropriate hemisphere instead of being
deleted.
If the argument buff is not supplied, a buffer is
calculated by taking the mean between the 10th and 90th percentile
of horizontal distances between occurrence points.
If getDynamicAlphaHull() cannot satisfy the provided conditions,
the occurrences are buffered and then a minimum convex hull is
drawn around the buffer polygons.
A SpatVector
library(terra) # Create sample raster r <- rast(ncol=10, nrow=10) values(r) <- 1:100 # Create test occurrences set.seed(0) longitude <- sample(-50:50, size = 20, replace = FALSE) set.seed(0) latitude <- sample(-30:30, size = 20, replace = FALSE) occurrences <- as.data.frame(cbind(longitude,latitude)) # Here's the function result <- marineBackground(occs = occurrences, buff = 100000, fraction = .9, partCount = 2, clipToOcean = FALSE)library(terra) # Create sample raster r <- rast(ncol=10, nrow=10) values(r) <- 1:100 # Create test occurrences set.seed(0) longitude <- sample(-50:50, size = 20, replace = FALSE) set.seed(0) latitude <- sample(-30:30, size = 20, replace = FALSE) occurrences <- as.data.frame(cbind(longitude,latitude)) # Here's the function result <- marineBackground(occs = occurrences, buff = 100000, fraction = .9, partCount = 2, clipToOcean = FALSE)
Uses MaxEnt from predicts package to test multiple models with
different feature class combinations, regularization multipliers, and a user
supplied partitioning scheme for training and testing. The function outputs
model objects, model results, as well as prediction SpatRasters of each model
projected back onto geographic space using the supplied environmental SpatRaster
stacks. Note: you will need to install rJava to successfully run this function.
maxent_3D( maxent_df, wanted_fc, wanted_rm, wanted_partition = NULL, projection_layers, occs, depth_list )maxent_3D( maxent_df, wanted_fc, wanted_rm, wanted_partition = NULL, projection_layers, occs, depth_list )
maxent_df |
'data.frame' where the first column is a vector of presences named "p" containing 1's and 0's. Each row represents a cell in the spatRaster volume with an x, y, z coordinate, and 1's are presences while 0's are absences, or background points. Other columns are environmental variable values extracted at the occurrence and background points, and should have the same names as the names of the environmental layers in the projection_layers list of SpatRaster stacks. |
wanted_fc |
a character vector giving what feature class combinations should be tried. "L" refers to linear, "Q" refers to quadratic, "H" refers to hinge, and "P" refers to product. Should be in the format c("L", "Q", "LQ") etc. |
wanted_rm |
regularization multipliers to be tried in format c(1:4) |
wanted_partition |
optional, should be the output of 'partition_3D'. if no partition is supplied, all points will be used for training |
projection_layers |
|
occs |
|
depth_list |
vector of depths corresponding to depth slices of list elements of projection_layers. Should be positive and go from shallowest depth to deepest depth |
The names of the projection_layers should be the same as the column names of the environmental variables in maxent_df. The number of models output will be the number of feature classes multiplied by the number of regularization multipliers. For example, a wanted_fc of c("L", "Q", "P") and a wanted_rm of c(1:3) will output 9 total models.
An object of class list with four components:
$models, a list containing each model object produced.
$results, a data.frame where each row is a model corresponding to the list element
in $models and $predictions. If there was no partition supplied for training and
testing, each column will report the feature class, regularization multiplier, AUC,
total coefficients, nonzero coefficients, AICc, and delta AICc. if a partition is
used, it will report the average of these statistics across all partitions.
$predictions, a list of spatRaster stacks where each list element is a model
corresponding to the rows of $results and elements of $models projected onto the
supplied projection_layers. Each layer in the stack is a depth slice.
$partition_results, a list object of the same length as the number of partition
groups containing a data.frame of results for each model for each partition. Only
produced if a partition is supplied.
library(dplyr) library(predicts) library(terra) # creating list of spatraster stacks where each element is a depth slice r1_d1 <- rast(ncol = 100, nrow = 100) set.seed(0) values(r1_d1) <- as.numeric(sample(c(1:100), size = 1000, replace = TRUE)) r2_d1 <- rast(ncol = 100, nrow = 100) set.seed(0) values(r2_d1) <- as.numeric(sample(c(1:1000), size = 1000, replace = FALSE)) r1_d2 <- r1_d1 values(r1_d2) <- as.numeric(values(r1_d1)+10) r2_d2 <- r2_d1 values(r2_d2) <- as.numeric(values(r2_d1)+10) d1 <- c(r1_d1, r2_d1) names(d1) <- c("valsr1", "valsr2") d2 <- c(r1_d2, r2_d2) names(d2) <- c("valsr1", "valsr2") envlist <- list(d1, d2) # creating occs and bgs set.seed(0) occs <- sample(c(1:nrow(crds(envlist[[1]][[1]]))), size = 50, replace = FALSE) bgs <- sample(c(1:nrow(crds(envlist[[1]][[1]]))), size = 500, replace = FALSE) occ_indices <- sample(c(1:nrow(crds(envlist[[1]][[1]]))), size = 50, replace = FALSE) bg_indices <- sample(c(1:nrow(crds(envlist[[1]][[1]]))), size = 500, replace = FALSE) occs_d1 <- crds(envlist[[1]][[1]])[occ_indices[1:25],] occs_d2 <- crds(envlist[[2]][[1]])[occ_indices[26:50],] bg_d1 <- crds(envlist[[1]][[1]])[bg_indices[1:250],] bg_d2 <- crds(envlist[[2]][[1]])[bg_indices[251:500],] # extracting at occs and bgs occ_valsr1_d1 <- extract(envlist[[1]][[1]], occs_d1) occ_valsr1_d2 <- extract(envlist[[2]][[1]], occs_d2) occ_valsr2_d1 <- extract(envlist[[1]][[2]], occs_d1) occ_valsr2_d2 <- extract(envlist[[2]][[2]], occs_d2) occ_valsr1 <- rbind(occ_valsr1_d1, occ_valsr1_d2) occ_valsr2 <- rbind(occ_valsr2_d1, occ_valsr2_d2) bg_valsr1_d1 <- extract(envlist[[1]][[1]], bg_d1) bg_valsr1_d2 <- extract(envlist[[2]][[1]], bg_d2) bg_valsr2_d1 <- extract(envlist[[1]][[2]], bg_d1) bg_valsr2_d2 <- extract(envlist[[2]][[2]], bg_d2) bg_valsr1 <- rbind(bg_valsr1_d1, bg_valsr1_d2) bg_valsr2 <- rbind(bg_valsr2_d1, bg_valsr2_d2) valsr1 <- rbind(occ_valsr1, bg_valsr1) valsr2 <- rbind(occ_valsr2, bg_valsr2) p1 <- rep(1, times = 50) p0 <- rep(0, times = 500) p <- c(p1, p0) maxdf <- data.frame(p, valsr1, valsr2) coords <- rbind(occs_d1, occs_d2) colnames(coords) <- c("longitude", "latitude") depth_vector <- c(rep(1, times = 25), rep(2, times = 25)) # Use data.frame instead of cbind so $depth is completely valid inside maxent_3D occs_dataframe <- data.frame(coords, depth = depth_vector) # Pass the clean data.frame to the function if(requireNamespace("rJava", quietly = TRUE)){ result <- maxent_3D(maxent_df = maxdf, wanted_fc = c("L", "Q"), wanted_rm = c(1:2), projection_layers = envlist, occs = occs_dataframe, depth_list = c(1,2)) }library(dplyr) library(predicts) library(terra) # creating list of spatraster stacks where each element is a depth slice r1_d1 <- rast(ncol = 100, nrow = 100) set.seed(0) values(r1_d1) <- as.numeric(sample(c(1:100), size = 1000, replace = TRUE)) r2_d1 <- rast(ncol = 100, nrow = 100) set.seed(0) values(r2_d1) <- as.numeric(sample(c(1:1000), size = 1000, replace = FALSE)) r1_d2 <- r1_d1 values(r1_d2) <- as.numeric(values(r1_d1)+10) r2_d2 <- r2_d1 values(r2_d2) <- as.numeric(values(r2_d1)+10) d1 <- c(r1_d1, r2_d1) names(d1) <- c("valsr1", "valsr2") d2 <- c(r1_d2, r2_d2) names(d2) <- c("valsr1", "valsr2") envlist <- list(d1, d2) # creating occs and bgs set.seed(0) occs <- sample(c(1:nrow(crds(envlist[[1]][[1]]))), size = 50, replace = FALSE) bgs <- sample(c(1:nrow(crds(envlist[[1]][[1]]))), size = 500, replace = FALSE) occ_indices <- sample(c(1:nrow(crds(envlist[[1]][[1]]))), size = 50, replace = FALSE) bg_indices <- sample(c(1:nrow(crds(envlist[[1]][[1]]))), size = 500, replace = FALSE) occs_d1 <- crds(envlist[[1]][[1]])[occ_indices[1:25],] occs_d2 <- crds(envlist[[2]][[1]])[occ_indices[26:50],] bg_d1 <- crds(envlist[[1]][[1]])[bg_indices[1:250],] bg_d2 <- crds(envlist[[2]][[1]])[bg_indices[251:500],] # extracting at occs and bgs occ_valsr1_d1 <- extract(envlist[[1]][[1]], occs_d1) occ_valsr1_d2 <- extract(envlist[[2]][[1]], occs_d2) occ_valsr2_d1 <- extract(envlist[[1]][[2]], occs_d1) occ_valsr2_d2 <- extract(envlist[[2]][[2]], occs_d2) occ_valsr1 <- rbind(occ_valsr1_d1, occ_valsr1_d2) occ_valsr2 <- rbind(occ_valsr2_d1, occ_valsr2_d2) bg_valsr1_d1 <- extract(envlist[[1]][[1]], bg_d1) bg_valsr1_d2 <- extract(envlist[[2]][[1]], bg_d2) bg_valsr2_d1 <- extract(envlist[[1]][[2]], bg_d1) bg_valsr2_d2 <- extract(envlist[[2]][[2]], bg_d2) bg_valsr1 <- rbind(bg_valsr1_d1, bg_valsr1_d2) bg_valsr2 <- rbind(bg_valsr2_d1, bg_valsr2_d2) valsr1 <- rbind(occ_valsr1, bg_valsr1) valsr2 <- rbind(occ_valsr2, bg_valsr2) p1 <- rep(1, times = 50) p0 <- rep(0, times = 500) p <- c(p1, p0) maxdf <- data.frame(p, valsr1, valsr2) coords <- rbind(occs_d1, occs_d2) colnames(coords) <- c("longitude", "latitude") depth_vector <- c(rep(1, times = 25), rep(2, times = 25)) # Use data.frame instead of cbind so $depth is completely valid inside maxent_3D occs_dataframe <- data.frame(coords, depth = depth_vector) # Pass the clean data.frame to the function if(requireNamespace("rJava", quietly = TRUE)){ result <- maxent_3D(maxent_df = maxdf, wanted_fc = c("L", "Q"), wanted_rm = c(1:2), projection_layers = envlist, occs = occs_dataframe, depth_list = c(1,2)) }
Calculates multivariate environmental similarity surface based on model calibration and projection data
MESS3D(calibration, projection)MESS3D(calibration, projection)
calibration |
A |
projection |
A named |
MESS3D is a wrapper around MESS from the modEvA
package. It calculates MESS for each depth layer. Negative values
indicate areas of extrapolation which should be interpreted with
caution (see Elith et al, 2010 in MEE).
A SpatRaster vector with MESS scores for each
voxel; layer names correspond to layer names of first
SpatRaster vector in projection list.
The calibration dataset should include both presences and background/pseudoabsence points used to calibrate an ecological niche model.
Elith J, Kearney M, and Phillips S. 2010. The art of modelling range-shifting species. Methods in Ecology and Evolution, 1, 330-342.
library(terra) library(dplyr) # Create sample rasterBricks r1 <- rast(ncol=10, nrow=10) values(r1) <- 1:100 r2 <- rast(ncol=10, nrow=10) values(r2) <- c(rep(20, times = 50), rep(60, times = 50)) r3 <- rast(ncol=10, nrow=10) values(r3) <- 8 envBrick1 <- c(r1, r2, r3) names(envBrick1) <- c(0, 10, 30) r1 <- rast(ncol=10, nrow=10) values(r1) <- 100:1 r2 <- rast(ncol=10, nrow=10) values(r2) <- c(rep(10, times = 50), rep(20, times = 50)) r3 <- rast(ncol=10, nrow=10) values(r3) <- c(rep(c(10,20,30,25), times = 25)) envBrick2 <- c(r1, r2, r3) names(envBrick2) <- c(0, 10, 30) rastList <- list("temperature" = envBrick1, "salinity" = envBrick2) # Create test reference set set.seed(0) longitude <- sample(ext(envBrick1)[1]:ext(envBrick1)[2], size = 10, replace = FALSE) set.seed(0) latitude <- sample(ext(envBrick1)[3]:ext(envBrick1)[4], size = 10, replace = FALSE) set.seed(0) depth <- sample(0:35, size = 10, replace = TRUE) occurrences <- as.data.frame(cbind(longitude,latitude,depth)) # Sample data at occurrences to characterize calibration region cal_temp <- xyzSample(occurrences, rastList$temperature) cal_sal <- xyzSample(occurrences, rastList$salinity) calibration <- data.frame( temperature = cal_temp, salinity = cal_sal) # Run the function messStack <- MESS3D(calibration = calibration, projection = rastList) plot(messStack)library(terra) library(dplyr) # Create sample rasterBricks r1 <- rast(ncol=10, nrow=10) values(r1) <- 1:100 r2 <- rast(ncol=10, nrow=10) values(r2) <- c(rep(20, times = 50), rep(60, times = 50)) r3 <- rast(ncol=10, nrow=10) values(r3) <- 8 envBrick1 <- c(r1, r2, r3) names(envBrick1) <- c(0, 10, 30) r1 <- rast(ncol=10, nrow=10) values(r1) <- 100:1 r2 <- rast(ncol=10, nrow=10) values(r2) <- c(rep(10, times = 50), rep(20, times = 50)) r3 <- rast(ncol=10, nrow=10) values(r3) <- c(rep(c(10,20,30,25), times = 25)) envBrick2 <- c(r1, r2, r3) names(envBrick2) <- c(0, 10, 30) rastList <- list("temperature" = envBrick1, "salinity" = envBrick2) # Create test reference set set.seed(0) longitude <- sample(ext(envBrick1)[1]:ext(envBrick1)[2], size = 10, replace = FALSE) set.seed(0) latitude <- sample(ext(envBrick1)[3]:ext(envBrick1)[4], size = 10, replace = FALSE) set.seed(0) depth <- sample(0:35, size = 10, replace = TRUE) occurrences <- as.data.frame(cbind(longitude,latitude,depth)) # Sample data at occurrences to characterize calibration region cal_temp <- xyzSample(occurrences, rastList$temperature) cal_sal <- xyzSample(occurrences, rastList$salinity) calibration <- data.frame( temperature = cal_temp, salinity = cal_sal) # Run the function messStack <- MESS3D(calibration = calibration, projection = rastList) plot(messStack)
Samples in 2D at resolution of raster
mSampling2D(occs, rasterTemplate, mShp, verbose = TRUE)mSampling2D(occs, rasterTemplate, mShp, verbose = TRUE)
occs |
A dataframe with at least two columns named "longitude" and "latitude", or that can be coerced into this format. |
rasterTemplate |
A |
mShp |
A shapefile defining the area from which background points should be sampled. |
verbose |
|
This function is designed to sample background points
for distributional modeling in two dimensions. The returned
data.frame contains all points from across the designated
background. It is up to the user to determine how to
appropriately sample from those background points.
A data.frame with 2D coordinates of points
for background sampling.
library(terra) # Create sample raster r <- rast(ncol=10, nrow=10) values(r) <- 1:100 # Create test occurrences set.seed(0) longitude <- sample(ext(r)[1]:ext(r)[2], size = 10, replace = FALSE) set.seed(0) latitude <- sample(ext(r)[3]:ext(r)[4], size = 10, replace = FALSE) occurrences <- data.frame(longitude,latitude) # Generate background sampling buffer buffPts <- vect(occurrences, c("longitude", "latitude")) crs(buffPts) <- crs(r) mShp <- aggregate(buffer(buffPts, width = 1000000)) # Here's the function result <- mSampling2D(occs = occurrences, rasterTemplate = r, mShp = mShp)library(terra) # Create sample raster r <- rast(ncol=10, nrow=10) values(r) <- 1:100 # Create test occurrences set.seed(0) longitude <- sample(ext(r)[1]:ext(r)[2], size = 10, replace = FALSE) set.seed(0) latitude <- sample(ext(r)[3]:ext(r)[4], size = 10, replace = FALSE) occurrences <- data.frame(longitude,latitude) # Generate background sampling buffer buffPts <- vect(occurrences, c("longitude", "latitude")) crs(buffPts) <- crs(r) mShp <- aggregate(buffer(buffPts, width = 1000000)) # Here's the function result <- mSampling2D(occs = occurrences, rasterTemplate = r, mShp = mShp)
Samples XYZ coordinates from a shapefile from maximum to minimum occurrence depth at XYZ resolution of envBrick.
mSampling3D(occs, envBrick, mShp, depthLimit = "all", verbose = TRUE)mSampling3D(occs, envBrick, mShp, depthLimit = "all", verbose = TRUE)
occs |
A |
envBrick |
A |
mShp |
A shapefile defining the area from which background points should be sampled. |
depthLimit |
An argument controlling the depth
extent of sampling. Refer to |
verbose |
|
This function is designed to sample background points for
distributional modeling in three dimensions. If a voxel (3D pixel)
in the SpatRaster vector intersects with an occurrence from
occs, it is removed. Note that this function returns points
representing every voxel in the background area within the
specified depth range. It is up to the user to downsample from
these data as necessary, depending on the model type being used.
depthLimit argument options:
occs Samples background from the full depth extent of occs.
all Samples background from the full depth extent of envBrick.
A vector of length 2 with maximum and minimum depth values from
which to sample.
A data.frame with 3D coordinates of points for background
sampling.
library(terra) # Create test raster r1 <- rast(ncol=10, nrow=10) values(r1) <- 1:100 r2 <- rast(ncol=10, nrow=10) values(r2) <- c(rep(20, times = 50), rep(60, times = 50)) r3 <- rast(ncol=10, nrow=10) values(r3) <- 8 envBrick <- c(r1, r2, r3) names(envBrick) <- c(0, 10, 30) # Create test occurrences set.seed(0) longitude <- sample(ext(envBrick)[1]:ext(envBrick)[2], size = 10, replace = FALSE) set.seed(0) latitude <- sample(ext(envBrick)[3]:ext(envBrick)[4], size = 10, replace = FALSE) set.seed(0) depth <- sample(0:35, size = 10, replace = TRUE) occurrences <- data.frame(longitude,latitude,depth) # Generate background sampling buffer buffPts <- vect(occurrences, c("longitude", "latitude")) crs(buffPts) <- crs(envBrick) mShp <- aggregate(buffer(buffPts, width = 1000000)) # Here's the function occSample3d <- mSampling3D(occs = occurrences, envBrick = envBrick, mShp = mShp, depthLimit = "occs")library(terra) # Create test raster r1 <- rast(ncol=10, nrow=10) values(r1) <- 1:100 r2 <- rast(ncol=10, nrow=10) values(r2) <- c(rep(20, times = 50), rep(60, times = 50)) r3 <- rast(ncol=10, nrow=10) values(r3) <- 8 envBrick <- c(r1, r2, r3) names(envBrick) <- c(0, 10, 30) # Create test occurrences set.seed(0) longitude <- sample(ext(envBrick)[1]:ext(envBrick)[2], size = 10, replace = FALSE) set.seed(0) latitude <- sample(ext(envBrick)[3]:ext(envBrick)[4], size = 10, replace = FALSE) set.seed(0) depth <- sample(0:35, size = 10, replace = TRUE) occurrences <- data.frame(longitude,latitude,depth) # Generate background sampling buffer buffPts <- vect(occurrences, c("longitude", "latitude")) crs(buffPts) <- crs(envBrick) mShp <- aggregate(buffer(buffPts, width = 1000000)) # Here's the function occSample3d <- mSampling3D(occs = occurrences, envBrick = envBrick, mShp = mShp, depthLimit = "occs")
A convenient wrapper around ggplot
to generate a formatted plot of a single raster.
oneRasterPlot( rast, land = NA, landCol = "black", scaleRange = NA, graticule = TRUE, title = "A Raster", verbose = TRUE, ... )oneRasterPlot( rast, land = NA, landCol = "black", scaleRange = NA, graticule = TRUE, title = "A Raster", verbose = TRUE, ... )
rast |
A single |
land |
An optional coastline polygon shapefile
of types |
landCol |
Color for land on map. |
scaleRange |
Optional numeric vector containing
maximum and minimum values for color scale. Helpful
when making multiple plots for comparison. Defaults
to minimum and maximum of input |
graticule |
|
title |
A title for the plot. |
verbose |
|
... |
Additional optional arguments to pass to
|
A plot of mapping the values of the input raster layer
library(terra) rast <- rast(ncol=10, nrow=10) values(rast) <- seq(0,99, 1) oneRasterPlot(rast = rast)library(terra) rast <- rast(ncol=10, nrow=10) values(rast) <- seq(0,99, 1) oneRasterPlot(rast = rast)
Creates partition schemes in 3D for model training and testing.
Can create block or kfold partitions, with output as a vector or list
partition_3D( maxent_df, coord_df, which_partition = "k.fold", kfolds = NULL, orientation = "lon_lat", return_format = "vector", ensure_all_folds = TRUE, max_attempts = 100, na_strategy = "NA" )partition_3D( maxent_df, coord_df, which_partition = "k.fold", kfolds = NULL, orientation = "lon_lat", return_format = "vector", ensure_all_folds = TRUE, max_attempts = 100, na_strategy = "NA" )
maxent_df |
Data frame with column 'p' (1=presence, 0=absence). |
coord_df |
Data frame with 'longitude','latitude','depth' aligned with maxent_df. |
which_partition |
"k.fold" (default) or "block". |
kfolds |
Integer >= 2 for k.fold. |
orientation |
For "block": "lon_lat" (default) or "lat_lon". |
return_format |
"vector" (default) or "list". |
ensure_all_folds |
(k.fold) Ensure all folds appear among presences with known depth (default TRUE). |
max_attempts |
(k.fold) Retry cap when ensure_all_folds=TRUE (default 100). |
na_strategy |
(k.fold) How to handle NA depth rows: "NA" (default) leaves them NA; "random" assigns a random fold. |
maxent_df and coord_df should be the same dataframes to be provided
to maxent_3D() for model production.
The spatial block partition, similarly to a traditional 2D partition,
will separate the occurrences into 4 groups of equal (or about equal) size,
with two groups making up two blocks on the upper half of the depth distribution
and two groups on the lower half. Background points are assigned to spatial groups
according to how they fall into the spatial partitions delimited by the occurrences.
Note that this means the number of background points in each partition may not be as
close to equal as the occurrences.
If return_format="vector": integer vector (1..k or 1..4), length == nrow(maxent_df). Unassignable rows get NA. If "list": list(occ_partitions, bg_partitions).
# create test dataframe occ <- rep(1, times = 10) bg <- rep(0, times = 1000) env1 <- sample(c(1:100), size = 1010, replace = TRUE) env2 <- sample(c(1:1000), size = 1010, replace = TRUE) p <- c(occ, bg) testdf <- data.frame(p, env1, env2) # create test coord data r <- terra::rast(ncol = 100, nrow = 100) set.seed(0) longitude <- sample(terra::ext(r)[1]:terra::ext(r)[2], size = 1010, replace = TRUE) set.seed(0) latitude <- sample(terra::ext(r)[3]:terra::ext(r)[4], size = 1010, replace = TRUE) depth <- sample(c(0, 5, 10, 15, 20, 25, 30, 35, 40, 45), size = 1010, replace = TRUE) test_coords <- data.frame(longitude, latitude, depth) # Here's the function result_kfold <- partition_3D(maxent_df = testdf, coord_df = test_coords, which_partition = 'k.fold', kfolds = 3) result_block <- partition_3D(maxent_df = testdf, coord_df = test_coords, which_partition = 'block', orientation = 'lat_lon')# create test dataframe occ <- rep(1, times = 10) bg <- rep(0, times = 1000) env1 <- sample(c(1:100), size = 1010, replace = TRUE) env2 <- sample(c(1:1000), size = 1010, replace = TRUE) p <- c(occ, bg) testdf <- data.frame(p, env1, env2) # create test coord data r <- terra::rast(ncol = 100, nrow = 100) set.seed(0) longitude <- sample(terra::ext(r)[1]:terra::ext(r)[2], size = 1010, replace = TRUE) set.seed(0) latitude <- sample(terra::ext(r)[3]:terra::ext(r)[4], size = 1010, replace = TRUE) depth <- sample(c(0, 5, 10, 15, 20, 25, 30, 35, 40, 45), size = 1010, replace = TRUE) test_coords <- data.frame(longitude, latitude, depth) # Here's the function result_kfold <- partition_3D(maxent_df = testdf, coord_df = test_coords, which_partition = 'k.fold', kfolds = 3) result_block <- partition_3D(maxent_df = testdf, coord_df = test_coords, which_partition = 'block', orientation = 'lat_lon')
This script plots a semitransparent layer of suitable habitat for each depth layer. The redder the color, the shallower the layer, the bluer, the deeper. The more saturated the color, the more layers with suitable habitat.
plotLayers( rast, land = NA, landCol = "black", title = NULL, graticule = TRUE, ... )plotLayers( rast, land = NA, landCol = "black", title = NULL, graticule = TRUE, ... )
rast |
A |
land |
An optional coastline polygon shapefile
of types |
landCol |
Color for land on map. |
title |
A title for the plot. If not title is
supplied, the title "Suitability from (MINIMUM
DEPTH) to (MAXIMUM DEPTH)" is inferred from
names of |
graticule |
Do you want a grid of lon/lat lines? |
... |
Additional optional arguments. |
A plot of class recordedplot
Only include the depth layers that you actually want to plot.
library(terra) rast1 <- rast(ncol=10, nrow=10) values(rast1) <- rep(0:1, 50) rast2 <- rast(ncol=10, nrow=10) values(rast2) <- c(rep(0, 50), rep(1,50)) rast3 <- rast(ncol=10, nrow=10) values(rast3) <- rep(c(1,0,0,1), 25) distBrick <- c(rast1, rast2, rast3) plotLayers(distBrick)library(terra) rast1 <- rast(ncol=10, nrow=10) values(rast1) <- rep(0:1, 50) rast2 <- rast(ncol=10, nrow=10) values(rast2) <- c(rep(0, 50), rep(1,50)) rast3 <- rast(ncol=10, nrow=10) values(rast3) <- rep(c(1,0,0,1), 25) distBrick <- c(rast1, rast2, rast3) plotLayers(distBrick)
A convenient wrapper around ggplot
to generate formatted plots comparing two sets of
occurrence point plots.
pointCompMap( occs1, occs2, spName, land = NA, occs1Col = "#bd0026", occs2Col = "#fd8d3c", agreeCol = "black", occs1Name = "Set 1", occs2Name = "Set 2", landCol = "gray", waterCol = "steelblue", ptSize = 1, verbose = TRUE, ... )pointCompMap( occs1, occs2, spName, land = NA, occs1Col = "#bd0026", occs2Col = "#fd8d3c", agreeCol = "black", occs1Name = "Set 1", occs2Name = "Set 2", landCol = "gray", waterCol = "steelblue", ptSize = 1, verbose = TRUE, ... )
occs1 |
A |
occs2 |
A |
spName |
A character string with the species name to be used in the plot title. |
land |
An optional coastline polygon shapefile
of types |
occs1Col |
Color for occurrence points on map |
occs2Col |
Color for occurrence points on map |
agreeCol |
Color for occurrence points shared
between |
occs1Name |
An optional name for the first set
of occurrences, which will be color-coded to
|
occs2Name |
An optional name for the first set
of occurrences, which will be color-coded to
|
landCol |
Color for land on map |
waterCol |
Color for water on map |
ptSize |
|
verbose |
|
... |
Additional optional arguments to pass to
|
A ggplot plot object.
The x and y column names of occs1 and occs2
must match.
set.seed(5) occs <- data.frame(cbind(decimalLatitude = sample(seq(7,35), 24), decimalLongitude = sample(seq(-97, -70), 24))) set.seed(0) occs1 <- occs[sample(1:nrow(occs), size = 12, replace = FALSE),] set.seed(10) occs2 <- occs[sample(1:nrow(occs), size = 12, replace = FALSE),] pointCompMap(occs1 = occs1, occs2 = occs2, occs1Col = "red", occs2Col = "orange", agreeCol = "purple", occs1Name = "2D", occs2Name = "3D", waterCol = "steelblue", spName = "Steindachneria argentea", ptSize = 2, verbose = FALSE)set.seed(5) occs <- data.frame(cbind(decimalLatitude = sample(seq(7,35), 24), decimalLongitude = sample(seq(-97, -70), 24))) set.seed(0) occs1 <- occs[sample(1:nrow(occs), size = 12, replace = FALSE),] set.seed(10) occs2 <- occs[sample(1:nrow(occs), size = 12, replace = FALSE),] pointCompMap(occs1 = occs1, occs2 = occs2, occs1Col = "red", occs2Col = "orange", agreeCol = "purple", occs1Name = "2D", occs2Name = "3D", waterCol = "steelblue", spName = "Steindachneria argentea", ptSize = 2, verbose = FALSE)
A convenient wrapper around ggplot to generate formatted occurrence point plots.
pointMap( occs, spName, land = NA, ptCol = "#bd0026", landCol = "gray", waterCol = "steelblue", ptSize = 1, verbose = TRUE, ... )pointMap( occs, spName, land = NA, ptCol = "#bd0026", landCol = "gray", waterCol = "steelblue", ptSize = 1, verbose = TRUE, ... )
occs |
A |
spName |
A character string with the species name to be used in the plot title. |
land |
An optional coastline polygon shapefile
of types |
ptCol |
Color for occurrence points on map |
landCol |
Color for land on map |
waterCol |
Color for water on map |
ptSize |
|
verbose |
|
... |
Additional optional arguments to pass to
|
A ggplot plot object.
occs <- read.csv(system.file("extdata/Steindachneria_argentea.csv", package='voluModel')) spName <- "Steindachneria argentea" pointMap(occs = occs, spName = spName, land = rnaturalearth::ne_countries(scale = "small", returnclass = "sf")[1])occs <- read.csv(system.file("extdata/Steindachneria_argentea.csv", package='voluModel')) spName <- "Steindachneria argentea" pointMap(occs = occs, spName = spName, land = rnaturalearth::ne_countries(scale = "small", returnclass = "sf")[1])
A convenient wrapper around terra::plot
to generate formatted plots comparing two rasters.
This is used in the context of voluModel to
overlay semi-transparent distributions (coded as 1)
in two different RasterLayers.
rasterComp( rast1 = NULL, rast2 = NULL, col1 = "#1b9e777F", col2 = "#7570b37F", rast1Name = "Set 1", rast2Name = "Set 2", land = NA, landCol = "black", title = "A Raster Comparison", graticule = TRUE, ... )rasterComp( rast1 = NULL, rast2 = NULL, col1 = "#1b9e777F", col2 = "#7570b37F", rast1Name = "Set 1", rast2Name = "Set 2", land = NA, landCol = "black", title = "A Raster Comparison", graticule = TRUE, ... )
rast1 |
A single |
rast2 |
A single |
col1 |
Color for |
col2 |
Color for |
rast1Name |
An optional name for the first set
of occurrences, which will be color-coded to
|
rast2Name |
An optional name for the first set
of occurrences, which will be color-coded to
|
land |
An optional coastline polygon shapefile
of types |
landCol |
Color for land on map. |
title |
A title for the plot. |
graticule |
Do you want a grid of lon/lat lines? |
... |
Additional optional arguments to pass to
|
A plot of class recordedplot overlaying mapped,
semitransparent extents of the input rasters
The extents of rast1 and rast2
must match.
library(terra) rast1 <- rast(ncol=10, nrow=10) values(rast1) <- rep(0:1, 50) rast2 <- rast(ncol=10, nrow=10) values(rast2) <- c(rep(0, 50), rep(1,50)) rasterComp(rast1 = rast1, rast2 = rast2)library(terra) rast1 <- rast(ncol=10, nrow=10) values(rast1) <- rep(0:1, 50) rast2 <- rast(ncol=10, nrow=10) values(rast2) <- c(rep(0, 50), rep(1,50)) rasterComp(rast1 = rast1, rast2 = rast2)
Uses thin plate spline regression from
fields package to smooth raster values.
smoothRaster(inputRaster, fast = FALSE, ...)smoothRaster(inputRaster, fast = FALSE, ...)
inputRaster |
An object of class |
fast |
A logical operator. Setting to |
... |
For any additional arguments passed to |
Original raster is smoothed using a thin
plate spline. This may be desirable in cases where
the user has a reasonable expectation of spatial autocorrelation,
but observes putative measurement errors in a raster. The user has
the option of choosing fastTps to speed calculation,
but be advised that this is only an approximation
of a true thin plate spline.
An object of class SpatRaster
library(terra) library(fields) # Create sample raster r <- rast(ncol=100, nrow=100) values(r) <- 1:10000 # Introduce a "bubble" values(r)[720:725] <- 9999 plot(r) # Smooth bubble with smoothRaster fastSmooth <- smoothRaster(r, fast = TRUE, aRange = 10.0) plot(fastSmooth)library(terra) library(fields) # Create sample raster r <- rast(ncol=100, nrow=100) values(r) <- 1:10000 # Introduce a "bubble" values(r)[720:725] <- 9999 plot(r) # Smooth bubble with smoothRaster fastSmooth <- smoothRaster(r, fast = TRUE, aRange = 10.0) plot(fastSmooth)
Creates 3D presence/absence layers by setting all values in a suitability layer above a certain threshold to 1 and all values below that threshold to 0. The threshold value is determined by the sensitivity given by the user, where the sensitivity is the percentage of occurrences to be counted as present in the resulting presence/absence layer. For example, with a sensitivity of 0.9 or 90%, The threshold is the suitability value where 90% of the occurrence points fall on grid cells with suitability scores above that value.
The function can try multiple different sensitivity levels, and for each will output a presence/absence spatRaster stack, as well as the specificity (proportion of psuedoabsences correctly considered absent), the true skill statistic (TSS), which is a measure of how well the threshold balances commission and omission error, and the suitability value of the threshold.
The user can also downweight sensitivity when calculating TSS with the "weights" parameter. This may be important to do as 3D niche models often have a higher ratio of psuedoabsences to occurrences than 2D models.
threshold_3D( predicted_layers, thresholding_vals, maxent_df, coord_df, weights = NULL )threshold_3D( predicted_layers, thresholding_vals, maxent_df, coord_df, weights = NULL )
predicted_layers |
A spatRaster stack of suitability layers where each layer corresponds to a depth slice |
thresholding_vals |
A vector of sensitivity levels for creating different thresholds and presence/absence layers. If one wanted to test sensitivities of 90%, and 95%, this would be input as c(0.9, 0.95) |
maxent_df |
'data.frame' where first column is a vector of presences named "p" containing 1's and 0's. Each row represents a cell in the spatRaster volume with an x, y, z coordinate, and 1's are presences while 0's are absences, or background points. Other columns are environmental variable values extracted at the occurrence and background points. |
coord_df |
A dataframe containing the longitude, latitude, and depth for each cell in maxent_df, named "longitude," "latitude," and "depth" |
weights |
a numeric giving what the sensitivity should be downweighted by. If no value is given, TSS will be calculated according to its original formula |
a 'list' with two components:
$threshold_layers, a spatRaster stack of presence absence rasters thresholded at the sensitivity with the highest TSS
$tss_results, a 'data.frame' containing the specificity, TSS, and suitability score for each sensitivity given in thresholding_vals
Allouche O, Tsoar A, and Kadmon R. 2006. Assessing the accuracy of species distribution models: prevalence, kappa and the true skill statistic (TSS). Journal of Applied Ecology, 43, 1223-1232.
library(terra) library(dplyr) # creating list of spatRaster stacks where each element is a depth slice r1_d1 <- rast(ncol = 100, nrow = 100) set.seed(0) values(r1_d1) <- sample(c(1:100), size = 1000, replace = TRUE) r2_d1 <- rast(ncol = 100, nrow = 100) set.seed(0) values(r2_d1) <- sample(c(1:1000), size = 1000, replace = TRUE) r1_d2 <- r1_d1 values(r1_d2) <- values(r1_d1)+10 r2_d2 <- r2_d1 values(r2_d2) <- values(r2_d1)+10 d1 <- c(r1_d1, r2_d1) names(d1) <- c("valsr1", "valsr2") d2 <- c(r1_d2, r2_d2) names(d2) <- c("valsr1", "valsr2") envlist <- list(d1, d2) # creating occs and bgs set.seed(0) occs <- sample(c(1:nrow(crds(envlist[[1]][[1]]))), size = 50, replace = FALSE) bgs <- sample(c(1:nrow(crds(envlist[[1]][[1]]))), size = 500, replace = FALSE) occs_d1 <- crds(envlist[[1]][[1]])[occs[1:25],] occs_d2 <- crds(envlist[[2]][[1]])[occs[26:50],] bg_d1 <- crds(envlist[[1]][[1]])[bgs[1:250],] bg_d2 <- crds(envlist[[1]][[1]])[bgs[251:500],] # extracting at occs and bgs occ_valsr1_d1 <- extract(envlist[[1]][[1]], occs_d1) occ_valsr1_d2 <- extract(envlist[[2]][[1]], occs_d2) occ_valsr2_d1 <- extract(envlist[[1]][[2]], occs_d1) occ_valsr2_d2 <- extract(envlist[[2]][[2]], occs_d2) occ_valsr1 <- rbind(occ_valsr1_d1, occ_valsr1_d2) occ_valsr2 <- rbind(occ_valsr2_d1, occ_valsr2_d2) bg_valsr1_d1 <- extract(envlist[[1]][[1]], bg_d1) bg_valsr1_d2 <- extract(envlist[[2]][[1]], bg_d2) bg_valsr2_d1 <- extract(envlist[[1]][[2]], bg_d1) bg_valsr2_d2 <- extract(envlist[[2]][[2]], bg_d2) bg_valsr1 <- rbind(bg_valsr1_d1, bg_valsr1_d2) bg_valsr2 <- rbind(bg_valsr2_d1, bg_valsr2_d2) valsr1 <- rbind(occ_valsr1, bg_valsr1) valsr2 <- rbind(occ_valsr2, bg_valsr2) p1 <- rep(1, times = 50) p0 <- rep(0, times = 500) p <- c(p1, p0) maxdf <- data.frame(p, valsr1, valsr2) # creating coord_df coord_df <- rbind(occs_d1, occs_d2, bg_d1, bg_d2) z <- c(rep(1, times = 25), rep(2, times = 25), rep(1, times = 250), rep(2, times = 250)) coord_df <- cbind(coord_df, z) colnames(coord_df) <- c("longitude", "latitude", "depth") # creating suitability rasters suitd1 <- d1[[1]] values(suitd1) <- runif(min = 0, max = 1, n = length(values(suitd1))) suitd2 <- d2[[1]] values(suitd2) <- runif(min = 0, max = 1, n = length(values(suitd2))) suit <- c(suitd1, suitd2) # here's the function result <- threshold_3D(predicted_layers = suit, thresholding_vals = c(0.9, 0.95), maxent_df = maxdf, coord_df = coord_df, weights = 2/3)library(terra) library(dplyr) # creating list of spatRaster stacks where each element is a depth slice r1_d1 <- rast(ncol = 100, nrow = 100) set.seed(0) values(r1_d1) <- sample(c(1:100), size = 1000, replace = TRUE) r2_d1 <- rast(ncol = 100, nrow = 100) set.seed(0) values(r2_d1) <- sample(c(1:1000), size = 1000, replace = TRUE) r1_d2 <- r1_d1 values(r1_d2) <- values(r1_d1)+10 r2_d2 <- r2_d1 values(r2_d2) <- values(r2_d1)+10 d1 <- c(r1_d1, r2_d1) names(d1) <- c("valsr1", "valsr2") d2 <- c(r1_d2, r2_d2) names(d2) <- c("valsr1", "valsr2") envlist <- list(d1, d2) # creating occs and bgs set.seed(0) occs <- sample(c(1:nrow(crds(envlist[[1]][[1]]))), size = 50, replace = FALSE) bgs <- sample(c(1:nrow(crds(envlist[[1]][[1]]))), size = 500, replace = FALSE) occs_d1 <- crds(envlist[[1]][[1]])[occs[1:25],] occs_d2 <- crds(envlist[[2]][[1]])[occs[26:50],] bg_d1 <- crds(envlist[[1]][[1]])[bgs[1:250],] bg_d2 <- crds(envlist[[1]][[1]])[bgs[251:500],] # extracting at occs and bgs occ_valsr1_d1 <- extract(envlist[[1]][[1]], occs_d1) occ_valsr1_d2 <- extract(envlist[[2]][[1]], occs_d2) occ_valsr2_d1 <- extract(envlist[[1]][[2]], occs_d1) occ_valsr2_d2 <- extract(envlist[[2]][[2]], occs_d2) occ_valsr1 <- rbind(occ_valsr1_d1, occ_valsr1_d2) occ_valsr2 <- rbind(occ_valsr2_d1, occ_valsr2_d2) bg_valsr1_d1 <- extract(envlist[[1]][[1]], bg_d1) bg_valsr1_d2 <- extract(envlist[[2]][[1]], bg_d2) bg_valsr2_d1 <- extract(envlist[[1]][[2]], bg_d1) bg_valsr2_d2 <- extract(envlist[[2]][[2]], bg_d2) bg_valsr1 <- rbind(bg_valsr1_d1, bg_valsr1_d2) bg_valsr2 <- rbind(bg_valsr2_d1, bg_valsr2_d2) valsr1 <- rbind(occ_valsr1, bg_valsr1) valsr2 <- rbind(occ_valsr2, bg_valsr2) p1 <- rep(1, times = 50) p0 <- rep(0, times = 500) p <- c(p1, p0) maxdf <- data.frame(p, valsr1, valsr2) # creating coord_df coord_df <- rbind(occs_d1, occs_d2, bg_d1, bg_d2) z <- c(rep(1, times = 25), rep(2, times = 25), rep(1, times = 250), rep(2, times = 250)) coord_df <- cbind(coord_df, z) colnames(coord_df) <- c("longitude", "latitude", "depth") # creating suitability rasters suitd1 <- d1[[1]] values(suitd1) <- runif(min = 0, max = 1, n = length(values(suitd1))) suitd2 <- d2[[1]] values(suitd2) <- runif(min = 0, max = 1, n = length(values(suitd2))) suit <- c(suitd1, suitd2) # here's the function result <- threshold_3D(predicted_layers = suit, thresholding_vals = c(0.9, 0.95), maxent_df = maxdf, coord_df = coord_df, weights = 2/3)
Plots cell values along a vertical transect
transectPlot( rast = NULL, sampleAxis = "lon", axisValue = NA, scaleRange = NA, plotLegend = TRUE, depthLim = as.numeric(max(names(rast))), transRange = c(-90, 90), transTicks = 20, verbose = FALSE, ... )transectPlot( rast = NULL, sampleAxis = "lon", axisValue = NA, scaleRange = NA, plotLegend = TRUE, depthLim = as.numeric(max(names(rast))), transRange = c(-90, 90), transTicks = 20, verbose = FALSE, ... )
rast |
A multilayer |
sampleAxis |
Specifies whether a latitudinal ("lat") or longitudinal ("long") transect is desired. |
axisValue |
Numeric value specifying transect postion. |
scaleRange |
A numeric vector of length 2, specifying the range that should be used for the plot color scale. |
plotLegend |
|
depthLim |
A single vector of class |
transRange |
A |
transTicks |
|
verbose |
|
... |
Additional optional arguments to pass to |
A ggplot showing a vertical slice through the SpatRaster.
Only unprojected SpatRaster files are supported.
library(terra) rast1 <- rast(ncol=10, nrow=10) values(rast1) <- rep(0:3, 50) rast2 <- rast(ncol=10, nrow=10) values(rast2) <- c(rep(0, 50), rep(1,25), rep(2,25)) rast3 <- rast(ncol=10, nrow=10) values(rast3) <- rep(c(1,3,2,1), 25) distBrick <- c(rast1, rast2, rast3) names(distBrick) <- c(0:2) transectPlot(distBrick, depthLim = 3)library(terra) rast1 <- rast(ncol=10, nrow=10) values(rast1) <- rep(0:3, 50) rast2 <- rast(ncol=10, nrow=10) values(rast2) <- c(rep(0, 50), rep(1,25), rep(2,25)) rast3 <- rast(ncol=10, nrow=10) values(rast3) <- rep(c(1,3,2,1), 25) distBrick <- c(rast1, rast2, rast3) names(distBrick) <- c(0:2) transectPlot(distBrick, depthLim = 3)
SpatRaster vector using 3D coordinatesGets values at x,y,z occurrences from a given 3D environmental variable brick
xyzSample(occs, envBrick, verbose = TRUE)xyzSample(occs, envBrick, verbose = TRUE)
occs |
A |
envBrick |
A |
verbose |
|
The SpatRaster vector object should
have numeric names that correspond with the beginning
depth of a particular depth slice. For example, one
might have three layers, one from 0 to 10m, one from
10 to 30m, and one from 30 to 100m. You would name the
layers in this brick names(envBrick) <- c(0, 10, 30.
xyzSample identifies the layer name that is closest
to the depth layer value at a particular X, Y
coordinate, and samples the environmental value at that
3D coordinate.
Vector of environmental values equal in length
to number of rows of input occs data.frame.
library(terra) # Create test raster r1 <- rast(ncol=10, nrow=10) values(r1) <- 1:100 r2 <- rast(ncol=10, nrow=10) values(r2) <- c(rep(20, times = 50), rep(60, times = 50)) r3 <- rast(ncol=10, nrow=10) values(r3) <- 8 envBrick <- c(r1, r2, r3) names(envBrick) <- c(0, 10, 30) # Create test occurrences set.seed(0) longitude <- sample(ext(envBrick)[1]:ext(envBrick)[2], size = 10, replace = FALSE) set.seed(0) latitude <- sample(ext(envBrick)[3]:ext(envBrick)[4], size = 10, replace = FALSE) set.seed(0) depth <- sample(0:35, size = 10, replace = TRUE) occurrences <- as.data.frame(cbind(longitude,latitude,depth)) # Test function occSample3d <- xyzSample(occurrences, envBrick) # How to use occurrences$envtValue <- occSample3dlibrary(terra) # Create test raster r1 <- rast(ncol=10, nrow=10) values(r1) <- 1:100 r2 <- rast(ncol=10, nrow=10) values(r2) <- c(rep(20, times = 50), rep(60, times = 50)) r3 <- rast(ncol=10, nrow=10) values(r3) <- 8 envBrick <- c(r1, r2, r3) names(envBrick) <- c(0, 10, 30) # Create test occurrences set.seed(0) longitude <- sample(ext(envBrick)[1]:ext(envBrick)[2], size = 10, replace = FALSE) set.seed(0) latitude <- sample(ext(envBrick)[3]:ext(envBrick)[4], size = 10, replace = FALSE) set.seed(0) depth <- sample(0:35, size = 10, replace = TRUE) occurrences <- as.data.frame(cbind(longitude,latitude,depth)) # Test function occSample3d <- xyzSample(occurrences, envBrick) # How to use occurrences$envtValue <- occSample3d