\documentclass{article} \usepackage{graphicx} \usepackage[colorlinks=true,urlcolor=blue]{hyperref} \usepackage{color} \usepackage{Sweave} \begin{document} \section{Data import} First, we'll import the stations and meteo data from .csv files, and give temperature a decent name: \begin{Schunk} \begin{Sinput} > stations = read.csv("stations.csv") > lapply(stations, class) \end{Sinput} \begin{Soutput} $X01.stationkennung [1] "factor" $X02.station.name [1] "factor" $X03.messnetz [1] "factor" $X04.koord..x [1] "integer" $X05.koord..y [1] "integer" $X06.hoehe.ueber.meer..m. [1] "integer" $X07.min..datum [1] "factor" $X08.max..datum [1] "logical" $X09.kanton [1] "factor" $X10.klimat.region.nr [1] "integer" $X11.landkarte [1] "logical" \end{Soutput} \begin{Sinput} > slf_ams = read.csv("slf_ams.csv") > lapply(slf_ams, class) \end{Sinput} \begin{Soutput} $X01.stationsabk.rzung [1] "factor" $X02.standortnummer [1] "integer" $X03.datum.und.zeit [1] "factor" $X04.windgeschwindigkeit..m.s. [1] "numeric" $X05.windrichtung..grad. [1] "integer" $X06.gr.sste.sekundenb.e..m.s. [1] "numeric" $X07.lufttemperatur...c. [1] "numeric" $X08.rel..luftfeuchtigkeit.... [1] "integer" $X09.reflekt..strahlung..w.m2. [1] "integer" $X10.schneeh.he...cm. [1] "integer" $X11.schneeoberfl.chetemp....c. [1] "numeric" \end{Soutput} \begin{Sinput} > slf_ams$temp = as.numeric(slf_ams$X07.lufttemperatur...c.) \end{Sinput} \end{Schunk} \section{Sanity checks} Next, we'll check if the station names \begin{Schunk} \begin{Sinput} > levels(stations$X01.stationkennung) \end{Sinput} \begin{Soutput} [1] "DAV-3" "DAV-4" "DAV-5" "KLO-3" "PAR-2" "SLF-2" "WFJ-1" "WFJ-2" \end{Soutput} \end{Schunk} agree with the slf\_ams stations, defined by two fields: \begin{Schunk} \begin{Sinput} > unique(paste(slf_ams$X01.stationsabk.rzung, slf_ams$X02.standortnummer, + sep = "-")) \end{Sinput} \begin{Soutput} [1] "DAV-3" "DAV-4" "DAV-5" "KLO-3" "PAR-2" "SLF-2" "WFJ-1" "WFJ-2" \end{Soutput} \end{Schunk} Then, we create a common field, called station: \begin{Schunk} \begin{Sinput} > stations$station = stations$X01.stationkennung > slf_ams$station = paste(slf_ams$X01.stationsabk.rzung, slf_ams$X02.standortnummer, + sep = "-") \end{Sinput} \end{Schunk} \section{Time, monthly means} Now we'll work on the time field, which is now a factor that looks like this: \begin{Schunk} \begin{Sinput} > slf_ams$X03.datum.und.zeit[1] \end{Sinput} \begin{Soutput} [1] 2002.10.01 00:00 104973 Levels: 2002.10.01 00:00 2002.10.01 00:30 ... 2006.10.01 00:00 \end{Soutput} \end{Schunk} and we need to manually set the format of the fields (see ?strptime): \begin{Schunk} \begin{Sinput} > slf_ams$time = as.POSIXct(strptime(as.character(slf_ams$X03.datum.und.zeit), + "%Y.%m.%d %H:%M")) \end{Sinput} \end{Schunk} Make the selection for february and august data: \begin{Schunk} \begin{Sinput} > feb = format.POSIXct(slf_ams$time, "%m") == "02" > summary(feb) \end{Sinput} \begin{Soutput} Mode FALSE TRUE NA's logical 448706 37463 0 \end{Soutput} \begin{Sinput} > feb[1:10] \end{Sinput} \begin{Soutput} [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE \end{Soutput} \begin{Sinput} > aug = format.POSIXct(slf_ams$time, "%m") == "08" \end{Sinput} \end{Schunk} Make the selection for the 14:00 data; if we onlys select 14:00 we loose 2 out of 8 stations! \begin{Schunk} \begin{Sinput} > at1400 = format.POSIXct(slf_ams$time, "%H") == "14" & format.POSIXct(slf_ams$time, + "%M") == "00" > at1340 = format.POSIXct(slf_ams$time, "%H") == "13" & format.POSIXct(slf_ams$time, + "%M") == "40" > afternoon = at1400 | at1340 \end{Sinput} \end{Schunk} Select the tables (data.frames) on [rows, columns]: \begin{Schunk} \begin{Sinput} > temp_feb = na.omit(slf_ams[feb & afternoon, c("temp", "station")]) > temp_aug = na.omit(slf_ams[aug & afternoon, c("temp", "station")]) \end{Sinput} \end{Schunk} ... and compute station means: \begin{Schunk} \begin{Sinput} > T.feb = lapply(split(temp_feb$temp, as.factor(temp_feb$station)), + mean) > T.aug = lapply(split(temp_aug$temp, as.factor(temp_aug$station)), + mean) \end{Sinput} \end{Schunk} Here, {\tt split} splits the first argument according to groups in the second, and {\tt lapply} applies something (function in second argument) to each element of a list (1st arg.) The function unlist creates a simple vector: \begin{Schunk} \begin{Sinput} > out = data.frame(T.aug = unlist(T.aug), T.feb = unlist(T.feb)) > out \end{Sinput} \begin{Soutput} T.aug T.feb DAV-3 8.677419 -6.755357 DAV-4 9.556098 -6.666964 DAV-5 10.742157 -5.366667 KLO-3 9.579032 -5.471429 PAR-2 9.453659 -6.125893 SLF-2 15.146341 -1.681651 WFJ-1 7.157258 -9.573451 WFJ-2 7.616129 -7.438938 \end{Soutput} \end{Schunk} \section{Adding coordinates; matching tables} Now we will match table {\tt stations} with {\tt out} in order to get the coordinates and altitude of the temperature means: \begin{Schunk} \begin{Sinput} > m = match(row.names(out), as.character(stations$X01.stationkennung)) > m \end{Sinput} \begin{Soutput} [1] 5 3 4 6 7 8 2 1 \end{Soutput} \end{Schunk} Check that we match right; \verb|x[m]| puts x in order m; {\tt all.equal} checks equality of two things \begin{Schunk} \begin{Sinput} > all.equal(as.character(stations$X01.stationkennung[m]), row.names(out)) \end{Sinput} \begin{Soutput} [1] TRUE \end{Soutput} \begin{Sinput} > out$x = stations$X04.koord..x[m] > out$y = stations$X05.koord..y[m] > out$alt = stations$X06.hoehe.ueber.meer..m.[m] > out \end{Sinput} \begin{Soutput} T.aug T.feb x y alt DAV-3 8.677419 -6.755357 778300 184580 2450 DAV-4 9.556098 -6.666964 779125 184125 2330 DAV-5 10.742157 -5.366667 779530 184975 2300 KLO-3 9.579032 -5.471429 790100 190800 2310 PAR-2 9.453659 -6.125893 780430 191680 2290 SLF-2 15.146341 -1.681651 783800 187400 1560 WFJ-1 7.157258 -9.573451 780620 189650 2693 WFJ-2 7.616129 -7.438938 780850 189260 2540 \end{Soutput} \end{Schunk} \section{Regression, regression diagnostics} We can plot some regression diagnostics: \begin{Schunk} \begin{Sinput} > par(mfrow = c(2, 2)) > f1 = T.feb ~ alt > plot(f1, out) > abline(lm(f1, out)) > title(round(coefficients(lm(f1, out)), 6)) > f2 = T.aug ~ alt > plot(f2, out) > abline(lm(f2, out)) > title(round(coefficients(lm(f2, out)), 6)) > plot(T.feb ~ x, out) > plot(T.feb ~ y, out) \end{Sinput} \end{Schunk} \includegraphics{meteo-014} \section{Interpolating altitude, using idw} Next, we import the altitude points. Package {\tt rgdal} deals with raster and vector import; it loads package sp. \begin{Schunk} \begin{Sinput} > library(rgdal) > bm1197p = readOGR("../GISAufbaukurs_WS0708/data/1197-1198/DHM25/basis/point", + "bm1197p") \end{Sinput} \begin{Soutput} OGR data source with driver: ESRI Shapefile Source: "../GISAufbaukurs_WS0708/data/1197-1198/DHM25/basis/point", layer: "bm1197p" with 800 features and 4 fields Feature type: wkbPoint with 3 dimensions \end{Soutput} \begin{Sinput} > bm1198p = readOGR("../GISAufbaukurs_WS0708/data/1197-1198/DHM25/basis/point", + "bm1198p") \end{Sinput} \begin{Soutput} OGR data source with driver: ESRI Shapefile Source: "../GISAufbaukurs_WS0708/data/1197-1198/DHM25/basis/point", layer: "bm1198p" with 784 features and 4 fields Feature type: wkbPoint with 3 dimensions \end{Soutput} \begin{Sinput} > summary(bm1197p) \end{Sinput} \begin{Soutput} Object of class SpatialPointsDataFrame Coordinates: min max coords.x1 777510.9 794962.5 coords.x2 182009.4 193996.9 coords.x3 1179.0 3085.0 Is projected: NA proj4string : [NA] Number of points: 800 Data attributes: OBJECTID OBJECTVAL OBJECTORIG YEAROFCHAN Min. :8641614 Hoehenkote :799 LK25:800 Min. :1985 1st Qu.:8641814 Seebodenkote: 1 1st Qu.:1985 Median :8642014 Median :1985 Mean :8642014 Mean :1985 3rd Qu.:8642213 3rd Qu.:1985 Max. :8642413 Max. :1985 \end{Soutput} \begin{Sinput} > summary(bm1198p) \end{Sinput} \begin{Soutput} Object of class SpatialPointsDataFrame Coordinates: min max coords.x1 795017.1 812479.6 coords.x2 182009.4 193915.6 coords.x3 1298.0 3411.0 Is projected: NA proj4string : [NA] Number of points: 784 Data attributes: OBJECTID OBJECTVAL OBJECTORIG YEAROFCHAN Min. :8642414 Hoehenkote:784 LK25:784 Min. :1985 1st Qu.:8642610 1st Qu.:1985 Median :8642806 Median :1985 Mean :8642806 Mean :1985 3rd Qu.:8643001 3rd Qu.:1985 Max. :8643197 Max. :1985 \end{Soutput} \end{Schunk} The first two commands set the sets bpy.colors() colour ramp; then the two point data files are merged (rbind: row bind), and converted to a spatial data set: \begin{Schunk} \begin{Sinput} > library(lattice) > trellis.par.set(sp.theme()) > bmp = rbind(bm1197p, bm1198p) > bmp = as.data.frame(bmp) > names(bmp) = c(names(bmp)[1:4], "x", "y", "alt") > summary(bmp) \end{Sinput} \begin{Soutput} OBJECTID OBJECTVAL OBJECTORIG YEAROFCHAN Min. :8641614 Hoehenkote :1583 LK25:1584 Min. :1985 1st Qu.:8642010 Seebodenkote: 1 1st Qu.:1985 Median :8642406 Median :1985 Mean :8642406 Mean :1985 3rd Qu.:8642801 3rd Qu.:1985 Max. :8643197 Max. :1985 x y alt Min. :777511 Min. :182009 Min. :1179 1st Qu.:785796 1st Qu.:184869 1st Qu.:1884 Median :794789 Median :187925 Median :2314 Mean :794951 Mean :187871 Mean :2269 3rd Qu.:804052 3rd Qu.:190775 3rd Qu.:2642 Max. :812480 Max. :193997 Max. :3411 \end{Soutput} \begin{Sinput} > class(bmp) \end{Sinput} \begin{Soutput} [1] "data.frame" \end{Soutput} \begin{Sinput} > coordinates(bmp) = ~x + y > class(bmp) \end{Sinput} \begin{Soutput} [1] "SpatialPointsDataFrame" attr(,"package") [1] "sp" \end{Soutput} \begin{Sinput} > spplot(bmp["alt"]) > grd = makegrid(bmp, n = 1e+05) > names(grd) = c("x", "y") > coordinates(grd) = ~x + y > grd = as(grd, "SpatialPixels") > library(gstat) > grd.alt = idw(alt ~ 1, bmp, grd) \end{Sinput} \begin{Soutput} [inverse distance weighted interpolation] \end{Soutput} \begin{Sinput} > print(spplot(grd.alt["var1.pred"])) \end{Sinput} \end{Schunk} \includegraphics{meteo-016} \section{Spatial prediction with a regression model} \begin{Schunk} \begin{Sinput} > grd.alt$alt = grd.alt$var1.pred > grd.alt$T.feb = predict(lm(f1, out), grd.alt) > grd.alt$T.aug = predict(lm(f2, out), grd.alt) > print(spplot(grd.alt[c("T.feb", "T.aug")])) \end{Sinput} \end{Schunk} \begin{Schunk} \begin{Sinput} > grd.alt$T.feb.se = predict(lm(f1, out), grd.alt, se.fit = T)$se.fit > grd.alt$T.aug.se = predict(lm(f2, out), grd.alt, se.fit = T)$se.fit > print(spplot(grd.alt[c("T.feb.se", "T.aug.se")])) \end{Sinput} \end{Schunk} \includegraphics{meteo-018} \begin{Schunk} \begin{Sinput} > print(spplot(grd.alt[c("alt")], main = "altitude"), split = c(1, + 1, 1, 2), more = T) > print(spplot(grd.alt[c("T.feb", "T.aug")], main = "temperature"), + split = c(1, 2, 1, 2), more = F) \end{Sinput} \end{Schunk} \includegraphics{meteo-019} \begin{Schunk} \begin{Sinput} > print(spplot(grd.alt[c("alt")], main = "altitude"), split = c(1, + 1, 1, 2), more = T) > print(spplot(grd.alt[c("T.feb.se", "T.aug.se")], main = "standard errors"), + split = c(1, 2, 1, 2), more = F) \end{Sinput} \end{Schunk} \includegraphics{meteo-020} \section{Ordinary kriging of altitude} \begin{Schunk} \begin{Sinput} > library(gstat) > v = variogram(alt ~ 1, bmp) > v \end{Sinput} \begin{Soutput} np dist gamma dir.hor dir.ver id 1 5839 569.9613 19670.18 0 0 var1 2 17454 1271.5887 59863.63 0 0 var1 3 27096 2077.2693 98639.08 0 0 var1 4 35526 2889.7851 123895.60 0 0 var1 5 42545 3707.5985 143416.36 0 0 var1 6 47600 4526.5666 163425.41 0 0 var1 7 52173 5344.1203 181295.12 0 0 var1 8 55361 6163.5497 204598.89 0 0 var1 9 56802 6983.8903 227578.54 0 0 var1 10 58288 7804.0773 250701.92 0 0 var1 11 57726 8623.4404 272643.68 0 0 var1 12 55889 9443.8364 288598.30 0 0 var1 13 53647 10263.7473 297153.43 0 0 var1 14 50750 11085.3437 303038.14 0 0 var1 15 46832 11904.9439 297291.43 0 0 var1 \end{Soutput} \begin{Sinput} > plot(v) > v.fit = fit.variogram(v, vgm(1, "Exp", 3000)) > v.fit \end{Sinput} \begin{Soutput} model psill range 1 Exp 413491.0 8538.357 \end{Soutput} \begin{Sinput} > print(plot(v, v.fit)) \end{Sinput} \end{Schunk} \includegraphics{meteo-021} When executed, the following will break with an obscure error message: \begin{Schunk} \begin{Sinput} > grd.alt.kr = krige(alt ~ 1, bmp, grd, v.fit, nmax = 30) \end{Sinput} \end{Schunk} which is due to a duplicate observation. We can find the duplicate measurement(s) that cause(s) this error: \begin{Schunk} \begin{Sinput} > zerodist(bmp) \end{Sinput} \begin{Soutput} [,1] [,2] [1,] 1193 1198 \end{Soutput} \end{Schunk} and check if the attribute is duplicate as well: \begin{Schunk} \begin{Sinput} > bmp[c(1193, 1198), ] \end{Sinput} \begin{Soutput} coordinates OBJECTID OBJECTVAL OBJECTORIG YEAROFCHAN alt 1193 (796161, 187155) 8642806 Hoehenkote LK25 1985 2322 1198 (796161, 187155) 8642811 Hoehenkote LK25 1985 2322 \end{Soutput} \end{Schunk} Now, we can krige while de-selecting one of the duplicates: \begin{Schunk} \begin{Sinput} > grd.alt.kr = krige(alt ~ 1, bmp[-1193, ], grd, v.fit, nmax = 30) \end{Sinput} \begin{Soutput} [using ordinary kriging] \end{Soutput} \begin{Sinput} > print(spplot(grd.alt.kr[1])) \end{Sinput} \end{Schunk} \begin{Schunk} \begin{Sinput} > print(spplot(grd.alt.kr[1], cuts = 25, main = "kriged surface"), + split = c(1, 1, 1, 2), more = TRUE) > grd.alt.kr$se = sqrt(grd.alt.kr$var1.var) > print(spplot(grd.alt.kr["se"], cuts = 25, main = "kriging standard error surface"), + split = c(1, 2, 1, 2), more = FALSE) \end{Sinput} \end{Schunk} \includegraphics{meteo-026} \end{document}