Самый быстрый способ обрезать / обрезать сложный SpatialPolygonsDataFrame, используя другой многоугольник

Question

Какой самый быстрый способ обрезать (обрезать) сложный SpatialPolygonsDataFrame с @data, сохраненным в R, с использованием другого, возможно, сложного, SpatialPolygon? Я знаю два пути (показано под). Способ raster быстрее для менее сложных SpatialPolygonsDataFrames и возвращает SpatialPolygonsDataFrame, как в примере. Однако с большими SpatialPolygonsDataFrames он замедляется. Способ rgeos быстрее для больших объектов SpatialPolygonsDataFrame, но иногда дает сбой при очень сложных геометриях и не возвращает SpatialPolygonsDataFrames по умолчанию.

В последнее время я не обращал внимания на разработку ГИС в R, и вполне возможно, что сейчас есть более быстрые способы сделать это.

Примеры полигонов малы, чтобы учитывать компьютеры людей и пропускную способность. Рассмотрим «настоящие» полигоны 50–1000 МБ.

Настройка:

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

Большая часть геопространственной работы R теперь может быть выполнена с использованием пакета sf.

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

Похоже, что это может быть касание быстрее, чем растр и пр.

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

^{Создано в 2020-04-13 пакетом Представить (v0.3.0)}