2

A cold wave is defined if the minimum temperature at a meteorological station is below the normal temperature by 3 °C or more consecutively for 3 days or more. I have tried to calculate it multiple stations using following code

#Calculation for multistation
set.seed(123)
df <- data.frame("date"= seq(from = as.Date("1970-1-1"), to = as.Date("2000-12-31"), by = "day"),
                 "Station1" = runif(length(seq.Date(as.Date("1970-1-1"), as.Date("2000-12-31"), "days")), 10, 30),
                 "Station2" = runif(length(seq.Date(as.Date("1970-1-1"), as.Date("2000-12-31"), "days")), 11, 29),
                 "Station3" = runif(length(seq.Date(as.Date("1970-1-1"), as.Date("2000-12-31"), "days")), 9, 28))

head(df)

df$day <- format(df$date, format='%m-%d')

#Daily average (daily normal) calculation
df_summarise_all <- df %>% 
  as_tibble() %>% # for easier viewing 
  mutate(day = format(df$date, format='%m-%d')) %>% 
  group_by(day) %>%
  summarise_all(funs(mean)) %>% 
  pivot_longer(cols = -c(date, day), names_to = "variable", values_to = "value") 

#Coldwave event calculation
df %>%
  as_tibble() %>% # for easier viewing 
  pivot_longer(cols = -c(date, day), names_to = "Stations", values_to = "MinT") %>% 
  left_join(df_summarise_all %>% rename(mean_MinT = value), by = 'day') %>% 
  mutate(is_coldwave = zoo::rollapplyr(MinT < (mean_MinT - 3), 
                                       3, all,fill = NA))

As can be seen from the output, the joining of the daily normal and station minimum temperature is not correct. I have three questions

  1. How to correct the join?
  2. How to get the count of coldwaves per year for every station? and
  3. How to get overall coldwave count for every station?
UseR10085
  • 7,120
  • 3
  • 24
  • 54
  • By "normal temperature" are you referring to the `mean` temperature until that point in time? – Dunois Apr 13 '20 at 11:13
  • Normal temperature means long term average (generally 30 years or more). I have calculated that using the daily normal chunk of code. – UseR10085 Apr 13 '20 at 11:14
  • So that'd be the `mean` temperature for the entire dataset in this case? – Dunois Apr 13 '20 at 11:16
  • No, the daily normal has been calculated as for 1st row or value is mean of 1 Jan of all the years. Second row is the mean of 2nd Jan of all the years and so on. This way you will get 366 daily normal values. [See this link](https://stackoverflow.com/q/48386324/6123824) – UseR10085 Apr 13 '20 at 11:21

1 Answers1

1

You should join also by Stations and probably compare with mean_MinT - 3 for coldwave. 

library(dplyr)

df_out <- df %>%
             tidyr::pivot_longer(cols = -c(date, day), 
                 names_to = "Stations", values_to = "MinT") %>%
             left_join(df_summarise_all %>% rename(mean_MinT = value), 
                         by = c('day' = 'day', 'Stations' = 'variable')) %>%
             mutate(is_coldwave = zoo::rollapplyr(MinT < (mean_MinT - 3), 3, 
                     all,fill = NA))

To calculate yearly and total coldwave you can use the same method as in this answer. 

df_year <- df_out %>%
            group_by(Stations, year = format(date.x, "%Y")) %>%
           summarise(total_cold = with(rle(is_coldwave), 
                                        sum(values, na.rm = TRUE)))

and sum(df_year$total_cold) would give overall count.

Ronak Shah
  • 377,200
  • 20
  • 156
  • 213
  • Yea, it should be `mean_MinT - 3`. I have edited my question. It will be of great help if you answer my 2nd and 3rd questions also. – UseR10085 Apr 13 '20 at 11:29
  • @BappaDas Isn't that same as the previous answer I added link to? You have to use the exact same code. `df_year <- df_out %>% group_by(year = format(date, "%Y")) %>% summarise(total_heat = with(rle(is_heatwave), sum(values, na.rm = TRUE)))` – Ronak Shah Apr 13 '20 at 12:57
  • It is giving me ```Error: Column `year` must be length 33969 (the number of rows) or one, not 2```. – UseR10085 Apr 13 '20 at 13:11
  • Solved it, it should be `date.x`. It would be better if you add it here also. – UseR10085 Apr 13 '20 at 13:17
  • The problem with the updated answer is that it gives sum irrespective of stations. But I want stationwise sum of coldwaves. – UseR10085 Apr 13 '20 at 13:30
  • You need to add the column in `group_by`. – Ronak Shah Apr 13 '20 at 13:31
  • I have done it using ```df_out %>% group_by(year = format(date.x, "%Y"), Stations) %>% summarise(total_cold = with(rle(is_coldwave), sum(values, na.rm = TRUE))) %>% group_by(Stations) %>% summarise(Coldwave=sum(total_cold, na.rm = TRUE)) ```. Thank you very much for your help. – UseR10085 Apr 13 '20 at 13:33