Сохранить наибольшее кольцо многоугольного многоугольника - PullRequest
Новая
       22

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

0 голосов
/ 10 мая 2018

У меня есть 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 был только один полигон, один с площадью 5.5.

Возможно ли это?

EDIT

Калум. Спасибо, что ответили. Я никогда не использовал пакет 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)

1 Ответ

0 голосов
/ 10 мая 2018

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

Проблема в том, что, как указано, геометрия представляет собой один многоугольник с внешним внешним кольцом. Чтобы работать с геометрией отдельно, нам нужно разбить каждое кольцо на собственный многоугольник. В sf a 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 it, а затем используйте 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)

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

...