This document describes the calculation of climate trends at different wheather stations in the Bavarian Alps, that we published in our project Schnee Von Morgen. The results of our calculations are thus reproducible by knitting ./index.Rmd, which is also part of the archive including the relevant CSV-files.

Data description

The files in the input folder serve as input for the calculations described below:

The files in the output folder contain the results of the process after the calculations are finished. They are identical to the ones, which can be downloaded as zip-file in the info-layer of our application.

Calculations

First of all we load the csv-files that serve as input for the calculations.

climate <- read.csv("./input/klima.csv", na.strings=".", stringsAsFactors = FALSE)
stations <- read.csv("./input/stationIDs.csv", stringsAsFactors = FALSE)
trends <- read.csv("./input/trends.csv", stringsAsFactors = FALSE)

After having loaded necessary libraries in the background, we calculate the linear regression and the Mann-Kendall-Trend-Tests for each station and the period of maximum length.

perStation <- climate %>% 
  left_join(stations, by = c("id" = "key")) %>%
  filter(date<=trendEnde) %>%
  group_by(id) %>% 
  arrange(date) %>%
  do(
    tempera1 = lm(temperatur ~ seq(1961, 1960 + length(temperatur)), data= .),
    tempera2 = MannKendall(ts(.$temperatur, 1961, 1960 + length(.$temperatur), 1)),
    schneeh1 = lm(schneehoehe ~ seq(1961, 1960 + length(schneehoehe)), data= .),
    schneeh2 = MannKendall(ts(.$schneehoehe, 1961, 1960 + length(.$schneehoehe), 1)),  
    schneet1 = lm(schneetage ~ seq(1961, 1960 + length(schneetage)), data= .),
    schneet2 = MannKendall(ts(.$schneetage, 1961, 1960 + length(.$schneetage), 1))  
  ) 

Now we’re extracting the slope and the intercept from the linear models and the significance value of the 2-sided Mann-Kendall-Tests.

perStation <- perStation %>%  
  mutate(
    slope_tempera=summary(tempera1)$coeff[2],
    intercept_tempera=summary(tempera1)$coeff[1],
    S_tempera= tempera2$S,
    varS_tempera = tempera2$varS,
    p_value_tempera=tempera2$sl,
    sig_tempera=(1-p_value_tempera)*100
    #sig2=1-2*(1-pnorm(abs((S-1)/(sqrt(54*53*(2*54+5)/18)))))
  )

perStation <- perStation %>%  
  mutate(
    slope_schneeh=summary(schneeh1)$coeff[2],
    intercept_schneeh=summary(schneeh1)$coeff[1],
    S_schneeh= schneeh2$S,
    varS_schneeh = schneeh2$varS,
    p_value_schneeh=schneeh2$sl,
    sig_schneeh=(1-p_value_schneeh)*100
  )
  
perStation <- perStation %>%  
  mutate(
    slope_schneet=summary(schneet1)$coeff[2],
    intercept_schneet=summary(schneet1)$coeff[1],
    S_schneet= schneet2$S,
    varS_schneet = schneet2$varS,
    p_value_schneet=schneet2$sl,
    sig_schneet=(1-p_value_schneet)*100
  ) 

See the results e.g. of the calculated significance values in percent:

perStation %>%
  select(
    id, sig_tempera, sig_schneeh, sig_schneet
  )
## Source: local data frame [7 x 4]
## Groups: <by row>
## 
##      id sig_tempera sig_schneeh sig_schneet
##   (int)       (dbl)       (dbl)       (dbl)
## 1  1550    98.49881    91.24526    98.05491
## 2  3307    86.12456    72.73195    75.96270
## 3  3730    98.88045    98.88045    99.92794
## 4  4125    99.94023    94.04047    96.81707
## 5  5467    99.13749    99.75389    84.63989
## 6  5792    99.67115    93.12925    47.34111
## 7  5941    97.44182    99.44692    96.37140

Now, we prepare the data that is written to the output-files.

stations <- stations %>%
  left_join(perStation, by = c("key" = "id")) %>%
  mutate(
    temperaS = sig_tempera,
    schneehS = sig_schneeh,
    schneetS = sig_schneet,
    trendEnde = substring(trendEnde, 1, 4)
  ) %>%
  select(key,name,hoehe,trendEnde,temperaS,schneehS,schneetS)
  
trends <- trends %>%
  left_join(perStation, by = c("id" = "id")) %>%
  mutate(
    tempera = as.numeric(substring(date, 1, 4)) * slope_tempera + intercept_tempera,
    schneeh = as.numeric(substring(date, 1, 4)) * slope_schneeh + intercept_schneeh,
    schneet = as.numeric(substring(date, 1, 4)) * slope_schneet + intercept_schneet
  ) %>%
  select(
    id, date, tempera, schneeh, schneet
  )

trends
##      id       date    tempera    schneeh   schneet
## 1  5467 1961-01-01 -3.2537026 103.090243 177.12471
## 2  5467 2010-01-01 -1.8490696  48.578553 167.67529
## 3  1550 1961-01-01  0.3288514  18.031445 121.64310
## 4  1550 2014-01-01  1.6831482   9.222206  93.20875
## 5  3307 1961-01-01  1.0014342  18.494307 110.30541
## 6  3307 2014-01-01  1.9357054  12.157861  96.83336
## 7  5792 1961-01-01 -9.9212199 276.913153 181.28956
## 8  5792 2014-01-01 -8.5052880 219.753263 181.11785
## 9  3730 1961-01-01 -0.2533847  30.719244 138.39259
## 10 3730 2014-01-01  1.1484374  11.900981 102.01481
## 11 5941 1961-01-01 -0.5170069  47.187684 149.17940
## 12 5941 2002-01-01  0.6551565  20.238097 123.96346
## 13 4125 1961-01-01  1.2899504  10.691365  92.75778
## 14 4125 2005-01-01  3.7623961   4.426972  63.83882

Finally, we write the results to the CSV-files:

write.csv(trends, "./output/trends.csv", quote=FALSE, row.names=FALSE)
write.csv(stations, "./output/stationIDs.csv", quote=FALSE, row.names=FALSE)
file.copy("./input/klima.csv", "./output/klima.csv", overwrite=TRUE)
## [1] TRUE