Find start of river

Edit: Follow-up question here (= Find headwater polygons).

How can I determine the start of a river in PostGIS?

I have a river network (Multiline) and want to find the startpoints of the rivers. I can select the startpoints (rectangles in the graphic) using

SELECT
ST_StartPoint(ST_LineMerge(a.geom)) as stp, a.gid AS gid
FROM
spatial.stream AS a

and similarly the endpoints (crosses). But how can I find the start of the river - not of the line segments?

I tried something like (find the startpoints, that do not intersect with endpoints):

SELECT
ST_StartPoint(ST_LineMerge(a.geom)) as stp, a.gid AS gid
FROM
spatial.stream AS a, spatial.stream AS b
WHERE
ST_Disjoint(ST_StartPoint(ST_LineMerge(a.geom)), ST_EndPoint(ST_LineMerge(b.geom)))

But this takes forever:

"Nested Loop  (cost=0.00..1019343418.86 rows=525839789 width=323)"
"  Join Filter: st_disjoint(st_startpoint(st_linemerge(a.geom)), st_endpoint(st_linemerge(b.geom)))"
"  ->  Seq Scan on stream_typ b  (cost=0.00..4498.18 rows=39718 width=323)"
"  ->  Materialize  (cost=0.00..6364.77 rows=39718 width=319)"
"        ->  Seq Scan on stream_typ a  (cost=0.00..4498.18 rows=39718 width=319)"

Edit: I did not run it till it finished, so I even don't know if this query returns the desired result. But looking at the returned rows, this does not look correct (too many)).

I the end I want to find the polygons (grey in the graphic) where a/any river starts (the green polygons), but not those where the rivers just passes (the red polygon, no starts inside). But having just the startpoints would be good start!

Any idea how this could be done in PostGis? (also other open source solutions like GRASS, R etc. are welcome).

Update: My idea just to extract the startpoints was not concise enough :( E.g. consinder the following situation: The green polygons both have (true) startspoints inside and are those that I want. The red polygon has also startpoints inside, but the river flows through (so no headwater polygon). With John's solution below I get both. I only want the headwater streams. I am not sure if my idea with startpoints will lead to the solution (i make up a new question if desired). I though of two Where clauses:

• polygon contains a start point
• no flow through (= only 1 intersection of polygon with stream).

I tried this:

SELECT
polyg.*
FROM
polyg, start_points, stream
WHERE
st_contains(polyg.geom, start_points.geom)
AND ST_Npoints(ST_Intersection(poly.geom, stream.geom)) < 2

But it never finishes:( Perhabs there is a better way to query the headwater polygons ?

• Is the problem with what the query produces, or with its runtime? You have ST_LineMerge twice in the where clause, which is probably part of your problem, so I am wondering if you just want the query speed up, or if it is producing the wrong results. It isn't 100% clear from the question. – John Powell Jan 26 '15 at 16:22
• @JohnBarça: I clarified my question. I didn't run the query till it finished (aborted after 30 minutes or so), so I don't know if this query returns the desired results! – EDi Jan 26 '15 at 16:29
• OK, next question, where is the polygon to the right? – John Powell Jan 26 '15 at 16:30
• @JohnBarça: There is selected line segment in the graphic. This is a polygon where rivers start in. In the polygon right of it, the river just passes by (in the south of the polygon, no starts inside the polygon). Hope this is unterstandable, otherwise I will mark them... – EDi Jan 26 '15 at 16:34
• @JohnBarça I uploaded a colored picture with the desired polygons (green) and those that I don't want (red). – EDi Jan 26 '15 at 16:43

ST_StartPoint is the correct function to find single nodes at the start of a Linestring, however it does not work with MultiLinestring, so you will need to use ST_Dump to get the constituent Linestrings. If I have understood you question correctly, you then want all start points which are not also end points for more than one line, ie, points where two rivers join.

As an example, the following MultiLinestring has 3 line segments with start points that meet at the point (0, 0), and one that flows out of it.

WITH start_nodes as
(SELECT ST_StartPoint((ST_Dump(ST_Geomfromtext('MULTILINESTRING((10 10, 0 0), (5 5, 0 0),
(1 1, 0 0), (0 0, 20 20))'))).geom) as geom),
end_nodes as
(SELECT ST_EndPoint((ST_Dump(st_geomfromtext('MULTILINESTRING((10 10, 0 0), (5 5, 0 0),
(1 1, 0 0), (0 0, 20 20))'))).geom) as geom)
SELECT st_astext(sn.geom)
FROM start_nodes sn
WHERE sn.geom NOT IN (SELECT geom FROM end_nodes);

which returns

POINT(10 10)
POINT(5 5)
POINT(1 1)

excluding Point(0,0), as expected. The idea is to only include those start points that are not also end points.

In your case, the query could be written as,

WITH linestrings as
(SELECT (ST_DUMP(ST_LineMerge(geom))).geom as geom FROM spatial.stream),
start_points as
(SELECT ST_StartPoint(geom) as geom, from linestrings),
end_points as
(SELECT ST_EndPoint(geom) as geom from linestrings)
SELECT sp.geom from start_points sp
WHERE sp.geom not in (SELECT geom from end_points);

The idea it to first merge and split (st_dump) the lines making up the river, then grab the start and end points, and then finally select only those start points that are not also end points -- as this is where two, or more, rivers join.

Disclaimer: I have not tested this second query and I may have misunderstood the question or made a logical error, but I believe this is a good starting point.

• Thanks! The query runs in reasonable time and seems to give the correct results (though i did not test it extensively). – EDi Jan 28 '15 at 23:18
• Did you find a problem? – John Powell Jan 29 '15 at 7:14
• Thanks for you help! I am sorry, I wasn't concise enough with my question. I updated it. Basically I want only headwater polygons (that is) no polygon where the river passes. Your solution works perfectly. However I did not mention, that I can can starts also where the river passes. See my edits. So basically I need to add one more filter (no flow through polygon), but I couldn't come up with a fast solutiom (mine took the whole night and then there was a connection error :() – EDi Jan 29 '15 at 9:24
• I'll open a new question for this! – EDi Jan 29 '15 at 9:36
• Follow up question is here. – EDi Jan 29 '15 at 9:47

Have a look at http://blog.cleverelephant.ca/2010/07/network-walking-in-postgis.html?m=1 (Just the opposite "walking" direction)

• That's a very cool query, but answers a different question, which is building a network, not finding start points. – John Powell Jan 27 '15 at 10:35
• sure, but it gives you the complete linestring of the river so you just have to select the startpoint from this line – Thomas B Jan 27 '15 at 10:47