Introduction to distance-to

Overview

The distanceto package is designed to quickly sample distances from points features to other vector layers. Normally the approach for calculating distance to (something) involves generating distance surfaces using raster based approaches eg. raster::distance or gdal_proximity and subsequently point sampling these surfaces. Since raster based approaches are a costly method that frequently leads to memory issues or long and slow run times with high resolution data or large study sites, we have opted to compute these distances using vector based approaches. As a helper, there’s a decidedly low-res raster based approach for visually inspecting your region’s distance surface. But the workhorse is distance_to.

The distanceto package provides two functions:

Install

# Enable the robitalec universe
options(repos = c(
    robitalec = 'https://robitalec.r-universe.dev',
    CRAN = 'https://cloud.r-project.org'))

# Install distanceto
install.packages('distanceto')

distance_to

library(distanceto)
library(sf)
#> Linking to GEOS 3.9.1, GDAL 3.3.1, PROJ 8.0.1

Long-lat / unprojected coordinates

# Load nc data
nc <- st_read(system.file("shapes/sids.shp", package="spData"))
#> Reading layer `sids' from data source 
#>   `/home/alecr/R/x86_64-pc-linux-gnu-library/4.1/spData/shapes/sids.shp' 
#>   using driver `ESRI Shapefile'
#> Simple feature collection with 100 features and 22 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
#> CRS:           NA
st_crs(nc) <- "+proj=longlat +datum=NAD27"

# Set number of sampling points
npts <- 1e3

# Sample points in nc
ncpts <- st_sample(nc, npts)

# Select first 5 of nc
ncsub <- nc[1:5,]

# Measure distance from ncpts to first 5 of nc
dist <- distance_to(ncpts, ncsub, measure = 'geodesic')

# or add to ncpts
ncpts$dist <- dist

head(dist, 30)
#>  [1]      0.000 100251.919 111498.893   2973.749 119187.050  88176.616
#>  [7] 104028.802 142176.968 238416.936  73387.072  69765.270  51110.931
#> [13] 187526.880 166932.029 104379.414 109645.455  68473.130  46627.667
#> [19]  83015.046  60307.908  77483.893  40163.755  98120.287 140469.537
#> [25]  60011.006      0.000 175616.894  82308.304 138212.896  93276.728
hist(dist)

Projected coordinates

# Load seine data
data('seine', package = 'spData')

# Buffer seine by 1000 metres
bufseine <- st_buffer(seine, 1000)

# Set number of sampling points
npts <- 1e2

# Sample points within buffer seine
seinepts <- st_sample(bufseine, npts)

# Measure distance from seine points to seine
dist <- distance_to(seinepts, seine)

# or add to seine points
seinepts$dist <- dist

head(dist, 30)
#>  [1]  729.9492 1295.7741 1654.6630 1202.7673  433.3607  505.3997 1027.4409
#>  [8]  558.5800  713.9796  931.8196  440.2434 2055.2045 1153.4170 4546.5534
#> [15]  849.2587 2269.7334 1013.3373  955.0065 1134.6948  444.4242 1446.7029
#> [22]  693.0760  654.7551  932.4086  566.9460 1219.2643 1029.1893  337.6810
#> [29]  623.9310 1220.7289
hist(dist)

distance_raster

library(raster)
#> Loading required package: sp

rdist <- distance_raster(seine, 1e4)
#> Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
#> = prefer_proj): Discarded datum Unknown based on GRS80 ellipsoid in Proj4
#> definition

plot(rdist)