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 all about meeting two requirements that the team planning the marathon routes had:
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 identify tool function that shows the X, Y, Z and M coordinates of the nearest vertex in the 'derived' section, like this:
That point is 2.27km along the 5km fun run route and 13m above sea level. Note: In QGIS 3.4 this will also show the measure (distance) interpolated to exactly where you clicked.
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
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:
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.