4

I'm trying to find islands (type not water, completely surrounded by water) in a polygon table. The polygons are not overlapping.

Data:

CREATE TABLE islands (
   ID serial PRIMARY KEY,
   type VARCHAR (10),
   geom geometry);
INSERT INTO islands(type, geom)
   VALUES
   ('water',st_geomfromewkt('SRID=3006;POLYGON((510796 6429408,513939 6420063,526040 6421096,527763 6430183,517341 6436643,510796 6429408),(516954 6429623,516652 6425747,518203 6426436,519667 6428073,516954 6429623))')),
   ('land',st_geomfromewkt('SRID=3006;POLYGON((516954 6429623,516652 6425747,518203 6426436,519667 6428073,516954 6429623))')),
   ('land',st_geomfromewkt('SRID=3006;POLYGON((510796 6429408,517341 6436643,512561 6436427,508384 6434619,510796 6429408),(511011 6434231,512475 6435652,513207 6433025,511011 6434231))')),
   ('land',st_geomfromewkt('SRID=3006;POLYGON((511011 6434231,512475 6435652,513207 6433025,511011 6434231))'));

Only id 2 should be selected: enter image description here

I have tried query below based on Determine if a polygon is not enclosed by other polygons, but it is not selecting anything.

SELECT a.id
FROM islands a 
INNER JOIN islands b ON (ST_Intersects(a.geom, b.geom) AND a.id != b.id)
where a.type = 'water' and b.type <> 'water'
GROUP BY a.id
HAVING 1e-6 > 
  abs(ST_Length(ST_ExteriorRing(a.geom)) - 
  sum(ST_Length(ST_Intersection(ST_Exteriorring(a.geom), st_interiorringn(b.geom,1)))))
3

Assuming that you have a guaranteed topological equality between island and water polygon rings, using a simpler EXISTS filter should increase performance:

SELECT  a.id
FROM    islands AS a
WHERE   a.type <> 'water'
  AND   EXISTS (
    SELECT  1
    FROM    islands AS b,
            LATERAL ST_DumpRings(b.geom) AS dmp   -- rings to polygons
    WHERE   b.type = 'water'
      AND   dmp.path[1] > 0                       -- only interior rings
      AND   ST_Equals(a.geom, dmp.geom)
  )
;
  • ST_DumpRings is relatively cheap compared to constructing geometries, and enables finding matches more generically for all inner rings of a water body polygon
  • EXISTS will quit its execution once the first result row of its query is truthy
  • no additional aggregation (GROUP BY) required

@Cyril's answer will also work great with ST_DumpPoints, i.e. using

...
tblb AS (SELECT dmp.geom FROM islands, LATERAL ST_DumpRings(geom) AS dmp WHERE type = 'water' AND dmp.path[1] > 0)
... 

If topological equality is not guaranteed, you could check for ST_Within on the exterior ring polygon:

SELECT  a.id
FROM    islands AS a
WHERE   a.type <> 'water'
  AND   EXISTS (
    SELECT  1
    FROM    islands AS b,
            LATERAL ST_DumpRings(b.geom) AS dmp  -- rings to polygon
    WHERE   b.type = 'water'
      AND   dmp.path[1] = 0                      -- only exterior ring
      AND   ST_Within(a.geom, dmp.geom)
);
|improve this answer|||||
  • 1
    Nice! (My real data was multipolygons so I had to st_dump them into polygons first) – BERA Feb 17 at 9:49
  • 1
    @BERA Cool! Btw. if you expect many inner rings to be checked, you can add a.geom && b.geom to the inner querys WHERE filter. This should give a little boost even without indexes (that are likely not being used after the dump to singletons). – geozelot Feb 17 at 20:16
2

Here's another way to solve the question asked.

Run the script:

WITH                                                                          
tbla AS (SELECT * FROM islands WHERE type = 'land'),
tblb AS (SELECT ST_InteriorRingN(geom,1) geom FROM islands WHERE type = 'water')
SELECT ID, type, (a.geom) geom FROM tbla a INNER JOIN tblb b ON ST_Intersects(a.geom, b.geom);

Check the result, and if I misunderstood the question, let me know.

|improve this answer|||||
1

The query will work only for islands placed in the first interior ring of a water body, and for topological correct geometries. Your case fills both requeriments.

The only change is that you are referring as table a as if they were the land bodies, but you are filtering a.type = 'water', and the opposite for the table b.

SELECT a.id
FROM islands a 
INNER JOIN islands b ON (ST_Intersects(a.geom, b.geom) AND a.id != b.id)
where a.type = 'land' and b.type <> 'land'
GROUP BY a.id
HAVING 1e-6 > 
  abs(ST_Length(ST_ExteriorRing(a.geom)) - 
  sum(ST_Length(ST_Intersection(ST_Exteriorring(a.geom), st_interiorringn(b.geom,1)))));

Returns:

 id
----
  2
(1 row)
|improve this answer|||||

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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