Calculating minimum river width for a route within the river in PostGIS/Python

Question

I have a complex ancient river system containing islands and multiple channels. I also have line segments defining a vessel's path along the river:

I use PostGIS to hold the data, the river system is saved both as a multipolygon and a set of multilinestrings (for the outlines). The paths are simple lines.

For each given route (the points are sure to be in the river, but they are not in the middle of the river), I would like to calculate (or estimate) the minimum river width. I have been looking for approaches to this problem, but are unsure on what the best approach would be:

• My naive approach would be to take each point along the line and check the closest shapes.
• For each given point, I would have to find the shortest line that connects to the nearest shores. I have found some hints on this topic in the math section of StackExchange, but I have a multipolygon and that makes things complicated. *If I had this shortest line for each point, it would be trivial to calculate the minimum river width along this section.

How would I be able to do this in either PostGis and/or Python?

• I have thought about cutting the river into smaller segments somehow, mabye this would help.
• I might also try to create a buffered shape (using ST_buffer) on the line. Unfortunately, my line is not in the middle of the river, so this only helps partially. On the other hand, I might be able to identify shapes that touch/cut this new polygon and use this to simplify my problem...
• I guess, I could separate islands from the river itself (all shapes contained by another one are islands) - this might simplify things a bit - or not?

Any hints?

Answer 1

Difficult question, the holes are tricky.

One idea would be to start by getting the lines of the polygon. The difficulty is that you need to separate them to be sure to don't have both sides in the same line. For example on your image you can take a buffer big enough around the line between your two points to be sure to have what you need but not big enough to have the line on the bottom and the line on the top to be in the same linestring.

You would do something like ST_Intersection(ST_Boundary(poly), ST_Buffer(line, big_distance).

Then once you have your different linestrings of riverside, you can check the distance between the 2 closest lines, and they should be your two closest riverside. Once you have their id, you can get the distance between the 2 linestrings of riverside directly.

It can be a bit difficult notably if your buffer is not big enough to always have the line of a riverside completly inside, in that case both your closests could be from the same side... Also you could select the size of another branch of the river without knowing.

Another option would be, once you have all the riverside lines, not only check the distance between them, but also get the ST_Closest_Point for each line from each other, and make a line between both points to see if it cross your travel line. That way you can select the branch that you actually have crossed, and you also have in bonus the cross section where your river is the narrowest.

Anyway I think the first step you should take is to properly select the linestring for all the riversides, and distinguish between them.

Finally, you can also take a look at the ST_StraightSkeleton function, maybe there is something to do with it.

Answer 2

Edit:

I just realized you want the minimum river width - this answer can't help you with that, but I leave it here in case someone needs it.

---

Get the line length of the ST_ShortestLine between your routes and the ST_Boundary of the rivers (using their MultiPolygons):

SELECT
  rt.<id>,
  ST_Length(sl_geom::GEOGRAPHY) AS distance,
  sl_geom AS geom
FROM
  <routes> AS rt
  JOIN
  <river> AS rv
    ON ST_Intersects(rt.geom, rv.geom)
  CROSS JOIN LATERAL
    ST_ShortestLine(rt.geom, ST_Boundary(rv.geom)) AS sl_geom
;

Notes:

• only considers pairs of routes and rivers that actually ST_Intersects
• calculates the spheroidal ST_Length (in meter) by casting to a GEOGRAPHY; requires coordinates to be in a geographic reference system (e.g. EPSG:4326)
• returns the actual ST_ShortestLine as geom
  
  
• use a LATERAL expression for ST_ShortestLine to reuse its geometry in ST_Length

This will get you the closest distance to the shoreline for each individual LineString and for each intersecting MultiPolygon. Get the minimum of these distances via grouping if you want several of those LineStrings to be treated as a single route.

Answer 3

I found a solution, thanks to @robin loche's comment! I post my solution here in case others have similar questions.

The solution uses ST_ApproximateMedialAxis - it is part of the postgis_sfcgal extension of Postgis and creates aproximate medial axes within shapes, so i can be used to calculate the distance to shorelines.

Preparation

Make sure, you have the extension installed:

CREATE EXTENSION postgis_sfcgal;

I have three tables in my archaeological db:

• topology.water_body containing the water body/river system
• topology.recrivers reconstructed river travel paths (simple linepaths)
• topology.water_lines containing the river system's borders as linestrings.

The last table can be created from your river system shape using a SELECT...INTO query:

SELECT (dump_set).path[1] as id, (dump_set).geom as geom
into topology.water_lines
FROM (SELECT ST_Dump(ST_Boundary(geom)::geometry(MULTILINESTRING, 4326)) as dump_set
      from topology.water_body) as dump_results;

As you can see, I use EPSG:4326, you might want to adjust this if you use something else. So, for my example, I use the data from my question above. The line shown has the id 3.

The Solution

I had problems using ST_ApproximateMedialAxis, because my colleagues from archaeology created an insanely complex river system with touching rings... The solution was to extract only a part of the river system and check this smaller part. So, I get the bounding box of my path, expand it a bit and clip that:

SELECT clipped.geom as axis
INTO topology.clipped_water_body
from (SELECT (ST_Dump(ST_ClipByBox2D(geom, (SELECT ST_Expand(Box2D(geom), 0.006) AS arr
                                            FROM topology.recrivers
                                            WHERE id = 3)))).geom AS geom
      FROM topology.water_body) as clipped
WHERE ST_contains(clipped.geom, (SELECT geom FROM topology.recrivers WHERE id = 3)) = true;

Remember, the id of my shape is 3, thus the query. I also excluded all geometries that do not contain my path (in my example, this would be the case close to the northeast corner). This will create a clipped subpart of my insane river shape - and there will be only one shape left (no touching rings, hopefully):

A second query will create the approximated ideal medians in my single geometry:

SELECT ST_ApproximateMedialAxis(axis) as axis
INTO topology.ideal_paths
FROM topology.clipped_water_body

Of course, I can skip the creation of the clipped_water_body table, but I wanted to make clear what I am doing. This will create a web of lines, as seen in this closeup:

Now, here it becomes interesting. The medial axis is designed to be in the middle of each river arm, so the closest shore point of any point along the axes leads to half the river width! I just need a sample. I can do this by taking the points (using ST_DumpPoints) of the path I want to take through the river (or rather that of the vessel travelling the path). I can create the closest points to the median axes:

SELECT ST_ClosestPoint((ST_DumpPoints(geom)).geom, (SELECT axis FROM topology.ideal_paths)) as geom
INTO topology.closest_points
FROM topology.recrivers
where id = 3;

This will yield a number of sample points on the median axes:

We are almost done now! We now have a couple of samples we can check. The rest is a bit of SQL magic:

SELECT 2 * MIN(ST_Distance(cp.geom::geography, st_closestpoint(wl.geom, cp.geom)::geography)) as min_width
FROM topology.closest_points as cp,
     topology.water_lines as wl

Here we compare the distance of each of our water lines (i.e. the shorelines) to each of our sample points. We then take the minimum of all distances and multiply by 2 - this yields be the minimum river width from our sample points - or at least a good approximation.

Summary

Using ST_ApproximateMedialAxis and your path, you can create sample points to measure the distance to the closest shore line. Naturally, this is an approximation only.

Possible Enhancements

You could enhance this method by creating more points on your original path (using ST_LineInterpolatePoints). This will yield better results and might catch some weird edge cases.

This will be done by changing the closest points-calculation into the following two statements:

-- create enough sample points for a good approximation
SELECT (ST_Dump(ST_LineInterpolatePoints(geom, 0.02))).geom
INTO topology.sample_points
FROM topology.recrivers
where id = 3;

-- closest points in the ideal axis
SELECT ST_ClosestPoint(geom, (SELECT axis FROM topology.ideal_paths)) as geom
INTO topology.closest_points
FROM topology.sample_points;

This created 50 samples, as shown here:

Feedback welcome!