mettebs

Predicting with R - anyone?

Discussion created by mettebs on Apr 8, 2014
Latest reply on Apr 9, 2014 by jevans02
I hope someone outthere can give me some ideas as to what might be wrong. I'm using R for predicting a resulting raster from several rasters used as predictors. One of these rasters is causing me problems. If I use it in the prediction, R will crash. If I leave it out and only run with the 16 other predicting rasters everything is fine. I've tried reclassifying it and cutting it up into smaller chunks. When cut into three parts, two of them will run and one will crash. All predictor rasters are ESRI Grids, the problematic one is a categorical raster and it looks fine in ArcGIS.

All rasters have the same projection, cell size and extent.

Please if anyone have an idea or suggestion I could try, I'd be very grateful. It's a little off topic, but I might be lucky that someone have tried this before. My code is below, but it runs fine without that one particular raster, so...

TIA,
Mette

library(Cubist)
library(caret)
library(class)

setwd("D:/dtemp/Texture/ModellingMarch2014")
data<-read.table(file="clay_silt_fsand_Adata.txt", header=T, sep=";")
str(data)####data structure

###########give factor data####
data$georeg<-as.factor(data$georeg)
data$landscape<-as.factor(data$landscape)
data$wetlands<-as.factor(data$wetlands)
data$landuse<-as.factor(data$landuse)
data$geology<-as.factor(data$geology)
data$fgjord<-as.factor(data$fgjord)


names(data)

nc<-ncol(data)
data_predictors<-data[,c(6:nc)]


####make cubist model####
ctrl<- trainControl(method = "LGOCV",repeats=10, number=10,
                    summaryFunction = defaultSummary, selectionFunction = "best")
cubist_model<-train(data[,c(6:nc)],data[,2],method="cubist",trControl=ctrl)


###########package for read rasters###
library(raster)
library(rgdal)

#########read raster file###
aspect<-raster("D:/dtemp/Texture/Tiles/Tile1/aspect")
centallav0<-raster("D:/dtemp/Texture/Tiles/Tile1/centallav0")
dirinsola<-raster("D:/dtemp/Texture/Tiles/Tile1/dirinsola")
elevation<-raster("D:/dtemp/Texture/Tiles/Tile1/elevation")
fgjord<-raster("D:/dtemp/Texture/Tiles/Tile1/fgjord")
flowaccu<-raster("D:/dtemp/Texture/Tiles/Tile1/flowaccu")
geology<-raster("D:/dtemp/Texture/Tiles/Tile1/geology")
georeg<-raster("D:/dtemp/Texture/Tiles/Tile1/georeg")
landscape<-raster("D:/dtemp/Texture/Tiles/Tile1/landscape")
landuse<-raster("D:/dtemp/Texture/Tiles/Tile1/landuse")
midslppos<-raster("D:/dtemp/Texture/Tiles/Tile1/midslppos")
mrvbf<-raster("D:/dtemp/Texture/Tiles/Tile1/mrvbf")
sagawi<-raster("D:/dtemp/Texture/Tiles/Tile1/sagawi")
slpdeg<-raster("D:/dtemp/Texture/Tiles/Tile1/slpdeg")
valldepth<-raster("D:/dtemp/Texture/Tiles/Tile1/valldepth")
vdistchn<-raster("D:/dtemp/Texture/Tiles/Tile1/vdistchn")
wetlands<-raster("D:/dtemp/Texture/Tiles/Tile1/wetlands")

predictors<-c(aspect,centallav0,dirinsola,elevation,fgjord,flowaccu,geology,
              georeg,landscape,landuse,midslppos,
              mrvbf,sagawi,slpdeg,valldepth,vdistchn,wetlands)

var.stack <- stack(predictors)###all input should have same extent and cell size

#Cubist_result ###print the result
capture.output(summary(cubist_model), file = "summary(cubist_model_clayA_tile1).txt")###write results out to text file
system.time(predict(var.stack, cubist_model,
                    filename="D:/dtemp/Texture/ModellingMarch2014/clayA_tile1.img",
                    fun=predict, ext=NULL, const=NULL, index=1, na.rm=TRUE, progress="text",
                    format="HFA", overwrite=TRUE))

Outcomes