library(dplyr)
library(tidyr)
library(ggplot2)
library(ggtext)

theme_set(theme_minimal())

Introduction

Contemplating about the weather, I wondered if I could find out the “most unusual” and “most ideal” years regarding air temperature in Germany, i.e. if I could identify the years in which the daily temperature deviated the most and the least from the expected seasonal temperature. So I decided to look into historical climate data, created an extremely simplified seasonal temperature model and then investigated the deviations from that model. Although it’s all quite simple, this little exploration gives some insights into how and why we can use a linear model for such data.

You can find the source code repository for this project on GitHub. The code is licensed under Apache License 2.0. See LICENSE file.

Data

I retrieved the historical climate data for a weather station in Berlin-Dahlem (a bit outside the city) from 1950 to now from the German Meteorological Service (Deutscher Wetterdienst – DWD). The data come as delimited files with semicolon as column separator. Historical data until 2022 and present data from 2022 to now come as separate files.

raw_hist <- read.delim('data/produkt_klima_tag_19500101_20221231_00403.txt', sep = ';')
head(raw_hist)
##   STATIONS_ID MESS_DATUM QN_3   FX   FM QN_4  RSK RSKF  SDK SHK_TAG  NM VPM
## 1         403   19500101 -999 -999 -999    5  2.2    7 -999       0 5.0 4.0
## 2         403   19500102 -999 -999 -999    5 12.6    8 -999       0 8.0 6.1
## 3         403   19500103 -999 -999 -999    5  0.5    1 -999       0 5.0 6.5
## 4         403   19500104 -999 -999 -999    5  0.5    7 -999       0 7.7 5.2
## 5         403   19500105 -999 -999 -999    5 10.3    7 -999       0 8.0 4.0
## 6         403   19500106 -999 -999 -999    5  7.2    8 -999      12 7.3 5.6
##       PM  TMK UPM  TXK  TNK  TGK eor
## 1 1025.6 -3.2  83 -1.1 -4.9 -6.3 eor
## 2 1005.6  1.0  95  2.2 -3.7 -5.3 eor
## 3  996.6  2.8  86  3.9  1.7 -1.4 eor
## 4  999.5 -0.1  85  2.1 -0.9 -2.3 eor
## 5 1001.1 -2.8  79 -0.9 -3.3 -5.2 eor
## 6  997.5  2.6  79  5.0 -4.0 -4.0 eor
raw_pres1 <- read.delim('data/produkt_klima_tag_20221107_20240509_00403.txt', sep = ';')
head(raw_pres1)
##   STATIONS_ID MESS_DATUM QN_3   FX   FM QN_4 RSK RSKF SDK SHK_TAG  NM  VPM
## 1         403   20221107 -999 -999 -999   10 0.0    6 4.5       0 6.2  9.6
## 2         403   20221108 -999 -999 -999   10 0.2    6 7.5       0 6.0 10.4
## 3         403   20221109 -999 -999 -999   10 1.0    6 3.7       0 6.6 11.4
## 4         403   20221110 -999 -999 -999   10 0.0    0 6.1       0 5.1 10.2
## 5         403   20221111 -999 -999 -999   10 0.0    0 1.9       0 6.3  9.6
## 6         403   20221112 -999 -999 -999   10 0.0    0 7.3       0 4.0  8.8
##       PM  TMK UPM  TXK TNK  TGK eor
## 1 1002.9 10.7  75 15.0 6.4  5.1 eor
## 2 1002.7 12.1  75 16.9 7.9  4.2 eor
## 3 1001.5 11.8  83 15.0 9.0  5.1 eor
## 4 1012.6 11.7  74 14.3 8.6  5.8 eor
## 5 1020.1  8.6  87 12.8 4.0  0.6 eor
## 6 1022.8  6.4  92 13.8 1.8 -0.9 eor
raw_pres2 <- read.delim('data/produkt_klima_tag_20230203_20240805_00403.txt', sep = ';')
head(raw_pres2)
##   STATIONS_ID MESS_DATUM QN_3   FX   FM QN_4 RSK RSKF SDK SHK_TAG  NM VPM
## 1         403   20230203 -999 -999 -999   10 5.1    8 1.8       0 7.4 7.8
## 2         403   20230204 -999 -999 -999   10 0.0    0 5.4       0 3.8 5.2
## 3         403   20230205 -999 -999 -999   10 0.1    7 3.8       0 4.8 4.1
## 4         403   20230206 -999 -999 -999   10 0.8    7 0.0       0 5.4 4.8
## 5         403   20230207 -999 -999 -999   10 0.0    0 8.9       0 1.3 4.0
## 6         403   20230208 -999 -999 -999   10 0.0    0 8.9       0 0.4 4.0
##       PM  TMK UPM TXK  TNK   TGK eor
## 1 1004.3  5.5  87 8.6  1.4   0.4 eor
## 2 1024.5  1.3  76 4.1 -2.1  -5.3 eor
## 3 1028.5 -2.4  82 1.3 -5.3  -7.3 eor
## 4 1029.7 -2.8  95 0.2 -5.9  -9.4 eor
## 5 1032.2 -3.5  85 2.5 -8.0 -10.6 eor
## 6 1031.0 -2.5  80 2.4 -6.5  -8.6 eor

After reading in the files, we merge them, select only the necessary variables, transform the dates and remove duplicates (since the historical and the present data both contain observations from 2022) to generate our final measurements dataset meas:

meas <- bind_rows(raw_hist, raw_pres1, raw_pres2) |>
    select(date = MESS_DATUM, temp = TMK) |>   # TMK is day-time average temperature in °C
    mutate(date = as.POSIXct(strptime(date, "%Y%m%d")),
           year = as.integer(as.numeric(format(date, "%Y"))),
           day = as.integer(as.numeric(format(date, "%j")))) |> # day of the year as decimal number from 1 to 366
    distinct(date, .keep_all = TRUE)   # remove duplicates
rm(raw_hist, raw_pres1, raw_pres2)         # don't need the raw data any more
stopifnot(all(count(meas, date)$n == 1))   # make sure there are no duplicates
head(meas)
##         date temp year day
## 1 1950-01-01 -3.2 1950   1
## 2 1950-01-02  1.0 1950   2
## 3 1950-01-03  2.8 1950   3
## 4 1950-01-04 -0.1 1950   4
## 5 1950-01-05 -2.8 1950   5
## 6 1950-01-06  2.6 1950   6

Visual analysis

Let’s visualize the time series with a simple plot. I will also add a smoothed curve showing an overall trend, which indicates a nearly linear increase in average annual temperature by about 2°C since the 1950’s. I’ll later come back to that. We can also see the typical seasonal changes.

ggplot(meas, aes(date, temp)) +
    geom_line() +
    geom_smooth(method = "gam") +
    labs(title = "Daily day-time average temperature in Berlin-Dahlem over time",
         x = "",
         y = "Temperature in °C")