PostGIS for Australian data: GDA2020, MGA zones and distances

Australian location data has a quirk most PostGIS tutorials skip: the continent moves. The Australian plate moves north-east by about 7 centimetres a year, so the country has had two modern datums in a generation, GDA94 and GDA2020, and the same street corner has slightly different coordinates in each. Get the datum wrong and your points are off by a metre or two. Get the type wrong and your "2 km radius" is measured in degrees. Here's how to set up PostGIS for Australian data so neither happens.

Three datums you will meet

A datum fixes coordinates to the earth. In Australia you will see three, each with an EPSG code that PostGIS uses as its SRID:

GDA2020 and GDA94 differ by roughly 1.5 to 1.8 metres depending on where you are. GDA2020 and WGS 84 are within a metre or so of each other today, and PostGIS treats a conversion between them as no change at all. For finding the nearest clinic or checking whether an address is in a delivery zone, that is fine. For surveying, cadastral or engineering work it is not, and you need the transformations your state's surveyor-general publishes. PostGIS also treats GDA94 to WGS 84 (4283 to 4326) as no change, so a GDA94 point “converted” to 4326 is still roughly 1.5 to 1.8 metres out: always transform GDA94 data to 7844.

All of these come built in: PostGIS ships the EPSG definitions in its spatial_ref_sys table, including GDA2020 and every MGA2020 zone.

Step 1: a column that measures in metres

PostGIS has two spatial types. geometry does flat-plane maths in the units of its coordinate system, which for latitude and longitude means degrees. geography does its maths on the curved earth and answers in metres. For points stored as latitude and longitude, geography is almost always what you want:

create table sites (
  id bigserial primary key,
  name text not null,
  location geography(Point, 7844) not null
);
create index sites_location_idx on sites using gist (location);

insert into sites (name, location) values
  ('Opera House', ST_SetSRID(ST_MakePoint(151.2153, -33.8568), 7844)::geography),
  ('Town Hall',   ST_SetSRID(ST_MakePoint(151.2065, -33.8732), 7844)::geography),
  ('Parramatta',  ST_SetSRID(ST_MakePoint(151.0036, -33.8150), 7844)::geography);

Note the order: ST_MakePoint takes longitude first, then latitude. Swap them and PostGIS does not raise an error: for geography it prints a notice, quietly forces the latitude into range and stores a point in the middle of the Atlantic. Check a row or two on a map after your first import.

Step 2: radius and nearest-first queries

Everything within 2 km of Sydney's GPO, using the index:

select name from sites
 where ST_DWithin(location,
                  ST_SetSRID(ST_MakePoint(151.2093, -33.8688), 7844)::geography,
                  2000);

The third argument is in metres because the column is geography. Use ST_DWithin rather than ST_Distance(...) < 2000: ST_DWithin can use the GiST index, while the ST_Distance comparison makes PostgreSQL calculate the distance to every row.

For the closest few, order by the <-> distance operator, which also walks the index:

select name,
       round(ST_Distance(location, ST_SetSRID(ST_MakePoint(151.2093, -33.8688), 7844)::geography)) as metres
  from sites
 order by location <-> ST_SetSRID(ST_MakePoint(151.2093, -33.8688), 7844)::geography
 limit 5;

--   name       | metres
--  Town Hall   |    553
--  Opera House |   1442
--  Parramatta  |  19952

The degrees trap

Store the same points as geometry and ask for a distance, and you get an answer in degrees, which is meaningless as a distance:

select ST_Distance(ST_SetSRID(ST_MakePoint(151.2093, -33.8688), 7844),
                   ST_SetSRID(ST_MakePoint(144.9631, -37.8136), 7844));
-- 7.39   (degrees, Sydney to Melbourne)

select ST_Distance(ST_SetSRID(ST_MakePoint(151.2093, -33.8688), 7844)::geography,
                   ST_SetSRID(ST_MakePoint(144.9631, -37.8136), 7844)::geography) / 1000;
-- 713.9  (kilometres)

A degree of longitude is also shorter in Hobart than in Darwin, so a "within 0.02 degrees" filter covers a different area in every city. If an existing table is geometry in degrees, cast to geography for distance queries, or convert the column.

Step 3: bringing GDA94 data forward

If you load an older dataset, check its metadata for the datum. A GDA94 column converts in place:

alter table old_sites
  alter column geom type geometry(Point, 7844)
  using ST_Transform(geom, 7844);

We measured the shift on PostGIS 3.6: a point at the Sydney Opera House moves 1.49 metres between GDA94 and GDA2020. Mixing the two datums without converting is the classic source of "the pin is on the wrong side of the road".

Step 4: areas, buffers and boundaries in MGA2020

geography handles distances, radius searches and areas, but a smaller set of functions. When you need planar operations in metres, such as buffering a site boundary or clipping polygons, transform into the Map Grid of Australia (MGA2020) zone that covers your data:

-- suburbs is the boundary table defined in the next section
select name, round(ST_Area(ST_Transform(boundary, 7856)) / 10000) as hectares
  from suburbs;

Each zone is six degrees of longitude wide. A query that spans several zones, or the whole country, should stay in geography instead.

Which suburb is this point in?

Boundary data such as the ABS's ASGS regions comes as polygons. Store them as geometry(MultiPolygon, 7844) with a GiST index, and join points to them:

create table suburbs (
  name text primary key,
  boundary geometry(MultiPolygon, 7844) not null
);
create index suburbs_boundary_idx on suburbs using gist (boundary);

select s.name as suburb, count(*) as sites
  from sites x
  join suburbs s on ST_Contains(s.boundary, x.location::geometry)
 group by s.name;

Point-in-polygon is a flat-plane test, so it uses geometry; the cast from the geography point keeps the SRID, so both sides are GDA2020.

Where the coordinates come from

For street addresses, G-NAF, the national address file, publishes a geocode for most addresses in GDA2020, so its coordinates go straight into a geography(Point, 7844) column. WattleAddr, our address autocomplete and verification API, returns G-NAF geocodes in GDA2020 for that reason. For points from phones and browsers, the coordinates are WGS 84; store them as 4326, or as 7844 if a metre doesn't matter to you, and keep a note of which you chose.

PostGIS on WattleDB

Every WattleDB database comes with PostGIS 3.6 installed, so there is nothing to request. It lives in a wattledb_ext schema that is on your search path, so ST_DWithin and the geography type resolve without a schema prefix, and every role in your database can use them. The coordinate tables take part of the roughly 16 MB that the catalogue and extensions use before you create a table. Raster, topology and pgRouting are available on request. The details are on the extensions page.

Your location data stays in a database hosted in Sydney by an Australian-owned company, with backups in Melbourne and point-in-time recovery within your plan’s recovery window. Like every query in this post, all of it is plain PostgreSQL and PostGIS, so it works wherever PostGIS 3 runs.