RFoptions(seed=0) ## *ANY* simulation will have the random seed 0; set
## RFoptions(seed=NA) to make them all random again
################################################################
## ##
## a geostatistical analysis that demonstrates ##
## features of the package `RandomFields' ##
## ##
################################################################
data(soil)
str(soil)
soil <- RFspatialPointsDataFrame(
coords = soil[, c("x.coord", "y.coord")],
data = soil[, c("moisture", "NO3.N", "Total.N",
"NH4.N", "DOC", "N20N")],
RFparams=list(vdim=6, n=1)
)
data <- soil["moisture"]
if (!interactive()) {
data(soil)
idx <- 1:10
data <- RFspatialPointsDataFrame(
coords = soil[idx, c("x.coord", "y.coord")],
data = data.frame(moisture=soil[["moisture"]][idx]),
RFparams=list(vdim=1, n=1)
)["moisture"]
}
## plot the data first
colour <- rainbow(100)
plot(soil["moisture"], col=colour)
## fit by eye
RFgui.model <- RFgui(soil["moisture"])
if (!interactive()) RFgui.model <- RMexp()
## fit by ML
model <- ~1 + RMplus(RMwhittle(scale=NA, var=NA, nu=NA), RMnugget(var=NA))
fit <- RFfit(model, data=data)
plot(fit, fit.method=c("ml", "plain", "sqrt.nr", "sd.inv"),
model = RFgui.model, col=1:8)
## Kriging ...
x <- seq(min(data@coords[, 1]), max(data@coords[, 1]),
l=if (interactive()) 121 else 4)
k <- RFinterpolate(fit["ml"], x=x, y=x, grid=TRUE, data=data)
plot(x=k, col=colour)
plot(x=k, y=data, col=colour)
## what is the probability that at no point of the
## grid given by x and y the moisture is greater than 24 percent?
cs <- RFsimulate(model=fit["ml"], x=x, y=x, data=data,
n=if (interactive()) 50 else 3)
plot(cs, col=colour)
plot(cs, y=data, col=colour)
Print(mean(apply(as.matrix(cs@data) <= 24, 2, all))) ## about 40 percent
##...
FinalizeExample()
Run the code above in your browser using DataLab