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.
The files in the input folder serve as input for the calculations described below:
./input/klima.csv contains the average values of temperature, snowheight and snowdays in the winter half-year since 1961. The numbers trace back to the daily values measured by the DWD../input/trends.csv contains the beginning and end of the periods, for which the trend and its significance is calculated. The longest period with no more than 2 missing values per winter half-year since 1961 is chosen../input/stationIDs.csv gives information about the wheather stations at which these indicators were measured.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.
./output/klima.csv is identical to the input file above../output/trends.csv contains the values of the indicators at the beginning / end of each period under consideration - according to the calculated linear regression../output/stationIDs.csv ressambles the input file, complemented with the significance values of each indicator according to the Mann-Kendall-Trends.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.37140Now, 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.83882Finally, 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