2

クラス「ポリゴン」の11589個の空間オブジェクトを持つSpatialPolygonsDataFrameがあります。これらのオブジェクトの 10699 は、正確に 1 つのポリゴンで構成されています。ただし、残りの空間オブジェクトは複数のポリゴン (2 ～ 22) で構成されています。

のオブジェクトが複数のポリゴンで構成されている場合、次の 3 つのシナリオが考えられます。

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

1 に答える 1

3

私はあなたの質問を誤解しているかもしれませんが、これを必要以上に難しくしている可能性があります.

(注: 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(...)、穴を埋めます。以前にこの問題が発生したことがあり (この質問を参照)、応答がないことに基づいて、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))
于 2014-02-23T02:02:10.327 に答える