根据多边形快速过滤掉一个sf点的网格

Question

我想在美国制作网格（在 x- 和 y- 坐标的数据框意义上），或者美国地区，丢弃原始网格中超出美国边界的任何点。我有一些代码似乎可以工作，但是当我将网格增量缩小到 1 公里 (1e3) 左右时，它的速度非常慢。如何才能更快地做到这一点？也许有一种方法可以在没有 lapply 或循环的情况下构建我需要的简单特征集合，或者这可以通过 MULTIPOINT 而不是 POINT 的简单特征集合来完成。

library(sf)

crs.us.atlas = 2163 # https://epsg.io/2163

us = read_sf(paste0("/vsizip/",
    "/tmp/us.zip")) # from: https://www2.census.gov/geo/tiger/GENZ2019/shp/cb_2019_us_state_500k.zip
# Filter down to the lower 48 + DC.
us = us[us$STUSPS %in% setdiff(c(state.abb, "DC"), c("AK", "HI")),]
us = st_transform(us, crs = crs.us.atlas)
l = as.list(st_bbox(us))

increment = 1e5

g = expand.grid(
    x = seq(l$xmin, l$xmax, by = increment),
    y = seq(l$ymin, l$ymax, by = increment))

message("Running the slow part")
print(system.time(g <- g[0 < sapply(FUN = length, st_intersects(
    st_as_sfc(crs = crs.us.atlas, lapply(1 : nrow(g), function(i)
        st_point(as.numeric(g[i, c("x", "y")])))),
    us)),]))

print(nrow(g))

Answer 1

我会按如下方式解决问题。首先，加载包。 tmap仅用于地图，您可以轻松忽略

library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
library(sfheaders)
library(tmap)

下载并读入数据

temp_zip <- tempfile(fileext = ".zip")
download.file(
  "https://www2.census.gov/geo/tiger/GENZ2019/shp/cb_2019_us_state_500k.zip", 
  destfile = temp_zip
)
unzip(temp_zip, exdir = tempdir())
us <- st_read(file.path(tempdir(), "cb_2019_us_state_500k.shp"))
#> Reading layer `cb_2019_us_state_500k' from data source `C:\Users\Utente\AppData\Local\Temp\RtmpCakscO\cb_2019_us_state_500k.shp' using driver `ESRI Shapefile'
#> Simple feature collection with 56 features and 9 fields
#> geometry type:  MULTIPOLYGON
#> dimension:      XY
#> bbox:           xmin: -179.1489 ymin: -14.5487 xmax: 179.7785 ymax: 71.36516
#> geographic CRS: NAD83

向下过滤到较低的 48 + DC。

us = us[us$STUSPS %in% setdiff(c(state.abb, "DC"), c("AK", "HI")),]
us = st_transform(us, crs = 2163)

定义增量、网格和sfc对象。关键部分是创建和sfc_point对象，用于后续操作

increment = 1e5
g = expand.grid(
  x = seq(st_bbox(us)$xmin, st_bbox(us)$xmax, by = increment),
  y = seq(st_bbox(us)$ymin, st_bbox(us)$ymax, by = increment)
)
g_sfc <- sfc_point(as.matrix(g)) %>% 
  st_set_crs(2163)

查找包含在美国的点的 ID

my_ids <- unlist(st_contains(us, g_sfc))

可视化结果

tm_shape(g_sfc) + 
  tm_dots(col = "grey", size = 0.05) + 
tm_shape(g_sfc[my_ids]) + 
  tm_dots(col = "darkred", size = 0.05) + 
tm_shape(us) + 
  tm_borders(lwd = 1.3)

重复 1e3（但我不会添加任何情节，因为那将近 1300 万点）

increment = 1e3
g = expand.grid(
  x = seq(st_bbox(us)$xmin, st_bbox(us)$xmax, by = increment),
  y = seq(st_bbox(us)$ymin, st_bbox(us)$ymax, by = increment)
)

生成数据大约需要 20 秒

system.time({
  g_sfc <- sfc_point(as.matrix(g)) %>% 
    st_set_crs(2163)
})
#>    user  system elapsed 
#>   16.29    0.92   17.27

并在 80 秒内找到美国境内点的 ID。

system.time({
  my_ids <- unlist(st_contains(us, g_sfc))
})
#>    user  system elapsed 
#>   67.75    8.32   80.86

^{由 reprex package (v0.3.0)}

于 2021 年 1 月 13 日创建

如果您需要更高效的东西，我建议您 polyclip。

Answer 2

或者，您可以使用 fasterize 包以所需的分辨率围绕形状文件创建光栅网格，然后使用 raster::rasterToPoints 函数提取网格坐标。

这几乎可以立即收集 xy 位置。转换回 sf 对象大约需要 10 秒左右。

library(sf)
library(fasterize)
library(raster)

crs.us.atlas = 2163

# https://www2.census.gov/geo/tiger/GENZ2019/shp/cb_2019_us_state_500k.zip
us = read_sf(list.files(pattern='.shp$'))
us = us[us$STUSPS %in% setdiff(c(state.abb, "DC"), c("AK", "HI")),]
us = st_transform(us, crs = crs.us.atlas)

increment = 1e3

# create the empty raster to for fasterize base
usrast = raster::raster(us, resolution=rep(increment, 2))

# rasterize the shapefile
r = fasterize::fasterize(us, usrast)

# extract raster cell coordinates
coords = rasterToPoints(r)

# convert coordinates to sf
shp = coords %>%
    as.data.frame() %>%
    st_as_sf(crs=crs.us.atlas, coords=c('x', 'y'))

Answer 3

sf 有一个函数 st_make_grid 可以帮助解决这个问题（在旧版本的 sf 中，它甚至会自动裁剪成多边形），但奇怪的是，它是 quite slow 在撰写本文时。

通过使用 st_as_sf 将 g 转换为简单的点特征集合，我可以在没有额外依赖的情况下获得合理的性能：

g = st_as_sf(crs = crs.us.atlas, coords = c(1, 2), g)

然后，随着@agila使用unlist，慢的部分变成了

print(system.time(g <- g[unlist(st_intersects(us, g)),]))

increment = 1e3 在高端服务器上需要 3 分钟。

根据多边形快速过滤掉一个sf点的网格

Quickly filter down a grid of sf points according to a polygon

optimization

r

sf