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:
- collapes multi white space to a single one
- Reformat the linebreak
- Insert a " before the first test
- Insert " between words
- Remove a white space infront of the names
- Insert a " at the end of the line
- Remove the wrong inserted " in the header line
- 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)