在R中连接多边形

#1 楼

以下解决方案基于Roger Bivand在R-sig-Geo上的帖子。我以他的示例为例，用俄勒冈州的一些人口普查数据替换了德国的shapefile，您可以从此处下载（从“俄勒冈州的县和人口普查数据”中获取所有shapefile组件）。

让我们开始加载所需的程序包和将shapefile导入到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)
 

接下来，您需要一些分组变量才能汇总数据。在我们的示例中，分组仅基于单个县的坐标。参见下图，黑色边框表示原始多边形，而红色边框表示由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)
 

到目前为止，一切都很好。但是，执行unionSpatialPolygons时，与原始shapefile的子区域有关的数据属性（例如，人口密度，面积等）会丢失。我猜您也想汇总与shapefile相关联的普查数据，因此您需要一个中间步骤。

首先必须将多边形转换为数据框以执行聚合。现在，让我们取数据属性列六到八（“ AREA”，“ POP1990”，“ POP1997”），并根据上述ID应用函数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)
的SpatialPolygonsDataFrame，您将获得广义多边形和从上述汇总聚合步骤中获得的普查数据。 
 
 
 oregon.union 
 
 
  

评论
指向俄勒冈州shapefile的链接似乎已死，我很难从俄勒冈州加载其他shapefile



 –彼得·潘（Peter Pan）

 20/07/17在14:22 

#2 楼

这是使用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)