gbm現在、多項式モデルをラスターブリックに予測することはできないようです。ただし、比較的小さなラスター グリッドの場合は、これを回避する簡単な方法があることを確認してください。これについては以下で説明します。しかし、ここでのプロセスは非常に遅く、大規模なラスター、多くのクラス (私の場合は植生群集)、および予測変数を扱う場合、課題がないわけではありません。以下の情報が、同じ課題に直面する人の役に立てば幸いです。
以下では、多項 gbm モデルと 20 の予測変数を使用して、36 の植生群落の発生確率を予測しようとしています。私の研究領域は、213,000,000 ピクセルの 30x30m ラスター グリッドですが、以下のコードは、プロセスの開発/テストに使用した 1221 セルの小さなグリッドに関連しています。
> require (gbm)
> require (raster)
> require (rgdal)
> load("gbmmodel_p20.Rda")
> print(gbmmodel)
gbm(formula = as.formula(Nclustal_1 ~ tcd_coast_disa_f + tce_raddq_f +
tce_radwq_f + tct_temp_minwin_f + tct_tempdq_f + tcw_clim_etaaann_f +
tcw_precipseas_f + tcw_precipwq_f + tcw_rain1mm_f + tdd_strmdstge6_i +
tlf_logre10_f + tlf_rough0500_f + trs_land_pfc_2008 + trs88_sspr_g_50p +
trs88_ssum_b_50p + trs88_ssum_d_50p + tsp_bd200_f + tsp_cly200a_f +
tsp_ph200_f + tsp_tn060a_f), distribution = "multinomial",
data = gbmdata, n.trees = 2500, interaction.depth = 2, n.minobsinnode = 3,
shrinkage = 0.003, bag.fraction = 0.75, train.fraction = 1,
cv.folds = 8, keep.data = TRUE, verbose = TRUE, class.stratify.cv = TRUE,
n.cores = 8)
A gradient boosted model with multinomial loss function.2500 iterations were performed.
The best cross-validation iteration was 2500.
There were 20 predictors of which 20 had non-zero influence.
次のように、予測変数をラスター スタックにスタックしました。
> img.files <- list.files("/mnt/scratch/mcilwea/R/TSG/inmodel20_test",
pattern='\\.img$', full.names=TRUE)
> rasStack <- stack(img.files)
> NAvalue(rasStack) <- -9999
> projection(rasStack)
"+proj=longlat +ellps=GRS80 +towgs84=-16.237,3.51,9.939,1.4157e-006,2.1477e-006,1.3429e-006,1.91e-007 +no_defs"
rasStack の名前が上のモデルの名前と同じであることを確認することが重要です。
> names(rasStack)
[1] "tcd_coast_disa_f" "tce_raddq_f" "tce_radwq_f"
[4] "tct_tempdq_f" "tct_temp_minwin_f" "tcw_clim_etaaann_f"
[7] "tcw_precipseas_f" "tcw_precipwq_f" "tcw_rain1mm_f"
[10] "tdd_strmdstge6_i" "tlf_logre10_f" "tlf_rough0500_f"
[13] "trs88_sspr_g_50p" "trs88_ssum_b_50p" "trs88_ssum_d_50p"
[16] "trs_land_pfc_2008" "tsp_bd200_f" "tsp_cly200a_f"
[19] "tsp_ph200_f" "tsp_tn060a_f"
predict.gbm を実行する前に、最適な反復モデルを呼び出します
> best.iter <- gbm.perf(gbmmodel, method = "cv", plot.it = TRUE)
グリッド セルを一連の空間ポイントに変換することで、1221 個のセルを含むテスト エリアの一連のラスター出力画像を作成できます (以下を参照)。
points<-raster(img.files[1])
points.df <- as.data.frame(rasterToPoints(points))
coordinates(points.df) <- ~x+y
plot(points.df)
coords <- coordinates(points.df)
rasterOut <- extract(rasStack, coords)
outTable<- as.data.frame(cbind(coords, rasterOut))
outTable[1:1,1:22]
x y tcd_coast_disa_f tce_raddq_f tce_radwq_f tct_temp_minwin_f tct_tempdq_f
149.1269 -35.6457 1.052329 10.82778 23.63533 -0.9852222 5.928154
tcw_clim_etaaann_f tcw_precipseas_f tcw_precipwq_f tcw_rain1mm_f tdd_strmdstge6_i tlf_logre10_f
600 13.93321 179.9841 80.2064 491 1.945529
tlf_rough0500_f trs_land_pfc_2008 trs88_sspr_g_50p trs88_ssum_b_50p trs88_ssum_d_50p tsp_bd200_f
15.6701 0 0.38 0.09000003 0.55 1.590021
tsp_cly200a_f tsp_ph200_f tsp_tn060a_f
33.33834 5.648166 0.03193555
predict.gbm モデルを実行するには
predtable <- as.data.frame(predict.gbm(gbmmodel, outTable, n.trees=best.iter, type="response"))
predout <- cbind(coords,predtable)
predout[1:1,1:38]
x y e24.2500 e26.2500 e59.2500 g152.2500 g157.2500 g94.2500 m31.2500
149.1269 -35.6457 0.001286283 0.0006473167 0.002043077 0.4973372 8.686316e-05 0.0006710651 0.01067058
m36.2500 m68.2500 MU11.2500 MU45.2500 OTHER.2500 p14.2500 p15.2500 p17.2500
0.004314056 0.007128109 0.0005012718 0.0006254022 0.1727706 0.1411112 0.0009099294 0.0002520156
p19.2500 p20.2500 p22.2500 p220.2500 p23.2500 p24.2500 p27.2500 p338.2500
0.003205936 0.002534798 0.0001474091 0.001214219 0.008455798 0.01701965 0.001879607 0.002238932
p420.2500 p520.2500 p54.2500 p9.2500 u118.2500 u179.2500 u21.2500 u22.2500
0.001456685 0.00108458 0.0003695966 0.02501649 0.0005977814 0.01711885 0.0558054 0.002357498
u23.2500 u27.2500 u28.2500 u78.2500 Unit5.2500
0.00040357 0.001422519 0.0002764237 0.01699094 4.835942e-05
write.csv(predout, "Predout.csv", row.names=TRUE)
次の方法で、predtable から 36 個の新しいラスター イメージのセットに発生確率値を書き込むことができます。
names <- names(predtable)
for (i in 1:length(names)) {
SpatialPointspredTable <- SpatialPointsDataFrame (coords=coords, data=predtable[i])
gridded(SpatialPointspredTable)=TRUE
rasValues <- raster(SpatialPointspredTable)
projection(rasValues) <- "+proj=longlat +ellps=GRS80 +towgs84=-16.237,3.51,9.939,1.4157e-006,2.1477e-006,1.3429e-006,1.91e-007 +no_defs"
plot(rasValues)
writeRaster(rasValues, filename=names[i], format="HFA", overwrite=TRUE)
}
これにより、必要な出力が得られます-ただし、データフレームを予測する必要はありません-ラスターブリックに直接予測することができれば、プロセスははるかに高速で効率的になります。
私が走れば
predict(rasStack,
gbmmodel,
n.trees=best.iter,
filename="multiclass_BRT_20p_test_idrisi",
format="IDRISI",
na.rm=FALSE,
type="response",
overwrite=TRUE,
progress="text",
cores=8)
出力は、予測したい最初の植生コミュニティを表すラスター グリッドです。
|=========================================================| 100%
class : RasterLayer
dimensions : 33, 37, 1221 (nrow, ncol, ncell)
resolution : 0.0002777778, 0.0002777778 (x, y)
extent : 149.1268, 149.1371, -35.65473, -35.64556 (xmin, xmax, ymin, ymax)
coord. ref. : NA
data source : /mnt/scratch/mcilwea/R/TSG/multiclass_BRT_20p_test_idrisi.rdc
names : layer
values : 3.762369e-06, 0.9337785 (min, max)
IDRISI ファイル形式はマルチバンド イメージをサポートしていないため、ミックスに index=1:36 を追加してマルチバンド ラスターブリックを出力として生成することはできません。これを実行しようとすると、format="GTiff" または "HFA" (または rgdal を必要とするその他の形式) を設定すると、次のエラー メッセージが表示
されます。 1、オフセット = オフ) : ラスター IO 中のエラー"
ただし、format="raster" を設定すると、rasterbrick 出力を取得できますが、idrisi 画像 (predict.gbm モデルの最初の出力列) 以外のデータを読み書きすることはできません。
「警告メッセージ: In .rasterFromRasterFile(grdfile, band = band, objecttype, ...) : 値ファイルのサイズがセルの数と一致しません (指定されたデータ型)」
predrast <- predict(object=rasStack,
model=gbmmodel,
n.trees=best.iter,
filename="multi_test",
fun=predict.gbm,
format="raster",
index=1:5,
bandorder="BIL",
ext=extent(rasStack[[1:20]]),
na.rm=FALSE,
type="response",
datatype="FLT4S",
overwrite=TRUE,
progress="text",
cores=8)
|=====================================================================100%
predrast
class : RasterBrick
dimensions : 33, 37, 1221, 5 (nrow, ncol, ncell, nlayers)
resolution : 0.0002777778, 0.0002777778 (x, y)
extent : 149.1268, 149.1371, -35.65473, -35.64556 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +ellps=GRS80 +towgs84=-16.237,3.51,9.939,1.4157e-006,2.1477e-006,1.3429e-006,1.91e-007 +no_defs
data source : C:\Data\FINAL_TSG\test\multi_test.grd
names : layer.1, layer.2, layer.3, layer.4, layer.5
min values : 3.762369e-06, 3.762369e-06, 3.762369e-06, 3.762369e-06, 3.762369e-06
max values : 0.9337785, 0.9337785, 0.9337785, 0.9337785, 0.9337785
上記のラスターブリックを個々のラスター画像のセットに変換しようとすると
writeRaster(predrast, filename="multi_test.img", format="HFA", bylayer=TRUE, suffix="numbers", overwrite=TRUE)
どの画像も意味をなさない。
また、マルチバンド CDF イメージとして書き込もうとすると、rgdal エラーとは異なる一連の警告メッセージが表示されることも少し不可解です。
| 0%
Error in ncdf::put.var.ncdf(nc, x@title, v, start = c(1, start, lstart), :
put.var.ncdf: error: you asked to write 11988 values, but the passed data array only has 11840 entries!
|======== | 25%
Error in ncdf::put.var.ncdf(nc, x@title, v, start = c(1, start, lstart), :
put.var.ncdf: error: you asked to write 11988 values, but the passed data array only has 11840 entries!
|================== | 50%
Error in ncdf::put.var.ncdf(nc, x@title, v, start = c(1, start, lstart), :
put.var.ncdf: error: you asked to write 11988 values, but the passed data array only has 11840 entries!
|=============================================== | 75%
Error in ncdf::put.var.ncdf(nc, x@title, v, start = c(1, start, lstart), :
put.var.ncdf: error: you asked to write 7992 values, but the passed data array only has 7955 entries!
|=============================================================| 100%
ここで、何が起こっているのかわかりませんか??
上記の問題に遭遇することなく、rasterbrick に直接予測できるようにするために、gbm パッケージの作成者と協力できる方法を知っている人がいれば素晴らしいことです。
完全なラスター グリッドで使用したコードを誰かが知りたい場合は、下にコメントを残してください。喜んで提供します。
乾杯アレン
sessionInfo()
R version 3.1.2 (2014-10-31)
Platform: x86_64-w64-mingw32/x64 (64-bit)
locale:
[1] LC_COLLATE=English_Australia.1252 LC_CTYPE=English_Australia.1252 LC_MONETARY=English_Australia.1252
[4] LC_NUMERIC=C LC_TIME=English_Australia.1252
attached base packages:
[1] parallel splines stats graphics grDevices utils datasets methods base
other attached packages:
[1] ncdf_1.6.8 rgdal_0.9-1 gbm_2.1 lattice_0.20-30 survival_2.37-7 raster_2.3-24 sp_1.0-17
loaded via a namespace (and not attached):
[1] grid_3.1.2 tools_3.1.2
# Traceback error for
Error in rgdal::putRasterData(x@file@transient, v, band = 1, offset = off) :
Failure during raster IO
> traceback()
7: .Call("RGDAL_PutRasterData", raster, rasterData, as.integer(offset),
PACKAGE = "rgdal")
6: rgdal::putRasterData(x@file@transient, v, band = 1, offset = off)
5: writeValues(predrast, predv, tr$row[i])
4: writeValues(predrast, predv, tr$row[i])
3: .local(object, ...)
2: predict(object = rasStack, model = gbmmodel, n.trees = best.iter,
filename = "multi_img", format = "HFA", na.rm = FALSE, type = "response",
datatype = "FLT4S", overwrite = TRUE, progress = "text")
1: predict(object = rasStack, model = gbmmodel, n.trees = best.iter,
filename = "multi_img", format = "HFA", na.rm = FALSE, type = "response",
datatype = "FLT4S", overwrite = TRUE, progress = "text")