使用 R 从随机森林模型生成预测栅格？

Question

我将随机森林模型拟合到来自 R 中的测试站点的表格数据，现在想使用对应于相同预测变量（例如，坡度、海拔、pH）的栅格数据生成一个显示预测概率值的栅格该模型。

RF 模型旨在使用不同的环境和地球物理数据预测 0/1 二进制变量SITE_NONSITE。

#random forest model
set.seed(321)
rf1 <- randomForest(formula=SITE_NONSITE ~., data=dcc.s.dummy, ntree=500, mtry=10)

dcc.s.dummy 包含以下数据：

str(dcc.s.dummy)
'data.frame':   7899 obs. of  25 variables:
 $ COST_DIST_ECOTONE        : num  -0.232 0.176 -0.443 -0.478 -0.305 ...
 $ COST_DIST_HEA            : num  -0.233 -0.659 -1.055 -0.999 -0.455 ...
 $ COST_DIST_MEDSTR         : num  0.74388 0.63933 0.55964 0.50768 0.00993 ...
 $ COST_DIST_RIV_COAST      : num  0.59 0.63 0.621 0.639 0.617 ...
 $ DEM30_ASP_RE_2           : num  0 0 0 0 1 0 0 0 0 0 ...
 $ DEM30_ASP_RE_3           : num  0 1 0 0 0 0 0 0 1 0 ...
 $ DEM30_ASP_RE_4           : num  1 0 0 0 0 0 0 1 0 0 ...
 $ DEM30_ASP_RE_5           : num  0 0 1 1 0 1 1 0 0 1 ...
 $ DEM30_M                  : num  0.916 0.72 0.499 0.54 1.114 ...
 $ DEM30_SLOPE              : num  0.2063 0.4631 -0.6445 -0.0512 -0.8235 ...
 $ LOC_REL_RE               : num  -0.489 -0.476 -0.476 -0.459 -0.661 ...
 $ LOC_SD_SLOPE             : num  -0.118 -0.135 -0.316 -0.367 -0.57 ...
 $ SSURGO_ESRI_DRAINAGE_RE_2: num  0 0 0 0 0 0 0 0 0 0 ...
 $ SSURGO_ESRI_DRAINAGE_RE_3: num  1 1 1 1 1 1 1 1 1 1 ...
 $ SSURGO_ESRI_DRAINAGE_RE_4: num  0 0 0 0 0 0 0 0 0 0 ...
 $ SSURGO_ESRI_DRAINAGE_RE_5: num  0 0 0 0 0 0 0 0 0 0 ...
 $ SSURGO_ESRI_DRAINAGE_RE_6: num  0 0 0 0 0 0 0 0 0 0 ...
 $ SSURGO_ESRI_EROSION_RE_2 : num  0 0 0 0 0 1 1 0 0 1 ...
 $ SSURGO_ESRI_EROSION_RE_3 : num  1 1 1 0 1 0 0 1 1 0 ...
 $ SSURGO_ESRI_EROSION_RE_4 : num  0 0 0 0 0 0 0 0 0 0 ...
 $ SSURGO_ESRI_LOC_DIV      : num  -0.328 -0.188 -0.157 -0.213 -0.652 ...
 $ SSURGO_ESRI_NATIVEVEG_2  : num  1 1 1 0 1 0 0 1 1 1 ...
 $ SSURGO_ESRI_NATIVEVEG_3  : num  0 0 0 0 0 1 1 0 0 0 ...
 $ SSURGO_PH                : num  0.813 0.059 1.529 2.32 -1.298 ...
 $ SITE_NONSITE             : Factor w/ 2 levels "0","1": 2 2 2 2 2 1 1 2 2 2

然后，我在整个研究区域中获取与这些相同预测变量相对应的栅格，并将它们组合成一个栅格堆栈：

#plot model predictions
COST_DIST_ECOTONE <- raster("cost_dist_ecotone_s.tif.tif")
COST_DIST_HEA <- raster("cost_dist_hea_s.tif.tif")
COST_DIST_MEDSTR <- raster("cost_dist_medstr_s.tif.tif")
COST_DIST_RIV_COAST <- raster("cost_dist_riv_coast_s.tif.tif")
DEM30_ASP_RE_2 <- raster("dem30_asp_rel_2.tif.tif")
DEM30_ASP_RE_3 <- raster("dem30_asp_rel_3.tif.tif")
DEM30_ASP_RE_4 <- raster("dem30_asp_rel_4.tif.tif")
DEM30_ASP_RE_5 <- raster("dem30_asp_rel_5.tif.tif")
DEM30_M <- raster("dem30_m_s.tif.tif")
DEM30_SLOPE <- raster("dem30_slope_s.tif.tif")
LOC_REL_RE <- raster("loc_rel_re_s.tif.tif")
LOC_SD_SLOPE <- raster("loc_sd_slope_s.tif.tif")
SSURGO_ESRI_DRAINAGE_RE_2 <- raster("SSURGO_ESRI_drainage_reclass_nulfill_2.tif.tif")
SSURGO_ESRI_DRAINAGE_RE_3 <- raster("SSURGO_ESRI_drainage_reclass_nulfill_3.tif.tif")
SSURGO_ESRI_DRAINAGE_RE_4 <- raster("SSURGO_ESRI_drainage_reclass_nulfill_4.tif.tif")
SSURGO_ESRI_DRAINAGE_RE_5 <- raster("SSURGO_ESRI_drainage_reclass_nulfill_5.tif.tif")
SSURGO_ESRI_DRAINAGE_RE_6 <- raster("SSURGO_ESRI_drainage_reclass_nulfill_6.tif.tif")
SSURGO_ESRI_EROSION_RE_2 <- raster("SSURGO_ESRI_erosion_reclass_nulfilll_2.tif.tif")
SSURGO_ESRI_EROSION_RE_3 <- raster("SSURGO_ESRI_erosion_reclass_nulfilll_3.tif.tif")
SSURGO_ESRI_EROSION_RE_4 <- raster("SSURGO_ESRI_erosion_reclass_nulfilll_4.tif.tif")
SSURGO_ESRI_LOC_DIV <- raster("SSURGO_ESRI_loc_div_s.tif.tif")
SSURGO_ESRI_NATIVEVEG_2 <- raster("SSURGO_ESRI_nativeveg_nullfill_2.tif.tif")
SSURGO_ESRI_NATIVEVEG_3 <- raster("SSURGO_ESRI_nativeveg_nullfill_3.tif.tif")
SSURGO_PH <- raster("SSURGO_pH_nullfill_s.tif.tif")

ApPl_stack <- stack(COST_DIST_ECOTONE, COST_DIST_HEA, COST_DIST_MEDSTR, COST_DIST_RIV_COAST, DEM30_ASP_RE_2, DEM30_ASP_RE_3, DEM30_ASP_RE_4, DEM30_ASP_RE_5, DEM30_M, DEM30_SLOPE, LOC_REL_RE, LOC_SD_SLOPE, SSURGO_ESRI_DRAINAGE_RE_2, SSURGO_ESRI_DRAINAGE_RE_3, SSURGO_ESRI_DRAINAGE_RE_4, SSURGO_ESRI_DRAINAGE_RE_5, SSURGO_ESRI_DRAINAGE_RE_6, SSURGO_ESRI_EROSION_RE_2, SSURGO_ESRI_EROSION_RE_3, SSURGO_ESRI_EROSION_RE_4, SSURGO_ESRI_LOC_DIV, SSURGO_ESRI_NATIVEVEG_2, SSURGO_ESRI_NATIVEVEG_3, SSURGO_PH)

但是，尝试在raster::predict() 中使用此光栅堆栈ApPl_stack 失败并出现以下错误：

ApPl_prob <- raster::predict(rf1, newdata=ApPl_stack, type="prob")

as.data.frame.default(x[[i]], optional = TRUE) 中的错误：不能 强制类“结构（“RasterLayer”，包=“光栅”）'到一个 数据帧

转换为数据框并改用它会生成此错误：

ApPl_df <- as.data.frame(ApPl_stack, xy=TRUE)
ApPl_prob <- raster::predict(rf1, newdata=ApPl_df, type="prob")

model.frame.default 中的错误（条款，新数据，na.action = na.omit）：
对象不是矩阵另外：警告消息：'newdata'有 658242 行但找到的变量有 754 行

在我的每个预测栅格中都有 658242 个像元和 754 行，这绝非巧合。我在这里想念什么？我觉得其中一个函数需要一个它没有得到的数据类型。

Answer 1

“对象名称”与图层名称无关，因此您需要设置这些名称以匹配用于拟合模型的 data.frame 中的名称。在大多数工作流程中，您会执行类似

的操作

f <- c("cost_dist_ecotone_s.tif.tif", "cost_dist_hea_s.tif.tif", "cost_dist_medstr_s.tif.tif")
s <- stack(f)
names(s) <- gsub(".tif.tif", "", f)

然后从 RasterStack 中提取值以适合您的模型 --- 在这种情况下，名称已经匹配。

但你犯的主要错误在这里

ApPl_prob <- raster::predict(rf1, newdata=ApPl_stack, type="prob")

第一个参数应该是 RasterStack：

ApPl_prob <- raster::predict(ApPl_stack, rf1, type="prob")

或者像你在回答中那样使用命名参数

raster::predict(model=rf1, object=ApPl_stack, type="prob")

Answer 2

在仔细检查了上面代码生成的所有对象的结构后，我发现了问题。无论出于何种原因，stack() 正在将栅格图层的名称改回其原始文件名，而不是我指定的对象名称。最初我没有注意到这个问题，因为plot(ApPl_stack) 显示了我希望出现的名称，即使它们实际上并没有反映在光栅堆栈的结构中。因此，提供给raster::predict() 的堆栈中的栅格名称与随机森林模型中的栅格名称不匹配。

添加一个额外的步骤来分配匹配的名称解决了这个问题：

names(ApPl_stack) <- c("COST_DIST_ECOTONE", "COST_DIST_HEA", "COST_DIST_MEDSTR", "COST_DIST_RIV_COAST", "DEM30_ASP_RE_2", "DEM30_ASP_RE_3", "DEM30_ASP_RE_4", "DEM30_ASP_RE_5", "DEM30_M", "DEM30_SLOPE", "LOC_REL_RE", "LOC_SD_SLOPE", "SSURGO_ESRI_DRAINAGE_RE_2", "SSURGO_ESRI_DRAINAGE_RE_3", "SSURGO_ESRI_DRAINAGE_RE_4", "SSURGO_ESRI_DRAINAGE_RE_5", "SSURGO_ESRI_DRAINAGE_RE_6", "SSURGO_ESRI_EROSION_RE_2", "SSURGO_ESRI_EROSION_RE_3", "SSURGO_ESRI_EROSION_RE_4", "SSURGO_ESRI_LOC_DIV", "SSURGO_ESRI_NATIVEVEG_2", "SSURGO_ESRI_NATIVEVEG_3", "SSURGO_PH")

然后我能够使用以下代码毫无问题地生成和绘制预测：

#plot predictions and save raster to file
ApPl_prob <- 1- raster::predict(model=rf1, object=ApPl_stack, type="prob")
palette <- matlab.like(20)
plot(ApPl_prob, col=palette)
writeRaster(ApPl_prob, "ApPl_prob", format='GTiff')