Unir polígonos en R

29

Me pregunto cómo unir polígonos espaciales con código R.

Estoy trabajando con datos del censo donde ciertas áreas cambian con el tiempo y deseo unir los polígonos y los datos correspondientes y simplemente informar sobre las áreas unidas. Estoy manteniendo una lista de polígonos que tienen cambios de censo a censo y que planeo fusionar. Me gustaría usar esta lista de nombres de áreas como una lista de búsqueda para aplicar a los datos del censo de diferentes años.

Me pregunto qué función R usar para fusionar polígonos seleccionados y datos respectivos. Lo busqué en Google pero simplemente me confundo con los resultados.

Geoconfundido
fuente
La respuesta a la mayoría de las operaciones de geometría, como disolución de polígonos, superposición, punto en polígono, intersección, unión, etc., es el paquete rgeos.
Spacedman
1
La Oficina del Censo de los Estados Unidos publica tablas para hacer esto para 1990-2000 y 2000-2010. Se pueden administrar con combinaciones de bases de datos , que se implementan mediante Rla mergefunción de.
whuber

Respuestas:

39

La siguiente solución se basa en una publicación de Roger Bivand en R-sig-Geo . Tomé su ejemplo reemplazando el archivo de forma alemán con algunos datos del censo de Oregón que puede descargar desde aquí (tome todos los componentes del archivo de forma de 'condados de Oregon y datos del censo').

Comencemos cargando los paquetes requeridos e importando el archivo shape en R.

# Required packages
libs <- c("rgdal", "maptools", "gridExtra")
lapply(libs, require, character.only = TRUE)

# Import Oregon census data
oregon <- readOGR(dsn = "path/to/data", layer = "orcounty")
oregon.coords <- coordinates(oregon)

A continuación, necesita alguna variable de agrupación para agregar los datos. En nuestro ejemplo, la agrupación se basa simplemente en las coordenadas de un solo condado. Vea la imagen a continuación, los bordes negros indican los polígonos originales, mientras que los bordes rojos representan polígonos agregados por oregon.id.

# Generate IDs for grouping
oregon.id <- cut(oregon.coords[,1], quantile(oregon.coords[,1]), include.lowest=TRUE)

# Merge polygons by ID
oregon.union <- unionSpatialPolygons(oregon, oregon.id)

# Plotting
plot(oregon)
plot(oregon.union, add = TRUE, border = "red", lwd = 2)

Original y agrupado shapefile de Oregon

Hasta aquí todo bien. Sin embargo, los atributos de datos relacionados con las subregiones del archivo de forma original (por ejemplo, densidad de población, área, etc.) se pierden cuando se realiza unionSpatialPolygons. Supongo que también le gustaría agregar sus datos censales asociados al archivo shape, por lo que necesitará un paso intermedio.

Primero debe convertir sus polígonos en un marco de datos para realizar la agregación. Ahora tomemos las columnas de atributos de datos de seis a ocho ("AREA", "POP1990", "POP1997") y agrúpelas de acuerdo con la función de aplicación de IDs anterior sum.

# Convert SpatialPolygons to data frame
oregon.df <- as(oregon, "data.frame")

# Aggregate and sum desired data attributes by ID list
oregon.df.agg <- aggregate(oregon.df[, 6:8], list(oregon.id), sum)
row.names(oregon.df.agg) <- as.character(oregon.df.agg$Group.1)

Finalmente, reconvierta su marco de datos de nuevo a un archivo de forma SpatialPolygonsDataFramepreviamente unificado oregon.uniony obtenga tanto polígonos generalizados como sus datos censales derivados del paso de agregación de resumen anterior.

# Reconvert data frame to SpatialPolygons
oregon.shp.agg <- SpatialPolygonsDataFrame(oregon.union, oregon.df.agg)

# Plotting
grid.arrange(spplot(oregon, "AREA", main = "Oregon: original county area"), 
             spplot(oregon.shp.agg, "AREA", main = "Oregon: aggregated county area"), ncol = 1)

Áreas de Oregon

fdetsch
fuente
10

Aquí hay una solución usando el paquete sf:

library(tidycensus)
library(dplyr)
library(sf)
library(ggplot2)

# get data from tindycensus for demonstration (note you need an API key, folow instructions here: https://walkerke.github.io/tidycensus/articles/basic-usage.html)
census <- tidycensus::get_acs(geography = "tract", variables = "B19013_001",
                           state = "TX", county = "Tarrant", geometry = TRUE) %>% 
  arrange(NAME)

# reduce dataset size
census <- census[1:8,]

# create grouping variable
group_1 <- census$GEOID[1:2]
group_2 <- census$GEOID[6:8]

census <- census %>% mutate(group = case_when(GEOID %in% group_1 ~ 'newgroup1',
                                              GEOID %in% group_2 ~ 'newgroup2',
                                              TRUE ~ GEOID))

# summarise by grouping variable (performs a union on grouped polygons and sums 'estimate')
census2 <- group_by(census, group) %>% 
  summarise(estimate = sum(estimate), do_union = TRUE)

# visualise using ggplot2 development version and facet by merged/unmerged datasets
plot_data <- rbind(census %>% select(group, estimate) %>%
                     mutate(facet = "unmerged"), 
                   census2 %>% mutate(facet = "merged"))

gp <- ggplot() + 
      geom_sf(data = plot_data, aes(fill = estimate), color = 'white') + 
      scale_fill_viridis_c() + 
      facet_wrap(~facet, ncol = 1)

ingrese la descripción de la imagen aquí

sebdalgarno
fuente
Pensé que solo agregaría una pequeña advertencia aquí, por si acaso: tenga cuidado de usar summarise()derivados con el do_unionargumento, ya que acabo de hacer algo así summarise_if(shapefile, predic.function, sum, na.rm = TRUE, do_union = TRUE), que terminó sumando un VERDADERO en cada celda (es decir, +1 para todas las operaciones). ¿Necesita investigar más para averiguar si eso es algo que debería informarse (al menos para una advertencia adicional) ...?
Stragu