data('swissRain')
swissAltitude = unwrap(swissAltitude)
swissRain = unwrap(swissRain)
swissRain$lograin = log(swissRain$rain)
swissRain[[names(swissAltitude)]] = extract(swissAltitude, swissRain, ID=FALSE)
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
swissBorder = unwrap(swissBorder)
swissKrige = krigeLgm(
data=swissRain,
formula = myTrend,
covariates = swissAltitude,
param=myParams,
grid = squareRaster(swissBorder, 40), expPred=TRUE)
plot(swissKrige[["predict"]], main="predicted rain")
plot(swissBorder, add=TRUE)
Run the code above in your browser using DataLab