我有一个 SpatialPolygonsDataFrame,其中包含 11589 个“多边形”类的空间对象。其中 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))
最佳答案
我可能误解了你的问题,但你可能会让这个问题变得比必要的更加困难。
(注意:我在 R 中处理 .zip 文件时遇到问题,因此我只是在操作系统中下载并解压缩它)。
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))
关于r - 绘制由多个多边形定义的空间区域,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/21962452/