Nothing
##' Generate polygon based on alpha hulls
##'
##' Generates an apha hull polygon, where the alpha parameter is determined by
##' the spatial distribution of the coordinates.
##'
##'
##' From a set of coordinates, this function will create an alpha hull with
##' \code{alpha = initialAlpha}, and will then increase \code{alpha} by
##' \code{alphaIncrement} until both the \code{fraction} and \code{partCount}
##' conditions are met.
##'
##' If the conditions cannot be satisfied, then a minimum convex hull is
##' returned.
##'
##' If \code{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.
##'
##' @param x dataframe of coordinates in decimal degrees, with a minimum of 3
##' rows.
##' @param fraction the minimum fraction of occurrences that must be included
##' in polygon.
##' @param partCount the maximum number of disjunct polygons that are allowed.
##' @param buff buffering distance in meters
##' @param initialAlpha the starting value for alpha
##' @param 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 \code{coordHeaders} is ignored.
##' @param clipToCoast Either "no" (no clipping), "terrestrial" (only
##' terrestrial part of range is kept) or "aquatic" (only non-terrestrial part
##' is clipped). See Details.
##' @param alphaIncrement the amount to increase alpha with each iteration
##' @param verbose prints the alpha value to the console, intended for
##' debugging.
##' @param alphaCap Max alpha value before function aborts and returns a
##' minimum convex hull.
##' @return a list with 2 elements: \item{hull}{ a sf polygon object }
##' \item{alpha}{ the alpha value that was found to satisfy the criteria. If a
##' convex hull was returned, this will list MCH. }
##' @author Pascal Title
##' @seealso Alpha hulls are created with \code{\link{ahull}}.
##' @examples
##'
##' 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)
##'
##' @export
getDynamicAlphaHull <- function(x, fraction = 0.95, partCount = 3, buff = 10000, initialAlpha = 3, coordHeaders = c('Longitude', 'Latitude'), clipToCoast = 'terrestrial', alphaIncrement = 1, verbose = FALSE, alphaCap = 400) {
# fraction = 0.95; partCount = 3; buff = 10000; initialAlpha = 3; coordHeaders = c('Longitude', 'Latitude'); clipToCoast = 'terrestrial'; alphaIncrement = 1; verbose = FALSE; alphaCap = 400
if (clipToCoast == FALSE) clipToCoast <- 'no'
clipToCoast <- match.arg(clipToCoast, c('no', 'terrestrial', 'aquatic'))
if (ncol(x) == 2) {
coordHeaders <- c(1,2)
}
#reduce to unique coordinates
x <- x[!duplicated(x[,coordHeaders]), coordHeaders]
x <- x[stats::complete.cases(x),]
if (nrow(x) < 3) {
stop('This function requires a minimum of 3 unique coordinates (after removal of duplicates).')
}
#Alpha hulls cannot be generated if first 3 points are linear.
while ((x[1,coordHeaders[1]] == x[2,coordHeaders[1]] & x[2, coordHeaders[1]] == x[3, coordHeaders[1]]) | (x[1,2] == x[2, coordHeaders[2]] & x[2, coordHeaders[2]] == x[3, coordHeaders[2]])) {
x <- x[sample(1:nrow(x),size = nrow(x)),]
}
#create spatial points object
x <- sf::st_as_sf(as.data.frame(x), coords = 1:2, crs = 4326)
if (nrow(x) < 3) {
stop('This function requires a minimum of 3 unique coordinates.')
}
#create alpha hull and calculate fraction of occurrences that fall within
#continue until fraction is reached
alpha <- initialAlpha
problem <- FALSE
if (verbose) cat('\talpha:', alpha, '\n')
hull <- try(alphahull::ahull(sf::st_coordinates(x), alpha = alpha), silent = TRUE)
# possible that ahull is detecting duplicate points despite previous checks
dropPt <- c()
while (inherits(hull, 'try-error') & any(grepl('duplicate points', hull))) {
ptDist <- sf::st_distance(x, x)
diag(ptDist) <- NA
units(ptDist) <- NULL
closest <- which(ptDist == min(ptDist, na.rm = TRUE), arr.ind = TRUE)
hull <- try(alphahull::ahull(sf::st_coordinates(x)[- closest[1,1]], alpha = alpha), silent = TRUE)
if (inherits(hull, 'ahull')) {
dropPt <- closest[1,1]
}
}
if (length(dropPt) > 0) {
x <- x[- dropPt,]
}
hull <- try(alphahull::ahull(sf::st_coordinates(x), alpha = alpha), silent = TRUE)
while (inherits(hull, 'try-error')) {
if (verbose) {cat('\talpha:', alpha, '\n')}
alpha <- alpha + alphaIncrement
hull <- try(alphahull::ahull(sf::st_coordinates(x), alpha = alpha), silent = TRUE)
if (alpha > alphaCap) {
problem <- TRUE
break
}
}
if (!problem) {
hull <- ah2sf(hull)
validityCheck <- function(hull) {
if (!is.null(hull) & !inherits(hull, 'try-error')) {
if (!all(sf::st_is_valid(hull))) {
TRUE
} else {
FALSE
}
} else {
FALSE
}
}
while (is.null(hull) | inherits(hull, 'try-error') | validityCheck(hull)) {
alpha <- alpha + alphaIncrement
if (verbose) {cat('\talpha:', alpha, '\n')}
hull <- try(ah2sf(alphahull::ahull(sf::st_coordinates(x), alpha = alpha)), silent = TRUE)
if (alpha > alphaCap) {
problem <- TRUE
break
}
}
#how many points are within hull?
pointWithin <- sf::st_intersects(x, hull)
alphaVal <- alpha
buffered <- FALSE
buff <- units::set_units(buff, 'm')
while (any(length(hull) > partCount, (sum(lengths(pointWithin)) / nrow(x)) < fraction, !all(sf::st_is_valid(hull)))) {
alpha <- alpha + alphaIncrement
if (verbose) {cat('\talpha:', alpha, '\n')}
hull <- try(alphahull::ahull(sf::st_coordinates(x), alpha = alpha), silent = TRUE)
while (inherits(hull, 'try-error') & alpha <= alphaCap) {
alpha <- alpha + alphaIncrement
hull <- try(alphahull::ahull(sf::st_coordinates(x),alpha = alpha), silent = TRUE)
}
if (!inherits(hull, 'try-error')) {
hull <- ah2sf(hull)
hull <- sf::st_transform(hull, crs = '+proj=eqearth')
if (all(sf::st_is_valid(hull))) {
hull <- sf::st_buffer(hull, dist = buff)
hull <- sf::st_transform(hull, crs = 4326)
buffered <- TRUE
pointWithin <- sf::st_intersects(x, hull)
}
}
alphaVal <- alpha
if (alpha > alphaCap) {
hull <- sf::st_convex_hull(x)
hull <- sf::st_transform(hull, crs = '+proj=eqearth')
hull <- sf::st_buffer(hull, dist = buff)
hull <- sf::st_transform(hull, crs = 4326)
buffered <- TRUE
alphaVal = 'MCH'
break
}
}
} else {
hull <- sf::st_convex_hull(x)
hull <- sf::st_transform(hull, crs = '+proj=eqearth')
hull <- sf::st_buffer(hull, dist = buff)
hull <- sf::st_transform(hull, crs = 4326)
buffered <- TRUE
alphaVal = 'MCH'
}
if (!buffered) {
hull <- sf::st_transform(hull, crs = '+proj=eqearth')
hull <- sf::st_buffer(hull, dist = buff)
hull <- sf::st_transform(hull, crs = 4326)
}
if (clipToCoast != 'no') {
world <- loadWorldMap()
if (clipToCoast == 'terrestrial') {
hull <- sf::st_intersection(hull, world)
} else {
hull <- sf::st_difference(hull, world)
}
}
return(list(hull, alpha = paste('alpha', alphaVal, sep = '')))
}
loadWorldMap <- function() {
# load medium-level resolution landmass vector map from natural-earth
## download if needed
worldpath <- paste0(find.package(package = 'rangeBuilder'), '/extdata/world.rds')
if (!file.exists(worldpath)) {
message('\n\tDownloading a world basemap for clipping. This only needs to happen once.\n\n')
world <- rnaturalearth::ne_download(scale = 'medium', type = 'land', category = 'physical', returnclass = 'sf')
world <- sf::st_geometry(world)
world <- sf::st_combine(world)
if (!dir.exists(paste0(find.package(package = 'rangeBuilder'), '/extdata'))) {
dir.create(paste0(find.package(package = 'rangeBuilder'), '/extdata'))
}
saveRDS(world, file = worldpath)
} else {
world <- readRDS(worldpath)
}
return(world)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.