r - 在大型数据集上加速 st_crop(sf 包)的方法

Question

我需要为 1 公顷的 ~ 4 mio 网格单元跨不同的 shapefile 提取信息。目前我在所有单元格的 for 循环中的每一层上使用 st_crop，但这会永远运行。我想加快使用“data.table”(DT)的方式按坐标裁剪 shapefile 的过程。让我们考虑下面的示例，我在其中寻找感兴趣区域中多边形边缘的范围:

require(sf)
require(data.table)
require(ggplot2)
require(tidyverse)

# load shapefile
nc = st_read(system.file("shape/nc.shp", package="sf"))


# Define a bounding-box that mimic a mowing-window or area of interest
bb <- st_bbox(c(xmin= -79, xmax=-78,ymin= 34.5, ymax= 35.5))


# Commute 'nc' into some sort of data.table object for fast subsetting, in preserving object's integrity (i.e. same id to all points of a given polygon)
nobs <- mapview::npts(nc,by_feature=T)
NC <- data.table::data.table(id=rep(1:nrow(nc),nobs),st_coordinates(nc)[,1:2])
head(NC)

# Compare cropping methods amon
library(microbenchmark)
x = runif(100)
test <- microbenchmark(
  crop_nc <- st_crop(nc,bb),
  crop_NC <- NC[X >= bb[1] & X < bb[3] & Y>= bb[2] & Y < bb[4]]
)  

print(test)
Unit: microseconds
  expr      min       lq      mean   median        uq       max neval cld
 crop_nc  5205.051 5675.807 6837.9472 5903.219 6829.0865 16046.654   100   b
 crop_NC   405.334  528.356  624.8398  576.996  656.9245  1295.361   100  a 
There were 50 or more warnings (use warnings() to see the first 50)

正如预期的那样，子集化的 DT 方式更快。现在让我们从我们的 DT 对象回到 sf 对象，如下所示:

crop_NC_sf <- st_as_sf(crop_NC,coords=c("X","Y"),crs=st_crs(nc))  %>% group_by(id)  %>%  summarise(i=mean(id)) %>% st_cast("POLYGON")

现在比较我们研究区域中包含的多边形的周长:

sum(st_length(crop_nc),na.rm=T)
1307555 [m]

sum(st_length(crop_NC_sf),na.rm=T)
2610959 [m]

显然效果不是很好......

问题:

• 还有其他方法可以加速 st_crop()

• 有没有一种方法可以从点重新创建多边形以保留相互连接的“原始”顺序点？

Answer 1

通过 st_intersects(或 st_crop)实现正确的交点比仅仅对一堆坐标进行子集化要多一些。您将在下面找到一个适用于您的示例的有效解决方案，该解决方案试图回答您的两个问题。

对于问题 1，我建议在裁剪之前识别与 bbox 相交的多边形，在许多情况下可以大大加快速度，尤其是当你有很多多边形并且 bbox 的面积与范围相比较小时的多边形面积。

关于问题 2，您可以使用 sfheaders 包，它提供了从 data.frames/data.tables 等创建 sf 对象的方法。

最后，我稍微修改了计时代码，以更好地反射(reflect)每个步骤需要做什么才能提供类似的输出，从而使比较更加公平。

data.table 方法的最终结果与 st_crop 方法不同，因为这种方法仅将坐标子集添加到 bbox 内部，但不插入bbox 在适当位置的坐标。此操作可能代价高昂，因为我们需要确定这些正确的位置。因此，您将在结果中获得任意形状，并且事情甚至可能会破裂，因为您可能会产生无效的几何形状。

一般来说，我会坚持使用 st_crop 方法，这可能会在以后为您省去麻烦。如果你感兴趣的区域很大，有很多多边形，而你使用的 bbox 相对于整个区域来说很小，我希望在进行实际相交之前识别与 bbox 相交的多边形会加快速度。

library(sf)
#> Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3
library(mapview)
library(data.table)
library(ggplot2)
library(tidyverse)
library(sfheaders)
library(microbenchmark)

# load shapefile
nc = st_geometry(st_read(system.file("shape/nc.shp", package="sf")))
#> Reading layer `nc' from data source `C:\Users\Tim\R\win-library\3.5\sf\shape\nc.shp' using driver `ESRI Shapefile'
#> Simple feature collection with 100 features and 14 fields
#> geometry type:  MULTIPOLYGON
#> dimension:      XY
#> bbox:           xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
#> epsg (SRID):    4267
#> proj4string:    +proj=longlat +datum=NAD27 +no_defs

# Define a bounding-box that mimic a mowing-window or area of interest
bb <- st_bbox(c(xmin= -79, xmax=-78,ymin= 34.5, ymax= 35.5))
bb_sfc = st_as_sfc(bb)
st_crs(bb_sfc) <- st_crs(nc)


# Commute 'nc' into some sort of data.table object for fast subsetting, in preserving object's integrity (i.e. same id to all points of a given polygon)
nobs <- mapview::npts(nc,by_feature=T)
NC <- data.table::data.table(id=rep(1:length(nc),nobs),st_coordinates(nc)[,1:2])
head(NC)
#>    id         X        Y
#> 1:  1 -81.47276 36.23436
#> 2:  1 -81.54084 36.27251
#> 3:  1 -81.56198 36.27359
#> 4:  1 -81.63306 36.34069
#> 5:  1 -81.74107 36.39178
#> 6:  1 -81.69828 36.47178

# Compare cropping methods amon
test <- microbenchmark(
    sf_way = {
        # identify subset of nc that intersects with bbox
        isec = unlist(st_intersects(bb_sfc, nc))
        crop_nc <- st_crop(nc[isec], bb_sfc)
    }, 
    dt_way = {
        crop_NC <- NC[X >= bb[1] & X < bb[3] & Y>= bb[2] & Y < bb[4]]
        crop_NC_sf <- sfheaders::sf_polygon(crop_NC, "X", "Y", polygon_id = "id")
        st_crs(crop_NC_sf) = st_crs(nc)
    }
)  

print(test)
#> Unit: milliseconds
#>    expr      min       lq     mean   median       uq      max neval cld
#>  sf_way 4.270801 4.492551 4.945424 4.828702 5.051951 8.666301   100   b
#>  dt_way 3.101400 3.442551 3.705610 3.593301 3.887202 8.270801   100  a

sum(st_length(crop_nc),na.rm=T)
#> 1307555 [m]
sum(st_length(crop_NC_sf),na.rm=T)
#> 975793 [m]

mapview(st_bbox(crop_nc)) + crop_NC_sf

^{由 reprex package (v0.3.0) 创建于 2019-11-29}