R links



Sweave is a utility to create a report from a single source, that is processed once to perform the full analysis and render the output. This way, confusion about which analysis happened where or when can be excluded. The following files are of interest:

Exercise R, day 1:

repeat the meteo analysis, but use average temperature data for May, 14:00 (and 14:40?). Show a plot of the monthly mean values against altitude with regression line; give a map with predicted average May temperatures, based on linear regression and interpolated altitude data. If there is time left, give in addition a map with the regression prediction standard errors. Send results with a brief report to me by email.

Exercises R/GRASS, day 2

Today we will look into combined use of R and GRASS. The topic is the uncertainty in DEM-derived slopes, when the DEM is subject to statistical error (uncertainty), e.g. resulting from interpolation.
  1. restart the thin client (put it off by keeping the on-button pressed, and then on again), and log in
  2. Copy the grass data base from /home/epebe_01/grassdata/Schweiz/AnLe to your area, and start grass using this as data base
  3. within grass, on the shell prompt, start R
  4. install R packages gstat and spgrass6, and load them
  5. run the script below
    bm = readVECT6("bm1198_97p")
    rast = readRAST6("bm1198_97pidw2")
    p4s = proj4string(bm)
    bm = as.data.frame(bm)
  6. explain what fundamentally changed when bm was converted to data.frame, and back to a spatial structure again.
  7. Find out what the standard deviation of a typical interpolation error is, by using cross validation on inverse distance interpolation:
    idw.cv = krige.cv(coords.x3~1, bm, nfold=100)
    idw.cv$obs = bm$coords.x3
    s = sd(idw.cv$residual)
  8. 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())
  9. export these simulated surfaces to GRASS (R function writeRAST6), compute a slope map for each, and compare these slope maps visually in GRASS.
  10. Use, instead, variograms and kriging for interpolation:
    fit.variogram(v, vgm(1, "Exp", 5000))
    v.fit = fit.variogram(v, vgm(1, "Exp", 5000))
    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())
  11. 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())
    pars = list(elevation="elevation.dem", slope="slope", aspect="aspect")
    slope = sim
    aspect = sim
    for (i in 1:10) {
        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())
  12. describe what the parseGRASS and execGRASS commands do.
  13. look at differences in slope, and in aspect, between the different simulations
  14. compare the maps showing the area where the slope is larger than 30 based on inverse distance interpolation and kriging interpolation
  15. 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).

Lecture slides

Slides on R, found here; associated files: the Sweave file and the resulting R script

R sources

Open source GIS sources

The partly outdated paper on the state of open source gis (2006)