R语言绘制空间地图时实现聚类区域分色填充的方法咨询

你需要的按聚类类别填充对应区域效果，核心是基于带聚类标签的点位生成空间分区（泰森多边形/插值面），将分区裁剪到研究区行政边界范围内后按类别填色，以下分别提供你给出的两个绘图框架的可运行修改代码。

ggplot2 版本实现 #

用sf包原生支持的泰森多边形工具生成分区，无需额外加载空间插值包，代码逻辑更简洁：

library(ggplot2)
library(sf)
library(dplyr)

# 读取CASTRO市行政边界（和原代码逻辑一致）
temp <- tempfile()
temp2 <- tempfile()
download.file("https://geoftp.ibge.gov.br/organizacao_do_territorio/malhas_territoriais/malhas_municipais/municipio_2015/UFs/PR/pr_municipios.zip",temp)
unzip(zipfile = temp, exdir = temp2)
shp <- sf::read_sf(temp2)
shp_subset <- shp[shp$NM_MUNICIP == "CASTRO",]

# 原点位数据集
Points_properties<-structure(list(Latitude = c(-24.781624, -24.775017, -24.769196, 
                                               -24.761741, -24.752019, -24.748008, -24.737312, -24.744718, -24.751996, 
                                               -24.724589, -24.8004, -24.796899, -24.795041, -24.780501, -24.763376, 
                                               -24.801715, -24.728005, -24.737845, -24.743485, -24.742601, -24.766422, 
                                               -24.767525, -24.775631, -24.792703, -24.790994, -24.787275, -24.795902, 
                                               -24.785587, -24.787558, -24.799524), Longitude = c(-49.937369, 
                                               -49.950576, -49.927608, -49.92762, -49.920608, -49.927707, -49.922095, 
                                               -49.915438, -49.910843, -49.899478, -49.901775, -49.89364, -49.925657, 
                                               -49.893193, -49.94081, -49.911967, -49.893358, -49.903904, -49.906435, 
                                               -49.927951, -49.939603, -49.941541, -49.94455, -49.929797, -49.92141, 
                                               -49.915141, -49.91042, -49.904772, -49.894034, -49.86651), 
                                               cluster = c("1", "1", "1", "1", "2", "2", "2", "2", "2", "2", "1", "1", "1", "1", "1", 
                                               "1", "2", "2", "2", "2", "1", "1", "1", "1", "1", "1", "1", "1", 
                                               "1", "1")), row.names = c(NA, -30L), class = "data.frame")

# 核心处理：点位转空间对象 -> 生成泰森多边形 -> 裁剪到行政边界 -> 关联聚类属性
points_sf <- st_as_sf(Points_properties, coords = c("Longitude", "Latitude"), crs = st_crs(shp_subset))
voronoi <- st_voronoi(st_union(points_sf)) %>% 
  st_collection_extract(type = "POLYGON") %>% 
  st_intersection(st_union(shp_subset))
voronoi_cluster <- st_join(st_sf(geometry = voronoi), points_sf, join = st_nearest_feature)

# 绘图
ggplot() +
  # 聚类区域填色层
  geom_sf(data = voronoi_cluster, aes(fill = cluster), alpha = 0.6) +
  # 行政边界轮廓
  geom_sf(data = shp_subset, fill = NA, linewidth = 0.8, color = "black") +
  # 原始点位叠加
  geom_sf(data = points_sf, aes(color = cluster), size = 2) +
  # 配色设置：聚类1为绿色，聚类2为蓝色，和参考效果一致
  scale_fill_manual(values = c("1" = "forestgreen", "2" = "steelblue")) +
  scale_color_manual(values = c("1" = "darkgreen", "2" = "darkblue")) +
  # 坐标范围和原代码保持一致
  coord_sf(xlim = c(min(Points_properties$Longitude)-0.1,
                    max(Points_properties$Longitude)+0.1),
           ylim = c(min(Points_properties$Latitude)-0.1,
                    max(Points_properties$Latitude)+0.1),
           expand = FALSE) +
  theme_void()

raster 版本实现 #

基于dismo包的泰森多边形工具生成分区，转栅格后掩膜裁剪，适配raster绘图框架：

library(rgdal)
library(sf)
library(raster)
library(dplyr)
library(dismo)

# 读取CASTRO市行政边界（和原代码逻辑一致）
temp <- tempfile()
temp2 <- tempfile()
download.file("https://geoftp.ibge.gov.br/organizacao_do_territorio/malhas_territoriais/malhas_municipais/municipio_2015/UFs/PR/pr_municipios.zip",temp)
unzip(zipfile = temp, exdir = temp2)
shp <- readOGR(temp2)
shp_subset <- shp[shp$NM_MUNICIP == "CASTRO",]

# 原点位数据集
Points_properties<-structure(list(Latitude = c(-24.781624, -24.775017, -24.769196, 
     -24.761741, -24.752019, -24.748008, -24.737312, -24.744718, -24.751996, 
     -24.724589, -24.8004, -24.796899, -24.795041, -24.780501, -24.763376, 
     -24.801715, -24.728005, -24.737845, -24.743485, -24.742601, -24.766422, 
     -24.767525, -24.775631, -24.792703, -24.790994, -24.787275, -24.795902, 
     -24.785587, -24.787558, -24.799524), Longitude = c(-49.937369, 
     -49.950576, -49.927608, -49.92762, -49.920608, -49.927707, -49.922095, 
     -49.915438, -49.910843, -49.899478, -49.901775, -49.89364, -49.925657, 
     -49.893193, -49.94081, -49.911967, -49.893358, -49.903904, -49.906435, 
     -49.927951, -49.939603, -49.941541, -49.94455, -49.929797, -49.92141, 
     -49.915141, -49.91042, -49.904772, -49.894034, -49.86651), cluster = c("1", "1", 
     "1", "1", "2", "2", "2", "2", "2", "2", "1", "1", "1", "1", "1", 
     "1", "2", "2", "2", "2", "1", "1", "1", "1", "1", "1", "1", "1", 
     "1", "1")), row.names = c(NA, -30L), class = c("tbl_df", "tbl", 
     "data.frame"))

# 核心处理：点位转空间对象 -> 生成泰森多边形 -> 转栅格 -> 按行政边界掩膜
coordinates(Points_properties) <- ~Longitude+Latitude
proj4string(Points_properties) <- proj4string(shp_subset)
vor_sp <- voronoi(Points_properties, shp_subset)
r <- raster(vor_sp, res = 0.001)
vor_rast <- rasterize(vor_sp, r, field = "cluster")
vor_rast_mask <- mask(vor_rast, shp_subset)

# 绘图
marsize <- .2
par(mar = rep(marsize,4))
z <- .4
shp_subset_cropped <- raster::crop(shp_subset,
                                   extent(c(
                                     (extent(shp_subset)@"xmin"*(1-z)+z*min(coordinates(Points_properties)[,1])),
                                     (extent(shp_subset)@"xmax"*(1-z)+z*max(coordinates(Points_properties)[,1])),
                                     (extent(shp_subset)@"ymin"*(1-z)+z*min(coordinates(Points_properties)[,2])),
                                     (extent(shp_subset)@"ymax"*(1-z)+z*max(coordinates(Points_properties)[,2]))
                                   )))
# 聚类区域填色层
plot(vor_rast_mask, col = c("1" = "forestgreen", "2" = "steelblue"), alpha = 0.6, legend = FALSE)
# 行政边界轮廓
plot(shp_subset_cropped, add = TRUE)
# 原始点位叠加
pp1 <- as.data.frame(Points_properties) %>% filter(cluster==1)
pp2 <- as.data.frame(Points_properties) %>% filter(cluster==2)
points(x = pp1$Longitude, y= pp1$Latitude, col = "darkgreen",pch = 16,cex=1)
points(x = pp2$Longitude, y= pp2$Latitude, col = "darkblue",pch = 16,cex=1)

补充说明 #

如果需要平滑过渡的分区边界而非泰森多边形的硬折线，可以将泰森多边形替换为IDW反距离加权插值或克里金插值，生成连续的聚类概率栅格后，按最大概率归属类别划分填色即可。

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