Introduction

A once-in-100-year flood disaster has inundated Townsville and many other sunshine’s state suburbs, threatened lives after a slow moving monsoon dumped record-breaking rain over the Townsville region for days recently.

Ninenews
source: 9News, 04 February 2019

Here are some statistics released by regarding the recent record-breaking flood:

Purpose of Analysis

It was reported that a ten-day stretch from Jan. 28 to Feb. 6, Townsville reported more than a year’s worth of rainfall. This analysis examined the rainfall data collected by the Australian Bureau of Meteorology (BOM) in Townsville, from 01 January 1841 to 18 Feburary 2019 to examine and find insights for the severity of the recent flood, in particular:

Statement 1: 1257 mm of rain fell at Townsville Airport in 10 days between January 28 and February 6, passing the previous record of 925.5 mm set in January 1953. Source: Townsville Airport, accessed on 18 February 2019.

Statement 2: This is a once-in-100 year flood in Townsville. Source: [ABC News] (https://www.abc.net.au/news/2019-02-04/townsvillle-flood-crisis-in-pictures/10777212)

Here is the heatmap animation. To create the heatmaps simply follow the codes provided.

heatmap.gif

Similar Analysis

Daniel Ohem has conducted an analysis of Townsville’s (i) annual rainfall and (ii) rainfall over a 10 day period [The Most Amount of Rain over a 10 Day Period on Record] https://www.r-bloggers.com/the-most-amount-of-rain-over-a-10-day-period-on-record/. My analysis focuses on daily rainfall analysis by creating heatmaps using the ggplot package.

Method

Step 1: I went to the [BOM website] (http://www.bom.gov.au/jsp/ncc/cdio/wData/wdata?p_nccObsCode=136&p_display_type=dailyDataFile&p_stn_num=032040&p_startYear=), downloaded all year daily rainfall (mm) data for the ‘Townsville Aero’ station (station number: 032040, file name: “IDCJAC0009_032040_1800_Data.csv”).
Step 2: I cleaned the data and created new features using the dplyr package.
Step 3: I created multiple heatmaps using [looping] (https://www.r-bloggers.com/how-to-write-the-first-for-loop-in-r/).

Step 1: Create Data in R Format

Firstly I loaded all the required R packages and define folder paths.

# load libraries
library("ggplot2") # This is R most famous graphic package
library("ggthemes") # This is an add-on package for ggplot2
library("dplyr") # This is my favourite data transforming tool
library("lubridate") # This package provide the 'year' function I need
library("zoo") # enable as.yearmon function for heat-mapping 

# Define folder path and load CSV file
getwd()

## [1] "C:/[Local Folder]/R Local Folder/weather/townsville"

savepath = "C:/[Local Folder]/R Local Folder/weather/townsville/RData"
setwd(savepath)

Step 2: Data Cleaning and Transformation

I created a new gg dataset based on the raw rainfall data downloaded from BOM.

# Load and prepare rainfall files
rain <- read.csv("C:/[Local Folder]/R Local Folder/weather/townsville/Input/IDCJAC0009_032040_1800_Data.csv") # load the data and name it as 'rain'
str(rain) # a quick check of data fields

## 'data.frame':    28539 obs. of  8 variables:
##  $ Product.code                                  : Factor w/ 1 level "IDCJAC0009": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Bureau.of.Meteorology.station.number          : int  32040 32040 32040 32040 32040 32040 32040 32040 32040 32040 ...
##  $ Year                                          : int  1941 1941 1941 1941 1941 1941 1941 1941 1941 1941 ...
##  $ Month                                         : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ Day                                           : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ Rainfall.amount..millimetres.                 : num  0 6.6 16.5 205.5 175 ...
##  $ Period.over.which.rainfall.was.measured..days.: int  NA 1 1 1 1 1 1 1 1 1 ...
##  $ Quality                                       : Factor w/ 3 levels "","N","Y": 3 3 3 3 3 3 3 3 3 3 ...

# I created a new dataset 'gg' based on the raw rainfall data
gg <- rain %>%
  rename(rain = Rainfall.amount..millimetres.)%>% # rename the rainfall variable (the original name is too long)
  filter(!is.na(rain)) %>% # exlcude any date with rainfall data missing
  mutate(date = as.Date(paste(Year, Month, Day, sep="-"), origin = "1904-01-01")) %>% # reformat the date field
  mutate(weekday = as.factor(weekdays(date, abbreviate = TRUE))) %>% # create a field to index the weekday, e.g. Mon, Tue, Wed...
  mutate(yearmonth = as.factor(as.yearmon(date))) %>% # create a combined field with both year and month
  mutate(month = format(date,format="%B")) %>% # reformat the month field
  mutate(year = as.numeric(Year)) %>% # create a new year field as numeric value
  mutate(week = week(date)) %>% # create a new field that tells the week number, i.e. from 1 to 52
  mutate(monthweek = ceiling(day(date)/7)) %>% # create a new field tells the number of week within a month, ie. from 1 to 5
  mutate(range = cut(year, # create a new variable: range to categorise year in 10-year durations
                   breaks=c(-Inf, 1939, 1949, 1959, 1969, 1979, 1989, 1999, 2009, Inf), 
                   labels=c("90", "80", "70", "60", "50", "40", "30","20","10"))) %>% 
  mutate(range = as.numeric(as.character(range))) %>% # change the range variable to numeric 
  select(-Year) # remove unwanted variable

I rearranged the year, month and weekday variable categories in my preferred order.

## Reorder weekdays and months in prederred orther
gg$weekdayf <- ordered(gg$weekday, levels=c("Mon", "Tue", "Wed", "Thu","Fri", "Sat", "Sun"))
gg$monthweekf <- ordered(gg$monthweek, levels=c(seq(5,1)))
gg$monthf <- ordered(gg$month,levels=c("January","February","March","April",
                                       "May","June","July","August",
                                       "September","October","November","December"))

I plotted the daily rainfall data on a histogram. It was found that the maxinum daily rainfall recorded in the history was 549mm - within one day. This is crazy because the annual rainfall for Townsville was about [1143mm] (http://www.bom.gov.au/qld/townsville/climate_Townsville.shtml). This histogram indicates Townsville has a fairly dry tropical climate, with daily rainfall less than 10mm in most of the days.

# Summary of daily rainfall 
summary(gg$rain) # 500mm rainfall a day? also known as the "Night Of Noah Floods"

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   0.000   0.000   3.134   0.000 548.800

## Histogram
ggplot(gg, aes(rain)) + 
  geom_histogram(breaks=seq(0, 500, by = 10),
                 alpha = .8) +
  theme_bw()

Step 3: Make Heat-maps

I created a sequence of heatmaps to visualise the daily rainfalls. Each heatmap covers 10-years of rainfall data. I also provide an animated version below.

# Prepare looping variables
r <- seq(80, 10, by = -10) # this is the same as the range I defined
ystart <- seq(1940, 2010, by = 10) # this defines the start year of each heatmap
yend <- ystart+9 # this defines the end-year of each heatmap
ft.base = 16 # this is a graphic parameters in defining font size in the heatmaps 

ptm <- proc.time() # set a timer for the loop (it tooks a while to run)
# i = length(ystart)
for (i in 1:length(r)){
  
  # For each loop, filter data by year range
  dt_s <- gg %>%
    filter(range == r[i])

  # Create a heatmap
  hm_s <- ggplot(dt_s, aes(weekdayf, monthweekf, fill = rain)) + 
    geom_tile(colour = "white") + 
    facet_grid(year~monthf) + 
    scale_fill_gradient2(low="white", mid="black", high="red", #colors in the scale
                         midpoint=300, 
                           breaks=seq(0,max(gg$rain),50), #breaks in the scale bar
                           limits=c(0, 550)) +
    labs(title = paste0("Daily Rainfall (mm) in Townsville Airport, QLD: ", ystart[i], " to ", yend[i]),
         subtitle = "Unfold the 10-years rainfall data in Townsville",
         caption = "Bureau of Meteorology, Australia",
         x = "Weekday", y = "Week Number within Month") +
    # guides(fill=guide_legend(title="Daily Rainfall (mm)")) + 
    theme_bw(base_size = ft.base) +
    theme(plot.title = element_text(size=ft.base*1.5, color = "darkblue"), 
          plot.subtitle = element_text(size=ft.base*1, color = "red"),
          axis.text.x = element_text(angle = 90, hjust = 1))
  
  print(hm_s)
}

proc.time() - ptm

##    user  system elapsed 
##   20.29    1.01   21.79

Conclusion

Statement 1: 1257 mm of rain fell at Townsville Airport in 10 days between January 28 and February 6, passing the previous record of 925.5 mm set in January 1953. SUPPORTED The recent flood broke the 10-days rainfall record in 1953.

Statement 2: This is a once-in-100 year flood in Townsville. Inconclusive It depends. There can be more than one ‘once-in-100 year flood’ if we consider the total daily rainfall recorded within one day, i.e. the Night of Noah on 10 JAN 1998 (https://www.facebook.com/townsvillebulletin/posts/where-were-you-on-the-night-of-noah/10155198111091180/).

“Unfold How Townsville’s one-in-100 year flood crisis with data”

Lok Yiu

27 FEB 2019