将 Shapefile 中的多边形与通用标签合并为：联合空间多边形答案

【问题标题】：Merging Polygons in Shape Files with Common Tag IDs: unionSpatialPolygons将 Shapefile 中的多边形与通用标签合并为：联合空间多边形
【发布时间】：2018-12-22 05:20:17
【问题描述】：

我正在尝试从形状文件中读取数据并将多边形与通用标签 ID 合并。

library(rgdal)
library(maptools)
if (!require(gpclib)) install.packages("gpclib", type="source")
gpclibPermit()

usa <- readOGR(dsn = "./path_to_data/", layer="the_name_of_shape_file")
usaIDs <- usa$segment_ID
isTRUE(gpclibPermitStatus())
usaUnion <- unionSpatialPolygons(usa, usaIDs)

当我尝试绘制合并的多边形时：

for(i in c(1:length(names(usaUnion)))){
  print(i)
  myPol <- usaUnion@polygons[[i]]@Polygons[[1]]@coords
  polygon(myPol, pch = 2, cex = 0.3, col = i)
}

所有合并的片段看起来都很好，除了密歇根周围的那些合并以一种非常奇怪的方式发生，因此这个特定片段的结果区域只给出一个小多边形，如下所示。

i = 10
usaUnion@polygons[[i]]@Polygons[[1]]@coords

输出：

          [,1]     [,2] 
 [1,] -88.62533 48.03317
 [2,] -88.90155 47.96025
 [3,] -89.02862 47.85066
 [4,] -89.13988 47.82408
 [5,] -89.19292 47.84461
 [6,] -89.20179 47.88386
 [7,] -89.15610 47.93923
 [8,] -88.49753 48.17380
 [9,] -88.62533 48.03317

原来是一个北方的小岛：

我怀疑问题在于，由于某种原因，unionSpatialPolygons 函数不喜欢地理上分离的多边形 [密歇根州的左右两侧]，但我还没有找到解决方案。

这是您可以复制的link to input data。

【问题讨论】：

您能否通过制作一个带有几个多边形的玩具示例来实现此功能的重现，或者以某种方式提供您的 shapefile？
我得到了 shapefile，现在应该可以重现了。我用示例数据的链接更新了帖子。

标签： r plot rgdal maptools

【解决方案1】：

我认为问题不在于unionSpatialPolygons，而在于你的情节。具体来说，您只为每个 ID 绘制第一个 'sub-polygon'。运行以下命令来验证出了什么问题：

for(i in 1:length(names(usaUnion))){
  print(length(usaUnion@polygons[[i]]@Polygons))
}

对于每一个，你只拿了第一个。

我使用以下代码得到了正确的多边形连接/绘图：

library(rgdal)
library(maptools)
library(plyr)

usa <- readOGR(dsn = "INSERT_YOUR_PATH", layer="light_shape")
# remove NAs
usa <- usa[!is.na(usa$segment_ID), ]
usaIDs <- usa$segment_ID

#get unique colors
set.seed(666)
unique_colors <- sample(grDevices::colors()[grep('gr(a|e)y|white', grDevices::colors(), invert = T)], 15)

colors <- plyr::mapvalues(
  usaIDs, 
  from = as.numeric(sort(as.character(unique(usaIDs)))), #workaround to get correct color order
  to = unique_colors
)
plot(usa, col = colors, main = "Original Map")


usaUnion <- unionSpatialPolygons(usa, usaIDs)
plot(usaUnion, col = unique_colors, main = "Joined Polygons")

【讨论】：

【解决方案2】：

这是一个使用sf 绘制此图的示例，它突出了包与dplyr 和summarise 一起使用的能力如何使该操作极具表现力和简洁性。我filter 删除了缺失的 ID、group_by ID、summarise（默认情况下进行联合），并轻松使用geom_sf 进行绘图。

library(tidyverse)
library(sf)
# Substitute wherever you are reading the file from
light_shape <- read_sf(here::here("data", "light_shape.shp"))
light_shape %>%
  filter(!is.na(segment_ID)) %>% 
  group_by(segment_ID) %>%
  summarise() %>%
  ggplot() +
  geom_sf(aes(fill = factor(segment_ID)))

【讨论】：