时间戳

首页 > 解决方案 > 保持多环多边形的最大环

r - 保持多环多边形的最大环

问题描述

我有一个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 的多边形。

这可能吗？

编辑

Calum 你谢谢你的回答。我从未使用过 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)

标签： rgissp

解决方案

这是使用sf包的一种方法。不幸的是，我不太熟悉sp，但也许这将是切换的动力！

问题在于，给定的几何形状是具有外部外环的单个多边形。要单独处理几何图形，我们需要将每个环拆分为自己的多边形。在sf中，aPOLYGON是一个矩阵列表，每个矩阵对应一个环中的点；aMULTIPOLYGON是这些POLYGON矩阵列表的列表。因此，要转换，我们执行以下操作：

  1. （将 转换SpatialPolygons为sfcwith st_as_sfc）
  2. map(list)跨越POLYGON将每个环包装在一个列表中，然后将列表列表转换为st_multipolygon. sfc请注意，如果有更多几何图形，我还会跨列映射，尽管这里只有一个。
  3. 将此列表转换回sfc对象，然后再转换为sf对象。这sfc只是具有其他属性（如坐标系）的几何列表，sf它sfc仅构成数据框中的一列。这是必需的，因为我们想知道每个几何图形的面积。
  4. 使用奇妙的功能st_cast将其分割MULTIPOLYGON成单独的多边形。我们添加了 arowid以便我们知道哪些 newPOLYGON属于 original 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)

由reprex 包(v0.2.0) 于 2018 年 5 月 9 日创建。

推荐阅读