Joindre des polygones dans R

29

Je me demande comment joindre des polygones spatiaux en utilisant le code R?

Je travaille avec des données de recensement où certaines zones changent au fil du temps et je souhaite joindre les polygones et les données correspondantes et simplement signaler les zones jointes. Je tiens à jour une liste de polygones qui ont changé de recensement en recensement et que je prévois de fusionner. J'aimerais utiliser cette liste de noms de régions comme liste de recherche à appliquer aux données du recensement de différentes années.

Je me demande quelle fonction R utiliser pour fusionner les polygones sélectionnés et les données respectives. Je l'ai googlé, mais je suis tout simplement confus par les résultats.

Geoconfused
la source
La réponse à la plupart des opérations de géométrie comme la dissolution de polygones, la superposition, le point dans le polygone, l'intersection, l'union, etc., est le package rgeos.
Spacedman
1
Le US Census Bureau publie des tableaux à cet effet pour 1990-2000 et 2000-2010. Ils peuvent être gérés avec des jointures de base de données , qui sont implémentées par Rla mergefonction de.
whuber

Réponses:

39

La solution suivante est basée sur un article de Roger Bivand sur R-sig-Geo . J'ai pris son exemple en remplaçant le fichier de formes allemand par des données de recensement de l'Oregon que vous pouvez télécharger à partir d' ici (prenez toutes les composantes du fichier de formes des «comtés de l'Oregon et données de recensement»).

Commençons par charger les packages requis et importer le fichier de formes dans 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)

Ensuite, vous avez besoin d'une variable de regroupement afin d'agréger les données. Dans notre exemple, le regroupement est simplement basé sur les coordonnées d'un seul comté. Voir l'image ci-dessous, les bordures noires indiquent les polygones d'origine, tandis que les bordures rouges représentent les polygones agrégés par 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)

Fichier de formes Oregon original et groupé

Jusqu'ici tout va bien. Cependant, les attributs de données liés aux sous-régions du fichier de formes d'origine (par exemple la densité de population, la superficie, etc.) sont perdus lors de l'exécution unionSpatialPolygons. Je suppose que vous souhaitez également agréger vos données de recensement associées au fichier de formes, vous aurez donc besoin d'une étape intermédiaire.

Vous devez d'abord convertir vos polygones en une trame de données afin d'effectuer l'agrégation. Prenons maintenant les colonnes d'attributs de données six à huit ("AREA", "POP1990", "POP1997") et agrégons-les en fonction de la fonction d'application des ID ci-dessus 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)

Enfin, reconvertissez votre trame de données en SpatialPolygonsDataFramefournissant le fichier de formes précédemment unifié oregon.unionet vous obtenez à la fois des polygones généralisés et vos données de recensement dérivées de l'étape d'agrégation de résumé ci-dessus.

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

Régions de l'Oregon

fdetsch
la source
10

Voici une solution utilisant le package 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)

entrez la description de l'image ici

sebdalgarno
la source
Je pensais que j'ajouterais juste un petit avertissement ici, juste au cas où: méfiez-vous de l'utilisation de summarise()dérivés avec l' do_unionargument, car je viens de faire quelque chose comme summarise_if(shapefile, predic.function, sum, na.rm = TRUE, do_union = TRUE), qui finit également par additionner un VRAI dans chaque cellule (c'est-à-dire +1 pour toutes les opérations). Besoin d'enquêter davantage pour savoir si c'est quelque chose qui doit être signalé (au moins pour un avertissement supplémentaire) ...?
stragu