Vorbemerkungen

Dieses Dokument beschreibt die Vorprozessierung und explorative Analyse des Datensatzes, der Grundlage des auf srf.ch veröffentlichten Artikel Diesen Effekt hat der Klimawandel auf Ihren Wohnort ist.

SRF Data legt Wert darauf, dass die Datenvorprozessierung und -Analyse nachvollzogen und überprüft werden kann. SRF Data glaubt an das Prinzip offener Daten, aber auch offener und nachvollziehbarer Methoden. Zum anderen soll es Dritten ermöglicht werden, auf dieser Vorarbeit aufzubauen und damit weitere Auswertungen oder Applikationen zu generieren.

Die Endprodukte des vorliegenden Scripts, neben der vorliegenden explorativen Analyse, sind JSON-Dateien zu den Klimaszenarien für jede Gemeinde und ein JSON-File mit allen Metainformationen zu den Gemeinden (Datenbeschreibung siehe unten):

  • output/climate_projections_bfsId.json: Das Skript generiert 2212 Dateien, für jede Gemeinde eine Datei. Die Dateien tragen jeweils die Gemeindenummer des Bundesamts für Statistik (BfS) im Dateinamen. In den Dateien sind alle Klimavariablen enthalten, zwei simulierte Szenarien, alle drei Zeitperioden und alle drei Schätzungen.

  • output/municipalities.json: Metainformationen zu allen Gemeinden, Gemeindestand 1. Januar 2019.

R-Script & Daten

Die Vorprozessierung und Analyse wurde im Statistikprogramm R vorgenommen. Das zugrunde liegende Script sowie die prozessierten Daten können unter diesem Link heruntergeladen werden. Durch Ausführen von main.Rmd kann der hier beschriebene Prozess nachvollzogen und der für den Artikel verwendete Datensatz generiert werden. Dabei werden Daten aus dem Ordner input eingelesen und Ergebnisse in den Ordner output geschrieben.

SRF Data verwendet das rddj-template von Timo Grossenbacher als Grundlage für seine R-Scripts. Entstehen bei der Ausführung dieses Scripts Probleme, kann es helfen, die Anleitung von rddj-template zu studieren.

Debug-Informationen: This report was generated on 2019-11-30 10:31:44. R version: 3.5.3 on x86_64-apple-darwin15.6.0. For this report, CRAN packages as of 2019-03-01 were used.

GitHub

Der Code für die vorliegende Datenprozessierung ist auf https://github.com/srfdata/2019-11-auswirkungen-klimawandel zur freien Verwendung verfügbar.

Lizenz

Creative Commons Lizenzvertrag
2019-11-auswirkungen-klimawandel von SRF Data ist lizenziert unter einer Creative Commons Namensnennung - Weitergabe unter gleichen Bedingungen 4.0 International Lizenz.

Weitere Projekte

Code & Daten von SRF Data sind unter https://srfdata.github.io verfügbar.

Haftungsausschluss

Die veröffentlichten Informationen sind sorgfältig zusammengestellt, erheben aber keinen Anspruch auf Aktualität, Vollständigkeit oder Richtigkeit. Es wird keine Haftung übernommen für Schäden, die durch die Verwendung dieses Scripts oder der daraus gezogenen Informationen entstehen. Dies gilt ebenfalls für Inhalte Dritter, die über dieses Angebot zugänglich sind.

Datenbeschreibung

output/climate_projections_bfsId.json

Für jede, der 2212 Gemeinden wird eine JSON-Datei generiert mit allen nötigen Klimavariablen, Szenarien, Perioden und Schätzungen.

Attribut Typ Beschreibung
bfs_id Number Gemeindenummer
key String Klimavariable (FD, snowdays, ID, SD, HD, TN, tas und pr)
period String Zeitpunkt, den die Simulation abbiulden soll (1981-2010, 2035, 2060, 2085)
estimate String Schätzwerte für die Berechnungen, q5, q50 und q95
rcp String Unterschiedliche Klimaszenarios (obs, RCP2.6, RCP8.5)
season String Wert für die Jahreszeit, nur vorhanden bei den Variablen pr und tas

output/municipality.json

Attribut Typ Beschreibung
bfs_id Number Gemeindenummer
name String Gemeindename
altitude String Höhenlage der Gemeinde
region String Grossregion der Gemeinde
urban Boolean Grosstadt oder Agglomeration

Originalquelle

Die Daten des interaktiven Artikel stammen hauptsächlich von den «Klimaszenarien CH2018», die im November 2018 veröffentlicht wurden (CH2018 Project Team (2018): CH2018 - Climate Scenarios for Switzerland. National Centre for Climate Services. doi: 10.18751/Climate/Scenarios/CH2018/1.0). Grundlage bilden Simulationen mit insgesamt 21 verschiedenen Computermodellen, die an europäischen Forschungsinstitutionen – koordiniert durch das Projekt EURO-CORDEX – betrieben werden. Das Projekt EURO-CORDEX ist eine Initiative des WCRP.

Die regionalen Klimasimulationen berücksichtigen drei verschiedene Szenarien, je nach Entwicklung der Treibhausgasemissionen:

Die regionalen Klimasimulationen berücksichtigen drei verschiedene Szenarien, je nach Entwicklung der Treibhausgasemissionen:

  • Kein Klimaschutz (RCP8.5): Die klimawirksamen Emissionen nehmen stetig zu (hier verwendet als «pessimistisches Szenario»)

  • Konsequenter Klimaschutz (RCP2.6): Durch rasche und drastische Senkungen des Treibhausgasausstosses kann bis in 20 Jahren der Anstieg der Treibhausgase in der Atmosphäre gestoppt werden (hier verwendet als «optimistisches Szenario»). Damit lassen sich die Ziele des Pariser Klimaabkommens von 2015 wahrscheinlich erreichen und die globale Erwärmung auf zwei Grad Celsius gegenüber dem vorindustriellen Zustand begrenzen.

  • Begrenzter Klimaschutz (RCP4.5): Eine mittlere Entwicklung mit begrenztem Klimaschutz (hier nicht verwendet)

Die Klimaszenarien CH2018 beziehen sich jeweils auf einen Mittelwert der geschätzten klimatischen Verhältnisse über einen längeren Zeitraum. Wenn es im Text «2035» heisst, bezieht sich das auf die nahe Zukunft von 2020-2049. Ist die Rede von «2060», geht es um die Mitte des Jahrhunderts (2045-2074) und «2085» bezieht sich auf das Ende des Jahrhunderts (2070-2099). Als Referenzperiode gilt der Zeitraum von 1981 bis 2010. Werte von «heute» sind also gemittelte Messwerte über diesen Zeitraum.

Diese Projektionen der Klimamodelle sind jedoch nicht exakt, sondern streuen immer über einen gewissen Bereich. Der Median davon entspricht am ehesten dem absehbaren Wert und wird als «erwartetes» Ergebnis behandelt.

Gitternetz-Daten der Klimaszenarien 2018

Die Daten wurden im Dateiformat netCDF geliefert. Für jede Klimavariable gibt es einen Ordner.

input/FD_obs_CH2018

In diesem Ordner befinden sich die Dateien zu den Frosttagen.

input/HD_obs_CH2018

In diesem Ordner befinden sich die Dateien zu den Hitzetagen

input/ID_obs_CH2018

In diesem Ordner befinden sich die Dateien zu den Eistagen.

input/SD_obs_CH2018

In diesem Ordner befinden sich die Dateien zu den Sommertagen.

input/snowdays_obs_CH2018

In diesem Ordner befinden sich die Dateien zu den Schneetagen.

input/TN_obs_CH2018

In diesem Ordner befinden sich die Dateien zu den Tropennächten.

input/tas_obs_CH2018

In diesem Ordner befinden sich die Dateien zu den Tagesmittelwerten für die Temperatur nach Jahreszeit.

input/pr_obs_CH2018

In diesem Ordner befinden sich die Dateien zu den Tagesmittelwerten für Niederschlag nach Jahreszeit.

Klimatische Regionen

input/climate_regions

In diesem Ordner befinden sich die Dateien aufgeteilt und gemittelet nach Region. In jeder Datei sind alle Klimavariablen, Szenarien etc. vorhanden. Das Dateiformat ist ebenfalls netCDF. Die Daten stammen von MeteoSchweiz.

Gemeindegrenzen

input/gd-b-00.03-875-gg17

Die Shapedateien für die Gemeinden und Kantonen können hier runtergeladen werden: Generalisierte Gemeindegrenzen.

Grossregionen

input/Grossregionen Die Shapedatei zu den Grossregionen bildet die regionale Einteiung ab, die MeteoSchweiz für die Bewertung der Risiken und Chancen verwendet hat. Die Daten stammen von MeteoSchweiz.

Vorbereitungen

## [1] "package package:rmarkdown detached"
## Loading required package: knitr
## Loading required package: rstudioapi

Packages definieren

# from https://mran.revolutionanalytics.com/web/packages/checkpoint/vignettes/using-checkpoint-with-knitr.html
# if you don't need a package, remove it from here (commenting is probably not sufficient)
# tidyverse: see https://blog.rstudio.org/2016/09/15/tidyverse-1-0-0/
cat("
library(rstudioapi)
library(rgdal)
library(raster) # raster used for relief important: load before tidyverse
library(glue) # cooler string templating
library(tidyverse) # ggplot2, dplyr, tidyr, readr, purrr, tibble
library(magrittr) # pipes
library(readxl) # excel
library(scales) # scales for ggplot2
library(jsonlite) # json
library(lintr) # code linting
library(sf) # spatial data handling
library(ncdf4) # package to read netCDF files
library(ncdf4.helpers)
library(PCICt)
library(rmarkdown)",
file = "manifest.R")

Packages installieren

# if checkpoint is not yet installed, install it (for people using this
# system for the first time)
if (!require(checkpoint)) {
  if (!require(devtools)) {
    install.packages("devtools", repos = "http://cran.us.r-project.org")
    require(devtools)
  }
  devtools::install_github("RevolutionAnalytics/checkpoint",
                           ref = "v0.3.2", # could be adapted later,
                           # as of now (beginning of July 2017
                           # this is the current release on CRAN)
                           repos = "http://cran.us.r-project.org")
  require(checkpoint)
}
# nolint start
if (!dir.exists("~/.checkpoint")) {
  dir.create("~/.checkpoint")
}
# nolint end
# install packages for the specified CRAN snapshot date
checkpoint(snapshotDate = package_date,
           project = path_to_wd,
           verbose = T,
           scanForPackages = T,
           use.knitr = F,
           R.version = R_version)
rm(package_date)

Packages laden

source("manifest.R")
unlink("manifest.R")
sessionInfo()
## R version 3.5.3 Patched (2019-03-11 r76227)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Mojave 10.14.5
## 
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] de_CH.UTF-8/de_CH.UTF-8/de_CH.UTF-8/C/de_CH.UTF-8/de_CH.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] rmarkdown_1.11      PCICt_0.5-4.1       ncdf4.helpers_0.3-3
##  [4] ncdf4_1.16          sf_0.7-3            lintr_1.0.3        
##  [7] jsonlite_1.6        scales_1.0.0        readxl_1.3.0       
## [10] magrittr_1.5        forcats_0.4.0       stringr_1.4.0      
## [13] dplyr_0.8.3         purrr_0.3.0         readr_1.3.1        
## [16] tidyr_1.0.0         tibble_2.1.3        ggplot2_3.1.0      
## [19] tidyverse_1.2.1     glue_1.3.0          raster_2.8-19      
## [22] rgdal_1.3-9         sp_1.3-1            checkpoint_0.4.0   
## [25] rstudioapi_0.9.0    knitr_1.21         
## 
## loaded via a namespace (and not attached):
##  [1] tidyselect_0.2.5 xfun_0.5         haven_2.1.0      lattice_0.20-38 
##  [5] colorspace_1.4-0 vctrs_0.2.0      generics_0.0.2   htmltools_0.3.6 
##  [9] yaml_2.2.0       rlang_0.4.0      e1071_1.7-0.1    pillar_1.3.1    
## [13] DBI_1.0.0        withr_2.1.2      modelr_0.1.4     lifecycle_0.1.0 
## [17] plyr_1.8.4       munsell_0.5.0    gtable_0.2.0     cellranger_1.1.0
## [21] rvest_0.3.2      codetools_0.2-16 evaluate_0.13    rex_1.1.2       
## [25] class_7.3-15     broom_0.5.1      Rcpp_1.0.1       classInt_0.3-1  
## [29] backports_1.1.3  hms_0.4.2        digest_0.6.18    stringi_1.4.3   
## [33] grid_3.5.3       cli_1.0.1        tools_3.5.3      lazyeval_0.2.1  
## [37] crayon_1.3.4     pkgconfig_2.0.2  zeallot_0.1.0    xml2_1.2.0      
## [41] lubridate_1.7.4  assertthat_0.2.0 httr_1.4.0       R6_2.4.0        
## [45] units_0.6-2      nlme_3.1-137     compiler_3.5.3

Zusätzliche Scripts laden

# if you want to outsource logic to other script files, see README for 
# further information
knitr::read_chunk("scripts/my_script.R")
source("scripts/my_script.R")
my_function(5)
## [1] 5

Geographische Daten einlesen

Generalisierte Gemeindegrenzen und Grossregionen

# read municipal borders
municipality_geo <- read_sf(
  "input/gd-b-00.03-875-gg19/ggg_2019-LV95/shp/g2g19.shp",
  # set crs to lv95
  crs = 2056
)
## Warning: st_crs<- : replacing crs does not reproject data; use st_transform
## for that
# create new sf object with city centers
municipality_geo_point <- municipality_geo %>%
  st_drop_geometry() %>%
  st_as_sf(coords = c("E_CNTR", "N_CNTR"),  crs = 2056)

# shapefiles including urban areas
big_regions_shapefile <- read_sf(
  "input/Grossregionen/Grossraeume_neu.shp",
  crs = 21781
) %>%
  # transform from old ch1903 to new lv95
  st_transform(2056)
## Warning: st_crs<- : replacing crs does not reproject data; use st_transform
## for that

Daten zu den Klimaszenarien einlesen

In diesem Block werden die berechneten Klimaszenarien der Schweiz eingelesen.

Klimaregionen einlesen

Wir lesen die Daten zu den Klimaregionen ein und matchen diese dann ebenfalls mit den Gemeinden. Dadurch können wir die Gemeinden einer Klimaregion zuordnen.

climate_regions_geo_points <- c(
  "Alpen", 
  "AlpenS", 
  "Jura", 
  "Voralpen", 
  "Mittelland"
  ) %>%
  map_df(function(region) {
    # determine filepath
    filepath <- glue("input/climate_regions/{region}.nc")
    if (file.exists(filepath)) {
      # extract all values of the selected climate variable
      current_file <- nc_open(filepath)
      lon <- ncvar_get(current_file, varid = "lon")
      lat <- ncvar_get(current_file, varid = "lat")
      current_region <- ncvar_get(current_file, region)
      current_region_vec <- as.vector(current_region)
      lonlat <- as.matrix(expand.grid(lon, lat))
      lonlat_df <- data.frame(
        cbind(
          lonlat, 
          current_region_vec
          )
        ) %>%
        filter(
          current_region_vec > 0
          ) %>% 
          mutate(current_region_vec = case_when(
      current_region_vec == 1 ~ region
      ))
      nc_close(current_file)
      result <- tibble(
        region = lonlat_df$current_region,
        lon = lonlat_df$Var1,
        lat = lonlat_df$Var2
        )
    }
    return(result)
  }
  ) %>%
  st_as_sf(coords = c("lon", "lat"),  crs = 4326) %>%
  st_transform(2056)

# join the two datasets (matched and missing muncipalities)
climate_regions_all_municipalities <- st_join(
  municipality_geo_point %>% 
    select(GMDNR, GMDNAME, geometry),
  climate_regions_geo_points,
  join = st_nearest_feature
  ) %>% 
  st_drop_geometry() %>% 
  select(GMDNR, region)

Klimavariablen als Gitternetzpunkte einlesen

Wir lesen die Daten ein, machen einen Spatial-Join mit dem 2-km-Gitternetzpunkt, der einem Gemeindeortszentrum, am nächsten liegt.

lonlat_filename <- "input/FD_obs_CH2018/FD_RCP2.6_2035_q5.nc"
lonlat_file <- nc_open(lonlat_filename)
lon <- ncvar_get(lonlat_file, varid = "lon")
lat <- ncvar_get(lonlat_file, varid = "lat")
lonlat <- as.matrix(expand.grid(lon, lat))

# cleanup
rm(lonlat_filename, lonlat_file, lon, lat)

climate_projections_wide <- c(
  "FD", "HD", "ID", "SD", "snowdays", "TN", "pr", "tas"
  ) %>%
  map_dfc(function(climate_variable) {
    # map through different rcp scenarios
      c("RCP2.6", "RCP4.5", "RCP8.5") %>%
      map_dfc(function(rcp) {
        # map through different periods
         c("2035", "2060", "2085") %>%
          map_dfc(function(period) {
            # map through different estimates
            c("q5", "q50", "q95") %>%
              map_dfc(function(estimate) {
                 # determine folder
                folder <- glue("{climate_variable}_obs_CH2018/")
                if (climate_variable == "pr" || climate_variable == "tas") {
                  # map through different seasons
                  c("DJF", "MAM", "JJA", "SON") %>% 
                    map_dfc(function(season) {
                      # determine filename
                      filename <- glue(
                        "{climate_variable}_",
                        "{rcp}_",
                        "{period}_",
                        "{season}_",
                        "{estimate}.nc"
                        )
                      filepath <- glue("input/{folder}{filename}")
                      if (file.exists(filepath)) {
                        current_file <- nc_open(filepath)
                        # extract all values of the selected climate variable
                        current_climate_values <- ncvar_get(
                          current_file, 
                          climate_variable
                          )
                        current_climate_values_vec <- as.vector(
                          current_climate_values
                          )
                        nc_close(current_file)
                        col_name <- glue(
                          "{climate_variable}|",
                          "{rcp}|",
                          "{period}|",
                          "{estimate}|",
                          "{season}"
                          )
                        result <- tibble(
                          !!col_name := current_climate_values_vec
                          )
                       }
                    })
                  }
                  else {
                   # determine filename
                   filename <- glue(
                     "{climate_variable}_{rcp}_{period}_{estimate}.nc"
                     )
                   filepath <- glue("input/{folder}{filename}")
                   if (file.exists(filepath)) {
                       current_file <- nc_open(filepath)
                       # extract all values of the selected climate variable
                       current_climate_values <- ncvar_get(
                         current_file,
                         climate_variable
                         )
                       current_climate_values_vec <- as.vector(
                         current_climate_values)
                       nc_close(current_file)
                       col_name <- glue(
                         "{climate_variable}|{rcp}|{period}|{estimate}"
                         )
                       result <- tibble(
                         !!col_name := current_climate_values_vec
                       )
                   }
                  }
                   })
              })
          })
      }
    ) %>% 
  cbind(lonlat) %>%
  rename(
    lon = Var1,
    lat = Var2
    )

# for some weird reason we can't read the observational data in the same map
climate_observations_wide <- c(
  "FD", "HD", "ID", "SD", "snowdays", "TN", "pr", "tas"
  ) %>%
  map_dfc(function(climate_variable) {
    # map through different rcp scenarios
      c("obs") %>%
      map_dfc(function(rcp) {
        # map through different periods
         c("1981-2010") %>%
          map_dfc(function(period) {
            # map through different estimates
            c("mean", "yearly_mean") %>%
              map_dfc(function(estimate) {
                 # determine folder
                folder <- glue("{climate_variable}_obs_CH2018/")
                if (climate_variable == "pr" || climate_variable == "tas") {
                  # map through different seasons
                  c("DJF", "MAM", "JJA", "SON") %>% 
                    map_dfc(function(season) {
                      # determine filename
                      filename <- glue(
                        "{climate_variable}_",
                        "{rcp}_",
                        "{period}_",
                        "{season}_",
                        "{estimate}.nc"
                        )
                      filepath <- glue("input/{folder}{filename}")
                      if (file.exists(filepath)) {
                        current_file <- nc_open(filepath)
                        # extract all values of the selected climate variable
                        current_climate_values <- ncvar_get(
                          current_file, 
                          climate_variable
                          )
                        current_climate_values_vec <- as.vector(
                          current_climate_values
                          )
                        nc_close(current_file)
                        col_name <- glue(
                          "{climate_variable}|",
                          "{rcp}|",
                          "{period}|",
                          "{estimate}|",
                          "{season}"
                          )
                        result <- tibble(
                          !!col_name := current_climate_values_vec
                          )
                       }
                    })
                  }
                  else {
                   # determine filename
                   filename <- glue(
                     "{climate_variable}_{rcp}_{period}_{estimate}.nc"
                     )
                   filepath <- glue("input/{folder}{filename}")
                   if (file.exists(filepath)) {
                       current_file <- nc_open(filepath)
                       # extract all values of the selected climate variable
                       current_climate_values <- ncvar_get(
                         current_file, 
                         climate_variable
                         )
                       current_climate_values_vec <- as.vector(
                         current_climate_values
                         )
                       nc_close(current_file)
                       col_name <- glue(
                         "{climate_variable}|{rcp}|{period}|{estimate}"
                         )
                       result <- tibble(
                         !!col_name := current_climate_values_vec
                       )
                   }
                  }
                   })
              })
          })
      }
    ) %>% 
  cbind(lonlat) %>% 
  rename(
    lon = Var1, 
    lat = Var2
    )

# put the two data sets together (observations and projections)
climate_data_all_grid_points <- climate_projections_wide %>% 
  left_join(climate_observations_wide,
            by = c("lon" = "lon", "lat" = "lat")) %>% 
  group_by(lat, lon) %>% 
  drop_na()

# create geo point of lat lon from netcdf files
climate_values_geo_point <- climate_data_all_grid_points %>%
  st_as_sf(coords = c("lon", "lat"),  crs = 4326) %>%
  st_transform(2056)

# look for grid point that is closest (st_nearest feature) to the municipality
climate_values_all_municipalities <- st_join(
  municipality_geo_point %>% select(GMDNR, GMDNAME, geometry),
  climate_values_geo_point,
  join = st_nearest_feature
  ) %>% 
  st_drop_geometry() %>% 
  select(GMDNR, GMDNAME, everything())

# convert wide data frame into long data frame
climate_projections_by_municipality <- 
  climate_values_all_municipalities %>% 
  gather(string, value, -GMDNR, -GMDNAME) %>%
  # extract keys of climate variables
  tidyr::extract(
    col = "string",
    into = c("key", "rcp", "period", "estimate", "season"),
    regex = "(\\w+)\\|([\\w\\.]+)\\|([\\d-]+)\\|(\\w+)(?:\\|(\\w+))?",
    remove = TRUE
  ) %>%
  rename(bfs_id = GMDNR)

# join altitude information and climate regions
climate_projections_by_municipality %<>%
  left_join(
    municipality_geo %>% 
      st_drop_geometry() %>% 
      select(GMDNR, Z_CNTR), 
    by = c("bfs_id" = "GMDNR")
  ) %>% 
  mutate(
    altitude = case_when(
      Z_CNTR <= 800 ~ "0-800",
      Z_CNTR > 800 & Z_CNTR <= 1500 ~ "801-1500",
      Z_CNTR > 1500 ~ ">1500"
    )
  ) %>%
  select(-Z_CNTR) %>%
  # join climate regions from chunk above
  left_join(
    climate_regions_all_municipalities,
    by = c("bfs_id" = "GMDNR")
  )

#cleanup
rm(
  lonlat, 
  climate_projections_wide,
  climate_observations_wide,
  climate_data_all_grid_points, 
  climate_values_by_municipality,
  municipalities_with_no_match,
  climate_values_municipalities_no_match,
  climate_values_all_municipalities
  )

# save the data frame to an Rdata file
save(
  climate_projections_by_municipality, 
  file = "climate_projections_by_municipality.Rdata"
  )
save(
  climate_values_geo_point, 
  file = "climate_values_geo_point.Rdata"
  )
if (file.exists("climate_projections_by_municipality.Rdata")) {
  print("loading climate data")
  load("climate_projections_by_municipality.Rdata")
} else {
  print("please run the big chunk above")
}
## [1] "loading climate data"

Analyse

Zuerst werden unterschiedliche Funktionen erstellt, um dann in den Unterkapitel zu den einzelnen Variablen jeweils die Funktionen aufrufen zu können.

# function to translate keys to word describing the climate variable
translate_key_to_word <- function(key_id) {
  case_when(
    key_id == "FD" ~ "Frosttage",
    key_id == "HD" ~ "Hitzetage",
    key_id == "ID" ~ "Eistage",
    key_id == "SD" ~ "Sommertage",
    key_id == "snowdays" ~ "Schneetage",
    key_id == "TN" ~ "Tropennächte",
    key_id == "tas" ~ "Temperatur",
    key_id == "pr" ~ "Niederschlagsmenge"
  )
}

# display all municipalities in a hisogram
HistogramPlot <- function(data, key_id, period_val, rcp_val){
  colors <- c("#1cb0b5", "#e31f2b")
  names(colors) <- c("1981-2010", period_val)
  data %>%
    filter(key == key_id &
           period %in% c("1981-2010", period_val) &
           rcp %in% c("obs", rcp_val) &
           estimate %in% c("yearly_mean", "q50")) %>%
    group_by(period) %>%
    ggplot(aes(value, fill = period)) +
    geom_histogram(alpha = 0.4, position = "identity") +
    scale_fill_manual(values = colors) +
    theme_minimal() +
    labs(title = glue(
      "Wie hat sich die Anzahl 
      {translate_key_to_word(key_id)} im Vergleich 
      zum Mittelwert von 1981-2010 verändert?"
      ),
      subtitle = glue("Prognose für {rcp_val}"),
      x = glue(
          "Durchschnittliche Anzahl 
          {translate_key_to_word(key_id)} pro Jahr"
          ),
      y = "Anzahl Gemeinden",
      fill = "")
  }

# display all municipalities in a violin plot
ViolinPlot <- function(data, key_id, var, grouping_var){
  data %>%
    filter(key == key_id &
           (rcp == var | period == var) &
           estimate == "q50") %>%
    group_by(!!sym(grouping_var)) %>%
    ggplot(aes(x = !!sym(grouping_var), y = value)) +
    geom_violin() +
    theme_minimal() +
    labs(title = glue("Szenarien für {translate_key_to_word(key_id)}"),
        subtitle = ifelse(
          var %in% c("RCP2.6", "RCP4.5", "RCP8.5"), 
          glue("Prognose des {var}"), 
          glue("Klimaszenarien im Jahr {var}")
          ),
        x = "",
        y = glue(
          "Durchschnittliche Anzahl 
          {translate_key_to_word(key_id)} pro Jahr"
          ),
        fill = "")
  }
 
# find outlier municipalities
ExtremeMunicipalities <- function(data, key_id, period_val, rcp_val){
  data %>%
    filter(key == key_id &
           period %in% c("1981-2010", period_val) &
           rcp %in% c("obs", rcp_val) &
           estimate %in% c("yearly_mean", "q50")) %>%
    select(-rcp, -estimate, -season) %>%
    mutate(period = case_when(
      period == "1981-2010" ~ "now",
      period == period_val ~ "future"
    )) %>% 
    spread(key = period, value = value) %>% 
    mutate(
      delta = future - now
    ) %>%
    arrange(delta) %>%
    ungroup() %>%
    # Get first and last 10
    slice(1:10, (n() - 9):n())
}

# mean per season Plot
Seasons <- function(data, key_id, period_val, rcp_val){
  colors <- c("#1cb0b5", "#e31f2b")
  names(colors) <- c("1981-2010", period_val)
  
  data %>%
    filter(key == key_id &
           period %in% c("1981-2010", period_val) &
           rcp %in% c("obs", rcp_val) &
           estimate %in% c("mean", "q50")) %>%
    mutate(season = case_when(
        season == "MAM" ~ "Mär, Apr, Mai",
        season == "JJA" ~ "Jun, Jul, Aug",
        season == "SON" ~ "Sep, Okt, Nov",
        season == "DJF" ~ "Dez, Jan, Feb"
    ),
    season = factor(season, levels = c(
        "Mär, Apr, Mai",
        "Jun, Jul, Aug",
        "Sep, Okt, Nov",
        "Dez, Jan, Feb"))
    ) %>%
    group_by(period, season) %>%
    ggplot(aes(value, fill = period)) +
    geom_histogram(alpha = 0.4, position = "identity") +
    scale_fill_manual(values = colors) +
    theme_minimal() +
    facet_wrap(~ season, nrow = 1) +
    labs(title = glue(
      "Wie hat sich der saisonale Tagesmittelwerte der 
      {translate_key_to_word(key_id)} im
      Vergleich zum Mittelwert von 1981-2010 verändert?"
      ),
      subtitle = glue("Prognose des {rcp_val}"),
      x = glue("Tagesmittelwert 
               {translate_key_to_word(key_id)} 
               pro Jahreszeit"),
      y = "Anzahl Gemeinden",
      fill = "")
  }

# Plot which regions are sensitive to possible outcomes
Regions <- function(data, key_id, rcp_val){
  data %>%
    filter(key == key_id &
            rcp %in% c("obs", rcp_val) &
            estimate %in% c("yearly_mean", "q50")
            ) %>%
    ungroup() %>%
    mutate(
      region = factor(region, levels = c(
        "Jura",
        "Mittelland",
        "Voralpen",
        "Alpen",
        "AlpenS"
      )),
      altitude = factor(altitude, levels = c(
        "0-800",
        "801-1500",
        ">1500"))
        ) %>%
    group_by(period, region, altitude, rcp) %>%
    summarise(mean = mean(value)) %>%
    ggplot(aes(x = period,
               y = mean,
               fill = period)) +
    geom_bar(stat = "identity", 
               position = position_dodge2(reverse = TRUE)) +
    facet_grid(cols = vars(region), 
               rows = vars(altitude)) +
    scale_x_discrete(labels = NULL) +
    scale_fill_manual(values = c(
      "1981-2010" = "#ffd651",
      "2035" = "#f7a600",
      "2060" = "#ed7004",
      "2085" = "#ad3e14"
      )) +
    theme_minimal() +
    labs(title = glue(
      "Wie hat sich die durchschittliche Anzahl 
      {translate_key_to_word(key_id)} im
      Vergleich zum Mittelwert von 1981-2010 verändert?"
      ),
      subtitle = glue("Prognose des {rcp_val} nach Grossregion"),
      x = "",
      y = glue("Durchschnittliche Anzahl {translate_key_to_word(key_id)}"),
      fill = "")
  }

# Region and Season
translate_season <- function(season_val){
  case_when(
    season_val == "MAM" ~ "März, April, Mai",
    season_val == "JJA" ~ "Juni, Juli, August",
    season_val == "SON" ~ "September, Oktober, November",
    season_val == "DJF" ~ "Dezember, Januar, Februar")
}

RegionsSeason <- function(data, key_id, rcp_val, season_val){
  data %>%
    filter(key == key_id &
            rcp %in% c("obs", rcp_val) &
            estimate %in% c("mean", "q50") &
            season == season_val
            ) %>%
    ungroup() %>%
    mutate(
      region = factor(region, levels = c(
        "Jura",
        "Mittelland",
        "Voralpen",
        "Alpen",
        "AlpenS"
      )),
      altitude = factor(altitude, levels = c(
        "0-800",
        "801-1500",
        ">1500"))
        ) %>%
    group_by(period, region, altitude, rcp) %>%
    summarise(mean = mean(value)) %>%
    ggplot(aes(x = period,
               y = mean,
               fill = period)) +
    geom_bar(stat = "identity", 
               position = position_dodge2(reverse = TRUE)) +
    facet_grid(cols = vars(region), 
               rows = vars(altitude)) +
    scale_x_discrete(labels = NULL) +
    scale_fill_manual(values = c(
      "1981-2010" = "#ffd651",
      "2035" = "#f7a600",
      "2060" = "#ed7004",
      "2085" = "#ad3e14"
      )) +
    theme_minimal() +
    labs(title = glue(
      "Wie hat sich der Mittelwert der 
      {translate_key_to_word(key_id)} im Vergleich zum 
      Mittelwert von 1981-2010 verändert ({translate_season(season_val)})?"),
      subtitle = glue("Prognose des {rcp_val} nach Grossregion"),
      x = "",
      y = glue("Durchschnittliche Anzahl {translate_key_to_word(key_id)}"),
      fill = "")
  }

LineChartPerMunicipality <- function(data, municipality, key_id) {
  data %>%
    filter(
      GMDNAME == municipality,
      key == key_id,
      estimate %in% c("yearly_mean", "q50")) %>%
    ggplot(aes(x = period, y = value, group = rcp, color = rcp)) +
    geom_line() +
    geom_point() +
    theme_minimal() +
    labs(title = glue(
      "Wie hat sich der Mittelwert der 
      {translate_key_to_word(key_id)} im Vergleich zum 
      Mittelwert von 1981-2010 verändert?"
      ),
      subtitle = glue("Prognose für {municipality}"),
      x = "",
      y = glue("Durchschnittliche Anzahl {translate_key_to_word(key_id)}"),
      fill = "")
  }

# Small Multiple Maps
SmallMultipleMaps <- function(data, key_id, rcp_val, estimate_val) {
  municipality_geo %>%
  left_join(
    data %>%
      filter(
        key == key_id,
        rcp == rcp_val,
        estimate %in% c("yearly_mean", estimate_val)
      ),
    by = c("GMDNR" = "bfs_id")
  ) %>%
  ggplot(aes(
    fill = value
  )) +
    geom_sf(color = "transparent") +
    facet_wrap(~ period, ncol = 2) +
    theme(
      axis.line = element_blank(),
      axis.text.x = element_blank(),
      axis.text.y = element_blank(),
      axis.ticks = element_blank(),
      axis.title.x = element_blank(),
      axis.title.y = element_blank(),
      panel.background = element_blank(),
      panel.border = element_blank()
      ) +
    labs(
      title = glue(
        "{translate_key_to_word(key_id)}, {rcp_val}, {estimate_val}"
        ),
      fill = translate_key_to_word(key_id)
      ) +
    scale_fill_viridis_c()
  }
LineChartPerMunicipality(climate_projections_by_municipality, "Zürich", "FD")