# load bivalve occurrences to rasterise
library(terra)
data(bivalves)
# initialise Equal Earth projected coordinates
rWorld <- rast()
prj <- 'EPSG:8857'
rPrj <- project(rWorld, prj, res = 200000) # 200,000m is approximately 2 degrees
# coordinate column names for the current and target coordinate reference system
xyCartes <- c('paleolng','paleolat')
xyCell <- c('centroidX','centroidY')
# project occurrences and retrieve cell centroids in new coordinate system
llOccs <- vect(bivalves, geom = xyCartes, crs = 'epsg:4326')
prjOccs <- project(llOccs, prj)
cellIds <- cells(rPrj, prjOccs)[,'cell']
bivalves[, xyCell] <- xyFromCell(rPrj, cellIds)
# subsample 20 equal-area sites within 10-degree bands of absolute latitude
n <- 20
reps <- 100
set.seed(11)
bandAbs <- bandit(dat = bivalves, xy = xyCell,
iter = reps, nSite = n, output = 'full',
bin = 10, absLat = TRUE,
crs = prj
)
head(bandAbs[[1]]) # inspect first subsample
names(bandAbs)[1] # degree interval (absolute value) of first subsample
#> [1] "[10,20)"
unique(names(bandAbs)) # all intervals containing sufficient data
#> [1] "[10,20)" "[20,30)" "[30,40)" "[40,50)"
# note insufficient coverage to subsample at equator or above 50 degrees
# subsample 20-degree bands, where central band spans the equator
# (-10 S to 10 N latitude), as in Allen et al. (2020)
# (An alternative, finer-grain way to divide 180 degrees evenly into an
# odd number of bands would be to set 'bin' = 4.)
bandCent <- bandit(dat = bivalves, xy = xyCell,
iter = reps, nSite = n, output = 'full',
bin = 20, centr = TRUE, absLat = FALSE,
crs = prj
)
unique(names(bandCent)) # all intervals containing sufficient data
#> [1] "[-50,-30)" "[10,30)" "[30,50)"
Run the code above in your browser using DataLab