Last chance! 50% off unlimited learning
Sale ends in
Perform spatial prediction, producing a raster of predictions and conditional standard deviations.
krigeLgm(formula, data, grid, covariates = NULL,
param,
expPred = FALSE, nuggetInPrediction = TRUE,
mc.cores=getOption("mc.cores", 1L))
Either a model formula, or a data frame of linear covariates.
A SpatialPointsDataFrame containing the data to be interpolated
Either a raster
, or a single integer giving the
number of cells in the X direction which predictions will be made on. If the later
the predictions will be a raster of square cells covering the bounding box of data
.
The spatial covariates used in prediction, either a raster
stack or list of rasters.
A vector of named model parameters, as produced by likfitLgm
Should the predictions be exponentiated, defaults to FALSE
.
If TRUE
, predict new observations by adding the
nugget effect. The prediction variances will be adjusted accordingly, and the predictions
on the natural scale for logged or Box Cox transformed data will be affected.
Otherwise predict fitted values.
passed to mcmapply
if greater than 1.
A raster stack is returned with the following layers:
Estimated means from the fixed effects portion of the model
Predicted random effect
Conditional variance of predicted random effect (on the transformed scale if applicable)
Prediction of the response, sum of fixed and random effects. If exp.pred is TRUE, gives predictions on the exponentiated scale, and half of krige.var is added prior to exponentiating
If exp.pred=TRUE, the prediction of the logged process.
If a box cox transformation was used, the prediction of the process on the transformed scale.
A raster stack as above, though the random effect prediction is resampled to the same locations as the fixed effects.
the predictions and conditional variance of the random effects, on the same
raster as newdata
Given the model parameters and observed data, conditional means and variances of the spatial random field are computed.
# NOT RUN {
data('swissRain')
swissRain$lograin = log(swissRain$rain)
swissRain[[names(swissAltitude)]] = extract(swissAltitude, swissRain)
if(interactive() | Sys.info()['user'] =='patrick') {
swissFit = likfitLgm(data=swissRain,
formula=lograin~ CHE_alt,
param=c(range=46500, nugget=0.05,shape=1,
anisoAngleDegrees=35, anisoRatio=12),
paramToEstimate = c("range","nugget",
"anisoAngleDegrees", "anisoRatio")
)
myTrend = swissFit$model$formula
myParams = swissFit$param
dput(myParams)
# will give the following
} else {
myParams=structure(c(0.0951677553339405, 0.773131120360615, 49376.366554348,
1, 11.6731542330799, 0.649923845971017, 1, 2.26103865436198,
0.000146945159199578, 37.2378933790499), .Names = c("nugget",
"variance", "range", "shape", "anisoRatio", "anisoAngleRadians",
"boxcox", "(Intercept)", "CHE_alt", "anisoAngleDegrees"))
myTrend =lograin~ CHE_alt
}
# make sure krige can cope with missing values!
swissAltitude[1:50,1:50] = NA
swissKrige = krigeLgm(data=swissRain,
formula = myTrend,
covariates = swissAltitude,
param=myParams,
grid = 40, expPred=TRUE)
plot(swissKrige[["predict"]], main="predicted rain")
plot(swissBorder, add=TRUE)
# now with box cox and provide a raster for prediction, no covariates
if(interactive() | Sys.info()['user'] =='patrick') {
swissFit2 = likfitLgm(data=swissRain,
formula=rain~1,
param=c(range=52000, nugget=0.1,shape=1, boxcox=0.5,
anisoAngleDegrees=35, anisoRatio=8),
paramToEstimate = c("range","nugget",
"anisoAngleDegrees", "anisoRatio"),
parscale = c(range=5000,nugget=0.01,
anisoRatio=1,anisoAngleDegrees=5)
)
myTrend2 = swissFit2$trend
myParams2 = swissFit2$param
dput(myParams2)
} else {
myParams2=structure(c(0.865530531647866, 8.76993204385615, 54143.5826959284,
1, 7.36559089705556, 0.647158492167979, 0.5, 5.16254700135706,
37.0794502772753), .Names = c("nugget", "variance", "range",
"shape", "anisoRatio", "anisoAngleRadians", "boxcox", "(Intercept)",
"anisoAngleDegrees"))
myTrend2=rain~1
}
swissRaster = raster(extent(swissBorder), nrows=25, ncols=40)
swissKrige2 = krigeLgm(data=swissRain, formula = myTrend2,
param=myParams2,
grid = swissRaster)
plot(swissKrige2[["predict"]], main="predicted rain with box-cox")
plot(swissBorder, add=TRUE)
# }
Run the code above in your browser using DataLab