studio - superponer graficas en r

Trazar el área espacial definida por múltiples polígonos

Question

Tengo un SpatialPolygonsDataFrame con 11589 objetos espaciales de clase "polígonos". 10699 de esos objetos consiste en exactamente 1 polígono. Sin embargo, el resto de esos objetos espaciales constan de múltiples polígonos (2 a 22).

Si un objeto de consiste en múltiples polígonos, tres escenarios son posibles:

1) Los polígonos adicionales podrían describir un "agujero" en el área espacial descrita por el primer polígono. 2) Los polígonos adicionales también podrían describir áreas geográficas adicionales, es decir, la forma de la región es bastante compleja y se describe al juntar varias partes. 3) A menudo es una mezcla de ambos, 1) y 2).

Mi pregunta es: ¿cómo trazar ese objeto espacial que se basa en múltiples polígonos?

Pude extraer y trazar la información del primer polígono, pero no he descubierto cómo trazar todos los polígonos de un objeto espacial tan complejo a la vez.

A continuación encontrarás mi código. El problema es la decimoquinta línea.

# 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

Puedo malinterpretar tu pregunta, pero es posible que estés haciendo esto mucho más difícil de lo necesario.

(Nota: tuve problemas para tratar con el archivo .zip en R, así que acabo de descargarlo y descomprimirlo en el sistema operativo).

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

Tenga en cuenta que ggplot no trata los agujeros correctamente: geom_path(...) funciona bien, pero geom_polygon(...) llena los agujeros. He tenido este problema antes (vea esta pregunta ), y en función de la falta de respuesta puede ser un error en ggplot. Como no estás usando geom_polygon(...) , esto no te afecta ...

En su código anterior, reemplazaría la línea:

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

con:

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