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

Question

我不熟悉在 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)

但是，我没有得到两个多边形的组合 View (一个多边形在另一个多边形中的 View )。

如果有人对我如何编写代码将两个多边形区域合并为 R 中的单个多边形区域有一些建议，我将非常感激!

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

如果您需要将它们组合成一个 spatialPolygonDataFrame，请使用它。这里的细微差别是，如果两个层具有公共(public)多边形 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)

您会注意到，在这些图中，“圆”不是。这是因为我们使用的是经度/纬度，而您位于赤道以南很远的地方。数据需要重新投影。原来适合南非的 CRS 是 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()

现在圆圈是。