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 ?