Removing water from OpenStreetMap polygons?

Question

I have a full OSM import using non-customized osm2pgsql (slim). I need to remove all coastal water from the polygons requested by the user on the fly.

I have imported all water from osm water polygons shapes using shp2pgsql -I -s 3857 -g way water_polygons.shp water_polygons | psql -d gis, which creates a water_polygons table with 41,000+ rows.

It seems that ST_Difference is the PostGIS function to use, but I'm not sure of the best/most efficient approach.

1. Should I merge all the water rows into one row to speed up ST_Difference, or subtract each water polygon one by one?
2. I combine several polygons into one. Should I subtract water before combining the polygon, or after?

This is my current code to extract GeoJSON polygons without water subtraction. The polygons are grouped by their Wikidata ID.

SELECT id,
 ST_AsGeoJSON(ST_Transform(ST_ForceRHR(way), 4326)) as data FROM 
(
  SELECT id, ST_Multi(ST_Union(way)) AS way
  FROM (
    SELECT tags->'wikidata' AS id, (ST_Dump(way)).geom AS way
    FROM $1~
    WHERE tags ? 'wikidata' AND tags->'wikidata' IN ($2:csv)
    ) tbl1
  GROUP BY id
) tbl2

Answer 1

In Postgres, queries are a lot clearer and straightforward when you use CTEs (common table expressions), e.g. with my_cte as (select ... from some_table) ...

I'd go this way:

with
    all_polys as (select ...),
    water_polys as (select ... where "natural"='water' or  ...)
select ... 
   st_difference(all_polys.geometry, st_union(water_polys.geometry)) 
...
from all_polys, water_polys
where st_intersects(all_polys.geometry, water_polys.geometry)
group by all_polys.id```

And then you may also turn this last select into a CTE and do more operations on it, like turning it into GeoJSON. Also, you may do operations on CTEs that have geometries and IDs only, and then join it back with one of the original tables.