German Weather Service (DWD) Stations

Since summer 2017 the German weather service (Deutscher Wetterdienst, DWD) made it’s station data (due to a modified law) publicly available. Here is a simple script how to visualize the stations in Germany.

Since summer 2017 the German weather service (Deutscher Wetterdienst, DWD) made it’s station data (due to a modified law) publicly available. Here are two starting points for more information:

I used the list of stations with theirs geographic positions, to create a station density map with OpenSteetMap (OSM) in the background with ggplot in Rstats. Before mapping the station density, I transfomed the geographic coordinates to the equal distance UTM projection.

And I had to parse their text file with the unix commands sed and tr to transform it to a standard format readable by R:

  1. collapes multi white space to a single one
  2. Reformat the linebreak
  3. Insert a " before the first test
  4. Insert " between words
  5. Remove a white space infront of the names
  6. Insert a " at the end of the line
  7. Remove the wrong inserted " in the header line
  8. Remove the underline line.
tr -s " "  < KL_Tageswerte_Beschreibung_Stationen.txt | \
  tr -d '\r' | \
  sed -e 's/[A-ZÄÖÜ]/ "&/' | \
  sed -e 's/ [A-Za-zäöüÄÖÜ\-]* $/" "&/' | \
  sed -e 's/" " /" "/' | \
  sed -e 's/ $/"/' | \
  sed -e 's/^ "//' | \
  sed -e '/----/d'> KL_Tageswerte_Beschreibung_Stationen.csv
library(OpenStreetMap)
library(ggplot2)
library(viridisLite)
library(viridis)

## Read the table of all stations
data.dir = file.path(rprojroot::find_rstudio_root_file(), "data")
stations = read.table(file.path(data.dir, "KL_Tageswerte_Beschreibung_Stationen.csv"), header=TRUE)

bbox = c(min(stations$geoBreite) - 1,
         min(stations$geoLaenge) - 1,
         max(stations$geoBreite) + 1,
         max(stations$geoLaenge) + 1)

utm = rgdal::project(as.matrix(stations[c("geoLaenge", "geoBreite")]), proj = "+init=epsg:32632")
colnames(utm) <- c("Easting", "Northing")
stations = cbind(stations, utm)

map = openmap(bbox[1:2], bbox[3:4], type = "apple-iphoto", zoom=6)
map = openproj(map, projection = "+init=epsg:32632")
gg = autoplot(map)
gg = gg + coord_cartesian()
gg = gg + stat_binhex(data=stations, aes(Easting, Northing), alpha=0.7)
gg = gg + scale_fill_gradientn(colours = rev(viridis(255)))
print(gg)