- meteo.Rnw, the single Sweave source file; derived from this are:
- meteo.R the file with the extracted R sections, can be pasted straight into an R session for debugging purposes
- meteo.tex the latex file, to be processed by latex or pdflatex
- meteo.pdf the pdf file resulting from pdflatex'ing the latex file

- restart the thin client (put it off by keeping the on-button pressed, and then on again), and log in
- Copy the grass data base from /home/epebe_01/grassdata/Schweiz/AnLe to your area, and start grass using this as data base
- within grass, on the shell prompt, start R
- install R packages gstat and spgrass6, and load them
- run the script below
library(spgrass6) bm = readVECT6("bm1198_97p") rast = readRAST6("bm1198_97pidw2") summary(bm) p4s = proj4string(bm) bm = as.data.frame(bm) coordinates(bm)=~coords.x1+coords.x2 proj4string(bm)=p4s summary(bm)

- explain what fundamentally changed when bm was converted to data.frame, and back to a spatial structure again.
- Find out what the standard deviation of a typical interpolation
error is, by using cross validation on inverse distance interpolation:
library(gstat) idw.cv = krige.cv(coords.x3~1, bm, nfold=100) idw.cv$obs = bm$coords.x3 plot(var1.pred~obs,idw.cv) sd(idw.cv$residual) s = sd(idw.cv$residual)

- Generate white noise on a number of interpolated surfaces
inv = idw(coords.x3~1, bm, rast, nmax = 30) inv$sim1 = inv$var1.pred + s * rnorm(length(inv$var1.pred)) inv$sim2 = inv$var1.pred + s * rnorm(length(inv$var1.pred)) inv$sim3 = inv$var1.pred + s * rnorm(length(inv$var1.pred)) inv$sim4 = inv$var1.pred + s * rnorm(length(inv$var1.pred)) inv$sim5 = inv$var1.pred + s * rnorm(length(inv$var1.pred)) spplot(inv[3:7], col.regions=bpy.colors())

- export these simulated surfaces to GRASS (R function writeRAST6), compute a slope map for each, and compare these slope maps visually in GRASS.
- Use, instead, variograms and kriging for interpolation:
plot(variogram(coords.x3~1,bm)) v0=variogram(coords.x3~1,bm) fit.variogram(v, vgm(1, "Exp", 5000)) v.fit = fit.variogram(v, vgm(1, "Exp", 5000)) plot(v,v.fit) v = variogram(coords.x3~1,bm[-393,],width=50,cutoff=500) v.fit = fit.variogram(v, vgm(1, "Mat", 2000, kappa=1.5)) plot(v0, v.fit) plot(v, v.fit) kr = krige(coords.x3~1, bm[-393,], rast, v.fit, nmax=30) spplot(kr[1], col.regions=bpy.colors())

- carry out the following simulations
sim = krige(coords.x3~1, bm[-393,], rast, v.fit, nmax=30, nsim = 10) spplot(sim, col.regions=bpy.colors()) parseGRASS("r.slope.aspect") pars = list(elevation="elevation.dem", slope="slope", aspect="aspect") slope = sim aspect = sim for (i in 1:10) { print(i) writeRAST6(sim, "elevation.dem", i, overwrite = TRUE) execGRASS("r.slope.aspect", flags=c("overwrite"), parameters=pars) slope[[i]] = readRAST6("slope")[[1]] aspect[[i]] = readRAST6("aspect")[[1]] } spplot(slope, col.regions=bpy.colors()) spplot(aspect, col.regions=bpy.colors())

- describe what the parseGRASS and execGRASS commands do.
- look at differences in slope, and in aspect, between the different simulations
- compare the maps showing the area where the slope is larger than 30 based on inverse distance interpolation and kriging interpolation
- compare the areas with the area for which we have 95 percent confidence that the slope is smaller than 30, based on the simple (white noise) error model and the (correlated) kriging error model, and compare these areas with the interpolated ones (idw and kriging).

- CRAN - task views - spatial: CRAN spatial task view
- Generic R sources from the California soil research lab
- Grass + R on the spearfish data
- Tom Hengl's Spatial-analyst.net and geostat 2010 summer school
- Jon Fox' paper on R's future in the R journal