PostgreSQL: I have a PostGIS multipolygon table containing millions of features. It covers 95% of the country, however there are gaps between individual polygons. How can I fill these gaps with new polygons? Preserving and not changing the existing features.

See image below. I'd want to fill in that gap in the middle. And all the slithers and smaller polys on the right (potentially ignoring very small polys below an area of Xm2). Afterwards I can dice the shape up for better performance if it's too big.

enter image description here

note: these empty areas may very well extend all the way to the coast. I'd need to use the outline of the country to restrict this operation. I don't need new polygons appearing over the sea!


So the blue polys below are the shapes I'd like to create in this example

enter image description here

2 Answers 2


So my colleague came up with this solution. It works great for small areas, but I'd like to find a solution for millions of shapes across an entire country. I can see the ST_UNION here causing a blockage in this respect. The ST_BUFFER is just to close out thin slithers.

SELECT ST_DIFFERENCE(foo.geom, bar.geom)
FROM (SELECT ST_CONVEXHULL(ST_COLLECT(shape::geometry)) as geom FROM schema.polytable) as foo, 
(SELECT ST_BUFFER(ST_UNION(shape),0.5) as geom FROM schema.polytable) as bar


enter image description here

If anyone has suggestions for larger tables, I'm all ears.

UPDATE: I have found a solution for the entire country whereby I execute a similar process to that above, but using a gridded version of the country and iterating through each grid using ST_Intersect.

(optional) Before we start we may want to make a grid that does not extend beyond the country outline. So we'll take the entire 25x25km square grid table and a simple outline polygon of the country, then create a new table using SELECT (ST_DUMP(ST_INTERSECTION(a.geom,b.geom))).geom as geom to produce:

enter image description here

standard grid or country outline defined grid, we can then use:

SELECT ST_SUBDIVIDE(ST_DIFFERENCE(a.geom, b.geom)) as geom
(SELECT ST_BUFFER(ST_UNION(b.geom),0.5) as geom
 FROM schema.polytable b, schema.gridtable a 
 WHERE ST_INTERSECTS(b.geom,a.geom) AND a.grid_id = [use id number as a iteration variable here]) as b, schema.gridtable a 
WHERE a.grid_id = [use the same id number as a iteration variable here];

So slightly different from the previous SQL statement. There's no need for ST_CONVEXHULL this time as we're using a square grid to contain the output. Also again we use ST_BUFFER 0.5 to remove any thin inter-polygon slithers from the output. For better rendering and performance of the output, we use ST_SUBDIVIDE to divide up the resulting, and potentially huge, multipart polygon.

I need to put this in a python pipeline using the psycopg2 library, then I'll post the results on here. Testing on one grid (out of 500) takes 30secs. So could be 4 hours to run in total.

  • 1
    To scale this up/out, you could process batches of polygons using a rectangular grid defined over the data space. The constructed gap polygons can be clipped to grid cells. and optional unioned afterwards.
    – dr_jts
    Commented Jul 22, 2020 at 18:19
  • @dr_jts thanks, and this is actually what I've planned to do. Using 5kmx5km grid squares for the entire country and iterating through the grids by a grid id. This will be done in python (using psycopg2).
    – Theo F
    Commented Jul 22, 2020 at 18:30
  • 1
    Good stuff, would be interesting to hear how that works.
    – dr_jts
    Commented Jul 23, 2020 at 0:10
  • 3
    When looking to scale spatial processing queries up and/or out, gridding is usually the first technique to explore.
    – dr_jts
    Commented Jul 23, 2020 at 0:10
  • @dr_jts the pipeline is up and running. I'm using square 10x10km grids and iterating through them in a loop. I might have to reduce the grids to 5x5km as even at 10km the process can slow and even close server connection (Argh!) on more complex areas of land with many polygons.
    – Theo F
    Commented Jul 24, 2020 at 19:46

this code will fill the gaps and holes in the polygons. adjust according to your data

SELECT id, ST_Collect(ST_MakePolygon(geom)) As geom
    SELECT gid, ST_ExteriorRing((ST_Dump(geom)).geom) As geom
    FROM layer
    ) s
  • this fills in the gaps between the polygons, but it overlaps the existing polygons. I need the new fill shapes to respect the boundaries of the existing polygons. It also only fills in areas and holes of the SAME feature (like a donut poly). I'll update my OP to show an example.
    – Theo F
    Commented Jul 22, 2020 at 10:19
  • 1
    I'm not sure if I understand what you want then. Maybe post the data if you could
    – ziggy
    Commented Jul 22, 2020 at 12:48
  • hi ziggy, perhaps this image will help: i.ibb.co/7rgdy2Q/ziggy.png
    – Theo F
    Commented Jul 22, 2020 at 15:04

Your Answer

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

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