合併兩個多邊形區域爲一個多邊形區域中的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() 

現在圈。