将两个多边形区域合并为R中的单个多边形区域

我是 R 中处理空间数据和多边形的新手。

我有两个独立的形状文件，包含从 Google 地球中提取的两个多边形。因此，基本上第一个形状文件是一个位置（例如购物中心等），第二个形状文件是第一个位置周围三公里的半径。我将两个形状文件作为 SpatialPolygonsDataFrames 读入 R。我使用以下代码：

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 中将两个多边形区域合并为一个多边形区域进行编码有一些建议，我将非常感激！

如果您决心使用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()

如果您需要将它们组合成一个spatialPolygonDataFrame，请使用它。这里的细微差别是你不能使用spRbind(...)如果两个图层具有共同的多边形 ID。所以调用spChFIDs(...)更改其中一个多边形的 IDpoly2（圆圈）到“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)

您会注意到，在这些图上，“圆圈”不是。这是因为我们使用的是长/纬度，而您距离赤道以南很远。数据需要重新投影。事实证明，适合南非的 CRS 是utm-33 http://spatialreference.org/ref/epsg/32733/.

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()

现在圈子是。