neue Version mit leichter Veränderung

climate_indices/SDII.R

@@ -1,7 +1,8 @@
 #-------------------------------------------------------------------------------
 #Topic 4: Big data - DWD Stations with 40 years "Climate Data Baden-Württemberg"
 #
-#Indicator SDII: annual number of rain days
+# Simple daily intensity index: 
+# total amount of rainfall per year devided by the number of rain days per year
 #
 #M. Lorff, J. Opdenhoff, N. Veigel, Januar 2018
 #-------------------------------------------------------------------------------
@@ -16,10 +17,9 @@ library("tidyverse")
 #-------------------------------------------------------------------------------
 #Data
 #-------------------------------------------------------------------------------
-# The function requires the file "BW_Climate_1977_2016.txt" as downloaded from
-# the moodle course. Fill in the path to the original, unedited file.
+# The function requires the data from file "BW_Climate_1977_2016.txt" as downloaded from
+# the moodle course. Use e.g. read_tsv on the original, unedited file.
 
-path <- "... enter your location here.../BW_Climate_1977_2016.txt"
 
 # You need to specify your station id to run the function.
 
@@ -29,31 +29,31 @@ ID <-  2638
 #Function
 #-------------------------------------------------------------------------------
 
-SDII <-  function(path, ID) {
+SDII <-  function(data, ID) {
 
-data <- read_tsv(path)
 
 data <- 
     data %>%
   #extract years
     mutate(year = year(date)) %>%
     group_by(id, year) %>%
-  #count annual NA-values and days with more than 1 mm rainfall
+  #count annual NA-values, days with more than 1 mm rainfall and sum up the total rainfall per year
     summarise(count_na = sum(is.na(RSK)), count_precip = sum(RSK>=1, na.rm = T), total = sum(RSK, na.rm = T)) %>%
   #extract only your station ID
     filter(id == ID) %>%
   #Exclude years that have more than 60 NA-values
     mutate(rain_days = ifelse(count_na >= 60, NA, count_precip)) %>%
-    mutate(value = total/rain_days * 10) %>%
+    mutate(value_SDII = total/rain_days * 10) %>%
     ungroup()
 
-    mw <- mean(data$value, na.rm = T) #Index 1 average rain days over the whole period
+    mw <- mean(data$value_SDII, na.rm = T) #Index 1 average rain days over the whole period
 
-    avg1 <-  mean(data$value[which(data$year == 1977):which(data$year == 1996)], na.rm = T)
-    avg2 <-  mean(data$value[which(data$year == 1997):which(data$year == 2016)], na.rm = T)
+    #calculate mean of 20-year-serieses, BUT ONLY if NAs are less then 10
+    avg1 <-  ifelse(length(na.omit(data$value_SDII)[which(data$year == 1977):which(data$year == 1996)]) < 10, NA, mean(data$value_SDII[which(data$year == 1977):which(data$year == 1996)], na.rm = T))
+    avg2 <-  ifelse(length(na.omit(data$value_SDII)[which(data$year == 1997):which(data$year == 2016)]) < 10, NA, mean(data$value_SDII[which(data$year == 1997):which(data$year == 2016)], na.rm = T))   
 
-  #Exclude time series that has less than 10 years of data and calculate Index 2
-    diff <- ifelse(length(na.omit(data$value)) < 10, NA, avg2/avg1 * 100)
+    #calculate Index 2
+    diff <- (avg2/avg1) * 100
 
 return(list(rate_change =diff, total_mean = mw, data = data))
 
@@ -67,4 +67,4 @@ return(list(rate_change =diff, total_mean = mw, data = data))
 #Index 2 (rate_change) [%] and a data frame with the number of rain days for every
 #year at your station
 
-result <- SDII(path,ID)
+result <- SDII(raw_data,ID)