Thursday 13 October 2016

Check whether table has overlapping polygons, in PostGIS?


Using PostGIS, what is an efficient way to check whether a table with a geometry column has geometries that overlap each other?


This is my table:


  Column  |            Type             | Modifiers
---------+-----------------------------+-----------
sts | character(1) |
geom | geometry(MultiPolygon,5070) |

I guess I could see whether:



SUM(ST_Area(geom)) 

is the same as:


SUM(ST_Area(ST_Union(geom)))

But that doesn't feel like it's likely to be the most efficient way to do it, especially since my table has about 40k features.


Any ideas?


One comment on a related question suggests one could "search for overlaps by their bounding boxes" - is that a possibility?



Answer



In the situation where you only need to know whether a table contains any overlapping polygons, and you're not concerned with the size or abundance of overlaps, I recommend a query of the following form:



SELECT *
FROM my_table a
INNER JOIN my_table b ON
(a.geom && b.geom AND ST_Relate(a.geom, b.geom, '2********'))
WHERE a.ctid != b.ctid
LIMIT 1

Some pieces of this query to point out:



  • We're joining our input table to itself as a way of finding pairs of records in the table with a particular relationship.


  • We're including the && operator as part of our join condition. This operator returns true if the bounding boxes of its inputs intersect, which is a fast test that takes advantage of a spatial index.

  • The obscure part: we're also saying the two polygons involved in the join must have a relationship that satisfies the 2******** DE-9IM pattern. This is a way of checking whether the intersection of two polygons forms a polygon. (Confusingly, this is not necessarily what ST_Overlaps tells you, hence the use of the DE-9IM).

  • We prevent a row from being matched with itself by including a.ctid != b.ctid. The ctid is a unique, non-permanent row identifier that Postgres assigns to every row of every table. If your input data has a primary key, you could use a.id != b.id instead.

  • We add a LIMIT 1 at the end of the query so that Postgres can stop as soon as it's found a single overlap.


For this query to perform well, you'll need to have a spatial index on the geometry column. You can add one with CREATE INDEX ON my_table USING gist(geom) if you don't have one already.


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...