```{r, include = FALSE}
current_file <- knitr::current_input()
basename <- gsub(".Rmd$", "", current_file)

knitr::opts_chunk$set(
  fig.path = sprintf("images/%s/", basename),
  fig.width = 6,
  fig.height = 4,
  fig.align = "center",
  fig.retina = 3,
  echo = TRUE,
  warning = FALSE,
  message = FALSE,
  cache = TRUE,
  cache.path = "cache/"
)

if (!identical(getOption("bitmapType"), "cairo") && isTRUE(capabilities()[["cairo"]])) {
  options(bitmapType = "cairo")
}
```

```{r titleslide, child="assets/titleslide.Rmd"}
```

---

class: center middle bg-gray


.aim-box.tl.w-70[
Today you will:

- look at the ABS geographical boundaries for the 2021 census
- integrate data from different sources (census and election) to make exploratory inferences
]

--

.aim-box.tl.w-70[
Coding Perspective:

- Further expand our understanding of how to read and use spatial data in R
- Better understand how spatial data is organised 
- Learn how to intersect two spatial objections
- Practice re-projecting maps

]

---


# Recall <i class="fas fa-database"></i> 2021 Federal Election Data

.f6[
```{r aec-map, eval = FALSE}
library(tidyverse)
library(sf)
aec_map <- read_sf(here::here("data/vic-july-2021-esri/E_VIC21_region.shp")) #<<
votes <- read_csv("https://results.aec.gov.au/27966/Website/Downloads/HouseDopByDivisionDownload-27966.csv", skip = 1) #<<

electoral_winners = votes |> 
  mutate(DivisionNm = toupper(DivisionNm)) |>
  filter(Elected == "Y") |> 
  select(PartyAb, DivisionID, DivisionNm, Elected) |>
  distinct()

aec_map = aec_map |> 
  mutate(Elect_div = toupper(Elect_div))

winners_with_map = electoral_winners |> 
  left_join(aec_map, by = c("DivisionNm" = "Elect_div"))

aus_colours <- c(
  "ALP" = "#DE3533", "LNP" = "#ADD8E6", "KAP" = "#8B0000", "GVIC" = "#10C25B", "XEN" = "#ff6300",
  "LP" = "#0047AB", "NP" = "#0a9cca", "IND" = "#000000"
)

ggplot(winners_with_map) +
  geom_sf(aes(fill = PartyAb, geometry = geometry)) +
  scale_fill_manual(values = aus_colours)

#UPDATED CODE - previous variable assignments from past weeks 
winners_fix = winners_with_map
auscolors = aus_colours
```

]

---

# Recall <i class="fas fa-database"></i> 2021 Federal Election Data

.flex[
.w-50[
```{r aec-map, echo = FALSE}
```

]
.w-50[

There are two sources of data:

1. Electoral boundary
2. The votes for candidates in each electorate

]
]

---

# Recall <i class="fas fa-database"></i> 2021 ABS Census Data

* DataPacks <i class="fas fa-download"></i> https://datapacks.censusdata.abs.gov.au/datapacks/
* GeoPackages <i class="fas fa-download"></i> https://datapacks.censusdata.abs.gov.au/geopackages/

<img src="images/lecture-07/datapack-download.png" width = "70%"/>
---


class: transition

# ABS Census 2021

# GeoPackages

---

# GeoPackage

.blockquote[
A **GeoPackage** (GPKG) is an open, non-proprietary, platform-independent and standards-based data format for geographic information system implemented as a SQLite database container. Defined by the **Open Geospatial Consortium** (OGC) with the backing of the US military and published in 2014, GeoPackage has seen widespread support from various government, commercial, and open source organizations.

.right[
&mdash; Wikipedia
]
]

<br>

Recall: OGC also defines the WKT

---

# ABS GeoPackage (2021)

.info-box[
<i class="fas fa-download"></i> https://datapacks.censusdata.abs.gov.au/geopackages/

1. 2021
2. Victoria
3. Topic: Employment and Income or Table: G17 
4. GDA2020
]

* **Or use the [`strayr`](https://github.com/runapp-aus/strayr) package!** We'll use the one from the ABS website instead.

.f6.overflow-scroll.h-40[

```{r sa1-layers-2021}
# 2022
geopath_2021_G02 <- here::here("data/Geopackage_2021_G02_VIC_GDA2020/G02_VIC_GDA2020.gpkg")
st_layers(geopath_2021_G02)
```
]

---

# ABS GeoPackage (2016)

.info-box[
<i class="fas fa-download"></i> https://datapacks.censusdata.abs.gov.au/geopackages/

1. Victoria
2. Employment, Income and Unpaid Work (EIUW)
3. EIUW GeoPackage A
]

* **Or use the [`strayr`](https://github.com/runapp-aus/strayr) package!** We'll use the one from the ABS website instead.


.f6.overflow-scroll.h-40[

```{r sa1-layers-2016}
geopath_2016_eiuwa <- here::here("data/Geopackage_2016_EIUWA_for_VIC/census2016_eiuwa_vic_short.gpkg")
st_layers(geopath_2016_eiuwa)

#UPDATED CODE - old variable assignments from years
# geopath = geopath_2016_eiuwa 

```
]

---

# The Australian Statistical Geography Standard (ASGS)

.center[
<img src="images/lecture-08/acgs-structures.png" width="75%">

]

---

# The number of regions for each layer

.f5[
```{r}
st_layers(geopath_2021_G02) %>%
  # make it into a data.frame first
  tibble(!!!.) %>%
  # then you can the dplyr operations
  dplyr::arrange(features)
```

]

---

# <i class="fas fa-search-dollar"></i> Data in the layer

.f5.overflow-scroll.h-90[
```{r}
vicmap_ste_G02 <- read_sf(geopath_2021_G02, layer = "G02_STE_2021_VIC")
vicmap_ste_G02$geom
str(vicmap_ste_G02)

vicmap_ste = vicmap_ste_G02
```

]


---

name: ste

# State or Territory (STE)

.f6[
```{r map-ste}
vicmap_ste_G02 <- read_sf(geopath_2021_G02, layer = "G02_STE_2021_VIC")

ggplot(vicmap_ste_G02) +
  geom_sf(aes(geometry = geom, fill = Median_age_persons)) 

nrow(vicmap_ste)
```
]

---

# Breakout Session 

<center>
.aim-box.tl.w-70[
Try it yourself time:

- Download the geopackage for the 2021 census and variable G02
- Look at how that spatial data is organised 
- Look at the different layers

What are the differences between the regionalisations?  
What do you notice?  
Want a challenge visualise the 2016 data!  

]
</center>

---

name: gccsa

# Greater Capital City Statistical Areas (GCCSA)

* Each region with variable population

```{r map-gccsa}
vicmap_gccsa_G02 <- read_sf(geopath_2021_G02, layer = "G02_GCCSA_2021_VIC")
ggplot(vicmap_gccsa_G02) +
  geom_sf(aes(geometry = geom, fill = Median_age_persons))

nrow(vicmap_gccsa_G02)
```

---

name: sos

# Section of State (SOS)

* Major urban, other urban, bounded locally & rural balance

```{r map-sos}
vicmap_sos_G02 <- read_sf(geopath_2021_G02, layer = "G02_SOS_2021_VIC")
ggplot(vicmap_sos_G02) +
  geom_sf(aes(geometry = geom, fill = Median_age_persons))

nrow(vicmap_sos_G02)
```

---

name: ra

# Remoteness Areas (RA)

```{r map-ra}
vicmap_ra_G02 <- read_sf(geopath_2021_G02, layer = "G02_RA_2021_VIC")
ggplot(vicmap_ra_G02) +
  geom_sf(aes(geometry = geom, fill = Median_age_persons))

nrow(vicmap_ra_G02)
```

---

name: sosr

# Section of State Ranges (SOSR)

```{r map-sosr}
vicmap_sosr_G02 <- read_sf(geopath_2021_G02, layer = "G02_SOSR_2021_VIC")
ggplot(vicmap_sosr_G02) +
  geom_sf(aes(geometry = geom, fill = Median_age_persons))

nrow(vicmap_sosr_G02)
```

---

name: sa4

# Statistical Area Level 4 (SA4)

* Each region with population of 100,000 - 500,000

```{r map-sa4}
vicmap_sa4_G02 <- read_sf(geopath_2021_G02, layer = "G02_SA4_2021_VIC")
ggplot(vicmap_sa4_G02) +
  geom_sf(aes(geometry = geom, fill = Median_age_persons))

nrow(vicmap_sa4_G02)
```



---

name: sua

# Significant Urban Areas (SUA)

```{r map-sua}
vicmap_sua_G02 <- read_sf(geopath_2021_G02, layer = "G02_SUA_2021_VIC")
ggplot(vicmap_sua_G02) +
  geom_sf(aes(geometry = geom, fill = Median_age_persons))

nrow(vicmap_sua_G02)
```

---

name: ced

# Commonwealth Electoral Division (CED)


```{r map-ced}

vicmap_ced_G02 <- read_sf(geopath_2021_G02, layer = "G02_CED_2021_VIC")
ggplot(vicmap_ced_G02) +
  geom_sf(aes(geometry = geom, fill = Median_age_persons))

nrow(vicmap_ced_G02)
```
---

name: sa3

# Statistical Area Level 3 (SA3)

* Each region with population of 30,000 - 130,000

```{r map-sa3}
vicmap_sa3_G02 <- read_sf(geopath_2021_G02, layer = "G02_SA3_2021_VIC")
ggplot(vicmap_sa3_G02) +
  geom_sf(aes(geometry = geom, fill = Median_age_persons))

nrow(vicmap_sa3_G02)
```
---

name: lga

# Local Government Area (LGA)


```{r map-lga}
vicmap_lga_G02 <- read_sf(geopath_2021_G02, layer = "G02_SA3_2021_VIC")
ggplot(vicmap_lga_G02) +
  geom_sf(aes(geometry = geom, fill = Median_age_persons))

nrow(vicmap_lga_G02)
```

---

name: sed

# State Electoral Division (SED)

```{r map-sed}
vicmap_sed_G02 <- read_sf(geopath_2021_G02, layer = "G02_SED_2021_VIC")
ggplot(vicmap_sed_G02) +
  geom_sf(aes(geometry = geom, fill = Median_age_persons))

nrow(vicmap_sed_G02)
```
---

# Urban Centres and Localities (UCL)

```{r map-ucl}
vicmap_ucl_G02 <- read_sf(geopath_2021_G02, layer = "G02_UCL_2021_VIC")
ggplot(vicmap_ucl_G02) +
  geom_sf(aes(geometry = geom, fill = Median_age_persons))

nrow(vicmap_ucl_G02)
```

---

name: sa2

# Statistical Area Level 2 (SA2)

* Each region with populations in the range of 3,000-25,000

```{r map-sa2}
vicmap_sa2_G02 <- read_sf(geopath_2021_G02, layer = "G02_SA3_2021_VIC")
ggplot(vicmap_sa2_G02) +
  geom_sf(aes(geometry = geom, fill = Median_age_persons))

nrow(vicmap_sa2_G02)
```

---

name: poa

# Postal Areas (POA)

```{r map-poa}
vicmap_poa_G02 <- read_sf(geopath_2021_G02, layer = "G02_POA_2021_VIC")
ggplot(vicmap_poa_G02) +
  geom_sf(aes(geometry = geom, fill = Median_age_persons))

nrow(vicmap_poa_G02)
```

---

name: ssc

# State Area Localitites (SAL) (Formerly SSC)

```{r map-ssc}
vicmap_sal_G02 <- read_sf(geopath_2021_G02, layer = "G02_SAL_2021_VIC")
ggplot(vicmap_sal_G02) +
  geom_sf(aes(geometry = geom, fill = Median_age_persons))

nrow(vicmap_sal_G02)
```


---

name: sa1

# Statistical Area Level 1 (SA1)

* Each region with a population of range 200-800

```{r map-sa1, eval = TRUE}
vicmap_sa1_G02 <- read_sf(geopath_2021_G02, layer = "G02_SA1_2021_VIC")
ggplot(vicmap_sa1_G02) +
  geom_sf(aes(geometry = geom, fill = Median_age_persons))

nrow(vicmap_sa1_G02)
```

---

class: transition


# Electorate boundary <br>vs <br>Census boundary

🎯 Estimate a median age for an electorate 



---

# Comparing SED 2021 and electorates divisions 2022

.f4[See [here](https://wcd.numbat.space/lectures/lecture-05.html#35) for  `electoral_winners` data was.]

.f6[
```{r zoomed-map, fig.width = 8}
ggplot() +
  geom_sf(data = vicmap_sed_G02, 
          aes(geometry = geom, fill = Median_age_persons),
    alpha = 1, color = "white", size = 2) +
  geom_sf(data = winners_with_map, aes(geometry = geometry),
    fill = "transparent", color = "red", size = 2) +
  coord_sf(xlim = c(144.95, 145.24), ylim = c(-38.05, -37.85)) + 
  theme_bw()
```
]

---

# Closer look 🔬 Hotham electorate .circle.monash-bg-black.white[1]

.f5[

```{r district-map-no-eval, fig.show='hide', eval = FALSE, echo = TRUE}
electorate <- winners_with_map |>
  filter(DivisionNm == "HOTHAM") 

# Set projection to GDA1994 using EPSG:4283
st_crs(electorate$geometry,4283)
  
# Transform projection from GDA1994 to GDA2020 using EPSG:7844
electorate$geometry = st_transform(electorate$geometry, 7844) 

sed_intersect <- vicmap_sed_G02 |>
  filter(st_intersects(geom,
    electorate$geometry,
    sparse = FALSE
  )[, 1])

ggplot() +
  geom_sf(data = sed_intersect,
    aes(geometry = geom), color = "red", fill = "transparent") +
  geom_sf_text(data = sed_intersect,
    aes(label = SED_CODE_2021, geometry = geom), color = "red") +
  geom_sf(data = electorate, aes(geometry = geometry), fill = "transparent") +
  geom_sf_text(data = electorate, aes(geometry = geometry, label = DivisionNm))
```

] 

```{r district-map, fig.show='hide', eval = TRUE, echo = FALSE, warning=FALSE, message=FALSE}
electorate <- winners_with_map |>
  filter(DivisionNm == "HOTHAM") 

# Set projection to GDA1994 using EPSG:4283
st_crs(electorate$geometry, 4283)
  
# Transform projection from GDA1994 to GDA2020 using EPSG:7844
electorate$geometry = st_transform(electorate$geometry, 7844) 

sed_intersect <- vicmap_sed_G02 |>
  filter(st_intersects(geom,
    electorate$geometry,
    sparse = FALSE
  )[, 1])

sed_map <- ggplot() +
  geom_sf(
    data = sed_intersect,
    aes(geometry = geom),
    color = "red", fill = "transparent"
  ) +
  geom_sf_text(
    data = sed_intersect,
    aes(
      label = SED_CODE_2021,
      geometry = geom
    ),
    color = "red"
  ) +
  geom_sf(data = electorate, aes(geometry = geometry), fill = "transparent") +
  geom_sf_text(data = electorate, aes(geometry = geometry, label = DivisionNm))

```

---

# Closer look 🔬 Hotham electorate .circle.monash-bg-black.white[1]

There are `r nrow(sed_intersect)` SED regions that intersect with Hotham electorate.

```{r district-map-plot, echo = FALSE, fig.width = 7, fig.height = 7}
sed_map
```

---

# Closer look 🔬 Hotham electorate .circle.monash-bg-black.white[2]

.f5[

```{r district-map2, fig.show='hide', eval = TRUE}
sed_intersect2 <- sed_intersect |> #<<
  mutate(
    geometry = st_intersection(geom, electorate$geometry), #<<
    perc_area = 100 * st_area(geometry) / st_area(geom), #<<
    perc_area = as.numeric(perc_area)
  ) |> #<<
  filter(perc_area > 5) 
  
ggplot(sed_intersect2, aes(geometry = geometry)) +
  geom_sf(data = electorate) +
  geom_sf_text(
    data = electorate,
    aes(label = DivisionNm)
  ) +
  geom_sf(color = "red", aes(fill = Median_age_persons)) +
  geom_sf_text(
    aes(
      label = glue::glue("{SED_CODE_2021}
                        ({scales::comma(perc_area, 1)}%, {Median_age_persons})")
    ),
    color = "red"
  ) +
  theme(legend.position = "bottom")
```

]

---

# Closer look 🔬 Hotham electorate .circle.monash-bg-black.white[2]

* There are `r nrow(sed_intersect2)` SED areas with at least 5% intersection with the electoral area.
* **How would you characterise the median age for Hotham?**


```{r district-map2, echo = FALSE, fig.height = 6, fig.width = 6}
```

---

# Closer look 🔬 Hotham electorate .circle.monash-bg-black.white[3]

.flex[
.w-40[
```{r district-map2, echo = FALSE, fig.height = 6, fig.width = 6}
```
]
.w-60.f5[
**Strategy 1**

```{r}
sort(sed_intersect2$Median_age_persons)
```



]]

--

**Strategy 2**

```{r}
mean(sed_intersect2$Median_age_persons)
```



--


**Strategy 3**

```{r}
weighted.mean(
  sed_intersect2$Median_age_persons,
  sed_intersect2$perc_area
)
```


---

# Closer look 🔬 Hotham electorate .circle.monash-bg-black.white[4]

.flex[
.w-50[
.f5[
```{r sa1-intersect, eval = TRUE}
sa1_intersect <- vicmap_sa1_G02 %>% #<<
  filter(st_intersects(geom,
    electorate$geometry,
    sparse = FALSE
  )[, 1])

sa1_intersect2 <- sa1_intersect %>%
  mutate(
    geometry = st_intersection(geom, electorate$geometry),
    perc_area = 100 * st_area(geometry) / st_area(geom),
    perc_area = as.numeric(perc_area)
  ) %>% #<<
  filter(perc_area > 5)

ggplot(sa1_intersect) +
  geom_sf(color = "red", 
          aes( fill = Median_age_persons, geometry = geom)) +
  geom_sf(data = electorate, 
          color = "white", size = 2, fill = "transparent", 
          aes(geometry = geometry)) +
  theme(legend.position = "bottom")
```
]]

.w-50[
```{r sa1-intersect, echo = FALSE,fig.width = 6, fig.height = 6}
```
]]

---

# Closer look 🔬 Hotham electorate .circle.monash-bg-black.white[5]

.flex[
.w-40[
```{r sa1-intersect, echo = FALSE, fig.height = 6, fig.width = 6}
```
]
.w-60.f5[
**Strategy 1**

```{r}
fivenum(sa1_intersect2$Median_age_persons)
```



]]

--

**Strategy 2**

```{r}
mean(sa1_intersect2$Median_age_persons)
```



--


**Strategy 3**

```{r}
weighted.mean(sa1_intersect2$Median_age_persons, sa1_intersect2$perc_area)
```



--

**Strategy 4**

```{r sa1-region-histogram, fig.height=2}
ggplot(sa1_intersect2, aes(x = Median_age_persons)) +
  geom_histogram(binwidth = 1)
```

---

# Closer look 🕵️ Zero median age

(Hotham 2022)[https://www.aec.gov.au/profiles/vic/hotham.htm]

.flex[
.f6.w-50[
```{r strange-result}
sa1_intersect2 %>%
  filter(Median_age_persons == 0) %>%
  ggplot() +
  geom_sf() +
  geom_sf(
    data = electorate, color = "red",
    fill = "transparent",
    aes(geometry = geometry)
  )
```
]
.w-50[
<center>
<img src="images/lecture-08/2018-vic-hotham-detailed-map.jpg" width = "100%"/>
</center>

]
]

---

# Closer look 🔬 Hotham electorate .circle.monash-bg-black.white[6]

.flex[
.w-50.f5[

## Before

**Strategy 1**

```{r}
fivenum(sa1_intersect2$Median_age_persons)
```

**Strategy 2**

```{r}
mean(sa1_intersect2$Median_age_persons)
```

**Strategy 3**

```{r}
weighted.mean(sa1_intersect2$Median_age_persons, sa1_intersect2$perc_area)
```
]
.w-50.f5.pl3[

## After

```{r}
sa1_intersect3 <- sa1_intersect2 %>%
  filter(Median_age_persons != 0)
```





]]

--

**Strategy 1**

```{r}
fivenum(sa1_intersect3$Median_age_persons)
```



--

**Strategy 2**

```{r}
mean(sa1_intersect3$Median_age_persons)
```



--


**Strategy 3**

```{r}
weighted.mean(sa1_intersect3$Median_age_persons, sa1_intersect3$perc_area)
```



---

# Dorling Cartogram 

.flex[

.f6.w-50[

```{r zoomed-map-income-dorling}
sa1_intersect4 <- sa1_intersect %>%
  mutate(centroid = st_centroid(geom))

dorling_plot <- ggplot(sa1_intersect4) +
  geom_sf(
    data = electorate,
    aes(geometry = geometry), size = 4, fill = "grey60"
  ) +
  geom_sf(aes(geometry = centroid, color = Median_age_persons),
    size = 0.5, shape = 3
  ) +
  scale_color_viridis_c(name = "Median age", option = "magma")
```

]

.w-50[

```{r echo = FALSE}
dorling_plot
```

]
]

---

# Closer look 🔬 Hotham electorate .circle.monash-bg-black.white[7]

.f6[

```{r}
sa1_intersect5 <- sa1_intersect4 %>%
  filter(st_intersects(centroid, electorate$geometry, sparse = FALSE)[, 1],
    Median_age_persons != 0)
```

**Strategy 1**

```{r}
fivenum(sa1_intersect5$Median_age_persons)
```

**Strategy 2**

```{r}
mean(sa1_intersect5$Median_age_persons)
```

**Strategy 4**

```{r sa1-region-histogram2, fig.height=2}
ggplot(sa1_intersect5, aes(x = Median_age_persons)) +
  geom_histogram(binwidth = 1)
```

]

---

## Summary

<center>

.w-90[

<div class="idea-box">
<ul>
<li>There are many ways to characterise an electorate.</li>
<li>Estimates of median age of an electorate is more consistent using SA1 map data than SED map data.</li>
</ul>
</div>

]

--

.idea-box.tl.w-90[

* We looked at mapping the 2021 census boundaries and projected a summary of the census variable (i.e. median age) onto a 2022 electoral district

* Discovered some of the challenges with matching two different types of data

]

--

.info-box.tl.w-90[
Read [Forbes, Cook & Hyndman (2020) Spatial modelling of the two-party preferred vote in Australian federal elections: 2001–2016. *Australian & New Zealand Journal of Statisitcs*. ](https://onlinelibrary.wiley.com/doi/abs/10.1111/anzs.12292) for a more sophisticated approach to studying the census variables and election results together.
]

</center>

---

class: transition

#### Slides updated and maintained by Dr. Kate Saunders. Slides originally developed by Dr. Emi Tanaka

---

```{r endslide, child="assets/endslide.Rmd"}
```


