基于R语言tmap绘制边界内距数据点500米外区域的方法问询

我有一个包含数百个经纬度格式地理点的数据集，目前用tmap包的tm_dots绘制点图层，并叠加在tm_shape绘制的边界图层之上。现在需要实现：绘制出边界图层范围内、所有已绘制点500米范围之外的区域多边形。如果ggplot/ggmap等其他R绘图工具更适合，也可以采用。

# 加载所需包
library(tidyverse)
library(readxl)
library(maptools)
library(classInt)
library(RColorBrewer)
library(sf)
library(tmap)
library(scales)
library(tmaptools)
library(geodata)

# 读取边界多边形数据
shp_name <- "//ims.gov.uk//homedrive//users//JW2002//My Documents//Data//Demography, Mapping & Lookups//Shape Files//East of England//MSOA//Middle_Layer_Super_Output_Areas_December_2011_Generalised_Clipped_Boundaries_in_England_and_Wales.shp"

EofEMSOAs <- st_read(shp_name) %>% 
  st_as_sf()

# 读取 deprivation 数据（仅用于筛选英格兰东部的MSOA）
EofEMSOAsIMD <- read_excel("~/Data/Demography, Mapping & Lookups/IoD/National & EofE IoD 2019/National&IoD 2019 MSOAs.xlsx", 
                           sheet = "East of England MSOAs")

# 筛选英格兰东部的MSOA
EofEMSOAsCodeListOnly <- dplyr::pull(EofEMSOAsIMD, "Area Code")
EofEMSOAsCodeListOnly <- paste(EofEMSOAsCodeListOnly, collapse = '|')

EofEMSOAsFinalList <- EofEMSOAs[grep(EofEMSOAsCodeListOnly,  EofEMSOAs$msoa11cd),]

# 生成点数据
PointData <- read.table(textConnection("ID   Latitude   Longitude
A 52.9742585 0.5526301
B 52.972643 0.8495693
C 52.972643 0.8495693
D 51.46133804 0.36403501"), header=TRUE)

# 转换为空间点数据
PointDataPlotted = st_as_sf(PointData, coords = c('Longitude', 'Latitude'), crs = 4326)

# 创建点缓冲区（原代码此处dist设为5000，且未转换坐标系，单位错误）
PointDataPlotted2 <- PointDataPlotted %>%
  as.data.frame() %>%
  mutate(buffer = st_buffer(geometry, dist = 5000)) %>% 
  select(-geometry) %>% 
  st_as_sf()

# 创建边界合并多边形
union <- st_union(EofEMSOAsFinalList)

# 生成边界框
mask_union <- union %>% as_tibble() %>% 
  mutate(bbox = st_as_sfc(st_bbox(c(xmin = -5.5, xmax = 9, ymax = 51.5, ymin = 42), crs = st_crs(4326)))) %>% 
  st_as_sf() 

# 计算边界框与合并边界的差集（作为遮罩）
diff <- st_difference(mask_union$bbox, mask_union$geometry)

# 绘制地图
OutputMap <- 
  tm_shape(EofEMSOAsFinalList) +
  tm_fill(col = "red") +
  tm_shape(PointDataPlotted2)+
  tm_fill(col = "forestgreen") + 
  tm_shape(diff) +
  tm_fill(col = "white") +
  tm_shape(EofEMSOAsFinalList) +
  tm_borders(col = "white",
             lwd = 1,
             lty = "solid") +
  tm_add_legend(type = "symbol",
                labels = c("Restricted", "Public"),
                col = c("red", "forestgreen"),
                title = "Access type",
                size = 1.5,
                shape = 21)

核心思路是：用边界区域减去所有点的500米缓冲区合并后的区域，得到的就是边界内、点500米范围外的区域。需要注意的是，经纬度坐标系（EPSG:4326）下st_buffer的单位是度，不是米，所以必须先转换为投影坐标系（比如英国常用的EPSG:27700）。

步骤说明 #

1. 转换坐标系：将边界和点数据都转换为投影坐标系，确保缓冲区的单位是米。
2. 创建并合并缓冲区：给每个点创建500米缓冲区，然后合并为一个整体多边形，避免重复计算。
3. 计算差集：用边界区域减去合并后的缓冲区，得到目标区域。
4. 绘制结果：用tmap或ggplot2绘制最终区域。

修改后的代码示例 #

# 加载所需包
library(tidyverse)
library(sf)
library(tmap)

# ---------- 数据预处理（修正坐标系问题）----------
# 定义英国投影坐标系（EPSG:27700，单位为米）
uk_crs <- 27700

# 转换边界数据到投影坐标系
EofEMSOAsFinalList_proj <- EofEMSOAsFinalList %>% 
  st_transform(crs = uk_crs)

# 转换点数据到投影坐标系，并创建500米缓冲区，然后合并所有缓冲区
points_buffer_union <- PointDataPlotted %>% 
  st_transform(crs = uk_crs) %>% 
  st_buffer(dist = 500) %>%  # 这里dist单位是米，符合需求
  st_union() %>% 
  st_sfc() %>% 
  st_as_sf() %>% 
  rename(geometry = x)

# ---------- 计算目标区域：边界内、点500米外的区域 ----------
# 用边界减去合并后的缓冲区，得到目标区域
target_area <- st_difference(EofEMSOAsFinalList_proj, points_buffer_union)

# ---------- 绘制地图（用tmap）----------
# 转换回WGS84坐标系（可选，用于适配tmap默认显示）
target_area_wgs84 <- target_area %>% st_transform(crs = 4326)
EofEMSOAsFinalList_wgs84 <- EofEMSOAsFinalList %>% st_transform(crs = 4326)
points_buffer_union_wgs84 <- points_buffer_union %>% st_transform(crs = 4326)

OutputMap_final <- 
  # 绘制边界内的目标区域（点500米外的部分）
  tm_shape(target_area_wgs84) +
  tm_fill(col = "#f0f0f0", title = "区域类型") +
  # 绘制点的500米缓冲区
  tm_shape(points_buffer_union_wgs84) +
  tm_fill(col = "forestgreen", alpha = 0.6, title = "") +
  # 绘制边界线
  tm_shape(EofEMSOAsFinalList_wgs84) +
  tm_borders(col = "black", lwd = 0.8) +
  # 绘制原始点
  tm_shape(PointDataPlotted) +
  tm_dots(col = "red", size = 0.5, title = "") +
  # 添加图例
  tm_add_legend(type = "fill",
                labels = c("边界内、点500米外区域", "点500米缓冲区"),
                col = c("#f0f0f0", "forestgreen"),
                title = "区域类型") +
  tm_add_legend(type = "symbol",
                labels = c("原始点"),
                col = "red",
                shape = 20,
                size = 0.5)

# 显示地图
print(OutputMap_final)

# ---------- 用ggplot2绘制的备选方案 ----------
library(ggplot2)

ggplot() +
  # 绘制目标区域
  geom_sf(data = target_area_wgs84, fill = "#f0f0f0", color = NA) +
  # 绘制缓冲区
  geom_sf(data = points_buffer_union_wgs84, fill = "forestgreen", alpha = 0.6, color = NA) +
  # 绘制边界
  geom_sf(data = EofEMSOAsFinalList_wgs84, fill = NA, color = "black", linewidth = 0.8) +
  # 绘制点
  geom_sf(data = PointDataPlotted, color = "red", size = 2) +
  # 添加图例和主题
  labs(fill = "区域类型") +
  scale_fill_manual(values = c("#f0f0f0" = "边界内、点500米外区域", "forestgreen" = "点500米缓冲区")) +
  theme_minimal()

关键修正点 #

• 增加了坐标系转换：从WGS84（EPSG:4326）转换为英国投影坐标系（EPSG:27700），确保缓冲区的单位是米。
• 修正了缓冲区距离：把原代码的5000改为500，符合500米的需求。
• 合并了所有点的缓冲区：避免多个缓冲区重叠导致差集计算错误。
• 直接计算边界与缓冲区的差集：得到精确的目标区域。

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