Plot income per SA3.

diff --git a/AU-taxstats.R b/AU-taxstats.R
index d5e3c821b42bcfb0cc6a774adae27badc8def159..bda1448b4102219f7a2f3c9d98d03c4a8cf012a9 100644 (file)
--- a/AU-taxstats.R
+++ b/AU-taxstats.R
@@ -30,3 +30,27 @@ ggplot(taxstats.POA) +
   geom_sf(aes(fill = incomeperearningcapita, color = incomeperearningcapita)) +
   scale_fill_viridis("incomeperearningcapita") +
   scale_color_viridis("incomeperearningcapita")
   geom_sf(aes(fill = incomeperearningcapita, color = incomeperearningcapita)) +
   scale_fill_viridis("incomeperearningcapita") +
   scale_color_viridis("incomeperearningcapita")

+
+# Let's try by SA3
+if(!file.exists("data/1270055001_sa3_2016_aust_shape.zip")) 
+  download.file(url = "http://www.abs.gov.au/AUSSTATS/subscriber.nsf/log?openagent&1270055001_sa3_2016_aust_shape.zip&1270.0.55.001&Data%20Cubes&43942523105745CBCA257FED0013DB07&0&July%202016&12.07.2016&Latest", destfile = "data/1270055001_sa3_2016_aust_shape.zip")
+
+if(!file.exists("data/SA3_2016_AUST.shp")) 
+  unzip(zipfile = "data/1270055001_sa3_2016_aust_shape.zip", exdir = "data/")
+
+sa3 <- st_read(dsn = "data/", layer = "SA3_2016_AUST", stringsAsFactors = FALSE)
+taxstats.sa3 <- merge(x = taxstats, y = sa3, by.x = "Postcode", by.y = "POA_CODE16", all.y = TRUE)
+
+# Create a matrix of intersecting postal codes and SA3's
+
+POA_SAs <- st_intersects(x=sa3, y=POA, sparse=FALSE)
+taxstats.POA$incomeperearningcapita[is.na(taxstats.POA$incomeperearningcapita)] <- 0
+
+# Perform matrix multiplication to obtain the income metrix per SA3
+# Total income will be incorrect, as the POAs intersect with multiple SA3s
+sa3$incomeperearningcapita <- as.vector(POA_SAs %*% as.matrix(taxstats.POA$incomeperearningcapita))
+
+ggplot(sa3) +
+  geom_sf(aes(fill = incomeperearningcapita, color = incomeperearningcapita)) +
+  scale_fill_viridis("incomeperearningcapita") +
+  scale_color_viridis("incomeperearningcapita")