R语言气候数据Kriging插值异常问题排查

核心问题原因 #

1. 经纬度范围与网格映射完全颠倒
   站点数据中longitude范围在-16.0_{-15.8之间，`latitude`范围在-5.8}-5.6之间，但代码错误地将经度范围设为-5.8_{-5.6、纬度范围设为-16.1}-15.8，导致插值网格完全脱离站点的空间范围。Kriging在无有效空间关联时，会默认返回全局均值，也就是你看到的单一值。

2. 异常值干扰变异函数拟合
   数据中存在tempmax=1的极端异常点，远低于其他站点20-40℃的正常温度范围，会导致变异函数拟合偏差，甚至无法生成有效的空间结构模型，进一步触发全局均值插值。

3. 变异函数拟合失败
   空间范围不匹配加异常值干扰时，fit.variogram可能拟合出无效参数（如range=0），此时Kriging退化为简单的全局平均插值。

修正步骤与完整代码 #

步骤说明 #

• 纠正经度、纬度的范围对应关系，确保插值网格覆盖站点区域
• 移除或修正极端异常值
• 验证变异函数拟合参数的合理性

修正后的完整代码 #

# Load required packages
library(sp)
library(gstat)
library(raster)
library(ggplot2)

# 使用提供的示例数据（替代read.csv）
weather_data <- structure(list(latitude = c(-5.66761223, -5.72480242, -5.71591633, 
-5.74702767, -5.718042, -5.70906367, -5.72364083, -5.6800395, 
-5.680813, -5.68343283, -5.69915, -5.72150967, -5.72126, -5.7480005, 
-5.68805267, -5.685108, -5.71851317, -5.713897, -5.7104425, -5.720915, 
-5.6980405, -5.74182587, -5.68076091, -5.77276449, -5.6662859, 
-5.74445058), longitude = c(-15.9421786, -15.9560469, -15.9541762, 
-15.9777202, -15.923983, -15.9396147, -15.945323, -15.9431937, 
-15.9449277, -15.9484998, -15.96467, -15.9544655, -15.96869, 
-15.996367, -15.9782928, -15.945571, -15.9378142, -15.9792287, 
-15.9235108, -15.9277867, -15.9728712, -15.9987156, -15.9206547, 
-16.0006466, -15.9819352, -15.9575332), tempmax = c(24.3, 23.1, 
24.4, 24.6, 35.7, 38.1, 34.6, 28.2, 32.6, 35.8, 37.7, 37.1, 38.6, 
36.1, 37.6, 36.2, 36, 41.2, 39.4, 32.9, 22.2, 33.2, 1, 23.7, 
17.5, 23.5)), class = "data.frame", row.names = c(NA, -26L))

# 移除异常值（tempmax=1的站点）
weather_data <- weather_data[weather_data$tempmax != 1, ]

# 设置空间坐标（顺序为longitude+latitude，与数据字段对应）
coordinates(weather_data) <- ~longitude+latitude

# 生成并拟合变异函数，尝试多种模型确保拟合有效
variogram_model <- variogram(tempmax ~ 1, data = weather_data)
fitted_model <- fit.variogram(variogram_model, model = vgm(c("Sph", "Exp", "Gau")))
# 打印拟合结果，检查参数是否合理
print(fitted_model)

# 纠正空间范围：经度（longitude）-16.1到-15.8，纬度（latitude）-5.8到-5.6
lon_min <- -16.1
lon_max <- -15.8
lat_min <- -5.8
lat_max <- -5.6

# 创建插值网格：x对应longitude，y对应latitude
grd <- expand.grid(
  x = seq(lon_min, lon_max, by = 0.005),  # 调整分辨率平衡精度与效率
  y = seq(lat_min, lat_max, by = 0.005)
)
coordinates(grd) <- ~x + y
gridded(grd) <- TRUE

# 执行Kriging插值
kriging_result <- krige(tempmax ~ 1, locations = weather_data, newdata = grd, model = fitted_model)

# 转换为raster并可视化
kriging_raster <- raster(kriging_result)
kriging_df <- as.data.frame(rasterToPoints(kriging_raster))
colnames(kriging_df) <- c("longitude", "latitude", "temperature")

# 绘制插值地图，添加原始站点作为参考
ggplot() +
  geom_raster(data = kriging_df, aes(x = longitude, y = latitude, fill = temperature)) +
  scale_fill_gradient(low = "blue", high = "red", name = "Temperature (°C)") +
  geom_point(data = as.data.frame(weather_data), aes(x = longitude, y = latitude), color = "black", size = 1) +
  labs(title = "Interpolated Temperature Map",
       x = "Longitude", y = "Latitude") +
  theme_minimal()

结果验证 #

运行修正后的代码后，你会看到：

• 插值结果呈现合理的空间温度差异，不再是单一值
• 原始站点温度与插值区域趋势匹配
• 变异函数拟合参数（如range>0、sill接近样本方差）处于合理范围

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