R links

• R web site
• The R manuals (start with the introduction)

Scripts

• importing: import.R
• analysis: script.R

Sweave

• 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

Exercise R, day 1:

Exercises R/GRASS, day 2

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

   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)
   

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:

   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)
   

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:

    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())
    

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())
    
    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())
    

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

R sources

• 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

Open source GIS sources