在大型数据集上加速 st_crop（sf 包）的方法
Method(s) to speed up st_crop (sf package) on large datasets
我需要为 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()
有没有办法从点重新创建多边形以保持 'original' 顺序点相互连接？
通过 st_intersects（或 st_crop）进行适当的交叉比仅仅对一堆坐标进行子集化要多一些。您将在下面找到一个适用于您的示例的有效解决方案，该解决方案试图回答您的两个问题。
对于问题 1，我建议在裁剪之前识别与 bbox 相交的多边形在许多情况下可以加快速度，特别是如果你有很多多边形并且 bbox 的面积与多边形区域的范围。
关于问题 2，您可以使用包 sfheaders，它提供了从 data.frames/data.tables 等创建 sf 对象的方法。
最后，我稍微修改了计时代码，以更好地反映每个步骤需要做什么才能提供相似的输出，从而使比较更加公平。
data.table 方法的最终结果与 st_crop 方法不同，因为这种方法仅将坐标子集添加到 bbox 内部，但不会将 bbox 的坐标插入到合适的位置。此操作可能代价高昂，因为我们需要确定这些正确的位置。因此，您将在结果中获得任意形状，并且事情甚至可能会破裂，因为您可能会产生无效的几何形状。
一般来说，我会坚持使用 st_crop 方法，这可能会在以后为您省去麻烦。
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.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 创建
我找到了解决上述问题的方法。不要使用 'st_intersects'，而是应该使用 'st_intersection'（见此 post），因为它 lines/polygons 的 returns segments/area 与多边形相交。只需 4 分钟即可计算出与 4mio 1 公顷多边形相交的线的长度总和（与具有替代路径的年龄相比）。
library(sf)
#> Warning: le package 'sf' a été compilé avec la version R 3.5.3
#> Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3
library(microbenchmark)
#> Warning: le package 'microbenchmark' a été compilé avec la version R 3.5.3
# load shapefile
nc = st_geometry(st_read(system.file("shape/nc.shp", package="sf")))
#> Reading layer `nc' from data source `D:\Program Files\R\R-3.5.2\library\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)
# Compare cropping methods amon
TEST <- microbenchmark(
crop_way = {
# identify subset of nc that intersects with bbox
isec = unlist(st_intersects(bb_sfc, nc))
crop_nc <- st_crop(nc[isec], bb_sfc)
},
intersection_way = {
inter_nc = st_intersection(bb_sfc, nc)
}
)
print(TEST)
#> Unit: milliseconds
#> expr min lq mean median uq max
#> crop_way 4.445589 4.533909 5.024318 4.701732 4.825892 13.62889
#> intersection_way 2.687720 2.731952 3.029011 2.812023 2.904466 8.39283
#> neval cld
#> 100 b
#> 100 a
sum(st_length(crop_nc),na.rm=T)
#> 1307555 [m]
sum(st_length(inter_nc),na.rm=T)
#> 1307555 [m]
Created on 2019-12-04 by the reprex package (v0.3.0)
我需要为 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()
有没有办法从点重新创建多边形以保持 'original' 顺序点相互连接？
通过 st_intersects（或 st_crop）进行适当的交叉比仅仅对一堆坐标进行子集化要多一些。您将在下面找到一个适用于您的示例的有效解决方案，该解决方案试图回答您的两个问题。
对于问题 1，我建议在裁剪之前识别与 bbox 相交的多边形在许多情况下可以加快速度，特别是如果你有很多多边形并且 bbox 的面积与多边形区域的范围。
关于问题 2，您可以使用包 sfheaders，它提供了从 data.frames/data.tables 等创建 sf 对象的方法。
最后，我稍微修改了计时代码，以更好地反映每个步骤需要做什么才能提供相似的输出，从而使比较更加公平。
data.table 方法的最终结果与 st_crop 方法不同，因为这种方法仅将坐标子集添加到 bbox 内部，但不会将 bbox 的坐标插入到合适的位置。此操作可能代价高昂，因为我们需要确定这些正确的位置。因此，您将在结果中获得任意形状，并且事情甚至可能会破裂，因为您可能会产生无效的几何形状。
一般来说，我会坚持使用 st_crop 方法，这可能会在以后为您省去麻烦。
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.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 创建我找到了解决上述问题的方法。不要使用 'st_intersects'，而是应该使用 'st_intersection'（见此 post），因为它 lines/polygons 的 returns segments/area 与多边形相交。只需 4 分钟即可计算出与 4mio 1 公顷多边形相交的线的长度总和（与具有替代路径的年龄相比）。
library(sf)
#> Warning: le package 'sf' a été compilé avec la version R 3.5.3
#> Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3
library(microbenchmark)
#> Warning: le package 'microbenchmark' a été compilé avec la version R 3.5.3
# load shapefile
nc = st_geometry(st_read(system.file("shape/nc.shp", package="sf")))
#> Reading layer `nc' from data source `D:\Program Files\R\R-3.5.2\library\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)
# Compare cropping methods amon
TEST <- microbenchmark(
crop_way = {
# identify subset of nc that intersects with bbox
isec = unlist(st_intersects(bb_sfc, nc))
crop_nc <- st_crop(nc[isec], bb_sfc)
},
intersection_way = {
inter_nc = st_intersection(bb_sfc, nc)
}
)
print(TEST)
#> Unit: milliseconds
#> expr min lq mean median uq max
#> crop_way 4.445589 4.533909 5.024318 4.701732 4.825892 13.62889
#> intersection_way 2.687720 2.731952 3.029011 2.812023 2.904466 8.39283
#> neval cld
#> 100 b
#> 100 a
sum(st_length(crop_nc),na.rm=T)
#> 1307555 [m]
sum(st_length(inter_nc),na.rm=T)
#> 1307555 [m]
Created on 2019-12-04 by the reprex package (v0.3.0)