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).
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)
SELECT polyg.* FROM polyg, start_points, stream WHERE 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.