Edzer Pebesma April 2, 2014
load_ext rpy2.ipython
%%R
library(sp)
data(meuse)
class(meuse)
[1] "data.frame"
%%R
names(meuse)
[1] "x" "y" "cadmium" "copper" "lead" "zinc" "elev" [8] "dist" "om" "ffreq" "soil" "lime" "landuse" "dist.m"
%%R
coordinates(meuse) = ~x+y
class(meuse)
[1] "SpatialPointsDataFrame" attr(,"package") [1] "sp"
%%R
bubble(meuse, "zinc", col=c("#00ff0088", "#00ff0088"), main = "zinc concentrations (ppm)")
%%R
data(meuse.grid)
coordinates(meuse.grid) = ~x+y
%%R
gridded(meuse.grid) = TRUE
image(meuse.grid["dist"])
title("distance to river (red = 0)")
library(gstat)
zinc.idw = krige(zinc~1, meuse, meuse.grid)
[inverse distance weighted interpolation]
%%R
spplot(zinc.idw["var1.pred"], main = "zinc inverse distance weighted interpolations")
%%R
plot(log(zinc)~sqrt(dist), meuse)
abline(lm(log(zinc)~sqrt(dist), meuse))
%%R
lzn.vgm = variogram(log(zinc)~1, meuse)
lzn.fit = fit.variogram(lzn.vgm, model = vgm(1, "Sph", 900, 1))
plot(lzn.vgm, lzn.fit)
%%R
lznr.vgm = variogram(log(zinc)~sqrt(dist), meuse)
lznr.fit = fit.variogram(lznr.vgm, model = vgm(1, "Exp", 300, 1))
plot(lznr.vgm, lznr.fit)
%%R
lzn.kriged = krige(log(zinc)~1, meuse, meuse.grid, model = lzn.fit)
spplot(lzn.kriged["var1.pred"])
[using ordinary kriging]
%%R
lzn.condsim = krige(log(zinc)~1, meuse, meuse.grid, model = lzn.fit, nmax = 30, nsim = 4)
spplot(lzn.condsim, main = "three conditional simulations")
drawing 4 GLS realisations of beta... [using conditional Gaussian simulation]
%%R
lzn.condsim2 = krige(log(zinc)~sqrt(dist), meuse, meuse.grid, model = lznr.fit, nmax = 30, nsim = 4)
spplot(lzn.condsim2, main = "three UK conditional simulations")
drawing 4 GLS realisations of beta... [using conditional Gaussian simulation]
%%R
lzn.dir = variogram(log(zinc)~1, meuse, alpha = c(0, 45, 90, 135))
lzndir.fit = vgm(.59, "Sph", 1200, .05, anis = c(45, .4))
plot(lzn.dir, lzndir.fit, as.table = TRUE)
%%R
lznr.dir = variogram(log(zinc)~sqrt(dist), meuse, alpha = c(0, 45, 90, 135))
plot(lznr.dir, lznr.fit, as.table = TRUE)
%%R
vgm.map = variogram(log(zinc)~sqrt(dist), meuse, cutoff = 1500, width = 100,map = TRUE)
plot(vgm.map, threshold = 5)
%%R
g = gstat(NULL, "log(zn)", log(zinc)~sqrt(dist), meuse)
g = gstat(g, "log(cd)", log(cadmium)~sqrt(dist), meuse)
g = gstat(g, "log(pb)", log(lead)~sqrt(dist), meuse)
g = gstat(g, "log(cu)", log(copper)~sqrt(dist), meuse)
v = variogram(g)
g = gstat(g, model = vgm(1, "Exp", 300, 1), fill.all = TRUE)
g.fit = fit.lmc(v,g)
%%R
plot(v, g.fit)
%%R
vgm.map = variogram(g, cutoff = 1500, width = 100, map = TRUE)
plot(vgm.map, threshold = 5, col.regions = bpy.colors(), xlab = "", ylab = "")