使用另一个多边形裁剪/剪辑复杂 SpatialPolygonsDataFrame 的最快方法

Question

使用另一个可能很复杂的 SpatialPolygon 裁剪（裁剪）复杂的 SpatialPolygonsDataFrame 的最快方法是什么？我知道两种方法（如下所示）。 raster 方式对于不太复杂的 SpatialPolygonsDataFrames 更快，并返回一个 SpatialPolygonsDataFrame，如示例中所示。但是，大型 SpatialPolygonsDataFrames 会变慢。 rgeos 方式对于大型 SpatialPolygonsDataFrame 更快，但 it sometimes fails with very complex geometries 并且默认不返回 SpatialPolygonsDataFrames。

我最近没有关注 R 中的 GIS 开发，现在可能有更快的方法来做这件事。

示例多边形很小，以尊重人们的计算机和带宽。考虑 50-1000 MB 的“真实”多边形。

设置：

library(rnaturalearth)
library(sp)
library(raster)
library(rgeos)

dt <- rnaturalearth::ne_countries()
clip_boundary <- sp::SpatialPolygons(list(sp::Polygons(
list(sp::Polygon(data.frame(lon = c(0, 180, 180, 0), lat = c(40, 40, 80, 80)))), ID = 1)))

rgeos 方式：

system.time({
clipped.dt <- rgeos::gIntersection(dt, clip_boundary, byid = TRUE)
ids <- sapply(slot(clipped.dt, "polygons"), function(x) slot(x, "ID"))
ids <- sapply(strsplit(ids, " "), "[", 1) 

tmp.df <- data.frame(dt@data[ids,])
names(tmp.df) <- names(dt@data)

out <- sp::SpatialPolygonsDataFrame(clipped.dt, tmp.df, match.ID = FALSE)
})

#   user  system elapsed 
#  0.069   0.002   0.074 

class(out)
# [1] "SpatialPolygonsDataFrame"
# attr(,"package")
# [1] "sp"

raster 方式：

system.time({
  out <- raster::crop(dt, clip_boundary)
})

#   user  system elapsed 
#  0.042   0.001   0.043 

class(out)
# [1] "SpatialPolygonsDataFrame"
# attr(,"package")
# [1] "sp"

参考图（与问题无关）：

plot(out)

Answer 1

现在可以使用 sf 包完成 R 的大部分地理空间工作。

https://r-spatial.github.io/sf/index.html

看起来它可能比光栅和 sp 快一点。

library(rnaturalearth)
library(sp)
#library(raster)
#library(rgeos)
library(sf)
#> Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0
library(ggplot2)

# Original data 
dt <- rnaturalearth::ne_countries()
clip_boundary <- sp::SpatialPolygons(list(sp::Polygons(
  list(sp::Polygon(data.frame(lon = c(0, 180, 180, 0), lat = c(40, 40, 80, 80)))), ID = 1)))

# Original data as sf objects
dt_sf <- st_as_sf(dt)
boundary_sf <- st_as_sf(clip_boundary)

# clip 
clipped <- st_crop(dt_sf, boundary_sf)
#> although coordinates are longitude/latitude, st_intersection assumes that they are planar
#> Warning: attribute variables are assumed to be spatially constant throughout all
#> geometries

#plot
ggplot(clipped) + geom_sf()

system.time({
  out <- st_crop(dt_sf, boundary_sf)
})
#> although coordinates are longitude/latitude, st_intersection assumes that they are planar
#> Warning: attribute variables are assumed to be spatially constant throughout all
#> geometries
#>    user  system elapsed 
#>   0.019   0.000   0.018

system.time({
  out <- raster::crop(dt, clip_boundary)
})
#> Warning in raster::crop(dt, clip_boundary): non identical CRS
#>    user  system elapsed 
#>   0.526   0.000   0.525

system.time({
  clipped.dt <- rgeos::gIntersection(dt, clip_boundary, byid = TRUE)
  ids <- sapply(slot(clipped.dt, "polygons"), function(x) slot(x, "ID"))
  ids <- sapply(strsplit(ids, " "), "[", 1) 

  tmp.df <- data.frame(dt@data[ids,])
  names(tmp.df) <- names(dt@data)

  out <- sp::SpatialPolygonsDataFrame(clipped.dt, tmp.df, match.ID = FALSE)
})
#> Warning in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, drop_lower_td,
#> unaryUnion_if_byid_false, : spgeom1 and spgeom2 have different proj4 strings
#>    user  system elapsed 
#>   0.045   0.000   0.045

^{由reprex package (v0.3.0) 于 2020 年 4 月 13 日创建}