R语言中目标变量与辅助变量测点不同时的Co-kriging数据帧制备问题

别担心，第一次用gstat做协同克里金遇到位置不对齐的问题太正常了！我来一步步帮你搞定数据准备的事儿：

协同克里金本来就支持目标变量和辅助变量在不同位置观测，核心就是把所有观测点整合到一个空间数据集里，让每个位置对应的变量值（有观测的话）都正确记录，没观测的就留作NA就行。

步骤1：先确认坐标系统完全一致 #

这是最容易踩坑的第一步！你的目标变量数据框voi_gs和辅助变量的数据框（假设叫aux_gs，得包含x、y坐标和辅助变量列，比如z_aux）必须用一模一样的坐标系统，不然后续合并全白搭。用sf包来检查和转换很方便：

library(sf)
# 把目标变量数据转成sf空间对象，记得替换成你的坐标系EPSG码（比如32633）
voi_sf <- st_as_sf(voi_gs, coords = c("x", "y"), crs = 32633)
# 辅助变量数据也转成sf对象，CRS要和上面完全一致
aux_sf <- st_as_sf(aux_gs, coords = c("x", "y"), crs = 32633)

步骤2：整合所有空间观测点 #

我们要把两个数据集的所有位置合并，得到一个包含所有观测点的完整数据集，每个点对应的z_ar（目标）和z_aux（辅助）会在有观测的位置显示数值，没观测的位置就是NA：

# 先把目标变量的点和辅助变量的点做精确匹配合并，保留所有目标点
combined_sf <- st_join(voi_sf, aux_sf, join = st_equals, left = TRUE)
# 再反向合并一次，确保辅助变量的独特点也被包含进来
combined_sf <- st_join(combined_sf, aux_sf, join = st_equals, right = TRUE)
# 去重，避免出现重复的空间点
combined_sf <- combined_sf[!duplicated(st_geometry(combined_sf)), ]

如果你的坐标有微小误差（比如测量精度问题），没法精确匹配，可以换成st_nearest_join来匹配最近的点，但优先用st_equals保证精确性。

要是你习惯用传统的sp包，也可以这么做：

library(sp)
# 转换为SpatialPointsDataFrame
voi_sp <- SpatialPointsDataFrame(
  coords = voi_gs[, c("x", "y")], 
  data = voi_gs[, "z_ar", drop=FALSE], 
  proj4string = CRS("+init=epsg:32633") # 替换成你的坐标系
)
aux_sp <- SpatialPointsDataFrame(
  coords = aux_gs[, c("x", "y")], 
  data = aux_gs[, "z_aux", drop=FALSE], 
  proj4string = CRS("+init=epsg:32633")
)

# 合并两个空间对象
combined_sp <- spRbind(voi_sp, aux_sp)
# 去除重复点
combined_sp <- combined_sp[!duplicated(coordinates(combined_sp)), ]

步骤3：检查整理好的数据 #

现在combined_sf（或combined_sp）就是gstat需要的格式了！你可以用head(combined_sf)看看结果，应该是这样的：

x       y     z_ar     z_aux

8974 312216.6 530439.8 49.03470 NA
8283 312084.6 530559.8 57.15355 NA
...（辅助变量的观测点会有z_aux的值，z_ar列是NA）

步骤4：定义协同变异函数并跑Co-kriging #

数据准备好后，就可以定义变异函数模型然后执行预测了：

library(gstat)
# 先定义协同变异函数的初始模型：包含目标变量、辅助变量的变异函数，以及交叉变异函数
vgm_init <- vgm(list(
  vgm("Sph", psill = 10, range = 1000, nugget = 2),  # z_ar的变异函数
  vgm("Sph", psill = 8, range = 1200, nugget = 1),   # z_aux的变异函数
  vgm("Sph", psill = 5, range = 1000, nugget = 0)    # 交叉变异函数
))

# 用拟合好的协同变异函数模型
fit_vgm <- fit.lmc( ~ z_ar + z_aux, combined_sf, vgm_init)

# 假设你已经有了预测网格（比如叫grid_sf，是包含x,y的sf对象）
ck_pred <- krige(z_ar ~ z_aux, combined_sf, grid_sf, model = fit_vgm)

小提醒：交叉变异函数的参数（比如range）最好和两个变量的变异函数range接近，或者让fit.lmc自动帮你拟合最优参数。

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