################################################### ### chunk number 1: ################################################### stations = read.csv("stations.csv") # or use e.g. stations = read.csv(file.choose()) lapply(stations, class) slf_ams = read.csv("slf_ams.csv") lapply(slf_ams, class) ################################################### ### chunk number 2: ################################################### levels(stations$X01.stationkennung) ################################################### ### chunk number 3: ################################################### unique(paste(slf_ams$X01.stationsabk.rzung, slf_ams$X02.standortnummer, sep="-")) ################################################### ### chunk number 4: ################################################### stations$station = stations$X01.stationkennung slf_ams$station = paste(slf_ams$X01.stationsabk.rzung, slf_ams$X02.standortnummer, sep="-") ################################################### ### chunk number 5: ################################################### slf_ams$X03.datum.und.zeit[1] ################################################### ### chunk number 6: ################################################### slf_ams$time = as.POSIXct(strptime(as.character(slf_ams$X03.datum.und.zeit), "%Y.%m.%d %H:%M")) ################################################### ### chunk number 7: ################################################### feb = format.POSIXct(slf_ams$time, "%m") == "02" summary(feb) feb[1:10] aug = format.POSIXct(slf_ams$time, "%m") == "08" ################################################### ### chunk number 8: ################################################### at1400 = format.POSIXct(slf_ams$time, "%H") == "14" & format.POSIXct(slf_ams$time, "%M") == "00" # 14:00 at1340 = format.POSIXct(slf_ams$time, "%H") == "13" & format.POSIXct(slf_ams$time, "%M") == "40" # 13:40 afternoon = at1400 | at1340 ################################################### ### chunk number 9: ################################################### temp_feb = na.omit(slf_ams[feb & afternoon, c("temp", "station")]) temp_aug = na.omit(slf_ams[aug & afternoon, c("temp", "station")]) ################################################### ### chunk number 10: ################################################### 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) ################################################### ### chunk number 11: ################################################### out = data.frame(T.aug = unlist(T.aug), T.feb = unlist(T.feb)) out ################################################### ### chunk number 12: ################################################### m = match(row.names(out), as.character(stations$X01.stationkennung)) m ################################################### ### chunk number 13: ################################################### all.equal(as.character(stations$X01.stationkennung[m]), row.names(out)) out$x = stations$X04.koord..x[m] out$y = stations$X05.koord..y[m] out$alt = stations$X06.hoehe.ueber.meer..m.[m] out ################################################### ### chunk number 14: ################################################### par(mfrow = c(2,2)) # plot 2 x 2 plots 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) # should show nothing. plot(T.feb ~ y, out) # North-South gradient? ################################################### ### chunk number 15: ################################################### library(rgdal) bm1197p = readOGR("../GISAufbaukurs_WS0708/data/1197-1198/DHM25/basis/point", "bm1197p") bm1198p = readOGR("../GISAufbaukurs_WS0708/data/1197-1198/DHM25/basis/point", "bm1198p") summary(bm1197p) summary(bm1198p) ################################################### ### chunk number 16: ################################################### library(lattice) trellis.par.set(sp.theme()) # sets bpy.colors() ramp bmp = rbind(bm1197p, bm1198p) bmp = as.data.frame(bmp) # modify last three column names: names(bmp) = c(names(bmp)[1:4], "x", "y", "alt") summary(bmp) coordinates(bmp) = ~x+y spplot(bmp["alt"]) grd = makegrid(bmp, n = 100000) names(grd) = c("x", "y") coordinates(grd)=~x+y grd = as(grd, "SpatialPixels") library(gstat) # needed for idw; can also do variogram modelling and kriging grd.alt = idw(alt~1, bmp, grd) print(spplot(grd.alt["var1.pred"])) ################################################### ### chunk number 17: ################################################### # create grid with Temperature prediction: grd.alt$alt = grd.alt$var1.pred # make sure the predictor "alt" has the correct name 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")])) ################################################### ### chunk number 18: ################################################### # create grid with Temperature prediction standard errors: 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")])) ################################################### ### chunk number 19: ################################################### 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) ################################################### ### chunk number 20: ################################################### 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) ################################################### ### chunk number 21: ################################################### library(gstat) v = variogram(alt~1,bmp) v plot(v) v.fit = fit.variogram(v,vgm(1, "Exp", 3000)) v.fit print(plot(v, v.fit)) ################################################### ### chunk number 22: ################################################### grd.alt.kr = krige(alt~1, bmp, grd, v.fit, nmax = 30) # find the duplicate measurement(s) that cause(s) this error: zerodist(bmp) # check if the attribute is duplicate as well: bmp[c(1193,1198),] # krige while de-selecting one of the duplicates: grd.alt.kr = krige(alt~1, bmp[-1193,], grd, v.fit, nmax = 30) print(spplot(grd.alt.kr[1])) ################################################### ### chunk number 23: ################################################### # increase the number of colours to find artifacts in valleys: 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)