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.

 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
) tbl2
  • Using two sources seems error prone, especially if they haven't been created at the exact same time. Why not simply filtering by tags?(water, warerways,natural=water etc)
    – JGH
    Jan 21, 2018 at 16:23
  • @JGH In many cases OSM does not contain country/region areas without the coastal waters, or they are lines, not polygons. Doing the subtraction would ensure that regions don't have accidental coastal water in them. I don't think the time difference is a significant factor here. Jan 21, 2018 at 20:01

1 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:

    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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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