92577c60ba327ee07320311808efefc4ddfdc686
[R/project-au-taxstats.git] / AU-taxstats.R
1 install.packages("sf")
2 library(sf)
3 install.packages("devtools")
4 devtools::install_github("tidyverse/ggplot2") # needed for geom_sf
5 library(ggplot2)
6 library(viridis)
7 library(ggthemes)
8 library(animation) # for saveGIF
10 # Obtain the tax dataset if not available yet
11 if(!file.exists("data/taxstats2015individual06ataxablestatusstateterritorypostcode.csv")) 
12   download.file(url = "http://data.gov.au/dataset/5c99cfed-254d-40a6-af1c-47412b7de6fe/resource/90f7f4eb-2c44-4884-96c0-01060c820cfd/download/taxstats2015individual06ataxablestatusstateterritorypostcode.csv", destfile = "data/taxstats2015individual06ataxablestatusstateterritorypostcode.csv")
13 # http://data.gov.au/dataset/5c99cfed-254d-40a6-af1c-47412b7de6fe/resource/d3189e9d-533a-4893-b6a1-758781083418/download/taxstats2015individual06btaxablestatusstateterritorypostcode.csv
15 # Obtain shapefile with Australian postal codes if not available yet
16 if(!file.exists("data/1270055003_poa_2016_aust_shape.zip"))
17   download.file(url = "http://www.abs.gov.au/ausstats/subscriber.nsf/log?openagent&1270055003_poa_2016_aust_shape.zip&1270.0.55.003&Data%20Cubes&4FB811FA48EECA7ACA25802C001432D0&0&July%202016&13.09.2016&Latest", destfile = "data/1270055003_poa_2016_aust_shape.zip")
18 # Unzip it if not done already
19 if(!file.exists("data/POA_2016_AUST.shp")) 
20   unzip(zipfile = "data/1270055003_poa_2016_aust_shape.zip", exdir = "data/")
22 taxstats <- read.csv("data/taxstats2015individual06ataxablestatusstateterritorypostcode.csv", stringsAsFactors = FALSE)
23 taxstats <- dplyr::filter(taxstats, Taxable.status == "Taxable")
24 POA <- st_read(dsn = "data/", layer = "POA_2016_AUST", stringsAsFactors = FALSE)
26 taxstats.POA <- merge(x = taxstats, y = POA, by.x = "Postcode", by.y = "POA_CODE16", all.y = TRUE)
28 taxstats.POA$incomeperearningcapita <- taxstats.POA$`Total.Income.or.Loss..` / taxstats.POA$Total.Income.or.Loss.no. 
29 # Postal codes turn out not to be too interesting, as they're way more granular around 
30 # big cities - making the high income postal codes invisible on the chart below
31 ggplot(taxstats.POA) +
32   geom_sf(aes(fill = incomeperearningcapita, color = incomeperearningcapita)) +
33   scale_fill_viridis("incomeperearningcapita") +
34   scale_color_viridis("incomeperearningcapita")
36 # Let's try by SA3
37 if(!file.exists("data/1270055001_sa3_2016_aust_shape.zip")) 
38   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")
40 if(!file.exists("data/SA3_2016_AUST.shp")) 
41   unzip(zipfile = "data/1270055001_sa3_2016_aust_shape.zip", exdir = "data/")
43 sa3 <- st_read(dsn = "data/", layer = "SA3_2016_AUST", stringsAsFactors = FALSE)
45 # Create a matrix of intersecting postal codes and SA3's
47 POA_SAs <- st_intersects(x=sa3, y=POA, sparse=FALSE)
48 taxstats.POA$`Total.Income.or.Loss..`[is.na(taxstats.POA$`Total.Income.or.Loss..`)] <- 0
49 taxstats.POA$`Total.Income.or.Loss.no.`[is.na(taxstats.POA$`Total.Income.or.Loss.no.`)] <- 0
50 # Perform matrix multiplication to obtain the income metrix per SA3
51 # Total income will be incorrect, as the POAs intersect with multiple SA3s
52 sa3$TotalIncome <- as.vector(POA_SAs %*% as.matrix(taxstats.POA$`Total.Income.or.Loss..`))
53 sa3$TotalIncomeEarners <- as.vector(POA_SAs %*% as.matrix(taxstats.POA$`Total.Income.or.Loss.no.`))
54 sa3$incomeperearningcapita <- (sa3$TotalIncome / sa3$TotalIncomeEarners)/1000
56 # As SA3s are still to narrow around cities compared to in the country,
57 # let's simply look at Melbourne
59 cities = c("Perth","Melbourne","Sydney","Adelaide","Brisbane")
61 # Create a plot for each of these cities. This is wrapped in a function
62 # for use by saveGIF
64 plots <- function() {lapply(cities, function(x){
65  plot <- ggplot(dplyr::filter(sa3, data.table::`%like%`(GCC_NAME16, x) )) +
66   geom_sf(aes(fill = incomeperearningcapita, color = incomeperearningcapita)) +
67   scale_fill_viridis(name = "",
68                      limits = c(min(sa3$incomeperearningcapita, na.rm = TRUE),max(sa3$incomeperearningcapita, na.rm = TRUE))) +
69   scale_color_viridis(name = "",
70                                 limits = c(min(sa3$incomeperearningcapita, na.rm = TRUE),max(sa3$incomeperearningcapita, na.rm = TRUE))) +
71   coord_sf(datum = NA) + # Work around https://github.com/tidyverse/ggplot2/issues/2071 to remove gridlines
72   labs(title = paste0(x," \nincome distribution"),
73        subtitle = "2014/15, in 1000s AUD",
74        caption = "\nSource: Australian Taxation Office") +
75   theme_economist() +
76   theme(legend.position = "bottom",
77         legend.text = element_text(angle = 45, hjust = 1, size = 8),
78         axis.text = element_blank(),
79         axis.ticks = element_blank())
80  print(plot)
81 })
82 }
84 saveGIF(plots(),movie.name = "AUCitiesIncomeDistribution.gif", interval = 2, loop = 2)