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

Question

我是新来的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中的单个多边形区域有一些建议，我会非常感激！

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，用这个。这里的细微差别在于如果两个图层具有公共的多边形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() 

现在圈。