I have been using Mapnik, Postgresql and Postgis to create topographic maps using OpenStreetMap (OSM) data. In OSM things like gates are often mapped as points located on or near a road. When rendering it is nice if the gate icon is across the road. For example, the US Forest Service shows a gate as a red dumbbell across a road. I want to replicate that.
This same technique is useful for aligning a waterfall icon along a stream or river.
Disclaimer: My use of SQL has been limited to simple single and double table selects and I’ve never written a stored database function before. So what I describe below may not be the best way to achieve the goal. Simply a way that I came up with when I couldn’t find an example to copy off the web.
Mapping is done by volunteers with varying levels of skill and training so we may find a gate mapped on a point other than the road. If on a road then the gate is likely to be at the beginning or end of a OSM “way”. The reason for that is the access conditions are likely to change at the gate and that would be a reason to “split” the road into separate ways in OSM.
What I was trying for was a function that could be used in a Postgresql select statement that would take a point and return a bearing or azimuth for which way the point was facing. For a gate it would allow it to be positioned across the road. The same logic would allow a waterfall icon to be aligned along the waterway it is associated with.
From the top level, we want stored database function that takes a point and returns a direction. If things go wrong (no nearby road, etc.) we will return a bearing/direction of zero so the rendering will have an automatic default provided and not need special case logic in the style XML.
CREATE OR REPLACE FUNCTION highway_poi_azimuth(highway_point geometry)
RETURNS double precision AS $$
DECLARE
angle double precision DEFAULT null;
BEGIN
--
-- Some Database magic happens here
--
IF (angle is null) THEN
angle = 0.0;
END IF;
RETURN angle;
END; $$
LANGUAGE PLPGSQL;
The first thing to do is find the road nearest our point. The “&&” operator must have some sort of overloading when using the Postgis extensions otherwise I don’t know how this works. I guess I really need to learn more about Postgresql and Postgis to understand the details. In any case, we expand the point to a bounding box of 100 meters then find all the highway ways/lines that intersect that box. We are interested in the way nearest the point, so we sort them and limit our selection to a single row:
SELECT degrees(ST_Azimuth(ST_LineInterpolatePoint(way, -- -- Some Postgis magic happens here -- INTO angle FROM planet_osm_line WHERE way && st_expand(highway_point, 100) AND (highway IS NOT null) ORDER BY ST_distance(way, highway_point) ASC LIMIT 1;
Looking through the various Postgis functions available, we can find a point on the way that is closest to our gate point (remember while it is likely the gate is on the way, it is not guaranteed). So:
ST_closestpoint(way, highway_point)
Will get us the point we actually want to work with. Now we need to figure out how to get the tangent to the way at that point. Or if not the tangent a fairly close estimate of it. Unfortunately there does not seem to be a Postgis function for that. There is a Postgis function that returns an azimuth between two points and a function that can convert the radians to degrees:
degrees(ST_Azimuth(point1, point2))
So if we can get two points to work with, we can estimate the tangent.
There is a function that returns how far a point is along a line as a value between 0 and 1. And there is another function that returns the length of the line. With those we can figure out how many meters our point of interest is along the line:
ST_LineLocatePoint(way, ST_closestpoint(way, highway_point)) * ST_Length(way)
ST_LineLocatePoint() returns a value between 0 and 1 for where a point is located on a line. Multiplying that by the length of the way returned by ST_Length() gives us a distance, in meters, along the line to our point.
We can estimate the tangent by taking points a small distance on either side and then use the ST_Azimuth() function to give us a bearing. A point 10 meters before the gate will be at:
ST_LineInterpolatePoint(way, GREATEST(((ST_LineLocatePoint(way, ST_closestpoint(way, highway_point)) * ST_Length(way))-10)/ST_Length(way), 0.0))
The ST_LineInterpolatePoint() function takes a value between 0 and 1, basically the inverse of ST_LineLocatePoint(). After subtracting a small distance, in our case 10 meters, we divide by the length of the line to normalize it. Since that might go negative we use the GREATEST() function to assure we are not below zero.
The same can be used to find a point 10 meters beyond the gate:
ST_LineInterpolatePoint(way, LEAST(((ST_LineLocatePoint(way, ST_closestpoint(way, highway_point)) * ST_Length(way))+10)/ST_Length(way), 1.0)))
Here we use the LEAST() function to limit our value to 1.0. Now that we have two points for ST_Azimuth() to use we can put it all together. Our function is:
CREATE OR REPLACE FUNCTION highway_poi_azimuth(highway_point geometry)
RETURNS double precision AS $$
DECLARE
angle double precision DEFAULT null;
BEGIN
-- Find highway that is nearest to the point of interest and figure out
-- the tangent at there.
-- So we find points 25 meters above and below the closest and
-- compute the azimuth from those.
SELECT degrees(ST_Azimuth(ST_LineInterpolatePoint(way, GREATEST(((ST_LineLocatePoint(way, ST_closestpoint(way, highway_point)) * ST_Length(way))-10)/ST_Length(way), 0.0)),
ST_LineInterpolatePoint(way, LEAST(((ST_LineLocatePoint(way, ST_closestpoint(way, highway_point)) * ST_Length(way))+10)/ST_Length(way), 1.0))))
INTO angle
FROM planet_osm_line
WHERE way && st_expand(highway_point, 100) AND (highway IS NOT null)
ORDER BY ST_distance(way, highway_point) ASC
LIMIT 1;
IF (angle is null) THEN
angle = 0.0;
END IF;
RETURN angle;
END; $$
LANGUAGE PLPGSQL;
The Mapnik layer that uses this has a SQL select statement that looks like:
(SELECT way,barrier,highway_poi_azimuth(way)::text as direction
FROM planet_osm_point WHERE ("barrier" is not null)) as foo
When we render that same area as the US Forest Service map we now get a properly aligned gate:
Here is another location that has several gates:
