Нахождение пространственных полигонов внутри пространственного полигона

Question

У меня есть многоугольник округа США, который представляет собой простую схему этого района.Это SpatialPolygonsDataFrame.

CRS("+proj=longlat +datum=WGS84"))

Далее у меня есть блок переписи штата (Техас) SpatialPolygonsDataFrame.

CRS("+proj=longlat +datum=WGS84")

Моя цель - найти группы блоков переписи, которые вписываются в этот округ США.

Мой подход кажется некорректным.

mysp <- raster::aggregate(rbind(SUBBLOCK_ALL.shp, CD32.shp))

Ошибка в as (x, "SpatialPolygons"): нет метода или значения по умолчанию для приведения «NULL» к «SpatialPolygons»

mysp <- rgeos::gIntersection(CD32.shp,SUBBLOCK_ALL.shp,byid = TRUE)

Ошибка в RGEOSBinTopoFunc (spgeom1, spgeom2, byid, id, drop_lower_td, unaryUnion_if_byid_false,: попытка получить слот "proj4string" из объекта базового класса ("NULL") снет слотов

mysp <- rgeos::gIntersection(CD32.shp,SUBBLOCK_ALL.shp,byid = TRUE)

Ошибка в RGEOSBinTopoFunc (spgeom1, spgeom2, byid, id, drop_lower_td, unaryUnion_if_byid_false,: попытка получить слот "proj4string" из объекта базового класса (""NULL ") без слотов

Мой вопрос: какой подход лучше?

Answer 1

Вот способ использования пакета sf и включенного набора данных nc для Северной Каролины. Это, очевидно, в другом масштабе, чем ваша проблема, но предпосылка должна быть такой же. Мы можем создать большой квадратный «район», охватывающий некоторые округа Северной Каролины, а затем использовать filter с st_within, чтобы оставить только те, которые находятся в пределах границы. Вы можете использовать другие предикаты, такие как st_intersects, в зависимости от того, что вы хотите сделать с теми, которые лежат на границе. Это также применимо к большему количеству районов, поскольку предикаты геометрии sf фактически сравнивают все геометрии в своих двух аргументах.

library(tidyverse)
library(sf)
#> Linking to GEOS 3.6.1, GDAL 2.2.3, proj.4 4.9.3
nc <- system.file("shape", "nc.shp", package = "sf") %>%
  read_sf() %>%
  st_transform(3857)

district <- st_polygon(list(matrix(
  data = c(-9000000, 4000000, -8500000, 4000000, -8500000, 4500000, -9000000, 4500000, -9000000, 4000000),
  ncol = 2,
  byrow = TRUE
  ))) %>%
  st_sfc(crs = 3857)

contained <- nc %>%
  filter(st_within(., district, sparse = FALSE))

plot(nc$geometry)
plot(district, add = TRUE)
plot(contained$geometry, add = TRUE, col = "red")

Создано в 2018-05-11 пакетом Представ. (v0.2.0).