Welcome to the exercises of module 3! Please try to solve the following exercises, and if you do not remember how to write the code for certain things, have a look again at the previous scripts. Programming is about learning by doing! So, do not hesitate to look for help in the internet, too. Remember that there is almost always more than just one way to solve a problem in R!

Part A: Weather station measurements

Exercise 01 (10 points)

Import the data of the weather station assigned to your group. Write a function that calculates the extreme weather variable assigned to your group, i.e. the function should create two new columns in your station record: One column that identifies whether an extreme event happened on the respective day, and another column that indicates the difference between the threshold and the measured value of the respective weather parameter (e.g. for group A, the weather parameter is PRCP). Apply your function to your station record, with the threshold assigned to your group (e.g. for group A, the threshold is 15 mm), and save the output from this function as a CSV file.

library(dplyr)
library(Kendall)
library(ggplot2)
library(stringr)
library(raster)
library(rgdal)

dir <- "C:/Users/Max/Desktop/IAMO_neu/eLearning/Module_3_Extreme_weather/solutions/"
areni_raw <- read.csv(paste0(dir, "data/Areni_edit.csv"), sep=",")
areni <- filter(areni_raw, PARA_SHORT == "PRCP")
areni <- dplyr::select(areni, c("YEAR", "MONTH", "DAY", "DATE", "DOY", "VALUE"))
names(areni)[6] <- "PRCP"

HeavyPrecipEvents <- function(table, threshold) {
  
  heavyprecipdays <- which(table$PRCP > threshold)
  
  newcolumnA  <- paste0("above", threshold)
  table[[newcolumnA]] <- 0
  table[[newcolumnA]][heavyprecipdays] <- 1
  
  newcolumnB  <- paste0("above", threshold, "_diff")
  table[[newcolumnB]] <- 0
  table[[newcolumnB]][heavyprecipdays] <- table$PRCP[heavyprecipdays] - threshold
  
  return(table)
}

result <- HeavyPrecipEvents(areni, 15)

head(result)
YEAR MONTH DAY DATE DOY PRCP above15 above15_diff
2002 9 1 2002-09-01 244 0 0 0
2002 9 2 2002-09-02 245 0 0 0
2002 9 3 2002-09-03 246 0 0 0
2002 9 4 2002-09-04 247 0 0 0
2002 9 5 2002-09-05 248 0 0 0
2002 9 6 2002-09-06 249 0 0 0 
write.csv(result, paste0(dir, "ex01.csv"), row.names = F)

Exercise 02 (10 points)

Calculate the 10th (groups C) or 90th (all other groups) percentile of the weather parameter that is used to calculate the extreme weather variable assigned to your group, for each month from January to December, for the entire period covered by your weather station. Save the results in a data frame and export it as a CSV file.

monthly_thresholds <- data.frame(month = c(1:12), threshold90 = NA)

for (i in 1:12)
{ areni_month <- areni$PRCP[which(areni$MONTH == i)]
  monthly_thresholds$threshold90[i] <- quantile(areni_month, probs = 0.9) }

monthly_thresholds
month threshold90
1 4.00
2 3.50
3 4.90
4 8.21
5 5.36
6 2.80
7 1.16
8 0.00
9 0.00
10 4.10
11 2.81
12 2.90 
write.csv(monthly_thresholds, paste0(dir, "ex02.csv"), row.names = F)

Exercise 03 (10 points)

Integrate the monthly percentile thresholds into your station record and re-assess the occurrence of extreme weather events (i.e. whether an event happened or not) based on these thresholds. Plot the number of occurrences of these events in a facetted barplot - “year” should be on the x-axis, each facet should represent one month, and each bar should represent the amount of days with extreme events during the respective month and year. Export your plot as PNG.

areni$threshold90 <- monthly_thresholds$threshold90[areni$MONTH]
areni$event90      <- areni$PRCP - areni$threshold90
areni$event90[areni$event90 <= 0] <- 0
areni$event90[areni$event90 > 0]  <- 1

areni_summary <- areni %>% group_by(YEAR, MONTH) %>% summarize (heavyprecipdays = sum(event90)) %>% as.data.frame()
## `summarise()` has grouped output by 'YEAR'. You can override using the `.groups` argument.
plot_ex3 <- ggplot(data=areni_summary, aes(x=YEAR, y=heavyprecipdays)) +
    geom_bar(stat="identity") +
    labs(x="", y="") +
    theme(axis.text.x=element_text(angle=45, hjust=1)) +
    ggtitle("Heavy Precipitation days at Areni station during different months, by year") +
    facet_wrap(~MONTH)

plot_ex3

png(paste0(dir, "ex03.png"), height=400, width=600)
  plot_ex3
dev.off()
## png 
##   2

Exercise 04 (10 points)

Use the result from exercise 03 to perform, for each month from January to December, the Mann-Kendall test to assess whether there is a significant change in extreme events over time. Save the results of this analysis as a CSV file.

mannkendall_results <- data.frame(month=c(1:12), change=NA)

for (i in 1:12)
{ areni_summary_temp <- filter(areni_summary, MONTH == i)
  teststat <- MannKendall(areni_summary_temp$heavyprecipdays)  
  if (teststat[2]$sl[1] < 0.05) { mannkendall_results$change[i] <- "significant" 
  } else { mannkendall_results$change[i] <- "not significant" }
}

mannkendall_results
month change
1 not significant
2 not significant
3 not significant
4 not significant
5 not significant
6 not significant
7 not significant
8 not significant
9 not significant
10 not significant
11 not significant
12 not significant 
write.csv(mannkendall_results, paste0(dir, "ex04.csv"), row.names = F)

Exercise 05 (10 points)

Write a function that calculates the occurrence of the extreme weather event assigned to your group for any given vector of percentile values (e.g. 10th, 11th and 12th percentile, etc.). The function should add as many columns to the station record as there are percentile values supplied to the function. Apply the function to your station record using five different thresholds, and save the resulting table as a CSV file.

extremePercentiles <- function(data, percentiles) {

  monthly_thresholds <- NULL

  for (i in 1:12)
  { data_month <- data$PRCP[which(data$MONTH == i)]
    monthly_thresholds <- rbind(monthly_thresholds, quantile(data_month, probs = percentiles)) 
  }

  monthly_thresholds <- as.data.frame(monthly_thresholds)
  monthly_thresholds$MONTH <- c(1:12)
  names(monthly_thresholds)[c(1:length(percentiles))] <- paste0("threshold_", names(monthly_thresholds)[c(1:length(percentiles))])

  data_merged <- merge(data, monthly_thresholds, by="MONTH")

  n_perc <- length(percentiles)

  for (n in 1:n_perc)
  { data_merged[,n+6+n_perc]                              <- data_merged$PRCP - data_merged[,n+6] 
    data_merged[,n+6+n_perc][data_merged[,n+6+n_perc]<=0] <- 0
    data_merged[,n+6+n_perc][data_merged[,n+6+n_perc]>0]  <- 1
  }

  names(data_merged)[(6+n_perc+1):length(data_merged)] <- paste0("events_", percentiles)

  return(data_merged)
}

areni <- areni[,c(1:6)]
result2 <- extremePercentiles(areni, seq(from=0.9, to=0.94, by=0.01))

head(result2)
MONTH YEAR DAY DATE DOY PRCP threshold_90% threshold_91% threshold_92% threshold_93% threshold_94% events_0.9 events_0.91 events_0.92 events_0.93 events_0.94
1 2003 29 2003-01-29 29 0.0 4 4.7 5 5.601 6.858 0 0 0 0 0
1 2005 22 2005-01-22 22 0.0 4 4.7 5 5.601 6.858 0 0 0 0 0
1 2006 19 2006-01-19 19 0.0 4 4.7 5 5.601 6.858 0 0 0 0 0
1 2003 30 2003-01-30 30 0.7 4 4.7 5 5.601 6.858 0 0 0 0 0
1 2003 31 2003-01-31 31 0.0 4 4.7 5 5.601 6.858 0 0 0 0 0
1 2004 27 2004-01-27 27 0.0 4 4.7 5 5.601 6.858 0 0 0 0 0 
write.csv(result2, paste0(dir, "ex05.csv"), row.names = F)

Part B: Future Climate Projections

Exercise 06 (10 points)

Import as RasterStacks the two NetCDFs that correspond to the climate model and future time period assigned to your group. There is one file for RCP 4.5, and one for RCP 8.5. Convert both stacks from degrees Kelvin to degrees Celsius. Display the first layer of each stack in a separate map and add the outline of Armenia. Add a title to both plots and export them as separate PNG files.

data_rcp45 <- stack(paste0(dir, "data/tasmax_day_MIROC5_rcp45_r1i1p1_EWEMBI_20810101-20901231_lat38.0to45.0lon39.0to52.0.nc4"))
## Loading required namespace: ncdf4
data_rcp85 <- stack(paste0(dir, "data/tasmax_day_MIROC5_rcp85_r1i1p1_EWEMBI_20810101-20901231_lat38.0to45.0lon39.0to52.0.nc4"))

data_rcp45 <- data_rcp45 - 273.15
data_rcp85 <- data_rcp85 - 273.15

armenia <- readOGR(paste0(dir, "shapefiles/gadm41_ARM_0.shp"), verbose=F)

png(paste0(dir, "ex06_A.png"), height=400, width=600)
  plot(data_rcp45[[1]], main="Max. Temp. from MIROC5 on Jan 01 2081, RCP 4.5")
  plot(armenia, add=TRUE)
dev.off()
## png 
##   2
  plot(data_rcp45[[1]], main="Max. Temp. from MIROC5 on Jan 01 2081, RCP 4.5")
  plot(armenia, add=TRUE)

png(paste0(dir, "ex06_B.png"), height=400, width=600)
  plot(data_rcp85[[1]], main="Max. Temp. from MIROC5 on Jan 01 2081, RCP 8.5")
  plot(armenia, add=TRUE)
dev.off()
## png 
##   2

Exercise 07 (10 points)

Calculate the amount of heat events in Armenia in June 2055 (groups A, C, E and G) or in June 2085 (groups B, D, and F) under RCP 4.5. The occurrence of a heat event should be defined as the case when a cell that overlaps Armenia exceeds 30 degrees Celsius.

## identify layers that correspond to June 2085
index <- which( (substring(names(data_rcp45), 2, 5) == "2085") & (substring(names(data_rcp45), 7, 8) == "06") )

## select only those layers that correspond to June 2085
data_rcp45_june2085 <- data_rcp45[[index]]

## change raster values below the threshold to 0, and above the threshold to 1
data_rcp45_june2085[data_rcp45_june2085 < 30]  <- 0
data_rcp45_june2085[data_rcp45_june2085 >= 30] <- 1

## rasterize the shapefile "armenia" based on the resolution of "data_rcp45_june2085" to create a template that yields a proper result in the *mask()* function.
armenia_rastered <- rasterize(armenia, data_rcp45_june2085, getCover=TRUE)
armenia_rastered[armenia_rastered>0]  <- 1
armenia_rastered[armenia_rastered==0] <- NA

## use "armenia_rastered" as template when applying the *mask()* function to the raster stack "data_rcp45_june2085"  
data_rcp45_june2085_masked <- mask(data_rcp45_june2085, armenia_rastered)

## count number of events in data_rcp45_june2085_masked
sum(na.omit(values(data_rcp45_june2085_masked)))
## [1] 259

Exercise 08 (10 points)

Write a function that calculates the yearly heat events (threshold: 30 degrees Celsius) in Armenia during summer, i.e. between June 21 and September 23. With this function, calculate the amount of heat events for each year in your RCP 8.5 data set, and export the results as a CSV file.

SummerHeatArmenia <- function(data, year, threshold) {
  
  index <- which( (substring(names(data), 2, 5) == year) & ( ( (substring(names(data), 7, 8) == "06") & (substring(names(data), 10, 11) %in% as.character(c(21:30))) ) | 
                    (substring(names(data), 7, 8) %in% c("07", "08")) | 
                    ( (substring(names(data), 7, 8) == "09") & (substring(names(data), 10, 11) %in% str_pad(c(1:23), 2, pad="0")) ) ))
  
  data_sel <- data[[index]]
  
  data_sel[data_sel < 30]  <- 0
  data_sel[data_sel >= 30] <- 1
  
  data_sel_masked <- mask(data_sel, armenia_rastered)

  return(sum(na.omit(values(data_sel_masked))))
}

summerheat <- data.frame(year=as.character(c(2081:2090)), events=NA)
for (i in 1:NROW(summerheat))
{ summerheat$events[i] <- SummerHeatArmenia(data_rcp85, summerheat$year[i], 30) }

head(summerheat)
year events
2081 1700
2082 1682
2083 1557
2084 1933
2085 2074
2086 1665 
write.csv(summerheat, paste0(dir, "ex08.csv"), row.names=F)

Exercise 09 (10 points)

Write a function that calculates low temperature events during any month and year, i.e. the amount of cells with maximum daily temperatures below a certain threshold. For each year covered by your RCP 4.5 data set, calculate the amount of low temperature events for two different months and three different thresholds. The result should be a table with 6 entries. Export this table as a CSV file.

LowTempEvents <- function(data, year, month, threshold) {
  
   index <- which((substring(names(data), 2, 5) == year) & (substring(names(data), 7, 8) == month))
   
   data_sel <- data[[index]]
   
   data_sel[data_sel < threshold]  <- 1
   data_sel[data_sel >= threshold] <- 0

   return(sum(na.omit(values(data_sel))))
}

results3 <- data.frame(year=rep(as.character(c(2081:2090)),each=6), month=as.character(str_pad(rep(c(1:2)), 2, pad="0")), threshold=rep(c(5,10,15), each=2), events = NA)

for (i in 1:NROW(results3))
{ results3$events[i] <- LowTempEvents(data_rcp85, results3$year[i], results3$month[i], results3$threshold[i]) }

head(results3)
year month threshold events
2081 01 5 3823
2081 02 5 2295
2081 01 10 8245
2081 02 10 5855
2081 01 15 10762
2081 02 15 9016 
write.csv(results3, paste0(dir, "ex09.csv"), row.names=F)

Exercise 10 (10 points)

Calculate the mean and standard deviation of the ten yearly amounts of extreme heat events during the whole year (threshold: 35 degrees Celsius) for both RCPs. Compare the results in a bar plot. The bar plot should contain one bar for each RCP, and each bar should have whiskers that show the standard deviation across the respective 10 yearly values. Export your plot as PNG.

extrheat_results <- data.frame(rcp = rep(c("RCP 4.5", "RCP 8.5"), each=10), year = c(2081:2090), events = NA)

ExtremeHeat <- function(data, year, threshold) {

   index <- which(substring(names(data), 2, 5) == as.character(year))
  
   data_sel <- data[[index]]
   
   data_sel[data_sel < threshold]  <- 0
   data_sel[data_sel >= threshold] <- 1

   return(sum(na.omit(values(data_sel))))
}

for ( i in 1:NROW(extrheat_results)) 
{ if (extrheat_results$rcp[i] == "RCP 4.5") { dataset <- data_rcp45 
  } else if (extrheat_results$rcp[i] == "RCP 8.5") { dataset <- data_rcp85  }
  
  extrheat_results$events[i] <- ExtremeHeat(dataset, extrheat_results$year[i], 35)
}

plotdata <- extrheat_results %>% group_by(rcp) %>% summarize (heat_mean = mean(events), heat_sd = sd(events)) %>% as.data.frame()

plot_ex10 <- ggplot(data=plotdata, aes(x=rcp, y=heat_mean)) +
    geom_bar(stat="identity") +
    geom_errorbar(aes(ymin=heat_mean-heat_sd, ymax=heat_mean+heat_sd)) +
    labs(y="raster cells exceeding 35 °C", x="") +
    ggtitle("Mean yearly heat events in the Caucasus according to MIROC5 from 2081 to 2090")

plot_ex10

png(paste0(dir, "ex010.png"), height=400, width=600)
  plot_ex10
dev.off()
## png 
##   2