如何提取GWR模型参数并应用到精细分辨率栅格？

问题背景 #

我拥有两个栅格图层：粗分辨率的ntl和精细分辨率的tirs。目标是提取地理加权回归(GWR)的**截距(intercept)和斜率(slope)**参数，并将其应用到精细分辨率栅格tirs上，实现类似线性回归(LM)的计算逻辑：GWR截距 + GWR斜率 × 精细分辨率栅格。

线性回归(LM)的实现示例 #

线性回归中参数为全局固定值，可直接计算，示例代码如下：

library(terra)
library(sp)

# 读取栅格
tirs = rast("path/tirs.tif") # 精细分辨率栅格
ntl = rast("path/ntl.tif")   # 粗分辨率栅格
    
# 填充空值
tirs = focal(tirs, 
             w = 9, 
             fun = mean, 
             na.policy = "only", 
             na.rm = TRUE)
  
# 高斯滤波
gf <- focalMat(tirs, 0.10*400, "Gauss", 11)
r_gf <- focal(tirs, w = gf, na.rm = TRUE)
  
# 重采样到粗分辨率
r_gf = resample(r_gf, ntl, method = "bilinear")

# 组合数据集
s = c(ntl, r_gf)
names(s) = c('ntl', 'r_gf')

# 构建线性模型
model <- lm(formula = ntl ~ tirs, data = s)

# 将系数应用到精细分辨率栅格
lm_pred = model$coefficients[1] + model$coefficients[2] * tirs

GWR的问题与现有代码 #

GWR的参数随空间变化，而非全局固定值，其系数摘要如下：

GWR系数估计值摘要：

Min.     1st Qu.      Median     3rd Qu.     Max.

Intercept -1632.61196   -55.79680   -15.99683    15.01596 1133.299

tirs20      -42.43020     0.43446     1.80026     3.75802   70.987

现有GWR运行代码：

library(GWmodel)
library(raster)

block.data = read.csv(file = "path/block.data00.csv")

# 提取坐标矩阵
x = as.data.frame(block.data$x)
y = as.data.frame(block.data$y)
sint = as.matrix(cbind(x, y))

# 转换为空间点数据框
coordinates(block.data) = c("x", "y")

# 模型公式
eq1 <- ntl ~ tirs

# 计算距离矩阵
dist = GWmodel::gw.dist(dp.locat = sint, focus = 0, longlat = FALSE)

# 带宽选择
abw = bw.gwr(eq1, 
       data = block.data, 
       approach = "AIC", 
       kernel = "tricube",
       adaptive = TRUE, 
       p = 2, 
       longlat = F, 
       dMat = dist,
       parallel.method = "omp",
       parallel.arg = "omp")

# 运行GWR模型
ab_gwr = gwr.basic(eq1, 
          data = block.data, 
          bw = abw, 
          kernel = "tricube",
          adaptive = TRUE, 
          p = 2,
          longlat = FALSE, 
          dMat = dist,
          F123.test = FALSE,
          cv = FALSE,
          parallel.method = "omp",
          parallel.arg = "omp")

ab_gwr

补充信息 #

需应用GWR系数的精细分辨率栅格参数：

tirs = raster(ncols=407, nrows=342, xmn=509600, xmx=550300, ymn=161800, ymx=196000, crs='+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs')

解决方案 #

步骤1：提取GWR的空间化参数 #

GWR模型结果ab_gwr的SDF对象包含每个采样点的坐标及对应截距、斜率参数，先将其转换为栅格格式：

# 提取GWR系数的空间数据框
gwr_coeffs = ab_gwr$SDF

# 将系数转换为栅格（匹配粗分辨率栅格的范围与分辨率）
intercept_rast = rasterFromXYZ(data.frame(gwr_coeffs@coords, intercept = gwr_coeffs$Intercept))
slope_rast = rasterFromXYZ(data.frame(gwr_coeffs@coords, slope = gwr_coeffs$tirs))

# 确保栅格CRS与目标精细栅格一致
crs(intercept_rast) = crs(tirs)
crs(slope_rast) = crs(tirs)

步骤2：将GWR参数重采样到精细分辨率 #

把粗分辨率的系数栅格重采样到与tirs相同的精细分辨率：

# 重采样截距和斜率栅格到精细分辨率
intercept_fine = resample(intercept_rast, tirs, method = "bilinear")
slope_fine = resample(slope_rast, tirs, method = "bilinear")

步骤3：应用GWR参数计算精细分辨率结果 #

按照截距 + 斜率 × 精细栅格的逻辑计算最终结果：

# 计算GWR预测的精细分辨率栅格
gwr_pred = intercept_fine + slope_fine * tirs

关键说明 #

• 若使用自适应带宽的GWR，参数仅在采样点存在，需通过空间插值（如双线性插值）将参数扩展到整个研究区的精细栅格。
• 确保所有栅格的坐标系一致，避免投影错误。

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