st_intersection к пространственному полигону данных не работает

Я пытаюсь получить пересечение двух шейп-файлов (участков переписи, которые попадают в границы определенных городских районов). Я могу успешно получить пересекающиеся объекты, однако, когда я пытаюсь преобразовать вывод sf_intersect в SpatialPolygonsDataframe, я получаю сообщение об ошибке:

"Ошибка в as_Spatial(from): преобразование из типа объекта sfc_GEOMETRY в sp не поддерживается"

Это мой код:

library(sf)
library(dplyr)
library(tigris)
library(sp)

#download shapefiles corresponding to metro areas 
metro_shapefiles<-core_based_statistical_areas(cb = FALSE, year = 2016)
#convert to sf and filter
metro_shapefiles<-st_as_sf(metro_shapefiles)%>%filter(GEOID==31080 )
#Data for California
census_tracts_california<-tracts(state="CA",year=2016)
census_tracts_california<-st_as_sf(census_tracts_california)

#INTERSECT AND CONVERT BACK TO SP
census_tracts_intersected1<-st_intersection(census_tracts_california,
                                            metro_shapefiles)

#back to spatial
census_tracts_intersected1<-as(census_tracts_intersected1,"Spatial")

2 ответа

Сообщение об ошибке говорит вам, что вы не можете преобразовать sfc_GEOMETRY к Spatial объект. Здесь нет sp эквивалентный объект.

В вашем результате пересечения у вас есть смесь геометрий (и, следовательно, вы получили sfc_GEOMETRY как ваша "геометрия"). Вы можете увидеть все геометрии здесь:

types <- vapply(sf::st_geometry(census_tracts_intersected1), function(x) {
    class(x)[2]
}, "")

unique(types)
# [1] "POLYGON"         "MULTILINESTRING" "MULTIPOLYGON"

Если вы хотите, вы можете извлечь каждый тип геометрии и конвертировать их в SP по отдельности:

lines <- census_tracts_intersected1[ grepl("*LINE", types), ]
polys <- census_tracts_intersected1[ grepl("*POLYGON", types), ]

spLines <- as(lines, "Spatial")
spPolys <- as(polys, "Spatial")

Дополнительная информация

Я упомянул в комментариях вы можете использовать st_join, Однако это может не дать вам желаемого результата. В sf В библиотеке есть геометрические двоичные предикаты, такие как ?st_intersects и геометрические операции, такие как ?st_intersection

Предикаты возвращают разреженную (по умолчанию) или плотную матрицу, сообщающую, с какой геометрией y пересекается каждая геометрия x. Если вы используете это в st_join, он возвратит (оригинальные) геометрии, которые пересекаются, а не разреженную матрицу.

В то время как операции (например, st_intersection) вычислит пересечение и вернет новую геометрию.

Пример использования

Предикаты (st_intersects) можно использовать внутри st_join и они вернут исходные геометрии, которые "пересекаются"

sf_join <- sf::st_join(census_tracts_california, metro_shapefiles, join = st_intersects)

В этом случае это дает один type объекта

types <- vapply(sf::st_geometry(sf_join), function(x) {
  class(x)[2]
}, "")

unique(types)
# [1] "MULTIPOLYGON"

## so you can convert to a Spatial object
spPoly <- as(sf_join, "Spatial")

Но вам нужно решить, если результат st_intersect это то, что вам нужно, или вам нужны новые геометрии, данные st_intersection,

дальнейшее чтение

  • информация о каждом объединении находится в блоге sf.

  • пространственные предикаты и примеры того, что делают различные операции, в Википедии (с некоторыми хорошими иллюстрациями)


Благодарим пользователя @lbussett за описание разницы между st_intersect а также st_intersection

Преобразование в объект не может обрабатывать смешанные линии и полигоны. После пересечения вы можете извлечь только многоугольники (и отбросить любые линии), используя:
st_collection_extract("POLYGON")

Ваш пример больше не работает, поэтому я создал новый пример, который пересекает два многоугольника с общей стороной. Это приводит к выходу пересечения многоугольника и линии. Со второй попытки я передал пересечение через функцию st_collection_extract() до успешного преобразования вSpatialобъект.

      library(sf)
library(dplyr)
library(sp)

#Create some boxes
BoxA <- st_polygon(list(cbind(c(0,0,2,2,0),c(0,2,2,0,0))))
BoxB <- st_polygon(list(cbind(c(1,1,3,3,1),c(1,3,3,1,1))))
BoxC <- st_polygon(list(cbind(c(2,2,4,4,2),c(0,2,2,0,0))))

#Create a funny shaped union to help demonstrate the intersection issue
BoxAB <- st_union(BoxA,BoxB)

plot(BoxAB)
plot(BoxC,add=TRUE,border="blue")

Пример полигонов для пересечения

      #Intersect of BoxAB with BoxC results in a line and a polygon
BoxIntersects<-st_intersection(BoxAB,BoxC)
plot(BoxIntersects)

Пересечение, состоящее из многоугольника и линии

      #back to spatial fails
SpatialVersionOfIntersects<-as(BoxIntersects,"Spatial")

Ошибка в .as_Spatial(from, cast, ID): преобразование из типа объектов sfc_GEOMETRY в sp не поддерживается.

      #Intersect again, but this time extract only the polygons
BoxIntersects<-st_intersection(BoxAB,BoxC) %>% st_collection_extract("POLYGON")

#back to spatial suceeds!
SpatialVersionOfIntersects<-as(BoxIntersects,"Spatial")
Другие вопросы по тегам