Łączenie wielokątów w R.

29

Zastanawiam się, jak połączyć wielokąty przestrzenne za pomocą kodu R?

Pracuję z danymi spisu, w których niektóre obszary zmieniają się w czasie i chcę dołączyć do wielokątów i odpowiednich danych, a po prostu raportować o połączonych obszarach. Utrzymuję listę wielokątów, które zmieniły spis na spis i które planuję scalić. Chciałbym wykorzystać tę listę nazw obszarów jako listę odnośników do danych spisu ludności z różnych lat.

Zastanawiam się, jakiej funkcji R użyć do scalenia wybranych wielokątów i odpowiednich danych. Poszukałem go, ale po prostu jestem zdezorientowany wynikami.

Geoconfused
źródło
Odpowiedzią na większość operacji geometrycznych, takich jak rozpuszczanie wielokątów, nakładanie, wielokąt punktowy, przecięcie, łączenie itp., Jest pakiet rgeos.
Spacedman
1
Biuro Spisu Powszechnego USA publikuje tabele, aby to zrobić dla lat 1990-2000 i 2000-2010. Mogą one być zarządzane w bazie łączy, które są realizowane przez R„s mergefunkcji.
whuber

Odpowiedzi:

39

Poniższe rozwiązanie jest oparte na poście Rogera Bivanda na temat R-sig-Geo . Wziąłem jego przykład, zastępując niemiecki plik kształtów niektórymi danymi spisu ludności z Oregonu, które można pobrać stąd (weź wszystkie składniki pliku kształtu z „hrabstw i danych spisu stanu Oregon”).

Zacznijmy od załadowania wymaganych pakietów i importowania pliku shapefile do 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)

Następnie potrzebujesz jakiejś zmiennej grupującej, aby agregować dane. W naszym przykładzie grupowanie opiera się po prostu na współrzędnych jednego hrabstwa. Zobacz obrazek poniżej, czarne obramowania wskazują oryginalne wielokąty, podczas gdy czerwone obramowania reprezentują wielokąty agregowane przez 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)

Oryginalny i zgrupowany plik kształtu Oregon

Jak na razie dobrze. Jednak atrybuty danych związane z podregionami pierwotnego pliku kształtu (np. Gęstość zaludnienia, obszar itp.) Gubią się podczas wykonywania unionSpatialPolygons. Wydaje mi się, że chciałbyś również zagregować swoje dane spisowe związane z plikiem kształtu, więc potrzebujesz pośredniego kroku.

Najpierw musisz przekonwertować wielokąty na ramkę danych, aby przeprowadzić agregację. Teraz weźmy kolumny atrybutów danych od szóstego do ósmego („OBSZAR”, „POP1990”, „POP1997”) i agregujemy je zgodnie z powyższymi funkcjami stosującymi identyfikatory 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)

Na koniec przekonwertuj ramkę danych z powrotem na SpatialPolygonsDataFramezapewniający poprzednio ujednolicony plik kształtu, oregon.uniona otrzymasz zarówno uogólnione wielokąty, jak i dane ze spisu pochodzące z powyższego kroku agregacji podsumowania.

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

Obszary Oregon

fdetsch
źródło
10

Oto rozwiązanie wykorzystujące pakiet 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)

wprowadź opis zdjęcia tutaj

sebdalgarno
źródło
Pomyślałem, że dodam tutaj małe ostrzeżenie, na wszelki wypadek: strzeż się używania summarise()pochodnych z do_unionargumentem, ponieważ właśnie zrobiłem coś takiego summarise_if(shapefile, predic.function, sum, na.rm = TRUE, do_union = TRUE), co zakończyło się sumowaniem PRAWDA w każdej komórce (tj. +1 dla wszystkich operacji). Potrzebujesz dowiedzieć się więcej, aby dowiedzieć się, czy należy to zgłosić (przynajmniej w celu uzyskania dodatkowego ostrzeżenia) ...?
stragu