I am trying to apply the generalized random tessellated sampling (GRTS) algorithm implemented in spsurvey to sample a dataset of 3 million points on a map. I am running into many vector memory exhausted errors and occasionally, on other attempts, R simply crashes without warning.
In this MWE with 1,000,000 points, you can reproduce the error and see that editing the R environment file doesn’t seem to address the problem. The objective is to sample down to 250,000 points or fewer depending on what the dataset supports given the minimum 3 km distance between points constraint. Is there another package that has the same funtionality but is more efficient than spsurvey? For example, it doesn’t appear that terra has a GRTS implementation, although often seems to be more efficient than sf.
library(usethis)
edit_r_environ()
# write following line in r environ file and restart:
# R_MAX_VSIZE=100Gb
library(spsurvey)
mwe <- data.frame(geometry = paste0("POINT (", 3000001:4000000, " ", -1500001:-2500000, ")"),
property = rnorm(1000000)) %>%
st_as_sf(wkt = "geometry", crs = 5070)
samples <- grts(mwe, n_base=75000, mindis=3000)
# Error: vector memory exhausted (limit reached?)
Update: turns out that 100 Gb is not nearly enough memory to sample this many points with this implementation of the GRTS algorithm. Sampling ~150,000 points took a little over 330 Gb of RAM. If there are more efficient ways of achieving the same objective, I would very much like to know!
2
I’m not familiar with grts sampling, but I did find an implementation using grtsdb https://inbo.github.io/grtsdb/index.html that might help you.
It seems to work, and relatively quickly on the sample data you provided. You’ll need to double check the sampling results and see if they meet your requirements. An (extremely) simple check was run in the code below.
library(dplyr) # for %>% pipe
library(sf)
library(grtsdb) # https://inbo.github.io/grtsdb/index.html
library(DBI)
library(RSQLite)
library(ggplot2)
# your example data from question
mwe <- data.frame(geometry = paste0("POINT (", 3000001:4000000, " ", -1500001:-2500000, ")"),
property = rnorm(1000000)) %>%
st_as_sf(wkt = "geometry", crs = 5070)
# write mwe to sqlite database
mwedb <- dbConnect(RSQLite::SQLite(), "mwedb.sqlite")
dbWriteTable(mwedb, 'mwe', mwe)
mwe_dbfile <- 'mwedb.sqlite'
db <- connect_db(mwe_dbfile)
# grtsdb requires a bounding box as an argument
my_bbox <- rbind(
c(3000001, 4000000),
c(-2500000, -1500001)
)
# use extract_sample from grtsdb
ex_sample <- extract_sample(grtsdb = db,
samplesize = 75000,
bbox = my_bbox,
cellsize = 3000)
# size & shape of ex_sample
dim(ex_sample)
#[1] 75000 3
head(ex_sample)
# x1c x2c ranking
#1 3705500 -2154500 1
#2 3828500 -1773500 2
#3 3216500 -2460500 3
#4 3054500 -2277500 7
#5 3555500 -1524500 10
#6 3261500 -2106500 11
# check distances from at least one point to all others,
# here we just use the first point in ex_sample
# to all others, the shortest is 3000
ex_sample_sf <- st_as_sf(ex_sample, coords = c("x1c", "x2c"))
near_points <- st_nearest_points(ex_sample_sf[1,],
ex_sample_sf[-1,]) %>%
st_as_sf() %>%
mutate(length = st_length(.)) %>%
arrange(length)
head(near_points, 2)
#Simple feature collection with 2 features and 1 field
#Geometry type: LINESTRING
#Dimension: XY
#Bounding box: xmin: 3705500 ymin: -2154500 xmax: 3708500 ymax: -2151500
#CRS: NA
# length x
#1 3000.000 LINESTRING (3705500 -215450...
#2 4242.641 LINESTRING (3705500 -215450...
# Plotting
ggplot() +
geom_sf(data = ex_sample_sf, size = .1, alpha = .1) +
geom_sf(data = near_points[1:2,], color = 'red') +
geom_sf(data = ex_sample_sf[1,], color = 'blue') +
ggtitle("all sampled points, first point in blue")
ggplot() +
geom_sf(data = ex_sample_sf, size = 1, color = 'black') +
geom_sf(data = near_points[1:2,], color = 'red') +
geom_sf(data = ex_sample_sf[1,], color = 'blue') +
coord_sf(xlim = c(3700000, 3710000), ylim = c(-2150000, -2160000)) +
ggtitle("zoom in on first point & 2 nearest")
Showing all 75,000 sampled points, first one in blue. Lines are drawn to its two nearest neighbors, but they don’t show at this level of detail. All points look pretty well spaced out.
Here zoomed in on the first point, with the lines to the two nearest neighbors visible. The closest is 3000 (units?)
Hope this helps.
2