R 中大型空间数据集的 GRTS

Question

我正在尝试应用

spsurvey

中实现的广义随机镶嵌采样 (GRTS) 算法对地图上 300 万个点的数据集进行采样。我遇到了许多

vector memory exhausted

错误，有时，在其他尝试中，R 只是在没有警告的情况下崩溃了。

在这个具有 1,000,000 点的 MWE 中，您可以重现错误并看到编辑 R 环境文件似乎并不能解决问题。目标是根据数据集支持的点采样数量减少到 250,000 个或更少，点之间的距离至少为 3 公里。是否有另一个具有相同功能但比

spsurvey

更高效的软件包？例如，

terra

似乎没有 GRTS 实现，尽管通常看起来比 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?)

更新：事实证明，100 Gb 的内存不足以通过 GRTS 算法的实现来采样这么多点。采样约 150,000 个点需要略多于 330 Gb 的 RAM。如果有更有效的方法来实现相同的目标，我非常想知道！

Answer 1

我不熟悉 grts 采样，但我确实找到了一个使用

gtrsdb

https://inbo.github.io/grtsdb/index.html 的实现，可能会对你有所帮助。

它似乎有效，并且对您提供的示例数据相对较快。您需要仔细检查采样结果，看看它们是否符合您的要求。下面的代码中运行了一个（极其）简单的检查。

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")

显示所有 75,000 个采样点，第一个为蓝色。线条被绘制到其两个最近的邻居，但它们不会以这种细节级别显示。所有点看起来都分布得很好。

这里放大了第一个点，可以看到到两个最近邻居的线。最接近的是 3000（单位？）

希望这有帮助。