在R中合并两个多边形区域为一个多边形区域

Question

我是第一次在R中使用空间数据和多边形。

我有两个不同的多边形形状文件，这些形状文件是从谷歌地球中提取的。因此，第一个形状文件是一个位置（例如购物中心等），而第二个形状文件是第一个位置周围三公里的半径。我将两个形状文件都读入R作为SpatialPolygonsDataFrames。我使用以下代码：

library(maptools)
library(sp)
library(spatstat)
options(digits=10) 

# Read polygon a

a <- readShapeSpatial(file.choose())
class(a)

spatstat.options(checkpolygons=FALSE)

r <- slot(a,"polygons")
r <- lapply(r, function(a) { SpatialPolygons(list(a)) })
windows <- lapply(r, as.owin)
Ploy_One <- tess(tiles=windows)

# Read polygon b

b <- readShapeSpatial(file.choose())
class(b)

spatstat.options(checkpolygons=FALSE)

s <- slot(b,"polygons")
s <- lapply(s, function(b) { SpatialPolygons(list(b)) })

windows <- lapply(s, as.owin)
Poly_Two <- tess(tiles=windows)

# Read polygon b

Combined_Region <- intersect.tess(Poly_One, Poly_Two)
plot(Combined_Region)

然而，我无法得到两个多边形的组合视图（即一个多边形在另一个内部的视图）。

如果有人能给我一些建议，告诉我如何在R中编写将两个多边形区域合并为单个多边形区域的代码，我将不胜感激！

- Carmen

---

你能否发布你的两个多边形shapefile的链接？ - jlhoward

谢谢！这是新链接 https://skydrive.live.com/redir?resid=7286AE33F47C4A63%21127。我在 spatstat 上找到了这些课程笔记（http://www.csiro.au/resources/pf16h - 第48页的 pdf 文件），它说每个多边形的顶点应该逆时针遍历外部边界和顺时针遍历内部边界（孔）。我尝试了这个方法，但是出现了一个错误，提示“窗口面积为负数；请检查所有多边形是否按正确方向遍历”。我可能做错了什么，如果成功了，我会告诉你的。非常感谢您的帮助！谢谢！ - Carmen

好的，我们正在取得进展。 "shapefile" 不是单个文件，而是具有相同名称和不同扩展名的文件集合，通常包括 .shp、.dbf、.shx、.prj 和 .shp.xml（并非总是所有这些）。因此，如果您的 "shapefile" 是 poly，那么实际文件将是 ploy.shp、poly.dbf、poly.shx 等。您需要将所有这些文件收集到一个 zip 文件中，并为 Polygon_One 和 Polygon_Two 发布该文件。 - jlhoward

抱歉！这是链接：https://skydrive.live.com/redir?resid=7286AE33F47C4A63%21131 - Carmen

成功了！使用以下代码 Z <- owin(poly = list(list(x = Poly_One_Long, y = Poly_One_Lat), list(x = Poly_Two_Long, y = Poly_Two_Lat))) ，其中 Poly_One 以逆时针方向绘制，Poly_Two 以顺时针方向绘制。Poly_Two 在 Poly_One 中形成一个孔洞。我没有将这些点作为 .shp 文件导入，而是作为简单的 .csv 文件导入。逆时针导入圆形多边形并不容易，我认为对于未来和更大规模的分析，该函数将不太有效，但现在我很高兴得到了一张漂亮的图片。如果您有任何进一步的建议，请分享，谢谢。 - Carmen

Answer 1

如果您已经承诺使用spatstat软件包，那么这篇文章可能对您没有太大帮助。如果不是，请继续阅读...

如果您只想将多边形绘制为单独的图层，请使用ggplot。

library(ggplot2)
library(sp)
library(maptools)

setwd("<directory with all your files...>")

poly1 <- readShapeSpatial("Polygon_One")
poly2 <- readShapeSpatial("Polygon_Two")
# plot polygons as separate layers,,,
poly1.df <- fortify(poly1)
poly2.df <- fortify(poly2)
ggplot() +
  geom_path(data=poly1, aes(x=long,y=lat, group=group))+
  geom_path(data=poly2, aes(x=long,y=lat, group=group), colour="red")+
  coord_fixed()

如果需要将它们合并为一个空间多边形数据框，则使用此方法。这里的细微差别是，如果两个图层具有相同的多边形 ID，则不能使用 spRbind(...)。所以调用 spChFIDs(...) 将 poly2（圆圈）中的一个多边形的 ID 改为 "R.3km"。

# combine polygons into a single shapefile
poly.combined <- spRbind(poly1,spChFIDs(poly2,"R.3km"))
# plot polygons using ggplot aesthetic mapping
poly.df <- fortify(poly.combined)
ggplot(poly.df) + 
  geom_path(aes(x=long, y=lat, group=group, color=group)) + 
  scale_color_discrete("",labels=c("Center", "3km Radius")) +
  coord_fixed()
# plot using plot(...) method for spatialObjects
plot(poly.combined)

您会注意到，在这些图中，"圆圈"并不是一个真正的圆形。这是因为我们使用了经纬度坐标系，而您所在的位置相对赤道较南。数据需要重新投影。结果表明，适用于南非的合适投影坐标系是UTM-33。

wgs.84 <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
utm.33 <- "+proj=utm +zone=33 +south +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
proj4string(poly.combined) <- CRS(wgs.84)
poly.utm33 <- spTransform(poly.combined,CRS(utm.33))
poly.df    <- fortify(poly.utm33)
ggplot(poly.df) + 
  geom_path(aes(x=long, y=lat, group=group, color=group)) + 
  scale_color_discrete("",labels=c("Center", "3km Radius")) +
  theme(axis.text=element_blank()) + labs(x=NULL,y=NULL) +
  coord_fixed()

现在圆圈是这样的。

- jlhoward