This is a follow up question to this question.

I have a river network (multiline) and some drainage polygons (see picture below). My goal is to select only the headwater polygons (green).

enter image description here

With John's solution I can easily extract the river start points (stars). However, I can have situations (red polygon) where I have startpoints in a polygon, but the polygon is not a headwater polygon, because it is flown though by the river. I only want the headwater polygons.

I tried to select them by counting the number of intersection between polygons and rivers (rationale: a headwater polygon should have only 1 intersection with the river)

    polyg, start_points, stream
    st_contains(polyg.geom, start_points.geom)
    AND ST_Npoints(ST_Intersection(poly.geom, stream.geom)) = 1

, where poylg are the poylgons, start_points from johns answer and stream is my river network.

However, this take forever and i did not run it:

"Nested Loop  (cost=0.00..20547115.26 rows=641247 width=3075)"
"  Join Filter: _st_contains(ezg.geom, start_points.geom)"
"  ->  Nested Loop  (cost=0.00..20264906.12 rows=327276 width=3075)"
"        Join Filter: (st_npoints(st_intersection(ezg.geom, rivers.geom)) = 1)"
"        ->  Seq Scan on ezg_2500km2_31467 ezg  (cost=0.00..2161.52 rows=1648 width=3075)"
"              Filter: ((st_area(geom) / 1000000::double precision) < 100::double precision)"
"        ->  Materialize  (cost=0.00..6364.77 rows=39718 width=318)"
"              ->  Seq Scan on stream_typ rivers  (cost=0.00..4498.18 rows=39718 width=318)"
"  ->  Index Scan using idx_river_starts on river_starts start_points  (cost=0.00..0.60 rows=1 width=32)"
"        Index Cond: (ezg.geom && geom)"

So my question is: How can I efficiently query headwater polygons?

Update: I added some sample data to my dropbox. Data is from south-west Germany. It's two shape files - one with streams and one with polygons.

  • So, just to be clear, you want the polygons that only contain start points, not the start points themselves. And the start points are defined as in your previous question (which I answered, and as far as I know), correctly? Jan 29 '15 at 11:56
  • Jupp, just the polygons that contain start points AND are not passed by a river / are only starts of the river. The red polygon above contains startpoints, but is NOT a headwater polygon as the river flows through it / does not start within the polygon...
    – EDi
    Jan 29 '15 at 12:16
  • So, you want the set of polygons that contain only those points that are river sources (from the previous question) and to exclude any where two rivers meet. Sorry, for all the questions, just want to be sure. Jan 29 '15 at 21:15
  • No, e.g. in the lower green polygon also two rivers meet. I want to exclude those polygons that have a river passing by (the river enters and leaves the polygon) and keep those with starts (and rivers leave only this polygon).
    – EDi
    Jan 29 '15 at 21:21
  • 1
    I don't know any PostGIS, so I can't help with direct code, however, in ArcGIS, I'd go along these lines: (1) intersect between lines and polygons into a point file. (2) delete (spatially) identical points. (3) add a numeric field to the point parameter with the value of 1 for every point. (4) spatial join the polygon onto the points and using the sum of the numeric field to indicate type of drainage. A sum of 1 means it is a headland. Higher than 1 means that there are more than 1 entrance or exit. Feb 11 '15 at 13:19

I believe the general outline (partly tested so far) is:

  1. Find the points representing stream sources, as in this answer.

  2. Intersect with the polygons table to get a count of source vertices by polygon.

  3. Use ST_DumpPoints in conjunction with group by geometry to get a count of each point. The idea being to get a count of how many rivers meet at a given point.

An example of such a query:

SELECT count(geom), ST_AsText(geom) as wkt
   (SELECT (ST_DumpPoints(foo.geom)).geom 
        ST_Collect(ST_MakeLine(ST_MakePoint(0,0), ST_MakePoint(10,10)),
                   ST_MakeLine(ST_MakePoint(0,0), ST_MakePoint(20,20))
        ) as geom
     ) foo 
 ) bar 
 GROUP BY geom; 

which returns:

count  |  wkt      
 2     | POINT(0 0)
 1     | POINT(10 10)
 1     | POINT(20 20)
  1. Run an intersects of 3 against the polygon table, to get a count (sum of vertices) of river junctions per polygon.

  2. Join the polygons from 2 on 4, rejecting those where the count (sum of vertices) of points at a junction is greater than the sum of the river sources, obtained by summing the sources by polygon from steps 1 and 2. If this condition holds, it means that at least one of the rivers meeting at a junction, originated outside of the polygon in question.

These can all be combined together in a largish sequence of CTEs, unless some tables were created from the steps involving points (and indexed).

I have no idea what the runtime of this will be on a full data set, having only tested part of this on a subset, but with a spatial index on the polygons table, there will be some assistance -- it is obviously not possible to apply an index to the points that emerge from ST_DumpPoints, so a full scan will be required there, although they should be in memory by then.

This is not being posted as a full answer, but as a work in progress, and a chance to find flaws in the logic. Working queries coming soon.


This is the query I came up with, which appears to work on a small subset of your data, but runs for hours on the full dataset.

CREATE TABLE good_polys as  
     rivers as 
       (SELECT (ST_DUMP(ST_LineMerge(geom))).geom as geom FROM streams),
     start_points as
       (SELECT ST_StartPoint(geom) as geom FROM rivers),
     end_points as 
        (SELECT ST_EndPoint(geom) as geom FROM rivers),
     junctions as 
        (SELECT (ST_DumpPoints(geom)).geom 
        FROM (SELECT geom FROM streams) s),
     source_polygons as 
            count(rivers.geom) as source_count, 
         FROM rivers, polygons
         WHERE st_intersects(polygons.geom, rivers.geom) 
         GROUP BY polygons.geom, polygons.gid),
     junction_polygons as 
            count(junctions.geom) as junction_count, 
         FROM junctions, polygons
         WHERE st_intersects(polygons.geom, junctions.geom) 
         GROUP BY polygons.geom, polygons.gid)
       junction_polygons jp, source_polygons sp 
    WHERE ST_Equals(jp.geom, sp.geom) 
    AND junction_count <= source_count;

EDIT 2. While this appears to produce correct answers on a small subset, the run time on the full dataset is horrible, presumably because the final query is doing n^2 comparisons and not using a spatial index. Probable solution would be to break query down and create tables from the initial points and point in polygon queries, which can then be spatially indexed before the final step.

  • The query is currently running on my desktop. I have no idea how long it will take or if it will be correct, though it looked reasonable from a small sample of your data. Do you have any idea how many of the polygons meet your criteria? Feb 17 '15 at 16:53
  • I will run the query on a server. I think that only a minor part of the polygons will meet the selection criteria...
    – EDi
    Feb 17 '15 at 18:08
  • That is what I found on a subset. I will post my query once it finishes Feb 17 '15 at 18:23
  • Simplification tomorrow. Feb 17 '15 at 23:41
  • Sorry, very busy today. I think the answer is to run the source query and the river junctions queries first, intersect with the polygons table to get counts per polygon, save these as tables, and then index them. Then run the final step, where the geometries are equal, and comparing the point count from the two tables. I am hoping this will then use an index rather than doing n² comparisons as present. Will post back later. Feb 18 '15 at 10:37

In pseudo code, this should work:

  • select all from polygons
  • (FULL OUTER?) join with points on polygon intersects points
  • (FULL OUTER?) join lines where polygon intersects lines
  • were line.riverid doesn't equal point.riverid
  • group by polygonid
  • count (pointid) > 0

I'm not really sure how to build the query, and I can't test it without a database to test on. It's a pretty crazy query, I think. But it should work!

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.