| Title: | Occurrence Filtering, Geographic Standardization and Generation of Species Range Polygons |
|---|---|
| Description: | Provides tools for filtering occurrence records, generating alpha-hull-derived range polygons and mapping species distributions. |
| Authors: | Pascal Title [aut, cre] |
| Maintainer: | Pascal Title <[email protected]> |
| License: | GPL (>= 3) |
| Version: | 2.2 |
| Built: | 2024-10-31 13:42:42 UTC |
| Source: | https://github.com/ptitle/rangebuilder |
Adds a legend to an existing raster plot, with some additional manual control
addRasterLegend( r, direction, side, location = "right", nTicks = 2, adj = NULL, shortFrac = 0.02, longFrac = 0.3, axisOffset = 0, border = TRUE, ramp = "terrain", isInteger = "auto", ncolors = 64, breaks = NULL, minmax = NULL, locs = NULL, cex.axis = 0.8, labelDist = 0.7, digits = 2, bigmark = "", ... )addRasterLegend( r, direction, side, location = "right", nTicks = 2, adj = NULL, shortFrac = 0.02, longFrac = 0.3, axisOffset = 0, border = TRUE, ramp = "terrain", isInteger = "auto", ncolors = 64, breaks = NULL, minmax = NULL, locs = NULL, cex.axis = 0.8, labelDist = 0.7, digits = 2, bigmark = "", ... )
r |
the rasterLayer object that has been plotted |
direction |
direction of color ramp. If omitted, then direction is automatically
inferred, otherwise can be specified as |
side |
side for tick marks, see |
location |
either a location name (see |
nTicks |
number of tick marks, besides min and max. |
adj |
if location is top, left, bottom or right, use this argument to adjust the location of the legend, defined in percent of the figure width. See Details for additional information. |
shortFrac |
Percent of the plot width range that will be used as the short dimention of the legend. Only applies to preset location options. |
longFrac |
Percent of the plot width range that will be used as the long dimention of the legend. Only applies to preset location options. |
axisOffset |
distance from color bar for labels, as a percent of the plot width. |
border |
logical, should the color legend have a black border |
ramp |
either a vector of color names for defining the color ramp, or "terrain" (default raster behavior) |
isInteger |
If |
ncolors |
grain size of color ramp |
breaks |
If a custom set of color breaks were used in plotting the raster, pass those color breaks here. This overrides the minmax option. |
minmax |
min and max values from which the color ramp will be derived. If left
as |
locs |
locations of tick marks, if |
cex.axis |
size of axis labels |
labelDist |
distance from axis to axis labels (passed to |
digits |
number of decimal places for labels |
bigmark |
character used to separate thousands and millions, passed
to |
... |
additional parameters to be passed to |
A number of predefined locations exist in this function to make it easy to
add a legend to a raster plot.
Preset locations are: topleft, topright, bottomleft,
bottomright, left, right, top and bottom.
If more fine-tuned control is desired, then a numeric vector of length 4 can be
supplied to location, specifying the min x, max x, min y and max y
values for the legend.
Additionally, the adj argument can be used to more intuitively adjust where
the legend is placed. adj is defined as a percentage of the figure width or
height, left to right, or bottom to top, respectively. For example, if the legend is at
the bottom, adj = 0.8 will place the legend 80
the figure, horizontally centered.
See examples.
Invisibly returns a list with the following components.
coords: 2-column matrix of xy coordinates for each color bin in the legend.
width: Coordinates for the short dimension of the legend.
pal: the color ramp
tickLocs: the tick mark locations in plotting units
labels: the values associated with those tick locations.
Pascal Title
Determines which country a given point falls in.
closestCountry(pt, crs = 4326)closestCountry(pt, crs = 4326)
pt |
longitude and latitude, as a numeric vector, 2-column table, or spatial points object. |
crs |
the CRS of the coordinate. If |
Based on a predetermined set of global points, this function finds the country of occurrence. This can be useful for checking the validity of a point by comparing the returned country to the country listed with the occurrence record. If a point falls close to the boundary between two countries, the names of the nearby countries are returned. This function will not be of much value if the point falls in the ocean, as it will return the country that is closest, regardless of how far away it is.
If one point is provided, a character vector is returned. If multiple points are provided, a list of character vectors is returned.
Pascal Title
#point near a country border closestCountry(c(-115.436, 32.657)) # testing different input options samp <- sample(1:nrow(crotalus), 10) xy <- crotalus[samp, c('decimallongitude', 'decimallatitude')] sfpts <- sf::st_as_sf(xy, coords = c('decimallongitude', 'decimallatitude'), crs = 4326) sfptsEA <- sf::st_transform(sfpts, crs = '+proj=eqearth') spPts <- as(sfpts, 'Spatial') closestCountry(xy) closestCountry(sfpts) closestCountry(sfptsEA) closestCountry(spPts)#point near a country border closestCountry(c(-115.436, 32.657)) # testing different input options samp <- sample(1:nrow(crotalus), 10) xy <- crotalus[samp, c('decimallongitude', 'decimallatitude')] sfpts <- sf::st_as_sf(xy, coords = c('decimallongitude', 'decimallatitude'), crs = 4326) sfptsEA <- sf::st_transform(sfpts, crs = '+proj=eqearth') spPts <- as(sfpts, 'Spatial') closestCountry(xy) closestCountry(sfpts) closestCountry(sfptsEA) closestCountry(spPts)
Calculates the potential error in coordinates due to lack of coordinate precision.
coordError(coords, nthreads = 1)coordError(coords, nthreads = 1)
coords |
longitude and latitude in decimal degrees, either as a long/lat vector, or as a 2-column table. Can be either as numeric or character format |
nthreads |
number of threads to use for parallelization of the function.
The R package |
This function assumes that the true precision of the coordinates is equivalent to the
greatest number of decimals in either the longitude or latitude that are not trailing
zeroes. In other words: (-130.45670, 45.53000) is interpreted as (-130.4567, 45.5300) (-130.20000, 45.50000) is interpreted as (-130.2, 45.5)
If we use (-130.45670, 45.53000) as an example, these coordinates are interpreted
as (-130.4567, 45.5300) and the greatest possible error is inferred as two
endpoints:
(-130.45670, 45.53000) and (-130.45679, 45.53009)
The distance between these two is then calculated and returned.
Returns a vector of coordinate error in meters.
Pascal Title
data(crotalus) xy <- crotalus[1:100, c('decimallongitude','decimallatitude')] coordError(xy)data(crotalus) xy <- crotalus[1:100, c('decimallongitude','decimallatitude')] coordError(xy)
Identifies occurrence records that do not occur on land.
filterByLand(coords, crs = 4326)filterByLand(coords, crs = 4326)
coords |
coordinates in the form of a 2 column numeric matrix, data.frame, numeric vector, or spatial points object (sf or sp). If spatial object, crs must be defined. |
crs |
crs of input coords. Ignored if input coords are spatial object. |
This function uses a rasterized version of the GSHHG (global self-consistent, hierarchical, high-resolution geography database, https://www.soest.hawaii.edu/pwessel/gshhg/), that has been buffered by 2 km.
returns a logical vector where TRUE means the point falls on
land.
Pascal Title
data(crotalus) #identify points that fall off land filterByLand(crotalus[,c('decimallongitude','decimallatitude')]) # testing different input options samp <- sample(1:nrow(crotalus), 10) xy <- crotalus[samp, c('decimallongitude', 'decimallatitude')] sfpts <- sf::st_as_sf(xy, coords = c('decimallongitude', 'decimallatitude'), crs = 4326) sfptsEA <- sf::st_transform(sfpts, crs = '+proj=eqearth') spPts <- as(sfpts, 'Spatial') filterByLand(xy) filterByLand(sfpts) filterByLand(sfptsEA) filterByLand(spPts)data(crotalus) #identify points that fall off land filterByLand(crotalus[,c('decimallongitude','decimallatitude')]) # testing different input options samp <- sample(1:nrow(crotalus), 10) xy <- crotalus[samp, c('decimallongitude', 'decimallatitude')] sfpts <- sf::st_as_sf(xy, coords = c('decimallongitude', 'decimallatitude'), crs = 4326) sfptsEA <- sf::st_transform(sfpts, crs = '+proj=eqearth') spPts <- as(sfpts, 'Spatial') filterByLand(xy) filterByLand(sfpts) filterByLand(sfptsEA) filterByLand(spPts)
Filter occurrence records by their proximity to each other.
filterByProximity(xy, dist, returnIndex = FALSE)filterByProximity(xy, dist, returnIndex = FALSE)
xy |
longitude and latitude in decimal degrees, either as a matrix, dataframe, or spatial points object. |
dist |
minimum allowed distance in km |
returnIndex |
if |
This function will discard coordinates that fall within a certain distance from other points.
If returnIndex = TRUE, returns a numeric vector of indices.
If returnIndex = FALSE, returns coordinates of the same class as the
input.
Pascal Title
data(crotalus) # within the first 100 points in the dataset, identify the set of points to # drop in order to have points no closer to each other than 20 km subset <- crotalus[1:100,] tooClose <- filterByProximity(xy= subset[ ,c('decimallongitude','decimallatitude')], dist=20, returnIndex = TRUE) plot(subset[ ,c('decimallongitude','decimallatitude')], pch=1, col='blue', cex=1.5) points(subset[tooClose, c('decimallongitude','decimallatitude')], pch=20, col='red') # testing different input options samp <- sample(1:nrow(crotalus), 100) xy <- crotalus[samp, c('decimallongitude', 'decimallatitude')] sfpts <- sf::st_as_sf(xy, coords = c('decimallongitude', 'decimallatitude'), crs = 4326) sfptsEA <- sf::st_transform(sfpts, crs = '+proj=eqearth') spPts <- as(sfpts, 'Spatial') filterByProximity(xy, dist=20, returnIndex = TRUE) filterByProximity(sfpts, dist=20, returnIndex = TRUE) filterByProximity(sfptsEA, dist=20, returnIndex = TRUE) filterByProximity(spPts, dist=20, returnIndex = TRUE)data(crotalus) # within the first 100 points in the dataset, identify the set of points to # drop in order to have points no closer to each other than 20 km subset <- crotalus[1:100,] tooClose <- filterByProximity(xy= subset[ ,c('decimallongitude','decimallatitude')], dist=20, returnIndex = TRUE) plot(subset[ ,c('decimallongitude','decimallatitude')], pch=1, col='blue', cex=1.5) points(subset[tooClose, c('decimallongitude','decimallatitude')], pch=20, col='red') # testing different input options samp <- sample(1:nrow(crotalus), 100) xy <- crotalus[samp, c('decimallongitude', 'decimallatitude')] sfpts <- sf::st_as_sf(xy, coords = c('decimallongitude', 'decimallatitude'), crs = 4326) sfptsEA <- sf::st_transform(sfpts, crs = '+proj=eqearth') spPts <- as(sfpts, 'Spatial') filterByProximity(xy, dist=20, returnIndex = TRUE) filterByProximity(sfpts, dist=20, returnIndex = TRUE) filterByProximity(sfptsEA, dist=20, returnIndex = TRUE) filterByProximity(spPts, dist=20, returnIndex = TRUE)
Checks for coordinate sign mistakes by checking all possibilities against country occupancy.
flipSign( coordVec, country, returnMultiple = FALSE, filterByLand = TRUE, crs = 4326 )flipSign( coordVec, country, returnMultiple = FALSE, filterByLand = TRUE, crs = 4326 )
coordVec |
numeric vector of length 2: longitude, latitude |
country |
the country that is associated with the record |
returnMultiple |
if multiple sign flips lead to the correct country,
return all options. If |
filterByLand |
if |
crs |
the crs of the coordinate. |
This function generates all possible coordinates with different signs, and
runs closestCountry on each, returning the coordinates that
lead to a country match. It ignores coordinate options that do not pass
filterByLand.
If a point falls close to the boundary between two countries, it is still considered a match.
list with 2 elements
matched |
logical: Was the country matched |
newcoords |
matrix of coordinates that were successful. |
Pascal Title
#correct coordinates flipSign(c(4.28, 39.98), country = 'Spain') #mistake in coordinate sign flipSign(c(115.436, 32.657), country = 'United States') #incorrect sign on both long and lat, but not possible to distinguish for longitude #except when we consider which alternative coords fall on land. flipSign(c(-4.28, -39.98), country = 'Spain', filterByLand = FALSE, returnMultiple = TRUE) flipSign(c(-4.28, -39.98), country = 'Spain', returnMultiple = TRUE) #coordinates are incorrect flipSign(c(4.28, 59.98), country = 'Spain')#correct coordinates flipSign(c(4.28, 39.98), country = 'Spain') #mistake in coordinate sign flipSign(c(115.436, 32.657), country = 'United States') #incorrect sign on both long and lat, but not possible to distinguish for longitude #except when we consider which alternative coords fall on land. flipSign(c(-4.28, -39.98), country = 'Spain', filterByLand = FALSE, returnMultiple = TRUE) flipSign(c(-4.28, -39.98), country = 'Spain', returnMultiple = TRUE) #coordinates are incorrect flipSign(c(4.28, 59.98), country = 'Spain')
Generates an apha hull polygon, where the alpha parameter is determined by the spatial distribution of the coordinates.
getDynamicAlphaHull( x, fraction = 0.95, partCount = 3, buff = 10000, initialAlpha = 3, coordHeaders = c("Longitude", "Latitude"), clipToCoast = "terrestrial", alphaIncrement = 1, verbose = FALSE, alphaCap = 400 )getDynamicAlphaHull( x, fraction = 0.95, partCount = 3, buff = 10000, initialAlpha = 3, coordHeaders = c("Longitude", "Latitude"), clipToCoast = "terrestrial", alphaIncrement = 1, verbose = FALSE, alphaCap = 400 )
x |
dataframe of coordinates in decimal degrees, with a minimum of 3 rows. |
fraction |
the minimum fraction of occurrences that must be included in polygon. |
partCount |
the maximum number of disjunct polygons that are allowed. |
buff |
buffering distance in meters |
initialAlpha |
the starting value for alpha |
coordHeaders |
the column names for the longitude and latitude
columns, respectively. If x has two columns, these are assumed to be
longitude and latitude, and |
clipToCoast |
Either "no" (no clipping), "terrestrial" (only terrestrial part of range is kept) or "aquatic" (only non-terrestrial part is clipped). See Details. |
alphaIncrement |
the amount to increase alpha with each iteration |
verbose |
prints the alpha value to the console, intended for debugging. |
alphaCap |
Max alpha value before function aborts and returns a minimum convex hull. |
From a set of coordinates, this function will create an alpha hull with
alpha = initialAlpha, and will then increase alpha by
alphaIncrement until both the fraction and partCount
conditions are met.
If the conditions cannot be satisfied, then a minimum convex hull is returned.
If clipToCoast is set to "terrestrial" or "aquatic", the resulting
polygon is clipped to the coastline, using a basemap from naturalearth.
The first time this function is run, this basemap will be downloaded.
Subsequent calls will use the downloaded map.
a list with 2 elements:
hull |
a sf polygon object |
alpha |
the alpha value that was found to satisfy the criteria. If a convex hull was returned, this will list MCH. |
Pascal Title
Alpha hulls are created with ahull.
data(crotalus) # create a polygon range for Crotalus atrox x <- crotalus[which(crotalus$genSp == 'Crotalus_atrox'),] x <- x[sample(1:nrow(x), 50),] range <- getDynamicAlphaHull(x, coordHeaders=c('decimallongitude','decimallatitude'), clipToCoast = 'no') plot(range[[1]], col=transparentColor('dark green', 0.5), border = NA) points(x[,c('decimallongitude','decimallatitude')], cex = 0.5, pch = 3) # to add a basic coastline, you can use the internal map # world <- rangeBuilder:::loadWorldMap() # plot(world, add = TRUE, lwd = 0.5)data(crotalus) # create a polygon range for Crotalus atrox x <- crotalus[which(crotalus$genSp == 'Crotalus_atrox'),] x <- x[sample(1:nrow(x), 50),] range <- getDynamicAlphaHull(x, coordHeaders=c('decimallongitude','decimallatitude'), clipToCoast = 'no') plot(range[[1]], col=transparentColor('dark green', 0.5), border = NA) points(x[,c('decimallongitude','decimallatitude')], cex = 0.5, pch = 3) # to add a basic coastline, you can use the internal map # world <- rangeBuilder:::loadWorldMap() # plot(world, add = TRUE, lwd = 0.5)
Given a list of SpatialPolygons or sf objects, return a bounding box object that encompasses all items.
getExtentOfList(shapes)getExtentOfList(shapes)
shapes |
a list of SpatialPolygons or simple features |
An object of class bbox.
Pascal Title
data(crotalus) # create some polygons, in this case convex hulls sp <- split(crotalus, crotalus$genSp) sp <- lapply(sp, function(x) x[,c('decimallongitude','decimallatitude')]) sp <- lapply(sp, function(x) x[chull(x),]) poly <- lapply(sp, function(x) sf::st_convex_hull(sf::st_combine(sf::st_as_sf(x, coords = 1:2, crs = 4326)))) getExtentOfList(poly)data(crotalus) # create some polygons, in this case convex hulls sp <- split(crotalus, crotalus$genSp) sp <- lapply(sp, function(x) x[,c('decimallongitude','decimallatitude')]) sp <- lapply(sp, function(x) x[chull(x),]) poly <- lapply(sp, function(x) sf::st_convex_hull(sf::st_combine(sf::st_as_sf(x, coords = 1:2, crs = 4326)))) getExtentOfList(poly)
Provides tools for filtering occurrence records, standardizing countries names, generating alpha-hull-derived range polygons and mapping species distributions.
Pascal Title <[email protected]>
Davis Rabosky, A.R., C.L. Cox, D.L. Rabosky, P.O. Title, I.A. Holmes, A. Feldman and J.A. McGuire. 2016. Coral snakes predict the evolution of mimicry across New World snakes. Nature Communications 7:11484.
Useful links:
Included datasets in rangeBuilder
The crotalus dataset is the result of a query for genus Crotalus on the VertNet
search portal (http://portal.vertnet.org/search), and has been thinned and lightly
filtered, to serve as an example dataset for this package.
Takes a list of polygons and creates a multi-layer SpatRaster.
rasterStackFromPolyList( polyList, resolution = 50000, retainSmallRanges = TRUE, extent = "auto" )rasterStackFromPolyList( polyList, resolution = 50000, retainSmallRanges = TRUE, extent = "auto" )
polyList |
a list of spatial polygon objects, named with taxon names. It is assumed that all items in last have same crs. |
resolution |
vertical and horizontal size of raster cell, in units of the polygons' projection |
retainSmallRanges |
boolean; should small ranged species be dropped or preserved. See details. |
extent |
if 'auto', then the maximal extent of the polygons will be used. If not auto, must be a numeric vector of length 4 with minLong, maxLong, minLat, maxLat. |
In the rasterization process, all cells for which the polygon covers the
midpoint are considered as present and receive a value of 1. If
retainSmallRanges = FALSE, then species whose ranges are so small
that no cell registers as present will be dropped. If
retainSmallRanges = TRUE, then the cells that the small polygon is
found in will be considered as present.
an object of class SpatRaster where all rasters contain
values of either NA or 1.
Pascal Title
## Not run: data(crotalus) library(sf) library(terra) # get 10 species occurrence sets uniqueSp <- split(crotalus, crotalus$genSp) uniqueSp <- lapply(uniqueSp, function(x) x[!duplicated(x[, c('decimallongitude', 'decimallatitude')]),]) uniqueSp <- names(uniqueSp[sapply(uniqueSp, nrow) > 5]) uniqueSp <- uniqueSp[1:10] # create range polygons ranges <- vector('list', length = length(uniqueSp)) for (i in 1:length(uniqueSp)) { x <- crotalus[which(crotalus$genSp == uniqueSp[i]),] ranges[[i]] <- getDynamicAlphaHull(x, coordHeaders = c('decimallongitude', 'decimallatitude'), clipToCoast = 'terrestrial') } # name the polygons names(ranges) <- uniqueSp # keep only the polygons ranges <- lapply(ranges, function(x) x[[1]]) # Create a SpatRaster with the extent inferred from the polygons, and a cell # resolution of 0.2 degrees. # cells with the presence of a species get a value of 1, NA if absent. rangeStack <- rasterStackFromPolyList(ranges, resolution = 0.2) # calculate species richness per cell, where cell values are counts of species richnessRaster <- app(rangeStack, fun=sum, na.rm = TRUE) # set values of 0 to NA richnessRaster[richnessRaster == 0] <- NA #plot ramp <- colorRampPalette(c('blue','yellow','red')) plot(richnessRaster, col=ramp(100)) # to add a basic coastline, you can use the internal map # world <- rangeBuilder:::loadWorldMap() # plot(world, add = TRUE, lwd = 0.5) ## End(Not run)## Not run: data(crotalus) library(sf) library(terra) # get 10 species occurrence sets uniqueSp <- split(crotalus, crotalus$genSp) uniqueSp <- lapply(uniqueSp, function(x) x[!duplicated(x[, c('decimallongitude', 'decimallatitude')]),]) uniqueSp <- names(uniqueSp[sapply(uniqueSp, nrow) > 5]) uniqueSp <- uniqueSp[1:10] # create range polygons ranges <- vector('list', length = length(uniqueSp)) for (i in 1:length(uniqueSp)) { x <- crotalus[which(crotalus$genSp == uniqueSp[i]),] ranges[[i]] <- getDynamicAlphaHull(x, coordHeaders = c('decimallongitude', 'decimallatitude'), clipToCoast = 'terrestrial') } # name the polygons names(ranges) <- uniqueSp # keep only the polygons ranges <- lapply(ranges, function(x) x[[1]]) # Create a SpatRaster with the extent inferred from the polygons, and a cell # resolution of 0.2 degrees. # cells with the presence of a species get a value of 1, NA if absent. rangeStack <- rasterStackFromPolyList(ranges, resolution = 0.2) # calculate species richness per cell, where cell values are counts of species richnessRaster <- app(rangeStack, fun=sum, na.rm = TRUE) # set values of 0 to NA richnessRaster[richnessRaster == 0] <- NA #plot ramp <- colorRampPalette(c('blue','yellow','red')) plot(richnessRaster, col=ramp(100)) # to add a basic coastline, you can use the internal map # world <- rangeBuilder:::loadWorldMap() # plot(world, add = TRUE, lwd = 0.5) ## End(Not run)
Standardizes country names to the list of countries used internally by this package.
standardizeCountry(country, fuzzyDist = 1, nthreads = 1, progressBar = TRUE)standardizeCountry(country, fuzzyDist = 1, nthreads = 1, progressBar = TRUE)
country |
character vector of country names or ISO codes |
fuzzyDist |
for fuzzy searching, the maximum string distance allowed for a match; if 0, fuzzy searching is disabled. |
nthreads |
number of threads to use for parallelization of the
function. The R package |
progressBar |
if |
This package interacts with data from the Global Invasive Species Database
(GISD), the Reptile Database, as well as global maps that were used to
generate the internal dataset used by closestCountry. Efforts
have been made to make country names consistent across these separate
datasets. This function can be used to convert the user's Country
field to the same standardized set.
Fuzzy matching uses the function adist.
Parallelization with nthreads becomes more time-efficient only if
the input vector is of multiple thousands of country names.
Character vector of the standardized country names. If no match
found, "" is returned.
Pascal Title
standardizeCountry(c("Russian Federation", "USA", "Plurinational State of Bolivia", "Brezil"))standardizeCountry(c("Russian Federation", "USA", "Plurinational State of Bolivia", "Brezil"))
Converts a named color and opacity and returns the proper RGB code.
transparentColor(namedColor, alpha = 0.8)transparentColor(namedColor, alpha = 0.8)
namedColor |
a color name |
alpha |
a transparency value between 0 and 1, where 0 is fully transparent |
Returns the transparent color in RGB format.
Pascal Title