Uniendo polígonos en R

Question

Estoy preguntándome cómo unir polígonos espaciales usando código R?

Estoy trabajando con datos censales donde ciertas áreas cambian con el tiempo y deseo unir los polígonos y los datos correspondientes para simplemente informar sobre las áreas unidas. Estoy manteniendo una lista de polígonos que han cambiado 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 datos censales de diferentes años.

Me pregunto qué función de R utilizar para fusionar polígonos seleccionados y sus respectivos datos. He buscado en Google pero simplemente me confundo con los resultados.

Answer 1

La siguiente solución se basa en una publicación de Roger Bivand en R-sig-Geo. Tomé su ejemplo reemplazando el shapefile alemán con algunos datos censales de Oregón que puedes descargar desde aquí (toma todos los componentes del shapefile de 'Condados y datos censales de Oregón').

Comencemos cargando los paquetes requeridos e importando el shapefile a R.

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

# Importar datos censales de Oregón
oregon <- readOGR(dsn = "path/to/data", layer = "orcounty")
oregon.coords <- coordinates(oregon)

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

# Generar IDs para la agrupación
oregon.id <- cut(oregon.coords[,1], quantile(oregon.coords[,1]), include.lowest=TRUE)

# Unir polígonos por ID
oregon.union <- unionSpatialPolygons(oregon, oregon.id)

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

Hasta ahora, todo va bien. Sin embargo, los atributos de los datos relacionados con las subregiones originales del shapefile (por ejemplo, densidad de población, área, etc.) se pierden al realizar unionSpatialPolygons. Supongo que también te gustaría agregar tus datos censales asociados al shapefile, por lo que necesitarás un paso intermedio.

Primero debes convertir tus polígonos a un dataframe para poder realizar la agregación. Ahora tomemos las columnas de atributos de datos seis a ocho ("AREA", "POP1990", "POP1997") y agrupémoslas según los IDs anteriores aplicando la función sum.

# Convertir SpatialPolygons a dataframe
oregon.df <- as(oregon, "data.frame")

# Agregar y sumar los atributos de datos deseados por lista de ID
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, reconvierne tu dataframe de nuevo a un SpatialPolygonsDataFrame proporcionando el shapefile unificado previamente oregon.union y obtendrás tanto polígonos generalizados como tus datos censales derivados del paso de agregación de resumen anterior.

# Reconvertir dataframe a SpatialPolygons
oregon.shp.agg <- SpatialPolygonsDataFrame(oregon.union, oregon.df.agg)

# Graficar
grid.arrange(spplot(oregon, "AREA", main = "Oregón: área original del condado"), 
             spplot(oregon.shp.agg, "AREA", main = "Oregón: área del condado agregada"), ncol = 1)

Answer 2

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

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

# obtener datos de tidycensus para la demostración (nota que necesitas una clave API, sigue las instrucciones aquí: 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)

# reducir el tamaño del conjunto de datos
census <- census[1:8,]

# crear variable de agrupamiento
grupo_1 <- census$GEOID[1:2]
grupo_2 <- census$GEOID[6:8]

census <- census %>% mutate(grupo = case_when(GEOID %in% grupo_1 ~ 'nuevogrupo1',
                                              GEOID %in% grupo_2 ~ 'nuevogrupo2',
                                              TRUE ~ GEOID))

# resumir por variable de agrupamiento (realiza una unión en polígonos agrupados y suma 'estimate')
census2 <- group_by(census, grupo) %>% 
  summarise(estimate = sum(estimate), do_union = TRUE)

# visualizar usando la versión de desarrollo de ggplot2 y facetar por conjuntos de datos fusionados/no fusionados
datos_grafico <- rbind(census %>% select(grupo, estimate) %>%
                     mutate(facet = "no fusionado"), 
                   census2 %>% mutate(facet = "fusionado"))

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