I'm trying to compute some statistics for OSM data using PostgreSQL 9.3.5 and PostGIS 2.1.4. I started with a small bavaria extract which I downloaded from Geofabrik. The db schema is the normal API 0.6 schema, the data was imported via the dump approach into Postgres (using the pgsnapshot_schema_0.6*.sql scripts which come with osmosis). ANALYZE VACUUM was also performed.
The only custom made thing I'm using is a polygon table which contains multipolygons for all administrative boundary relations. The polygon geometry was not simplified in any way.
What I'm now trying to achieve is counting all nodes which are inside of the admin=6 boundaries of bavaria. Here is my SQL query:
SELECT relpoly.id, count(node)
FROM bavaria.relpolygons relpoly, bavaria.nodes node
WHERE relpoly.tags @> '"boundary"=>"administrative","admin_level"=>"6"'::hstore
AND ST_Intersects(relpoly.geom, node.geom)
GROUP BY relpoly.id;
The runtime of this query is terrible because Postgres is doing a nested loop join and scans over all nodes for every admin=6 boundary. FYI, bavaria is divided into 98 admin=6 polygons and there are about 30 million nodes in the bavaria extract.
Is it possible to avoid this sub-optimal query execution and to tell Postgres that it should scan all nodes only once (e.g., by incrementing a counter for the corresponding polygon in the result set or by using hints)?
Edit:
1) a spatial index exists on the bavaria nodes:
CREATE INDEX idx_nodes_geom ON bavaria.nodes USING gist (geom);
2) the query plan look like this:
HashAggregate (cost=284908.49..284908.75 rows=26 width=103)
-> Nested Loop (cost=111.27..283900.80 rows=201537 width=103)
-> Bitmap Heap Scan on relpolygons relpoly (cost=4.48..102.29 rows=26 width=5886)
Recheck Cond: (tags @> '"boundary"=>"administrative", "admin_level"=>"6"'::hstore)
-> Bitmap Index Scan on relpolygons_geom_tags (cost=0.00..4.47 rows=26 width=0)
Index Cond: (tags @> '"boundary"=>"administrative", "admin_level"=>"6"'::hstore)
-> Bitmap Heap Scan on nodes node (cost=106.79..10905.50 rows=983 width=127)
Recheck Cond: (relpoly.geom && geom)
Filter: _st_intersects(relpoly.geom, geom)
-> Bitmap Index Scan on idx_nodes_geom (cost=0.00..106.55 rows=2950 width=0)
Index Cond: (relpoly.geom && geom)
3)
I created the following two indexes, but the query plan (and runtime) did not change
CREATE INDEX relpolygons_tags_boundary on bavaria.relpolygons( (tags->'boundary') );
CREATE INDEX relpolygons_tags_admin on bavaria.relpolygons( (tags->'admin_level') );
ANALYZE bavaria.relpolygons;
boundary
andadmin_level
) into extra columns on the table, and use those directly.