使用gstat执行IDW插值失败，CRS不匹配问题求解

问题描述 #

执行IDW插值时触发错误：

Error in .local(formula, locations, ...) : sf::st_crs(locations) == sf::st_crs(newdata) is not TRUE

使用的样本数据：

name,long,lat,waterLevel,elevation,depth
EM_01,18.553392,-34.07027,14.4,20.358,63.0
EM_27,18.574777,-34.068709,16.196,19.966,48.0
EM_29,18.613985,-34.053271,18.766,25.477,39.0
EM_20,18.654089,-34.045177,20.102,36.502,45.0
EM_23,18.643495,-34.037468,22.915,32.715,33.0
EM_13,18.637267,-34.029194,25.516,29.716,33.0

错误原因 #

你将栅格对象grid2km转为普通数据框grid2km.df后，数据框仅保留了X/Y坐标数值，丢失了空间参考（CRS）信息。idw函数要求输入的newdata必须是带CRS的空间对象（sf/Spatial类），无法识别普通数据框的空间属性，因此触发CRS不匹配报错。

修正方案 #

以下两种方法均可解决问题，核心是确保插值网格保留CRS信息：

方法1：将栅格转为sf空间对象 #

install.packages("sf")
install.packages("gstat")
install.packages("units")
install.packages("terra")

library(sf) 
library(gstat)
library(units)
library(terra) 

file = 'text.txt'

# 导入数据
cfaq <- read.csv(file, header = 1, sep = ',', dec = '.') 

# 转为sf对象并设置初始CRS
cfaq.sf <- st_as_sf(cfaq, coords=c("long", "lat"), crs = 4326) 

# 转换为本地投影（32734）
cfaq.sf <- st_transform(cfaq.sf, crs = 32734)

# 计算栅格边界
min_x <- floor(min(st_coordinates(cfaq.sf)[,1])/1000)*1000
max_x <- ceiling(max(st_coordinates(cfaq.sf)[,1])/1000)*1000
min_y <- floor(min(st_coordinates(cfaq.sf)[,2])/1000)*1000
max_y <- ceiling(max(st_coordinates(cfaq.sf)[,2])/1000)*1000

# 构建空栅格（用WKT格式CRS更稳定）
grid2km <- rast(
  xmin=min_x, xmax=max_x,
  ymin=min_y, ymax=max_y, 
  crs = st_crs(cfaq.sf)$wkt,
  resolution = 2000, 
  names="waterLevel"
)

# 将栅格转为带CRS的sf点对象
grid_sf <- st_as_sf(grid2km, na.rm = FALSE)

# 执行IDW插值
idw_result <- idw(log(waterLevel) ~ 1, cfaq.sf, grid_sf, idp = 0.5)

# 可选：将插值结果转回栅格，方便可视化
idw_rast <- rasterize(idw_result, grid2km, field = "var1.pred")

方法2：使用terra::interpolate直接插值 #

# 前面数据导入、投影转换步骤和方法1一致

# 构建空栅格
grid2km <- rast(
  xmin=min_x, xmax=max_x,
  ymin=min_y, ymax=max_y, 
  crs = st_crs(cfaq.sf)$wkt,
  resolution = 2000, 
  names="waterLevel"
)

# 创建IDW模型
idw_model <- gstat(formula = log(waterLevel) ~ 1, data = cfaq.sf, idp = 0.5)

# 直接在栅格上执行插值
idw_rast <- interpolate(grid2km, idw_model)

关键注意事项 #

• 不要随意将空间对象转为普通数据框，若必须转换，需手动将其转为带CRS的sf对象
• 优先使用WKT格式的CRS（st_crs(obj)$wkt），proj4string格式已逐渐被弃用
• 插值前可通过st_crs(cfaq.sf) == st_crs(grid_sf)检查样本点与网格的CRS是否一致

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