---
title: "Hot and Unequal: How Tree Canopy Divides Cumberland"
author: "Chandima Rathnayake - s4150663"
date: "2026-07-20"
output:
  flexdashboard::flex_dashboard:
    orientation: rows
    vertical_layout: scroll
    theme: cosmo
    source_code: embed
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE, warning = FALSE, message = FALSE,
                       out.width = "600px", fig.width = 6, fig.height = 4.2)

library(tidyverse)
library(readxl)
library(janitor)
library(plotly)
library(scales)
library(flexdashboard)

# Brand colours from THE CONVERSATION
tc_critical       <- "#D40404"  # signature red - high heat / high disadvantage
tc_accent         <- "#0F71F0"  # functional blue - low heat / low disadvantage
tc_text           <- "#242424"
tc_text_secondary <- "#5B5B5B"
tc_text_disabled  <- "#919191"
tc_bg_light       <- "#F4F4F4"
tc_green          <- "#2E7D32"  # design choice for vegetation cue, not a brand colour

theme_tc <- theme_minimal(base_family = "sans") +
  theme(
    plot.title       = element_text(color = tc_text, face = "bold", size = 14),
    plot.subtitle    = element_text(color = tc_text_secondary, size = 10.5),
    axis.title       = element_text(color = tc_text_secondary, size = 10),
    axis.text        = element_text(color = tc_text_secondary, size = 9),
    panel.grid.minor = element_blank(),
    panel.grid.major = element_line(color = tc_bg_light),
    legend.position  = "bottom",
    legend.title     = element_text(color = tc_text, size = 9),
    strip.text       = element_text(color = tc_text, face = "bold", size = 9),
    plot.background  = element_rect(fill = "white", color = NA)
  )
theme_set(theme_tc)
```

```{r data-load}
uhi_raw    <- read_csv("data/nsw-urban-heat-island-to-sa2-level-2016.csv") %>% clean_names()
canopy_raw <- read_csv("data/greater-sydney-region-tree-canopy-to-sa2-level-20220.csv") %>% clean_names()
seifa_raw  <- read_excel("data/Statistical Area Level 2, Indexes, SEIFA 2021.xlsx",
                          sheet = "Table 1", skip = 5)

colnames(seifa_raw) <- c("sa2_code", "sa2_name", "irsd_score", "irsd_decile",
                          "irsad_score", "irsad_decile", "ier_score", "ier_decile",
                          "ieo_score", "ieo_decile", "population")

seifa <- seifa_raw %>%
  filter(!is.na(sa2_code), !is.na(sa2_name)) %>%
  mutate(
    sa2_name   = str_trim(sa2_name),
    
    irsd_score = as.numeric(gsub(",", "", irsd_score)),
    population = as.numeric(gsub(",", "", population))
  )

sa2_name_lookup <- tibble(
  sa2_main16     = c(125031482, 125031485),
  sa2_name_fixed = c("Greystanes - Pemulwuy", "Merrylands - Holroyd")
)

uhi <- uhi_raw %>%
  left_join(sa2_name_lookup, by = "sa2_main16") %>%
  mutate(sa2_name = coalesce(str_trim(sa2_name21), sa2_name_fixed)) %>%
  select(sa2_main16, sa2_name, uhi_mean)

canopy <- canopy_raw %>%
  mutate(sa2_name = str_trim(sa2_name16)) %>%
  rename(canopy_pct = matches("^canopy")) %>%
  select(sa2_name, canopy_pct)

# core_all: all 16 Cumberland SA2s (temperature + canopy only)
# core_residential: 12 SA2s with a SEIFA score. 4 excluded - 3 have no score
# due to negligible residential population (ABS methodology); 1 (Greystanes - Pemulwuy) was split into two SA2s under the 2021 ASGS boundary edition
# SEIFA uses, vs. the 2016 edition the UHI/canopy data use.
core_all <- uhi %>%
  inner_join(canopy, by = "sa2_name") %>%
  left_join(seifa, by = "sa2_name")

core_residential <- core_all %>% filter(!is.na(irsd_decile))
```

Row {data-height=130}
-------------------------------------

### Peak Heat Intensity

```{r vb-heat}
hottest <- core_all %>% filter(uhi_mean == max(uhi_mean, na.rm = TRUE))
valueBox(paste0(round(hottest$uhi_mean, 1), "°C"),
         caption = paste("Hottest suburb:", hottest$sa2_name),
         icon = "fa-thermometer-full", color = tc_critical)
```

### Lowest Tree Canopy

```{r vb-canopy}
least_canopy <- core_all %>% filter(canopy_pct == min(canopy_pct, na.rm = TRUE))
valueBox(paste0(round(least_canopy$canopy_pct, 1), "%"),
         caption = paste("Lowest canopy:", least_canopy$sa2_name),
         icon = "fa-tree", color = tc_green)
```

### Most Disadvantaged Areas

```{r vb-seifa}
valueBox(sum(core_residential$irsd_decile == min(core_residential$irsd_decile, na.rm = TRUE)),
         caption = "of 12 residential suburbs sit in SEIFA decile 1 (most disadvantaged nationally)",
         icon = "fa-users", color = tc_accent)
```

Row {data-height=480}
-------------------------------------

### The Heat Gap {data-width=500}

*Every one of Cumberland's 16 statistical areas experiences a measurable Urban Heat Island effect — but the intensity varies by several degrees from one suburb to the next.*

```{r chart-1}
p1_data <- core_all %>%
  arrange(uhi_mean) %>%
  mutate(sa2_name = fct_reorder(sa2_name, uhi_mean))

p1 <- ggplot(p1_data, aes(x = sa2_name, y = uhi_mean,
                           text = paste0("<b>", sa2_name, "</b><br>",
                                         "UHI intensity: ", round(uhi_mean, 1), "°C"))) +
  geom_col(fill = tc_critical, alpha = 0.88) +
  coord_flip() +
  labs(title = "Where the Heat Hits Hardest",
       subtitle = "Urban Heat Island intensity across Cumberland's suburbs, Summer 2015-16",
       x = NULL, y = "UHI Intensity (°C)")

ggplotly(p1, tooltip = "text", width = 500, height = 420) %>% layout(margin = list(l = 200))
```

*Figure 1: Auburn - North runs nearly 3.5°C hotter than Rookwood Cemetery, the LGA's coolest area.*

### The Canopy Deficit {data-width=500}

*Tree canopy cover mirrors the heat pattern in reverse — with one striking exception.*

```{r chart-2}
p2_data <- core_all %>%
  mutate(sa2_name = factor(sa2_name, levels = levels(p1_data$sa2_name)))

p2 <- ggplot(p2_data, aes(x = sa2_name, y = canopy_pct,
                           text = paste0("<b>", sa2_name, "</b><br>",
                                         "Tree canopy: ", round(canopy_pct, 1), "%"))) +
  geom_col(fill = tc_green, alpha = 0.88) +
  coord_flip() +
  labs(title = "The Green Divide",
       subtitle = "2022 tree canopy cover (%) reveals stark spatial inequality",
       x = NULL, y = "Tree Canopy Cover (%)")

ggplotly(p2, tooltip = "text", width = 500, height = 420) %>% layout(margin = list(l = 200))
```

*Figure 2: Rookwood Cemetery breaks the pattern entirely — the coolest area in the LGA also has by far the highest canopy cover, at 24%.*

Row {data-height=520}
-------------------------------------

### Less Canopy, Higher Heat, Greater Disadvantage

*Suburbs with lower socioeconomic status (intense red points) tend to combine higher heat with lower canopy cover — bearing the compounded weight of the urban heat island effect. Point size shows population.*

```{r chart-3}
p3 <- ggplot(core_residential, aes(x = canopy_pct, y = uhi_mean,
                           color = irsd_score, size = population,
                           text = paste0("<b>", sa2_name, "</b><br>",
                                         "Canopy: ", round(canopy_pct, 1), "%<br>",
                                         "UHI: ", round(uhi_mean, 1), "°C<br>",
                                         "SEIFA IRSD score: ", irsd_score,
                                         " (lower = more disadvantaged)<br>",
                                         "Population: ", comma(population)))) +
  geom_point(alpha = 0.85) +
  scale_color_gradient(low = tc_critical, high = tc_accent,
                        name = "SEIFA IRSD Score\n(low = more disadvantaged)") +
  scale_size(range = c(3, 15), guide = "none") +
  labs(title = "When Postcode Determines Temperature",
       subtitle = "Lower canopy, higher heat, and greater disadvantage compound together",
       x = "Tree Canopy Cover (%)", y = "UHI Intensity (°C)")

ggplotly(p3, tooltip = "text", width = 600, height = 420)
```

Row {data-height=480}
-------------------------------------

### Disadvantage, Canopy and Heat Compound Together

*Grouping suburbs by disadvantage tier shows the pattern at a glance —  notably, none of Cumberland's residential areas fall into the "lower disadvantage" range found elsewhere in Sydney.*

```{r chart-4}
p4_data <- core_residential %>%
  mutate(
    disadvantage_tier = case_when(
      irsd_decile <= 3 ~ "Higher Disadvantage (Decile 1-3)",
      irsd_decile <= 7 ~ "Mid (Decile 4-7)",
      TRUE ~ "Lower Disadvantage (Decile 8-10)"
    ),
    disadvantage_tier = factor(disadvantage_tier,
                                levels = c("Higher Disadvantage (Decile 1-3)",
                                           "Mid (Decile 4-7)",
                                           "Lower Disadvantage (Decile 8-10)")),
    sa2_name = fct_reorder(sa2_name, canopy_pct)
  )

p4 <- ggplot(p4_data, aes(x = sa2_name, y = canopy_pct, fill = uhi_mean,
                           text = paste0("<b>", sa2_name, "</b><br>",
                                         "Canopy: ", round(canopy_pct, 1), "%<br>",
                                         "UHI: ", round(uhi_mean, 1), "°C"))) +
  geom_col() +
  facet_wrap(~disadvantage_tier, scales = "free_y") +
  coord_flip() +
  scale_fill_gradient(low = tc_accent, high = tc_critical, name = "UHI (°C)") +
  labs(title = "The Compounding Effect",
       subtitle = "Disadvantage tiers show which suburbs bear the greatest burden",
       x = NULL, y = "Tree Canopy Cover (%)")

ggplotly(p4, tooltip = "text", width = 600, height = 420)
```

Row {data-height=520}
-------------------------------------

### Where Should Cumberland Plant Next?

*The lower-left quadrant marks suburbs that are both more disadvantaged and have below-median canopy cover — the clearest priority targets for future greening investment.*

```{r chart-5}
med_irsd   <- median(core_residential$irsd_score, na.rm = TRUE)
med_canopy <- median(core_residential$canopy_pct, na.rm = TRUE)

p5 <- ggplot(core_residential, aes(x = irsd_score, y = canopy_pct,
                                    size = population, color = uhi_mean,
                                    text = paste0("<b>", sa2_name, "</b><br>",
                                                  "IRSD score: ", irsd_score, "<br>",
                                                  "Canopy: ", round(canopy_pct, 1), "%<br>",
                                                  "UHI: ", round(uhi_mean, 1), "°C<br>",
                                                  "Population: ", comma(population)))) +
  geom_vline(xintercept = med_irsd, linetype = "dashed", color = tc_text_disabled) +
  geom_hline(yintercept = med_canopy, linetype = "dashed", color = tc_text_disabled) +
  geom_point(alpha = 0.85) +
  scale_color_gradient(low = tc_accent, high = tc_critical, name = "UHI (°C)") +
  scale_size(range = c(3, 15), guide = "none") +
  labs(title = "Where Cumberland Should Plant Next",
       subtitle = "Priority suburbs in the lower-left quadrant need urgent greening investment",
       x = "SEIFA IRSD Score (lower = more disadvantaged)",
       y = "Tree Canopy Cover (%)")

ggplotly(p5, tooltip = "text", width = 600, height = 420)
```

Row {data-height=200}
-------------------------------------

### References & Acknowledgements

**Data sources:**

Cumberland City Council. (n.d.). *NSW Urban Heat Island to SA2 level 2016* [Data set]. Cumberland City Council Open Data Portal. https://data.cumberland.nsw.gov.au/explore/dataset/nsw-urban-heat-island-to-sa2-level-2016/

Cumberland City Council. (n.d.). *Greater Sydney Region Tree Canopy to SA2 level 2022* [Data set]. Cumberland City Council Open Data Portal. https://data.cumberland.nsw.gov.au/explore/dataset/greater-sydney-region-tree-canopy-to-sa2-level-20220/

Australian Bureau of Statistics. (2023). *Socio-Economic Indexes for Areas (SEIFA), Australia, 2021* [Data set]. https://www.abs.gov.au/statistics/people/people-and-communities/socio-economic-indexes-areas-seifa-australia/2021

**Methodology notes:** The UHI data is area-weighted and aggregated from the NSW Mesh Block-level dataset (Summer 2015-16). Tree canopy data, derived from 0.48m-resolution imagery captured by ArborCarbon Pty Ltd in the "Greater Sydney Region Tree Canopy to Modified Mesh Block 2022" dataset, is aggregated to the SA2 level. Four of Cumberland’s 16 SA2s are excluded from the disadvantage-linked charts: three due to negligible residential population under the ABS SEIFA methodology, and one (Greystanes - Pemulwuy) due to a boundary mismatch between the 2016 and 2021 SA2 editions. Brand colours were extracted from The Conversation’s live website CSS via browser DevTools inspection.

**GenAI acknowledgement:** Generative AI (ChatGPT) was utilised to assist with debugging R code syntax, and brainstorming the initial narrative structure for the data visualisations. All aspects of data sourcing, data cleaning, visualisation design, statistical analysis, and the final written text were completed independently by myself.