3

Using PostgreSQL and PostGIS, I have two tables of polygons, and I need to know which of the records in table p1 do not exist in table p2. I do not have additional information available (e.g. ids) - I need to make the determination based solely on the geometry.

A typical database solution is to use a left join

SELECT p1.geom from schema1.parcel p1
   JOIN schema2.parcel p2
    ON st_equals(p2.geom, p1.geom)
  WHERE
    p2.pid is null

This works fine for very small numbers of rows, but as the number of rows increases, it becomes very slow. Both geometry fields are indexed using gist. I have ~200K rows in a test set, and I let it run for a few hours and it didn't complete.

I noticed that using on p2.geom = p1.geom instead of st_equals performs much faster (a matter of seconds). I might be able to get away with using =, but the two are not equivalent and I would prefer to use st_equals.

Is there a generally-accepted best approach for this? Surely, it's a fairly common problem?

I'm wondering if there is a superior approach for writing such a query with st_equals? Alternatively, it has crossed my mind to do something along the lines of hashing each geometry, storing the hashes, and comparing those instead. Not sure if the expense of doing so is worth it...

Edit: it's probably pertinent to mention that both tables will (likely) have 5-10 million records, and I would expect a relatively small number of differences. Maybe 5-10%, but that's just a guess - I anticipate many matching geometries and few non-matching.

Note: I posted a related question on stackexchange that was more about the Postgres optimizer. To me, they are different enough to warrant a new question...

6
  • I'm a SQL newbie but I think my approach would be to first subset possible matches by using ST_Intersects (or whatever is appropriate for your geometries). Since they're already indexed, this should be really fast. Then look for exact matches within the candidates. At least that's how I'd approach it--I'm not good enough to quickly write an example.
    – Jon
    Nov 20, 2023 at 16:49
  • You must not use =. See gis.stackexchange.com/questions/223342/…
    – user30184
    Nov 20, 2023 at 17:29
  • 2
    @user30184 that behavior has changed with PostGIS 2.4.0 - it checks byte sequence equality now, while ~= checks for bounding box equality.
    – geozelot
    Nov 20, 2023 at 17:32
  • @user30184 Yes, = is not the same as st_equals, but one of the comments in that post points out that = is improved. I may be able to get away with the changed functionality...but I would still prefer st_equals blog.cleverelephant.ca/2017/09/postgis-operators.html Nov 20, 2023 at 17:34
  • @Jon - thanks for the suggestion. I added to my post some detail that I didn't think to include earlier. Basically, I expect that a subset using st_intersects would return nearly the entire dataset because the two tables are largely the same. So I like the idea, but I suspect (though haven't tested) that I'd face a similar issue because I'd still have a similar number of records Nov 20, 2023 at 17:39

1 Answer 1

4

The plain = implementation for the GEOMETRY/GEOGRAPHY types directly compare their encoded byte sequences, which naturally is much faster than determining geometric likeliness, especially for higher dimensional geometries like LineStrings and - even more so - Polygons.

The most idiomatic way to check for existence is an EXISTS anti-join:

SELECT
  p1.geom
FROM
  schema1.parcel AS p1
WHERE
  NOT EXISTS (
    SELECT
      1
    FROM
      schema2.parcel AS p2
    WHERE
      p1.geom && p2.geom
      AND
      ST_Equals(p1.geom, p2.geom)
  )
;

EXISTS short-circuits out of execution once the first row is added to its virtual result set, meaning that you benefit greatly from any early hits - which coincidences with your assumption of only a small fraction being non-existent¹.

I added an explicit call to the && bounding box overlap operator to coerce the planner into using an index lookup - I do not see the index being used with the actual ~= bounding box equality operator, nor ST_Equals alone in most cases, but I still suggest to at least test ~= on Polygons!
From your linked question on the motherboard it seems PostgreSQL's planner also greatly confuses the worthiness of the index on the Primary Key.

That being said, this will still need to loop over all rows in schema1.parcel and potentially apply ST_Equals on one or more index hits in schema2.parcel - but that cost is a fixed given in an anti-join, no matter the expression.


¹ this is a tricky assumption, though: a negative index lookup should be much faster to resolve than an early hit with an ST_Equals check applied on potentially complex Polygons - even more so when the index is looking for bounding box overlap. This, again, is the nature of an anti-join, so there is not much you can change about it.

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.