如何在德国邮编Shapefile中合并NA值区域，保留有数据区域精细粒度？

核心思路 #

将原数据集拆分为有数据的精细邮编区域和无数据的邮编区域两个子集，对无数据子集执行多边形合并，再将合并后的区域与有数据区域整合，最后用leaflet绘制地图。

具体步骤与代码实现 #

1. 加载rgeos包用于多边形合并操作
2. 拆分出有数据（einwohner2非NA）和无数据（einwohner2为NA）的两个子集
3. 对无数据子集执行多边形合并，减少绘图要素数量
4. 将合并后的多边形转换为与原数据格式一致的SpatialPolygonsDataFrame
5. 整合两个数据集并通过leaflet绘图，区分有数据和无数据区域

完整代码如下：

# German postcode shapes 

# Create temp files
temp <- tempfile()
temp2 <- tempfile()

# Download the zip file and save to 'temp' 
URL <- "https://downloads.suche-postleitzahl.org/v2/public/plz-5stellig.shp.zip"
download.file(URL, temp)

# Unzip the contents of the temp and save unzipped content in 'temp2'
unzip(zipfile = temp, exdir = temp2)

# Read shape file 
library(rgdal)
GER_postcode <- readOGR(temp2)

# Create subsample 
library(tidyverse)
library(rgeos) # 用于多边形合并操作

GER_postcode@data$einwohner2 <- ifelse(substr(GER_postcode@data$plz, 1, 1) %in% c("0", "1", "7"), GER_postcode@data$einwohner, NA)

# 拆分数据集：分离有数据和无数据的区域
has_data <- GER_postcode[!is.na(GER_postcode@data$einwohner2), ]
no_data <- GER_postcode[is.na(GER_postcode@data$einwohner2), ]

# 合并所有无数据区域的多边形
no_data_union <- gUnaryUnion(no_data)

# 将合并后的多边形转为SpatialPolygonsDataFrame，添加属性列并匹配投影
no_data_spdf <- SpatialPolygonsDataFrame(
  no_data_union,
  data = data.frame(einwohner2 = NA, row.names = "merged_no_data")
)
proj4string(no_data_spdf) <- proj4string(GER_postcode)

# 整合有数据区域和合并后的无数据区域
combined_data <- rbind(has_data, no_data_spdf)

# 绘制地图
library(leaflet)

# 为有数据区域设置颜色，无数据区域设为灰色
qpal <- colorBin("Reds", combined_data@data$einwohner2, bins=10, na.color = "#cccccc")

leaflet(combined_data) %>%
  addPolygons(
    stroke = TRUE,
    opacity = 1,
    fillOpacity = 0.5,
    smoothFactor = 0.5,
    color = "black",
    fillColor = ~qpal(einwohner2),
    weight = 1
  ) %>%
  addLegend(
    values = ~einwohner2,
    pal = qpal,
    title = "Population",
    na.label = "No Data"
  )

关键说明 #

• gUnaryUnion(no_data)会自动合并所有无数据的邮编多边形，若区域不连通则生成多个合并后的多边形，大幅降低绘图时的要素数量
• 合并后的区域必须转换为SpatialPolygonsDataFrame并匹配原数据的投影，才能与有数据区域成功整合
• 图例通过na.label明确标注无数据区域，提升地图可读性

内容的提问来源于stack exchange，提问作者Marco