Terra栅格化插值输出出现横竖线问题求助

问题描述 #

我尝试用包含经度、纬度和响应值（mud）的点数据生成连续插值表面，测试了Inverse Distance Weighting（IDW）和Ordinary Kriging两种方法，但最终输出栅格都出现规则分布的水平和垂直白线，无法得到连续表面。已对interp_rast对象使用resample()函数，没有改善。怀疑grd与out_df对象不匹配，但叠加放大后目视显示正常。

可复现代码 #

library(terra)
library(sp)
library(gstat)

extent <- c(-7.67151099570083, -2.34424614123549, 50.159278977608, 56.2717492369751)
blankRaster <- rast(ext(extent))
cell_size <- c(0.00117, 0.002075) * 10
res(blankRaster) <- cell_size

x_range <- c(ext(blankRaster)[1], ext(blankRaster)[2])
y_range <- c(ext(blankRaster)[3], ext(blankRaster)[4])

grd <- expand.grid(x = seq(from = x_range[1], to = x_range[2], by = cell_size[1]), y = seq(from = y_range[1], to = y_range[2], by = cell_size[2]))
coordinates(grd) <- ~ x+y
gridded(grd) <- T

lon <- c(-6.054578, -6.030297, -6.148603, -6.001423, -6.029923, -6.041965, -6.026100, -6.013715, -6.027235, -6.026288, -6.004180, -6.020842, -5.997938, -6.030913, -6.025397, -6.107727, -6.089913, -6.054938, -6.042965, -6.072297)

lat <- c(53.93057, 53.93222, 53.83308, 53.57506, 53.57339, 53.57612, 53.56305, 53.54610, 53.60224, 53.61145, 53.52641, 53.51799, 53.51425, 53.72778, 53.76036, 53.90055, 53.89335, 53.90751, 53.90092, 53.90038)

mud <- c(0.032, 0.039, 0.126, 0.146, 0.536, 0.225, 0.222, 0.408, 0.145, 0.112, 0.241, 0.031, 0.186, 0.340, 0.074, 0.162, 0.379, 0.147, 0.482, 0.220)

df <- data.frame(lon, lat, mud)

obs_rast <- rasterize(x = vect(df), y = blankRaster, field = "mud", fun = median)
obs_rast <- focal(obs_rast, w = 9, fun = median, na.policy="only", na.rm = T)
obs_rast <- resample(obs_rast, blankRaster)
out_df <- as.data.frame(obs_rast, xy = T)
names(out_df) <- c("lon", "lat", "mud")
coordinates(out_df) <- c("lon", "lat")

#IDW interpolation for each interp_df   
idw <- idw(formula = mud ~ 1, locations = out_df, newdata = grd)
idw_output <- as.data.frame(idw)
names(idw_output)[1:3] <- c("lon", "lat", "mud")
interp_output <- idw_output[, 1:3]

interp_rast <- terra::rasterize(x=vect(interp_output), y=blankRaster, field='mud', fun=median)
interp_rast <- resample(interp_rast, blankRaster)
plot(interp_rast)

输出效果 #

叠加放大验证代码 #

plot(interp_rast, xlim = c(-6.4, -5.8), ylim = c(53.2, 54.1))
plot(grd, add = T, pch = 20)
plot(out_df, add = T, pch = 20, col = "red")

叠加放大效果 #

问题根源与解决方法 #

问题根源 #

1. 多余的栅格化处理：先将原始点数据df栅格化为obs_rast再做focal处理，会生成大量含NA值的栅格点，用这些稀疏点插值时，NA值会导致插值结果出现规则空白线。
2. 栅格对齐冗余操作：手动构建grd后再通过rasterize转换为栅格，容易出现细微坐标对齐误差，加剧白线问题。

修正方案 #

直接使用原始点数据插值，简化栅格构建流程，避免引入NA值和对齐误差：

修正后的IDW插值代码

library(terra)
library(gstat)

# 定义目标栅格参数
extent <- c(-7.67151099570083, -2.34424614123549, 50.159278977608, 56.2717492369751)
cell_size <- c(0.00117, 0.002075) * 10
target_rast <- rast(ext(extent), res = cell_size)

# 原始点数据转为矢量对象
lon <- c(-6.054578, -6.030297, -6.148603, -6.001423, -6.029923, -6.041965, -6.026100, -6.013715, -6.027235, -6.026288, -6.004180, -6.020842, -5.997938, -6.030913, -6.025397, -6.107727, -6.089913, -6.054938, -6.042965, -6.072297)
lat <- c(53.93057, 53.93222, 53.83308, 53.57506, 53.57339, 53.57612, 53.56305, 53.54610, 53.60224, 53.61145, 53.52641, 53.51799, 53.51425, 53.72778, 53.76036, 53.90055, 53.89335, 53.90751, 53.90092, 53.90038)
mud <- c(0.032, 0.039, 0.126, 0.146, 0.536, 0.225, 0.222, 0.408, 0.145, 0.112, 0.241, 0.031, 0.186, 0.340, 0.074, 0.162, 0.379, 0.147, 0.482, 0.220)
df <- data.frame(lon, lat, mud)
df_vect <- vect(df, geom = c("lon", "lat"))

# 直接基于原始点做IDW插值
idw_result <- idw(mud ~ 1, locations = df_vect, newdata = target_rast)
# 转换为terra栅格对象
interp_rast <- rast(idw_result)

plot(interp_rast)

普通克里金（Ordinary Kriging）修正代码

# 拟合变异函数
vgm_model <- variogram(mud ~ 1, df_vect) %>% fit.variogram(model = vgm("Sph"))
# 克里金插值
krig_result <- krige(mud ~ 1, locations = df_vect, newdata = target_rast, model = vgm_model)
krig_rast <- rast(krig_result)

plot(krig_rast)

补充说明 #

如果确实需要对原始点做栅格化预处理，需先过滤掉NA值：

# 栅格化后过滤NA值
obs_rast <- rasterize(x = vect(df), y = target_rast, field = "mud", fun = median)
out_df <- na.omit(as.data.frame(obs_rast, xy = T))
coordinates(out_df) <- c("x", "y")
# 再用过滤后的点做插值
idw_result <- idw(mud ~ 1, locations = out_df, newdata = target_rast)

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