How to create a point distance marker layer along a line in PostGIS

Posted by: Gavin Fleming | in GRASS GIS, PostGIS, QGIS | 6 days, 14 hours ago | Comments

If you look at the Cape Town Marathon interactive web map you'll see kilometre markers along the race routes.  Creating these is not as straighforward as you might think. For the team planning the marathon and the other routes, it used to involve a painstaking manual process in Google Earth.

I had my eye on a solution involving linear referencing from the start and tried a couple of approaches before settling on the one I describe in detail below. I did find solutions using v.segment from GRASS and the QGIS LRS plugin but I felt these were too complex for a non-GIS specialist and they also involve several steps.

My final solution was a view in PostGIS that just sits there untended on the server and generates the points automatically every time someone views the web map. The planners can edit the routes without worrying about the markers. When they save a route update, the markers just work.

I first created an event table with the ids of the routes and the distances at which I wanted markers. This is a sample:

Marker events

Then because we're working in the local CRS for Cape Town, which does not come standard in GIS software, I had to add the the WGS19 CRS (called 'HBK_NO19' in QGIS) to PostGIS by running this query, giving it the arbitrary id of 40019:

​INSERT INTO spatial_ref_sys (srid,proj4text) 
VALUES (40019,'+proj=tmerc +lat_0=0 +lon_0=19 +k=1 +x_0=0 +y_0=0 +axis=enu +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs');

I then created a view like this for each route.

CREATE OR REPLACE VIEW markers_marathon AS
(SELECT ST_Transform(m.geom,40019) AS geom, mm.dist, ST_Length(geom::geography) AS length
FROM marathon m
     JOIN marker_events mm
ON m.route_id = mm.route)
SELECT row_number() over()::int4 AS id, ST_Lineinterpolatepoint(geom,dist/length) AS geom,
CASE WHEN round(dist/1000.0,1)::character varying LIKE '%.0' THEN (dist/1000) ELSE round(dist/1000.0,1) END AS dist
    FROM m19
    WHERE dist <= length;

The WITH subquery prepares a table that has the route in each row projected to the local CRS, along with the marker distance from the event table and the length of the route (very accurately calculated on the spheroid and taking into account elevation changes, since the geometry has Z values).

The main part of the query then interpolates a point along the route at each marker distance and does some formatting of the distance so it can be used as a label in the map. It also filters out any marker that might go beyond the length of the route (which would be meaningless to interpolate and would cause the query to fail).

I had to project the routes because they are stored unprojected (EPSG:4326) and I wanted to minimise distortion for the linear referencing, which operates only on geometries and not geographies. You'll notice the ST_Length() was done with a geography, which is a PostGIS type that forces functions applied to it to operate on the WGS84 spheroid rather than on a plane. A length calculated on the spheroid is more accurate than a length calculated in any planar projection and is further enhanced if there are Z values in the coordinates (ST_Length() on a 3D geography returns a 3D length).

If you have an elegant solution to the same problem I'd love to know about it... Using generate_series() like this instead of joining to an event table is pretty elegant.

Extra note:

How come there are just a few markers when zoomed out and then one every km when zoomed in? We did this with rule-based rendering in QGIS with scale dependency in each rule and this filter in the 'zoomed out rule': dist % 2 = 0, which filters out everything that has a remainder when divided by two (or five in the case of the marathon).

Current rating: 5


Template by Blacktie Mezzanine theme by CodingHouse