This is the second in a three part series on the behind-the-scenes GIS work that can go into planning a complex event, in this case the Cape Town Marathon.
Following on from creating point markers along a line in PostGIS this post delves into 4D coordinates, PostGIS functions and rasters in PostGIS.
What I show here is about meeting two requirements of the marathon route planning team:
- Automatically update the length of a route after each edit with the most accurate calculation possible
- Enable a QGIS user to query the distance along a line at any point by clicking on it
For the first one, we want to calculate the length along the spheroid, which is easy enough, but we can make it even more accurate by taking into account elevation as the route goes up and down hills.
For the second one, we want to take advantage of the QGIS 3.4 identify tool function that shows the X, Y, Z and M coordinates of the nearest vertex as well as the interpolated M and Z values in the 'derived' section, like this:
That point is 15.74km along the marathon route and 7.5m above sea level.
We are going to update the Z (elevation) and M (measure) values of every vertex every time the geometry is edited and saved.
For elevation we're going to use a raster DEM (digital elevation model) within the database. The City of Cape Town has an accessible, open data portal and I found a 10m DEM there. So I downloaded it and then loaded it into my PostGIS database:
raster2pgsql -C -I -M -F -Y -s 40019 -t 100x100 10m_BA.tif dem | psql -d mydb
The marathon and long trail routes in QGIS 3D
Then I had to make sure that each route's geometry field catered for 4 dimensions:
ALTER TABLE marathon ALTER COLUMN geom TYPE geometry(LinestringZM,4326) USING ST_Force4D(geom);
Next I created this function:
--you need a DEM called 'dem' in SRID WGS19 loaded in the database for this to work
--adapted from http://blog.mathieu-leplatre.info/drape-lines-on-a-dem-with-postgis.html
CREATE OR REPLACE FUNCTION update_ZM()
RETURNS trigger AS
EXECUTE format('WITH line AS
(SELECT ST_Transform($1.geom,40019) as geom),
(SELECT (ST_DumpPoints(geom)).geom AS geom FROM line),
(SELECT p.geom AS geom, ST_Value(d.rast, 1, p.geom) AS val
FROM %I.dem d, points2d p
WHERE ST_Intersects(d.rast, p.geom)),
(SELECT ST_SetSRID(ST_MakePoint(ST_X(geom), ST_Y(geom), val), 40019) AS geom FROM cells),
(SELECT ST_MakeLine(ST_Transform(geom,4326)) AS geom FROM points3d),
(SELECT ST_AddMeasure(geom,0,st_length(geom::geography)/1000) AS geom FROM line3d)
UPDATE %I.%I a
SET geom = line4d.geom
WHERE a.id = $1.id', TG_TABLE_SCHEMA, TG_TABLE_SCHEMA, TG_TABLE_NAME) USING NEW;
update_ZM() function does the following:
- projects the line to the local CRS
- extracts all the vertices as points
- intersects the points with the DEM
- adds the elevation from the DEM to the Z value of the point coordinates
- builds the points back into a line
- adds accurate 3D distances to the M value of the coordinates as described in more detail in the previous article
- replaces the previous version of the line with the one generated by the above process
The function is a dynamic trigger function, which means it will run when a trigger is fired and work on any table that has a line geometry.
The last step is to add the trigger to any line table you want to run this function on:
DROP TRIGGER update_zm ON marathon;
CREATE TRIGGER update_zm
AFTER INSERT OR UPDATE OF geom
FOR EACH ROW
WHEN (pg_trigger_depth() = 0)
EXECUTE PROCEDURE update_ZM();
This trigger runs the
update_ZM() function on each record that you add or edit. If you are editing in QGIS, this means as soon as you hit Save.