Turning R into a GIS – Mapping the weather in Germany

Posted on January 29, 2015 by Bigdata Doc in R bloggers | 0 Comments

[This article was first published on Big Data Doctor » R, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Nothing has gotten more attention in the visualization world like the map-based insights, or in other words, just plotting on a map different KPIs to allow for a playful discovery experience. I must admit, maps are cool, an awesome tool to “show-off” and to visually gain some insights.

But let’s be also clear about the limitations of map based charts:

You can compare locations based on a KPI, but you cannot quantify the difference between them
Color is difficult to understand and often leads to misinterpretation (e.g: what’s the meaning of red? more sales? or worst results?).
Color gradients are also pretty challenging for the human eye.
Zooming in/out results in insights detailing down/aggregation, but it’s difficult to establish a quantification between different granularity levels.

Anyways, R can be really useful to create high-quality maps… There are awesome packages like rMaps, where you have a set of controls available to make your maps interactive, rgooglemaps, maptools, etc.

In this post I’m going to plot weather KPIs for over 8K different postal codes (Postleitzahl or PLZ) in Germany. I’m going to shade the different areas according to their values -as you would expect –

We are going to follow these steps to visualize the temperature, the humidity and the snow fall for the entire German country:

Preparation of the required assets (PLZ coordinates, polygon lines, weather API key, etc)
Querying the weather API for each PLZ to retrieve the weather values
Map creation and PLZ data frame merging with the obtained weather information
Map display for the weather metrics and high-resolution picture saving

1- Assets preparation

We need to prepare a few assets… Everything freely accessible and just a mouse click away… Amazing, isn’t it?

The list of all PLZ with city name and the lat/long coordinates of a centroid (you can download this data from geonames)
The shapefiles for the PLZ to know how to draw them on a map (kindly made available out of the OpenStreetMaps at suche-postleitzahl.org)
A key for the weather API (you need to register at openweathermap.org, takes literally a second and they are not going to bother you with newsletters)

2-Downloading the weather data

Basically, it’s just a JSON call we can perform for each PLZ passing the lat/long coordinates to the openweather api’s endpoint. Each weather entry is then stored as a 1 row data frame we keep appending to the one holding all entries:

library(jsonlite)
#load the plz info you download from the geonames resource
plz.ort.de<-read.csv(file = "../plzgeo.csv")
weather.de<-NULL
for (i in 1:nrow(plz.ort.de))
{
  url<-paste0('http://api.openweathermap.org/data/2.5/weather?lat=',plz.ort.de[i,]$lat, '&lon=',plz.ort.de[i,]$lon,'&units=metric&APPID=PUT_YOUR_KEY_HERE')
  weather.entry<-jsonlite::fromJSON(url,simplifyMatrix = F,simplifyDataFrame = F,flatten = T)
  temperature<-weather.entry$main$temp
  humidity<-weather.entry$main$humidity
  wind.speed<-weather.entry$wind$speed
  wind.deg<-weather.entry$wind$deg
  snow<-weather.entry$snow$`3h`
  if (is.null(wind.speed)){ wind.speed<-NA}
  if (is.null(wind.deg)){ wind.deg<-NA}
  if (is.null(snow)){ snow<-NA}
  if (is.null(humidity)){ humidity<-NA}
  if (is.null(temperature)){ temperature<-NA}
  weather.de<-rbind(data.frame(plz=plz.ort.de[i,]$plz,temperature,humidity,wind.speed,wind.deg,snow),weather.de)  
#you might want to take your process to sleep for some milliseconds to give the API a breath
}

3-Map creation and PLZ-weather data frames merging

Using the rgal for the required spatial transformations. In this case, we use the EPSG 4839 for the German geography (see spTransform)

library(ggplot2)
library(rgdal)           # required for readOGR and spTransform
library(RColorBrewer)
 
setwd("[your_path]/plz-gebiete.shp/")
# read shapefile
wmap <- readOGR(dsn=".", layer="plz-gebiete")
 
map <- readOGR(dsn=".", layer="plz-gebiete")
map <- spTransform(map,CRS=CRS("+init=epsg:4839"))
map$region<-substr(map$plz, 1,1)
map.data <- data.frame(id=rownames(map@data), map@data)
map.data$cplz<- as.character(map.data$plz)
 
weather.de$cplz<- as.character(weather.de$plz)
#normalization to force all PLZs having 5 digits
weather.de$cplz<- ifelse(nchar(weather.de$cplz)<5, paste0("0",weather.de$cplz), weather.de$cplz)
map.data<-merge(map.data,weather.de,by=c("cplz"),all=T)
map.df   <- fortify(map)
map.df   <- merge(map.df,map.data,by="id", all=F)

4-Map display for the weather metrics and high-resolution picture saving

We just rely on the standard ggplot functionality to plot the weather metric we’d like to. To make it more readable, I facetted by region.

temperature<-ggplot(map.df, aes(x=long, y=lat, group=group))+
  geom_polygon(aes(fill=temperature))+
  facet_wrap(~region,scales = 'free') +
  geom_path(colour="lightgray", size=0.5)+
  scale_fill_gradient2(low ="blue", mid = "white", high = "green", 
                       midpoint = 0, space = "Lab", na.value = "lightgray", guide = "legend")+  theme(axis.text=element_blank())+
  theme(axis.text=element_text(size=12)) +
  theme(axis.title=element_text(size=14,face="bold")) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
  theme(panel.background = element_rect(fill = 'white')) +
  theme(panel.grid.major = element_line( color="snow2")) 
 
ggsave("../plz-temperature-300.png",  width=22.5, height=18.25, dpi=300)

5-(Bonus) Underlying map tiles

You probably feel like having a map as reference to see city names, roads, rivers and all that stuff in each PLZ. For that we can use ggmap, a really cool package for spatial visualization with Google Maps and OpenStreetMap.

library(plyr)
library(ggmap)
# reading the shapes
area <- readOGR(dsn=".", layer="plz-gebiete")
# using the normalized version with "0" for the later join
weather.de$plz<-weather.de$cplz
# from factor to character
area.df$plz<-as.character(area.df$plz)
area.df <- data.frame(id=rownames(area@data), area@data)
# merging weather and geographical information
area.extend<-merge(area.df, weather.de, by=c("plz"),all=F)
# building 
area.points <- fortify(area)
area.points<-merge(area.points, area.extend, by=c("id"),all=F)
d <- join(area.points, area.extend, by="id")
# region extraction
d$region<-substr(d$plz, 1,1)
bavaria<-subset(d, region=="8")
# google map tiles request... location is where you want your map centered at
google.map <- get_map(location="Friedberg", zoom =8, maptype = "terrain", color = "bw", scale=4)
ggmap(google.map) +
  geom_polygon(data=bavaria, aes(x=long, y=lat, group=group, fill=temperature), colour=NA, alpha=0.5) +
  scale_fill_gradient2(low ="blue", mid = "yellow", high = "green", 
                       midpoint = 0, space = "Lab", na.value = "lightgray", guide = "legend")+  theme(axis.text=element_blank())+  
  labs(fill="") +
  theme_nothing(legend=TRUE)
 
ggsave("../plz-temperature-Bavaria.png",  width=22.5, height=18.25, dpi=300)

The results speak for themselves!
Temperature in Germany

Temperature in the area around Munich only

Snow across Germany