繪製由多個多邊形定義的空間區域

Question

我有一個具有11589個「多邊形」類空間對象的SpatialPolygonsDataFrame。這些對象中的10699恰好由1個多邊形組成。但是，這些空間對象的其餘部分由多個多邊形（2到22）組成。

如果一個對象由多個多邊形的，三種情況都是可能的：

1）附加的多邊形可以描述在由第一多邊形描述的空間區域中的「洞」。 2）附加的多邊形還可以描述附加的地理區域，即該區域的形狀非常複雜，並通過將多個部分放在一起進行描述。 3）通常它是兩者的組合，1）和2）。

我的問題是：如何繪製這樣一個基於多個多邊形的空間物體？

我已經能夠提取和繪製第一個多邊形的信息，但我還沒有想出如何一次繪製這樣複雜的空間對象的所有多邊形。

下面你會發現我的代碼。問題是第15條最後一行。

# Load packages 
# --------------------------------------------------------------------------- 
library(maptools) 
library(rgdal) 
library(ggmap) 
library(rgeos) 


# Get data 
# --------------------------------------------------------------------------- 
# Download shape information from the internet 
URL <- "http://www.geodatenzentrum.de/auftrag1/archiv/vektor/vg250_ebenen/2012/vg250_2012-01-01.utm32s.shape.ebenen.zip" 
td <- tempdir() 
setwd(td) 
temp <- tempfile(fileext = ".zip") 
download.file(URL, temp) 
unzip(temp) 

# Get shape file 
shp <- file.path(tempdir(),"vg250_0101.utm32s.shape.ebenen/vg250_ebenen/vg250_gem.shp") 

# Read in shape file 
x <- readShapeSpatial(shp, proj4string = CRS("+init=epsg:25832")) 

# Transform the geocoding from UTM to Longitude/Latitude 
x <- spTransform(x, CRS("+proj=longlat +datum=WGS84")) 

# Extract relevant information 
att <- attributes(x) 
Poly <- att$polygons 


# Pick an geographic area which consists of multiple polygons 
# --------------------------------------------------------------------------- 
# Output a frequency table of areas with N polygons 
order.of.polygons.in.shp <- sapply(x@polygons, function(x) x@plotOrder) 
length.vector <- unlist(lapply(order.of.polygons.in.shp, length)) 
table(length.vector) 

# Get geographic area with the most polygons 
polygon.with.max.polygons <- which(length.vector==max(length.vector)) 
# Check polygon 
# x@polygons[polygon.with.max.polygons] 

# Get shape for the geographic area with the most polygons 
### HERE IS THE PROBLEM ### 
### ONLY information for the first polygon is extracted ### 
Poly.coords <- data.frame(slot(Poly[[polygon.with.max.polygons ]]@Polygons[[1]], "coords")) 

# Plot 
# --------------------------------------------------------------------------- 
## Calculate centroid for the first polygon of the specified area 
coordinates(Poly.coords) <- ~X1+X2 
proj4string(Poly.coords) <- CRS("+proj=longlat +datum=WGS84") 
center <- gCentroid(Poly.coords) 

# Download a map which is centered around this centroid 
al1 = get_map(location = c(lon=center@coords[1], lat=center@coords[2]), zoom = 10, maptype = 'roadmap') 

# Plot map 
ggmap(al1) + 
    geom_path(data=as.data.frame(Poly.coords), aes(x=X1, y=X2)) 

Answer 1

我可能會誤解你的問題，但有可能你做得比所需要的難得多。

（注意：我在R中處理.zip文件時遇到了問題，所以我只是將它下載並解壓縮到OS中）。

library(rgdal) 
library(ggplot2) 

setwd("< directory with shapefiles >") 
map <- readOGR(dsn=".", layer="vg250_gem", p4s="+init=epsg:25832") 
map <- spTransform(map, CRS("+proj=longlat +datum=WGS84")) 

nPolys <- sapply(map@polygons, function(x)length(x@Polygons)) 
region <- map[which(nPolys==max(nPolys)),] 
plot(region, col="lightgreen") 

# using ggplot... 
region.df <- fortify(region) 
ggplot(region.df, aes(x=long,y=lat,group=group))+ 
    geom_polygon(fill="lightgreen")+ 
    geom_path(colour="grey50")+ 
    coord_fixed() 

注意ggplot不處理得當孔：geom_path(...)工作正常，但geom_polygon(...)填充孔。我之前遇到過這個問題（請參閱this question），並且基於缺乏響應，它可能是ggplot中的一個錯誤。由於不使用geom_polygon(...)，這並不影響你...

在上面的代碼，你將取代行：

ggmap(al1) + geom_path(data=as.data.frame(Poly.coords), aes(x=X1, y=X2)) 

有：

ggmap(al1) + geom_path(data=region.df, aes(x=long,y=lat,group=group))