Data literacy with Python

GeoPandas and Map Shapefiles

Python and Pandas have some powerful packages to make use of map data. Map data can be imported for example as a shapefile and visualised along with data using GeoPandas.

These shapefiles are not images based on pixels, but use vectors to effectively redraw the postcode boundaries in whatever environment it is used. They provide much more definition than is needed to display on a screen on a web page, so there are tools in Python, namely topojson, which resample the vectors to suit a webpage, greatly simplifying the file and making it smaller and quicker to render.

Once again, the map/shapefile data for Postcodes in Australia is readily available from the ABS website, with a suitable Creative Commons Licence.

To reduce time and resources needed for this workflow, we have created a preprocessed shapefile named pc_sf_raw.shz, which is already a simplified version of the map file and can be loaded with the code below.

pc_sf_raw = geopandas.read_file('pc_sf_raw.shz')
pc_sf_raw.head(5)

See also the list of points to draw the postcode boundary/polygon for Bentley, Perth (6102).

pc_sf_raw[pc_sf_raw['POA_CODE21'] == '6102']['geometry'].get_coordinates()

Full citation for source of modified map/shapefile:

Australian Bureau of Statistics (2021) 'Non ABS Structures: Postal Areas - 2021 [https://www.abs.gov.au/statistics/standards/australian-statistical-geography-standard-asgs-edition-3/jul2021-jun2026/access-and-downloads/digital-boundary-files]' [Shapefile], Digital boundary files: Australian Statistical Geography Standard (ASGS) Edition 3, accessed 27th February 2024.

How was this file created?

The file was created using the code below, there is no need to execute the code in this workflow, it can take a while and depending on the prequantize/epsilon values chosen (which here determine the level of simplification for each postcode’s boundaries), can cause ‘out of memory’ type errors.

import topojson as tp
pc_sf_url = 'https://www.abs.gov.au/statistics/standards/australian-statistical-geography-standard-asgs-edition-3/jul2021-jun2026/access-and-downloads/digital-boundary-files/POA_2021_AUST_GDA2020_SHP.zip'
pc_sf_in = geopandas.read_file(pc_sf_url).to_crs(epsg=4283)
pc_sf_topo = tp.Topology(pc_sf_in.dropna(), prequantize=False, topology=True)
pc_sf_raw = pc_sf_topo.toposimplify(epsilon=0.006).to_gdf()
pc_sf_raw.to_file('pc_sf_raw.shz')

Demonstration 1 - GeoPandas Visualisation

Here we create a visualisation demonstration including a map of Australia with postcodes shown in colours reflecting the Index of Education and Occupation, and hovering over each postcode will display a label detailing the combined tax/seifa data from earlier in abbreviated form.

Though explaining the code for the visualisation is beyond the scope of this workflow, a powerful visualisation has been created with relatively little code. We can also see the effects of not ‘cleaning’ the data earlier. We can see areas missing data, where tax data was summarised into ‘Other’ categories, or there was a change in postcodes between the two data sources. The merge commands from earlier couldn’t find a match between the two datasets for these Postcodes and thus there is no corresponding data.

# Create a dataframe from earlier # Tax data combined workflow
tax_seifa = (tax2022_raw.query('~State.isin(["Unknown","Overseas"])')
            .assign(TaxableIncome_dollarspr = lambda x: round(x.TaxableIncome_dollars/x.Returns/1000,0))
            .assign(PrivateHealth_percentpp = lambda x: round(x.PrivateHealth_returns/x.Returns*100,0))
            .merge(seifa2021_raw, how="inner", on="Postcode"))

# Combine map shapefile and tax data
pc_sf = pc_sf_raw.merge(tax_seifa, left_on='POA_CODE21', right_on='Postcode', how="inner")

# Add a label to data which combines all of the tax data into a single abbreviated field
pc_sf["Details"] = "PCode:" + pc_sf["POA_CODE21"] + " Income:$" + pc_sf["TaxableIncome_dollarspr"].astype('str') + "K PrivHlth:" + pc_sf["PrivateHealth_percentpp"].astype('str') + "% IEO:" + pc_sf["ieo_percentile"].astype('str')

# Demo Visualisation 1
pc_sf.explore("ieo_percentile",
              cmap = "YlOrRd_r",
              scheme='quantiles',
              tiles = "cartodb voyager",
              tooltip = ["Details"],
              style_kwds = dict(color='black',weight=0.25),
              highlight_kwds = dict(color="black",weight=0.75,fillOpacity=0.75),
              legend_kwds = dict(caption="SEIFA Index of Education/Occupation percentile"))