ST_Subdivide degrade data

[sirsosp]

I use the ST_Subdivide function very often without problem but for a dataset, I get a completely unexpected degraded result.

[sirsosp]

input data

[sirsosp]

output data

[sirsosp]

here is the SQL Query to reproduce the issue :

DROP TABLE IF EXISTS data_olivier.test_for_stsubdivide_processed;

CREATE TABLE data_olivier.test_for_stsubdivide_processed AS
SELECT  
  ST_Multi(ST_Subdivide(geom,500))::geometry(MultiPolygon , 2154) AS geom,
  cod19niv4,
  cod19niv5
FROM data_olivier.test2_for_stsubdivide
;

[strk]

I'm afraid ST_Subdivide uses intersection heavily, which is an algorithm that does not avoid creating new vertices, which in turn is prone to degradation/collapses.

Ideally we'd have an algorithm to only split at vertices, such functionality is currently completely missing from PostGIS core (C implmentation).

[sirsosp]

Did you see the screenshot st_subdivide_bug.PNG in attached file ? It seems that there is a kind of snaptogrid.

[strk]

Kind of snaptogrid is what you get when floating precision operation can't complete due to some robustness issue (lines almost parallel that intersect, usually).

Even if the problem occurs on a very small portion of the geometry, precision reduction is applied to the whole geometry, incrementally, until the operation succeeds.

[ledocteur]

Hi, Perhaps our company could help and fund this feature, just let me know how contact. Thx !

Replying to strk:

Ideally we'd have an algorithm to only split at vertices, such functionality is currently completely missing from PostGIS core (C implmentation).

[strk]

I've ticketed #4683 for a ST_SafeSubdivide

[pramsey]

The new JTS overlay code is almost complete and will be ready to port shortly. This might provide the reliable intersection required.

[strk]

Great news. Trying the OverlayNG implementation not only returns a non-degraded output but also does so in 13 seconds rather than 9 minutes !

I'm attaching the resulting subdivided data. And an image comparing old and new generation overlay outputs in a random spot.

[strk]

Comparing old and new generation GEOS Overlay outputs

[strk]

Processing output with OverlayNG

[strk]

I've enabled debugging of BinaryOp in GEOS to find out what was going on and I found that in one of the recursive intersection cases all the heuristics fail until geometry precision reduction step is entered. The precision reduction starts with a grid of 1e-16 and grows the grid by 10 everytime the intersection fails with a topology exception. The biggest grid eventually succeeding is 1, which is visible in the output as a grid.

This also explains the time difference, as it basically takes 16 attempts only in the precision reduction step, not counting earlier attempts with snapping and common-bits removal...

More stats: precision reduction is engaged 26 times. Of these, 2 times we have to get down to 1 as a grid size (I guess this is due to using the same middle-line for the binary subdivision of the same input polygon). All the other 24 cases work at the very first 1e-16 grid size.

[strk]

Interesting enough (Martin, are you reading this?) the OverlayNG code also requires (in a single case) multiple iterations over the precision reduction loop, but succeeds earlier (succeeds at 1e-10).

[strk]

The OverlayNGSnapIfNeeded instead succeeds after the first snapping, never entering precision reduction

[mdavis]

@strk: good news that OverlayNGSnapIfNeeded works first time. The new snapping heuristic works very well. How does performance compare?

What is the final resolution of this? Are you working on a new ST_Subdivide that will use OverlayNG?

[strk]

Performance is great, but it's comparing to the performance of entering the GEOS heuristic so not sure it is a valid comparison.

I'm working on exposing an ST_Subdivide taking a gridSize parameter to opt-in for the OverlayNG (it's an opt-in because of the Z drop)

[komzpa]

Moving to GEOS milestone as it is expected that 3.9 will fix the issue.

[pramsey]

GEOS 3.9 is released