16 - use case: plot all rainfall values around a given point

16.1 Find meteo stations around a given point

m <- nearbyStations(49.211784, 9.812475, radius=30,
    res=c("daily","hourly"), var=c("precipitation","more_precip","kl"),
    mindate=as.Date("2016-05-30"), statname="Braunsbach catchment center")
# Remove duplicates. if kl and more_precip are both available, keep only more_precip:
library("berryFunctions")
m <- sortDF(m, "var")
m <- m[!duplicated(paste0(m$Stations_id, m$res)),]
m <- sortDF(m, "res")
m <- sortDF(m, "dist", decreasing=FALSE)
rownames(m) <- NULL
DT::datatable(m, options=list(pageLength=5, scrollX=TRUE))

Interactive map of just the meteo station locations:

library(leaflet)
m$col <- "red" ; m$col[1] <- "blue"
leaflet(m) %>% addTiles() %>%
  addCircles(lng=9.812475, lat=49.211784, radius=30e3) %>%
  addCircleMarkers(~geoLaenge, ~geoBreite, col=~col, popup=~Stationsname)

16.2 Download and process data

Download and process data for the stations, get the rainfall sums of a particular day (Braunsbach flood May 2016):

prec <- dataDWD(m$url, fread=TRUE, dir="../localdata")
names(prec) <- m$Stations_id[-1]
prec29 <- sapply(prec[m$res[-1]=="daily"], function(x)
         {
         if(nrow(x)==0) return(NA)
         col <- "RS"
         if(!col %in% colnames(x)) col <- "R1"
         if(!col %in% colnames(x)) col <- "RSK"
         x[x$MESS_DATUM==as.POSIXct(as.Date("2016-05-29")), col]
         })
prec29 <- data.frame(Stations_id=names(prec29), precsum=unname(prec29))
prec29 <- merge(prec29, m[m$res=="daily",c(1,4:7,14)], sort=FALSE)
head(prec29[,-7]) # don't show url column with long urls
Stations_idprecsumStationshoehegeoBreitegeoLaengeStationsname
2848105  45849.39.86Langenburg-Atzenrod
278772  35449.29.68Kupferzell-Rechbach
520682.239649.19.9 Vellberg-Kleinaltdorf
257594  42649.29.98Kirchberg/Jagst-Herboldshausen
341682.529449.39.81Mulfingen/Jagst
626076.638549.39.7 Ingelfingen-Stachenhausen

7 of the files contain no rows. readDWD warns about this (but the warnings are suppressed in this website).
One example is daily/more_precip/historical/tageswerte_RR_07495_20070114_20181231_hist.zip

16.3 Plot rainfall sum on map

For a quick look without a map, this works:

plot(geoBreite~geoLaenge, data=m, asp=1)
textField(prec29$geoLaenge, prec29$geoBreite, prec29$precsum, col=2)

But it’s nicer to have an actual map. If OSMscale installation fails, go to https://github.com/brry/OSMscale#installation

library(OSMscale)
map <- pointsMap(geoBreite,geoLaenge, data=m, type="osm", plot=FALSE)
pp <- projectPoints("geoBreite", "geoLaenge", data=prec29, to=map$tiles[[1]]$projection)
prec29 <- cbind(prec29,pp) ; rm(pp)
pointsMap(geoBreite,geoLaenge, data=m, map=map, scale=FALSE)
textField(prec29$x, prec29$y, round(prec29$precsum), font=2, cex=1.5)
scaleBar(map, cex=1.5, type="line", y=0.82)
title(main="Rainfall sum  2016-05-29  7AM-7AM  [mm]", line=-1)