通过st_intersects
(或st_crop
)进行适当的交集不只是对一组坐标进行分组。在下面,您将为您的示例找到一个可行的解决方案,试图回答您的两个问题。
对于问题1,我建议在裁剪之前识别与bbox相交的多边形在很多情况下可以使速度大大提高,尤其是如果您有很多多边形并且bbox的面积相对于范围而言较小多边形区域。
关于问题2,您可以使用包 sfheaders ,该包提供了从data.frames / data.tables等创建sf对象的方法。
最后,我对时序代码进行了一些调整,以更好地反映出提供相似输出所需的每个步骤,从而使比较更加公平。
data.table
方式的最终结果与st_crop
有所不同,因为这种方法仅将bbox内的坐标子集作为子集,而没有将bbox的坐标插入适当的位置。由于我们需要确定这些正确的位置,因此此操作的成本可能很高。因此,您将得到任意形状的结果,甚至可能会断裂,因为您可能会产生无效的几何形状。
通常,我会坚持使用st_crop
方法,这可能会在以后为您省去麻烦。如果您感兴趣的区域很大,有很多多边形,并且您使用的bbox相对于整个区域来说很小,那么我希望在进行实际的相交之前确定与bbox相交的多边形会大大加快速度。 / p>
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创建
,
我找到了克服上述问题的方法。不要使用“ st_intersects”,而应该使用“ st_intersection”(请参见此post),因为它会返回与多边形相交的线/多边形的线段/区域。仅用了4分钟即可计算出与1公顷多边形相交4mio的线的总和(与使用替代路径的年龄相比)。
library(sf)
#> Warning: le package 'sf' a été compilé avec la version R 3.5.3
#> Linking to GEOS 3.6.1,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,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)
本文链接:https://www.f2er.com/3006687.html