2

I'm trying to run a fuzzy geometry comparison on a large polygon dataset (>2 million rows), but am getting an error like so:

ERROR: ptarray_contains_point called on unclosed ring
SQL state: XX000

The problem lies with a few annoying invalid geometries being present in one of the tables, and this error is raised on the 'ST_DWithin' statement (see code excerpt below):

SELECT
    DISTINCT ON (t1.gid) t1.gid, t1.id, t2.gid, t2.id,
    (CASE WHEN ST_INTERSECTS(t1.geom, ST_MakeValid(t2.geom)) IS false THEN '9999.99'
          ELSE ST_HausdorffDistance(t1.geom, t2.geom)
    END) as hausdorff
FROM 
    (SELECT id, gid, geom FROM new_lyr WHERE gid = 12345) as t1,
    (SELECT id, gid, geom FROM old_lyr) as t2
WHERE
    ST_DWithin(t1.geom, t2.geom, 2)
ORDER BY
    t1.gid, hausdorff;

My question is: How can I somehow bypass the error, while still making use of the in-built PostGIS spatial index (i.e. the ST_DWithin statement). I have tried to put ST_MakeValid(t2.geom) into the WHERE clause, but this takes almost forever to run since it tries to make all the geometries in table t2 valid! Is there any other workaround way to resolve this, either in SQL, or would it be better just to clean up the geometry data in those tables and re-run the query?

  • 2
    Did you try ST_IsValid= true instead? – user30184 May 1 '17 at 7:31
  • in general fixing the broken geometries will always be best other wise you always have to work round them – Ian Turton May 1 '17 at 15:44
  • As @iant says, fix the broken geometries first, then everything will be fast, and the indexes will behave as expected. So, yes, to your final sentence. – John Powell May 1 '17 at 16:30
  • What about skipping those invalid geometries? How would I incorporate that into the SQL loop? – weiji14 May 1 '17 at 20:56
3

This should skip over invalid geometries.

SELECT
   DISTINCT ON (t1.gid) t1.gid, t1.id, t2.gid, t2.id,
(CASE WHEN ST_INTERSECTS(t1.geom, t2.geom) IS false THEN '9999.99'
      ELSE ST_HausdorffDistance(t1.geom, t2.geom)
END) as hausdorff
 FROM 
    (SELECT id, gid, geom FROM new_lyr WHERE gid = 12345) as t1,
   (SELECT id, gid, geom FROM old_lyr WHERE ST_IsValid(geom) ) as t2
 WHERE
      ST_DWithin(t1.geom, t2.geom, 2)
ORDER BY
    t1.gid, hausdorff;

You might want to compare teh performance with putting the ST_IsValid check in the outter where. I think both should work but the planner isn't guaranteed to run them in order.

SELECT
   DISTINCT ON (t1.gid) t1.gid, t1.id, t2.gid, t2.id,
(CASE WHEN ST_INTERSECTS(t1.geom, t2.geom) IS false THEN '9999.99'
      ELSE ST_HausdorffDistance(t1.geom, t2.geom)
END) as hausdorff
 FROM 
    (SELECT id, gid, geom FROM new_lyr WHERE gid = 12345) as t1,
   (SELECT id, gid, geom FROM old_lyr ) as t2
 WHERE
      ST_IsValid(t2.geom)  AND ST_DWithin(t1.geom, t2.geom, 2)
ORDER BY
    t1.gid, hausdorff;

Update to original answer to force the planner's hand.

Here I use ST_Expand && which is just a bounding box check so should be safe with invalid geometries as well. This will materialized the t2 table for those where bounding box is close enough to t2. Using && should be safe with even invalid geometries since invalid geometries should have a bounding box still. This may on downside produce a large worktable with false-positives.

Note the next part I use _ST_DWithin instead of ST_DWithin. The reason is that ST_DWithin is just syntactic sugar for ST_Expand(geom,2) && geom and _ST_DWithin, so no need to do the && again if you've already done it.

 WITH t1 AS (SELECT id, gid, geom FROM new_lyr WHERE gid = 12345),
      t2 AS (SELECT o.id, o.gid, o.geom FROM old_lyr AS o 
       INNER JOIN t1 ON  ST_Expand(t1.geom, 2) && o.geom
             WHERE ST_IsValid(o.geom) )
 SELECT
   DISTINCT ON (t1.gid) t1.gid, t1.id, t2.gid, t2.id,
(CASE WHEN ST_INTERSECTS(t1.geom, t2.geom) IS false THEN '9999.99'
      ELSE ST_HausdorffDistance(t1.geom, t2.geom)
END) as hausdorff
 FROM 
    t1,
    t2
 WHERE
      _ST_DWithin(t1.geom, t2.geom, 2)
ORDER BY
    t1.gid, hausdorff;
  • I am still getting the 'ptarray_contains_point called on unclosed ring' error with both solutions. I think this could be due to the ST_DWIthin function, because if I omit the outer WHERE ST_DWithin check the code runs but very very slowly. Might need to explicitly make sure the ST_IsValid code runs first before the ST_DWithin check? Not sure why your first solution fails though, considering it should be filtering out the invalid geometries already. – weiji14 May 2 '17 at 2:43
  • You would think, but PostgreSQL doesn't think that way. It rearranges order if it thinks the plan would be better. What you could do to force it is use a CTE. It might be slower, but it will force and order. I'll add as a 3rd option. – LR1234567 May 2 '17 at 3:22
  • Would it be possible to implement a spatial check that takes advantage of the spatial index, e.g. using the centroid or bounding box of the invalid geometry and checking if it is ST_DWithin a set parameter to the main t1.geom? – weiji14 May 2 '17 at 3:28
  • Look at my 3rd answer. the && does use the spatial index. ST_Centroid might not work with very invalid geometries which it seems you have here. and also you'd need a functional index on it anyway to use. – LR1234567 May 2 '17 at 3:35
  • Works like a charm now! The query itself took about 7 seconds on that one problematic test case which isn't so bad (since I have a python wrapper that only uses your advanced query when the one in the original question fails). I'll do a more thorough check on the rest of my polygons but I should be able to tweak the code if need be. Many thanks :) – weiji14 May 2 '17 at 3:51

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.