krige: Spatial prediction, or Kriging

Description Usage Arguments Details Value See Also Examples

Description

Perform spatial prediction, producing a raster of predictions and conditional standard deviations.

Usage

1
2
3
4
krigeLgm(formula, data, grid,  covariates = NULL,
	param,  
    expPred = FALSE, nuggetInPrediction = TRUE,
    mc.cores=getOption("mc.cores", 1L))

Arguments

formula

Either a model formula, or a data frame of linear covariates.

data

A SpatialPointsDataFrame containing the data to be interpolated

grid

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.

covariates

The spatial covariates used in prediction, either a raster stack or list of rasters.

param

A vector of named model parameters, as produced by likfitLgm

expPred

Should the predictions be exponentiated, defaults to FALSE.

nuggetInPrediction

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.

mc.cores

passed to mcmapply if greater than 1.

Details

Given the model parameters and observed data, conditional means and variances of the spatial random field are computed.

Value

A raster stack is returned with the following layers:

fixed

Estimated means from the fixed effects portion of the model

random

Predicted random effect

krige.var

Conditional variance of predicted random effect (on the transformed scale if applicable)

predict

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

predict.log

If exp.pred=TRUE, the prediction of the logged process.

predict.boxcox

If a box cox transformation was used, the prediction of the process on the transformed scale.

If the prediction locations are different for fixed and random effects (typically coarser for the random effects), a list with two raster stacks is returned.

prediction

A raster stack as above, though the random effect prediction is resampled to the same locations as the fixed effects.

random

the predictions and conditional variance of the random effects, on the same raster as newdata

See Also

lgm

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
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.0951770829449953, 0.77308786208928, 49379.3845752436, 
1, 11.673076577513, 0.649925237653982, 1, 2.26103881494066, 0.000146945995279231, 
37.2379731166102), .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)

geostatsp documentation built on July 14, 2018, 9:01 a.m.