强化 shapefile 时使用区域参数时出错

Question

我正在尝试使用 ggplot 绘制 shapefile (this one) 的多边形，以便可视化有关其中包含的区域的一些数据。下面的代码会产生错误：

library(rgdal)
library(tidyverse)
library(maptools)

hsa <- readOGR(dsn = "Data/hsa_bdry", layer = "HSA_Bdry") %>% fortify(region = "HSA93")

ggplot() + geom_polygon(data = hsa, aes(x = long, y = lat, group = group))

Error: TopologyException: Input geom 0 is invalid: Self-intersection at or near point -103.18495682693303 24.986200983871445 at -103.18495682693303 24.986200983871445

如何将 shapefile @data 槽中的数据附加到强化数据框而不会出现绘图错误？

到目前为止我尝试过的：

如果我从 fortify 中删除 region 参数，它会成功生成地图，但我需要区域 ID，以便我可以合并其余数据（以传递给 aes 参数在 ggplot 中）。

我尝试在调用 fortify 之前将 spTransform(CRS("+init=epsg:4269")) 添加到我的管道中（shapefile 的文档说它使用 North American Datum 1983，即 EPSG 4269），但随后出现错误

Error in spTransform(xSP, CRSobj, ...) : No transformation possible from NA reference system

Answer 1

如错误消息所示，shapefile 中有无效的 self-intersecting 点，导致 fortify 出现问题。您可以先从 rgeos 包中 pre-process 使用 gBuffer 的 shapefile：

# read in the shapefile as SpatialPolygonsDataFrame, but don't fortify it yet
hsa <- readOGR(dsn = "hsa_bdry", layer = "HSA_Bdry")

# check if there are invalid geometries & correct that
rgeos::gIsValid(hsa) #returns FALSE
hsa <- rgeos::gBuffer(hsa, byid = TRUE, width = 0)
rgeos::gIsValid(hsa) #returns TRUE

# fortify shouldn't encounter problems now.
hsa <- hsa %>% fortify(region = "HSA93")

这是使用基本主题选项的结果示例：

ggplot() + 
  geom_polygon(data = hsa, aes(x = long, y = lat, group = group),
               fill = "white", col = "black") +
  coord_map()

来自 GIS Stack Exchange 的

This post 也讨论了这个问题。

编辑：在 gBuffer 之前添加无效几何图形的视觉检查（可选）

如果您希望在更正之前检查无效的几何图形，您可以查看 shapefile 并放大报告问题的坐标。我通常更喜欢自己在 GIS 软件中执行此操作，但可以使用基本包的 plot 函数在 R 中绘制 SpatialPolygonsDataFrame。

# after loading in the shapefile but before any transformation
plot(hsa, col = "lightblue");axis(1);axis(2); title(main = "Full view")
points(x = -103.18495682693303, y = 24.986200983871445, 
       col = "red", cex = 5)

# toggle xlim & ylim till we get a good zoom
plot(hsa, col = "lightblue",
     xlim = c(-103.185, -103.18),
     ylim = c(24.98, 24.99)); axis(1); axis(2); title(main = "Zoom in")
points(x = -103.18495682693303, y = 24.986200983871445, 
       col = "red", cex = 5)

我们可以看出问题出在哪里了。三角形中两条线相交的点略微重叠。它很小，您自己可能不会注意到它，但是当我们尝试 fortify 它时它会引起问题。

执行 gBuffer 步骤后，再次查看 shapefile（除标题外代码没有变化）：

强化 shapefile 时使用区域参数时出错

Error using region argument when fortifying shapefile

r

shapefile

ggplot2

到目前为止我尝试过的：