如何将GAM模型预测结果正确映射至全球栅格文件？

问题背景 #

构建物种分布模型（SDM）时，针对物种A建立了广义相加模型（GAM），代码如下：

response<- SDM_data$presence
predictors <- c("bioclim.bio2","bioclim.bio3", "bioclim.bio8","bioclim.bio9","bioclim.bio13","bioclim.bio14","bioclim.bio15","bioclim.bio18","bioclim.bio19", "Elevation", "Evapotranspiration", "Slope","Roads")
formula <- as.formula(paste("presence ~", paste0("s(", predictors, ", k = 5)", collapse = "+")))

GAM_model <- gam(formula, data = SDM_data, family = binomial, select = TRUE)

其中presence为二元因子变量（1=存在，0=背景点），预测因子包含海拔、坡度及多种生物气候变量，数据覆盖全球范围。为避免过拟合，将平滑项节点数设为5。

后续尝试将预测结果转为栅格时，出现数据长度不匹配的警告：

In matrix(prediction, nrow = nrow(prediction_raster), ncol = ncol(prediction_raster),  :
data length [175565] is not a sub-multiple or multiple of the number of rows [1680]

最终生成的栅格映射完全不准确。

错误原因分析 #

1. 数据长度不匹配：目标栅格的总单元格数为nrow(prediction_raster)*ncol(prediction_raster)=1680*4308=7,237,440，但prediction的长度仅为175565，说明env_subset并非完整的全球栅格所有单元格数据，而是部分点数据。
2. 手动转矩阵的逻辑错误：假设prediction的顺序与栅格单元格的行优先顺序完全一致，但实际env_subset的点顺序大概率和栅格单元格顺序不匹配，强制转矩阵会导致预测值被错误填充，甚至循环重复数据。

解决方案 #

情况1：env是栅格栈/栅格砖（包含所有预测因子的栅格图层） #

直接使用predict函数对栅格对象预测，自动生成匹配的栅格结果，无需手动处理矩阵：

# 确保env是包含所有predictors的栅格栈/砖
prediction_raster <- predict(GAM_model, newdata = env, type = "response", na.action = "na.pass")
plot(prediction_raster)

情况2：env是带坐标的点数据框 #

如果env仅包含部分点的预测因子数据（如示例中的175565行），需要通过坐标匹配栅格单元格，将预测值填充到对应位置：

1. 确认env数据框包含坐标列（如x和y，对应经度、纬度）
2. 按以下步骤填充栅格：

# 创建目标空栅格
world_extent <- extent(-180, 179, -56, 84)
prediction_raster <- raster(world_extent, resolution = 0.08333333)

# 根据坐标获取每个点对应的栅格单元格索引
cell_indices <- cellFromXY(prediction_raster, env[, c("x", "y")])

# 将预测值填充到对应单元格
prediction_raster[cell_indices] <- prediction

# 绘制验证
plot(prediction_raster)

关键验证步骤 #

• 检查env的类型：用class(env)确认是RasterStack/RasterBrick还是data.frame
• 核对单元格索引：确保length(cell_indices)等于length(prediction)，避免出现NA索引
• 验证坐标范围：确认env的坐标范围与prediction_raster的范围一致，避免点落在栅格外导致索引为NA

内容的提问来源于stack exchange，提问作者Isma Soto Almena