Сохраните самое большое кольцо из многокольцевого многоугольника

У меня есть SpatialPolygons объект. Этот объект имеет одну особенность, объект Polygons, который сам состоит из нескольких объектов Polygon.

Я хотел бы создать подмножество объекта SpatialPolygons, чтобы объект Polygons имел только один объект Polygon, то есть Polygon с наибольшей площадью.

Я пробовал много разных подходов и не могу понять, как выполнить подмножество на уровне суб-полигонов.

Вот пример:

#create polygon
Sr1 = Polygon(cbind(c(2,4,4,1,2),c(2,3,5,4,2)))
Sr2 = Polygon(cbind(c(5,4,2,5),c(2,3,2,2)))
Srs = Polygons(list(Sr1,Sr2), "s1/2")

SpP = SpatialPolygons(list(Srs))

В этом примере SpP имеет единственный объект Polygons. Первый субполигон объекта Polygons имеет площадь 5,5, а второй - 1,5. Я хотел бы создать подмножество объекта SpatialPolygons так, чтобы у объекта Polygons был только один Polygon, с площадью 5.5.

Это возможно?

РЕДАКТИРОВАТЬ

Калум Спасибо за ответ. Никогда не пользовался пакетом sf, но выглядит круто. Тем более, что он так хорошо играет с dplyr!

Однако в итоге я нашел другое решение. Моя главная проблема заключалась в том, что я не понимал, как работает структура объекта SpatialPolygons. Объект SpatialPolygons содержит объекты 1 .. * Polygons, поэтому конструктор принимает список объектов Polygons. Объект Polygons содержит объекты 1 .. * Polygon, поэтому конструктор принимает список объектов Polygon. Имея это в виду, я использовал следующее решение.

plst <- SpP@polygons[[1]]@Polygons #get the list of Polygon objects
plst <- plst[which.max(sapply(plst,function(p) return(p@area)))] #filter to the max area
spoly2 <- SpatialPolygons(list(Polygons(plst,'id')))

График показывает, что это фильтрует кольца по наибольшей площади.

plot(SpP)
plot(spoly2,col=alpha('red',0.1),add=T)

r sp gis
person Ben Carlson    schedule 09.05.2018    source источник
comment
Есть только один объект? Вы создаете их все вручную вот так или это всего лишь пример? Пытаюсь понять, что лучше всего.   -  person Calum You    schedule 10.05.2018

Ответы (1)


arrow_upward
1
arrow_downward

Вот один из подходов, в котором используется пакет sf. К сожалению, я не слишком хорошо знаком с sp, но, возможно, это будет стимулом для перехода!

Проблема в том, что данная геометрия представляет собой один многоугольник с внешним внешним кольцом. Чтобы работать с геометрией отдельно, нам нужно разбить каждое кольцо на отдельный многоугольник. В sf POLYGON - это список матриц, каждая матрица соответствует точкам в кольце; MULTIPOLYGON - список тезисов POLYGON списки матриц. Итак, чтобы преобразовать, мы делаем следующее:

  1. (преобразовать SpatialPolygons в sfc с st_as_sfc)
  2. map(list) через POLYGON обертывание каждого кольца в списке, а затем преобразовать список списков с помощью st_multipolygon. Обратите внимание, что я также отображаю столбец sfc на случай, если геометрических фигур больше, хотя здесь у вас есть только одна.
  3. Преобразуйте этот список обратно в объект sfc, а затем объект sf. sfc - это просто список геометрий с дополнительными свойствами, такими как система координат, sf делает sfc только одним столбцом во фрейме данных. Это необходимо, потому что мы хотим знать площадь каждой геометрии.
  4. Используйте замечательную функцию st_cast, чтобы разделить MULTIPOLYGON на отдельные многоугольники. Мы добавили rowid, чтобы знать, какие новые POLYGON принадлежали исходному MULTIPOLYGON.
  5. Добавьте столбец площади с mutate и st_area, group_by, а затем используйте top_n для фильтрации только до самого большого многоугольника в группе. Обратите внимание, как красиво sf объекты работают с dplyr глаголами!

Сюжет показывает результат. Исходный многоугольник выделен красным цветом, самое большое его кольцо - синим.

library(tidyverse)
library(sf)
#> Linking to GEOS 3.6.1, GDAL 2.2.3, proj.4 4.9.3
library(sp)
Sr1 = Polygon(cbind(c(2,4,4,1,2),c(2,3,5,4,2)))
Sr2 = Polygon(cbind(c(5,4,2,5),c(2,3,2,2)))
Srs = Polygons(list(Sr1,Sr2), "s1/2")

SpP = SpatialPolygons(list(Srs))
sfc <- st_as_sfc(SpP)
sfc
#> Geometry set for 1 feature 
#> geometry type:  POLYGON
#> dimension:      XY
#> bbox:           xmin: 1 ymin: 2 xmax: 5 ymax: 5
#> epsg (SRID):    NA
#> proj4string:    NA
#> POLYGON ((2 2, 1 4, 4 5, 4 3, 2 2), (5 2, 4 3, ...

sf_obj <- sfc %>%
  map(
    .f = function(polygon) polygon %>% map(list) %>% st_multipolygon()
  ) %>%
  st_sfc() %>%
  st_sf() %>%
  rowid_to_column() %>%
  st_cast("POLYGON") %>% 
  mutate(area = st_area(geometry)) %>%
  group_by(rowid) %>%
  top_n(1, area)
#> Warning in st_cast.sf(., "POLYGON"): repeating attributes for all sub-
#> geometries for which they may not be constant

plot(sfc, col = "red")
plot(sf_obj$geometry, col = "blue", add = TRUE)

Создано 9 мая 2018 г. пакетом REPEX (v0.2.0).

person Calum You    schedule 09.05.2018