OpenStreetMap
This post should be in my glob, bu it’s down for the moment :(
For a long time now I’ve been thinking on a problem: OSM data sometimes contains riverbanks that have no centerline. This means that someone mapped (part of) the coasts of a river (or stream!), but didn’t care about adding a line that would mark its centerline.
But this should be computationally solvable, right? Well, it’s not that easy. See, for given any riverbank polygon in OSM’s database, you have 4 types of segments: those representing the right and left riverbanks (two types) and the flow-in and flow-out segments, which link the banks upstream and downstream. With a little bit of luck there will be only one flow-in and one flow-out segment, but there are no guarantees here.
One method could try and identify these segments, then draw a line starting in the middle of the flow-in segment, calculating the middle by traversing both banks at the same time, and finally connect to the middle for the flow-out segment. Identifying the segments by itself is hard, but it is also possible that the result is not optimal, leading to a jagged line. I didn’t try anything on those lines, but I could try some examples by hand…
Enter topology, the section of maths that deals with this kind of problems. The skeleton of a polygon is a group of lines that are equidistant to the borders of the polygon. One of the properties this set of lines provides is direction, which can be exploited to find the banks and try to apply the previous algorithm. But a skeleton has a lot of ‘branches’ that might confuse the algo. Going a little further, there’s the medial axis, which in most cases can be considered a simplified skeleton, without most of the skeleton branches.
Enter free software :) CGAL
is a library that can compute a lot of topological properties. PostGIS is clever
enough to leverage those algorithms and present, among others, the functions
ST_StraightSkeleton() and ST_ApproximateMedialAxis(). With these two and the
original polygon I plan to derive the centerline. But first an image that will
help explaining it:
The green ‘rectangle’ is the original riverbank polygon. The thin black line is the skeleton for it; the medium red line is the medial. Notice how the medial and the center of the skeleton coincide. Then we have the 4 branches forming a V shape with its vertex at each end of the medial and its other two ends coincide with the ends of the flow in and flow out segments!
So the algorithm is simple: start with the medial; from its ends, find the branches in the skeleton that form that V; using the other two ends of those Vs, calculate the point right between them, and extend the medial to those points. This only calculates a centerline. The next step would be to give it a direction. For that I will need to see if there are any nearby lines that could be part of the river (that’s what the centerline is for, to possibly extend existing rivers/centerlines), and use its direction to give it to the new centerline.
For the moment the algorithm only solves this simple case. A slightly more
complex case is not that trivial, as skeletons and medials are returned as a
MultiLineString with a line for each segment, so I will have to rebuild them
into LineStrings before processing.
I put all the code
online, of course :)
Besides a preloaded PostgreSQL+PostGIS database with OSM data, you’ll need
python3-sqlalchemy, geoalchemy, python3-fiona and python3-shapely. The
first two allows me to fetch the data from the db. Ah! by the way, you will need
a couple of views:
CREATE VIEW planet_osm_riverbank_skel AS SELECT osm_id, way, ST_StraightSkeleton (way) AS skel FROM planet_osm_polygon WHERE waterway = 'riverbank';
CREATE VIEW planet_osm_riverbank_medial AS SELECT osm_id, way, ST_ApproximateMedialAxis (way) AS medial FROM planet_osm_polygon WHERE waterway = 'riverbank';
Shapely allows me to manipulate the polygonal data, and fiona is used to save the results to a shapefile. This is the first time I ever use all of them (except SQLAlchemy), and it’s nice that it’s so easy to do all this in Python.
Discussion
Comment from imagico on 26 July 2016 at 20:12
Not sure if you do this as a helper tool for mapping or for data use.
For mapping this works nicely - have done this as well a few times for long and winded riverbank polygons because i was too lazy to draw the centerline by hand. Would be useful as a JOSM plugin by the way.
For data use as a preprocessing step for missing data? Forget it if you want to process things globally - unless you have a real lot of computing power to spare. Skeleton algorithms generally scale poorly with polygon complexity, generating skeletons for big polygons is slow.
Comment from Warin61 on 26 July 2016 at 23:48
While a mapper may take the ‘center line’ of a river to be a geometric mean of the river banks boats and ships will probably prefer to take it as the deepest part of the river for navigational safety purposes and if the river ever gets low then this will still be the rivers center line.
So .. include a source statement of your ‘center-line’. And refrain from changing one already entered unless you know the intention/source of that center line.
Comment from BushmanK on 27 July 2016 at 00:26
@Warin61, Since
waterway=riverand similar tags have never been described in documentation, as something, related to clearways, no single sane person will ever use this geometry for waterway navigation. Use of unofficial third-party data for navigation, which could affect safety, is often explicitly prohibited by national regulators. So, I wouldn’t worry about that.I’m guessing, there is no rivers with center line mapped according to bathymetry in OSM.
@imagico,
It would be great to have some common spatial functions in JOSM some day, similar to ones of PostGIS, in form of more or less visual interface.
Comment from Marcos Dione on 27 July 2016 at 09:29
@Warin61: Also, remember this is just a tool for generating the centerline. If I ever use it to import into OSM, it will surely be with a note of ‘check this actually makes sense’. @imagico’s idea of a Josm plugin is most probably the outcome.
Comment from joost schouppe on 28 July 2016 at 13:54
I recently tried the inverse of this same problem: trying to measure the width of a polygon. Problem is that it only works decently if your polygons are long (as opposed to square). Same would apply here. If for some reason, the mapper has cut the riverbank up into short and wide bits, then I would guess your process breaks down.
Comment from Marcos Dione on 28 July 2016 at 14:19
@joost: yes, in fact, so far the tool has no way to guess the original river’s flow direction (the riverbanks do not contain that info) and in any case with branches, ST_ApproximateMedialAxis() returns LineStrings with almost arbitrary direction. At least the segments the tool adds are consistent with the direction of the LineStrings they’re added to.
@imagico: re: exploding cases: I don’t need to find a huge riverbank to see the algo meltdown, with the last segment of the Rhône river is enough. I tell you, the skeleton and resulting medial are not pretty and bordering the useless. I’m playing with the idea of some polygon simplification at the end of the processing pipeline and polygon partitioning at the beginning.
Comment from pnorman on 30 July 2016 at 05:39
When I looked at this problem I concluded it would be far simpler to add the missing data then add a complex step to data transformations. The tools have gotten better since then, but I’d still rather fix the data once for everyone.