Sydney and NSW are having a lot of rain this year, thought to be due to a La Nina weather event. But this raise a few questions:

To answer this question, I initially looked for river level data in NSW, but river levels over time can be affected by development of the river basin. A cleaner set of data is rainfall data, and on a BOM website I found a set of daily rainfall data for the Sydney Botanic Gardens (next to Circular Quay) that has continuous data back to 1910. The link to the data source is shown below, and the ‘station number’ is 066006.

Data source - http://www.bom.gov.au/climate/data/index.shtml

This data was downloaded from the BOM website on 15 August 2022, and had data up to 1st August 2022. I imported this data into the RStudio data science package. Altogether there were 50,250 observations of daily rainfall in these gardens. This included 201 missing values - I replaced these missing values with the rainfall data for the previous day.

setwd("G:/My Drive/Training and development/R sandpit/Sydney rainfall analysis")
# read file
rainfall <- read.csv(file = 'Daily_rainfall_data_15Oct2022.csv')  # import rainfall data
rainfall %<>% rename(month_num = Month,
                     rainfall_mm = 6)  # rename 2 columns
rainfall %<>% mutate(Month = case_when(
            month_num == 1 ~ 'January',
            month_num == 2 ~ 'February',
            month_num == 3 ~ 'March',
            month_num == 4 ~ 'April',
            month_num == 5 ~ 'May',
            month_num == 6 ~ 'June',
            month_num == 7 ~ 'July',
            month_num == 8 ~ 'August',
            month_num == 9 ~ 'September',
            month_num == 10 ~ 'October',
            month_num == 11 ~ 'November',
            month_num == 12 ~ 'December',
            TRUE ~ 'error'),
            year_month = paste0(Year,'-',Month))

rainfall %<>% filter(Year > 1909) # limit to data after 2009 - which is the continuous data
rainfall$Month <- factor(rainfall$Month, levels = c('January', 'February', 'March', 'April', 'May', 'June', 'July', 'August', 'September', 'October', 'November','December'), ordered = TRUE) # create ordered factor variable for month
rainfall$rainfall_mm <- nafill(rainfall$rainfall_mm, type = 'locf')    # replace missing values with the rainfall for the earlier day
rainfall %<>% mutate(wet_day =as.integer(if_else(rainfall_mm > 0,1,0)))

rainfall %<>% group_by(year_month) %>% mutate(wet_day_count = cumsum(wet_day)) %>% 
              mutate(wet_day_sum = max(wet_day_count))

#head(rainfall)

Month by month variation

Looking at the data, the first thing I noted is that rainfall is, not surprisingly, affected by the month of the year, as shown below. In this diagram, the error bars show the 95% confidence intervals.

alpha <- 0.05
monthly_mean <- rainfall %>% 
    group_by(Month) %>% 
    summarize(mean = mean(rainfall_mm),
              lower_rain = mean(rainfall_mm) - qt(1- alpha/2, (n() - 1))*sd(rainfall_mm)/sqrt(n()),
              upper_rain = mean(rainfall_mm) + qt(1- alpha/2, (n() - 1))*sd(rainfall_mm)/sqrt(n()),
              wet_days = mean(wet_day_sum),
            lower_days = mean(wet_day_sum) - qt(1- alpha/2, (n() - 1))*sd(wet_day_sum)/sqrt(n()),
             upper_days = mean(wet_day_sum) + qt(1- alpha/2, (n() - 1))*sd(wet_day_sum)/sqrt(n()))
monthly_rainfall_mean <- mean(monthly_mean$mean)
monthly_mean$rainfall_correction = monthly_rainfall_mean/monthly_mean$mean # determine correction factor
monthly_rainfall_count <- mean(monthly_mean$wet_days)
monthly_mean$wet_day_correction = monthly_rainfall_count/monthly_mean$wet_days
month_mean2 <- monthly_mean %>% select(Month, rainfall_correction,wet_day_correction)
monthly_mean %>%
  ggplot(aes(x = Month, y = mean)) +
  geom_bar(stat = "identity", fill = "steelblue") +
  geom_errorbar(aes(ymin = lower_rain, ymax = upper_rain), width=.2,
                 position=position_dodge(.9)) + 
  labs(title = "Fig. 1 - Average daily rainfall in Sydney Botanic Gardens peaks \nin March and June",
              caption = "Data source: BOM",
              x = "Month", y = "Average daily \nrainfall in mm") +
  theme(axis.title.y=element_text(angle=0, vjust = 0.5, size = 10)) +
  scale_x_discrete(breaks = function(x){x[c(TRUE, FALSE)]})

We can also look at the number of wet days:

monthly_mean %>%
  ggplot(aes(x = Month, y = wet_days)) +
  geom_bar(stat = "identity", fill = "darkolivegreen3") +
  geom_errorbar(aes(ymin = lower_days, ymax = upper_days), width=.2,
                 position=position_dodge(.9)) + 
  labs(title = "Fig. 2 - Wet days in Sydney Botanic Gardens peaks in March",
              caption = "Data source: BOM",
              x = "Month", y = "Average number \nof wet days \nin month") +
  theme(axis.title.y=element_text(angle=0, vjust = 0.5, size = 10)) +
  scale_x_discrete(breaks = function(x){x[c(TRUE, FALSE)]})

Has 2022 had more wet days?

We can also apply the T-test test to the number of wet days in a month. The results shown below shows lower p-value of 0.0113, which shows lower probability that this variation is due to chance alone.

Monthly_average_wet_days_before_2022 <- monthly_sum$corrected_mean_wet_days[monthly_sum$Year < 2022]
Monthly_average_wet_days_2022 <- monthly_sum$corrected_mean_wet_days[monthly_sum$Year == 2022]

#shapiro.test(Monthly_average_rainfall_2022)
#shapiro.test(Monthly_average_rainfall_before_2022)
#ggqqplot(Monthly_average_rainfall_2022)
#ggqqplot(Monthly_average_rainfall_before_2022)
t.test(Monthly_average_wet_days_before_2022,Monthly_average_wet_days_2022, alternative = 'two.sided') # conduct two-sided T test on 2022 vs earlier data
## 
##  Welch Two Sample t-test
## 
## data:  Monthly_average_wet_days_before_2022 and Monthly_average_wet_days_2022
## t = -3.5894, df = 6.0667, p-value = 0.0113
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -12.06626  -2.29918
## sample estimates:
## mean of x mean of y 
##  9.760983 16.943704

 

How does 2022 compare to the next wettest years?

Apart from 2022, which were the next wettest years? The table below shows that the next wettest year was 1950, then 1990. This further suggests that, apart from 2022, there does not appear to be a recent trend towards wetter years.

# https://gt.rstudio.com/articles/intro-creating-gt-tables.html
top_ten <- rainfall_year %>% arrange(desc(mean_rainfall)) %>% slice(1:10) %>%
       mutate(mean_rainfall = round(mean_rainfall, digits = 1),
              mean_wet_days = round(mean_wet_days, digits = 1)) %>%
       select(-wet_days_colour, -rainfall_colour)

library(gt)
top_ten %>% gt() %>%
  tab_header(title = "Wettest known years for Sydney Botanical Gardens",
    subtitle = "Average rain per day in mm - Jan to July") %>%
  cols_label(
    mean_rainfall = "Average daily rainfall, mm",
    mean_wet_days = "Wet days per month"
  )
Wettest known years for Sydney Botanical Gardens
Average rain per day in mm - Jan to July
Year Average daily rainfall, mm Wet days per month
2022 10.1 17.9
1950 8.5 14.6
1990 7.5 12.5
1956 7.0 13.5
1913 6.6 11.3
1963 6.4 11.5
1955 6.2 13.4
1989 6.0 11.9
1919 5.7 10.4
1949 5.7 10.4

 

We could also look at the wettest single days, which are shown in the table below:

rainfall2 <- rainfall %>% arrange(desc(rainfall_mm))
top_ten_days <- rainfall2[1:10,]
top_ten_days  %<>%  mutate(Date = paste0(Day,' ',Month,' ',Year)) %>% 
       select(Date, rainfall_mm)  %>% select(2,3) 
top_ten_days <- top_ten_days[,c(2:3)]
top_ten_days %>% gt() %>%
  tab_header(title = "Wettest known days for Sydney Botanical Gardens") %>%
  cols_label(
    rainfall_mm = "Total rainfall, mm")
Wettest known days for Sydney Botanical Gardens
Date Total rainfall, mm
6 August 1986 340.0
28 March 1942 277.4
9 February 1992 264.0
9 November 1984 248.2
3 February 1990 238.0
1 May 1988 230.0
2 May 1953 226.1
11 March 1975 217.2
1 May 1955 193.0
11 February 1956 191.3

 

None of these ten wettest days have occurred in the last 10 years.

 

Is there a statistically significant long term trend?

We could look to see if there has a statistically significant trend line in the weather, again with data averaged on a month-by-month basis. The graph below shows the outcome of this analysis on average monthly rainfall. A separate regression analysis showed that the best fit line had a ‘p-value’ of 0.24, which means that we cannot consider this trend to be statistically significant. Note that this analysis was based on corrected mean monthly data, as discussed for the anova analysis.

monthly_sum %<>% mutate(Date = as.Date(paste0(1,'-',Month,'-',Year),format = "%d-%B-%Y"))
reg_model_a <- lm(corrected_mean_rainfall ~ Date, data = monthly_sum)  # Run regression model for trial data
#summary(reg_model_a)
monthly_sum %>%
  ggplot(aes(x = Date, y = corrected_mean_rainfall)) +
  geom_point(color = "grey68") +
   theme_bw() +
  geom_smooth(data = monthly_sum, aes(x = Date, y = corrected_mean_rainfall, color = "deepskyblue3"), formula = y ~ x, method = "lm", fill = "light blue", size = 1.5) +
 labs(title = "Fig 6. - No significant trend in average daily rainfall \nin Sydney Botanic Gardens",
              subtitle = "Best fit line is shown, but the slight gradient was not statistically significant. \nData was averaged on a month by basis basis",
              caption = "Data source: BOM",
              x = "", y = "Average daily \nrainfall in mm") +
  theme(axis.title.y=element_text(angle=0, vjust = 0.5, size = 10)) + 
  theme(legend.position="none")

However - when we looked at the number of wet days, there is a statistically significant trend, with the average number of wet days increasing from around 8 to around 11.

reg_model_b <- lm(corrected_mean_wet_days ~ Date, data = monthly_sum)  # Run regression model for trial data
#summary(reg_model_b)
monthly_sum %>%
  ggplot(aes(x = Date, y = corrected_mean_wet_days)) +
  geom_point(color = "grey68") +
   theme_bw() +
  geom_smooth(data = monthly_sum, aes(x = Date, y = corrected_mean_wet_days, color = "deepskyblue3"), formula = y ~ x, method = "lm", fill = "light blue", size = 1.5) +
 labs(title = "Fig 7. - It is raining more often at Sydney Botanic Gardens",
              subtitle = "Best fit line is shown -  trend was significant at 1% level.",
              caption = "Data source: BOM",
              x = "", y = "Number \nof wet days \nper month") +
  theme(axis.title.y=element_text(angle=0, vjust = 0.5, size = 10)) + 
  theme(legend.position="none")

 

Variability in rainfall

On a similar basis, we can also investigate whether the variability in rainfall in increasing. In the graph below, I have determined the variability, as defined by standard deviation of the average (corrected) monthly rainfall, for each five year block since 2010:

monthly_sum$fives <- fives[findInterval(monthly_sum$Year, fives)]
alpha <- 0.05
five_sumb <- monthly_sum %>% 
    group_by(fives) %>% 
    summarize(standard_deviation = sd(corrected_mean_rainfall)) %>%
             mutate(colour = if_else(fives == 2020,'firebrick1','lightskyblue'))

five_sumb %>%
  ggplot(aes(x = fives, y = standard_deviation)) +
  geom_bar(stat = "identity", fill = five_sumb$colour) +
  labs(title = "Fig 9. - Variability in monthly average rainfall in Sydney Botanic Gardens",
              subtitle = "Shows the standard deviation in average rainfall \non a month-by-month basis for each five year period.",
              caption = "Data source: BOM",
              x = "", y = "Standard deviation \n in monthly \nrainfall,mm") +
  theme(axis.title.y=element_text(angle=0, vjust = 0.5, size = 10), 
        plot.title = element_text(size = 12)) + 
   annotate("text", x = 1995, y = 5, label = "Red bar is 2020 onwards")

This graph shows that while the most recent period (2020 to July 2022) has greater variability, there has no been no strong trend prior to this. The high standard deviation for the recent period will be due to the very high rainfall for 2022.

 

Conclusions

In the introduction, I asked:

This data analysis shows that the 2022 is much wetter than any previous year in Sydney. In terms of average monthly rainfall there is no long-term trend, but there is a long term increase in the number of wet days per month.

Regardless, it will be worth watching to see if 2022 is the start of a longer term trend.