基于R语言ggplot2的气温空间插值地图制作技术问询

一、调试Voronoi插值代码 #

你的初始代码问题出在Voronoi多边形与温度数据的关联逻辑，以及空间对象转换步骤上。修改后可正常读取temperature字段，代码如下：

library(tidyverse)
library(sf)
library(spatstat)
library(maptools)
library(gstat) # 后续IDW、克里金插值依赖

# 读取比利时边界并统一坐标系
belgium_gadm <- st_read("gadm41_BEL.gpkg", layer = "ADM_ADM_0") %>% 
  st_transform(crs = 4326)

# 观测点数据转sf对象
temps_data <- data.frame(
  city = c("Brussels", "Antwerp", "Ghent"), 
  lat = c(50.8503, 51.2194, 51.0543), 
  lon = c(4.3517, 4.4024, 3.7174), 
  temperature = c(15, 16, 17)
) %>% 
  st_as_sf(coords = c("lon", "lat"), crs = 4326)

# 生成Voronoi图
temp_ppp <- as.ppp(temps_data, W = as.owin(st_bbox(belgium_gadm)))
voronoi_result <- dirichlet(temp_ppp)

# 转sf并关联温度（Voronoi多边形与原始观测点索引一一对应）
voronoi_sf <- as(voronoi_result, "sf") %>% 
  mutate(temperature = temps_data$temperature) %>% 
  st_intersection(belgium_gadm) # 裁剪到比利时边界内

# 绘图
ggplot() +
  geom_sf(data = belgium_gadm, fill = NA, color = "black") +
  geom_sf(data = voronoi_sf, aes(fill = temperature), color = "white") +
  coord_sf() +
  theme_minimal()

修改核心点：

• 统一所有空间对象为WGS84坐标系，避免空间运算报错
• 直接通过索引匹配温度数据，替代复杂的over函数关联
• 增加边界裁剪，消除超出比利时范围的Voronoi多边形

二、实现IDW插值 #

使用gstat包的idw函数完成反距离权重插值，代码如下：

# 创建比利时范围内的插值网格（控制cellsize调整精度）
belgium_grid <- st_make_grid(belgium_gadm, cellsize = 0.01, what = "centers") %>% 
  st_sf() %>% 
  st_intersection(belgium_gadm)

# 执行IDW插值（idp为距离权重参数，默认2）
idw_result <- idw(temperature ~ 1, temps_data, newdata = belgium_grid, idp = 2)

# 绘图
ggplot() +
  geom_sf(data = belgium_gadm, fill = NA, color = "black") +
  geom_sf(data = idw_result, aes(fill = var1.pred), color = NA) +
  geom_sf(data = temps_data, color = "red", size = 3) + # 叠加原始观测点
  coord_sf() +
  theme_minimal()

三、实现普通克里金插值 #

先拟合变异函数，再执行克里金插值，代码如下：

# 拟合变异函数（观测点较少时可尝试"Exp"指数模型）
variogram_model <- variogram(temperature ~ 1, data = temps_data)
fit_vgm <- fit.variogram(variogram_model, vgm("Sph")) # 球面模型

# 执行普通克里金插值
krige_result <- krige(temperature ~ 1, temps_data, newdata = belgium_grid, model = fit_vgm)

# 绘图
ggplot() +
  geom_sf(data = belgium_gadm, fill = NA, color = "black") +
  geom_sf(data = krige_result, aes(fill = var1.pred), color = NA) +
  geom_sf(data = temps_data, color = "red", size = 3) +
  coord_sf() +
  theme_minimal()

四、复刻整度/半度气温刻度 #

通过scale_fill系列函数设置刻度，匹配示例图样式，以Voronoi插值为例：

ggplot() +
  geom_sf(data = belgium_gadm, fill = NA, color = "black") +
  geom_sf(data = voronoi_sf, aes(fill = temperature), color = "white") +
  coord_sf() +
  scale_fill_continuous(
    name = "气温(℃)",
    breaks = seq(14, 18, 0.5), # 整度/半度间隔的刻度
    labels = seq(14, 18, 0.5),
    guide = guide_colorbar(ticks = TRUE, barwidth = 15)
  ) +
  theme_minimal()

连续插值（IDW/克里金）可改用scale_fill_steps实现分段着色：

ggplot() +
  geom_sf(data = belgium_gadm, fill = NA, color = "black") +
  geom_sf(data = idw_result, aes(fill = var1.pred), color = NA) +
  geom_sf(data = temps_data, color = "red", size = 3) +
  coord_sf() +
  scale_fill_steps(
    name = "气温(℃)",
    breaks = seq(14, 18, 0.5),
    labels = seq(14, 18, 0.5),
    show.limits = TRUE
  ) +
  theme_minimal()

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