False positive "edge crosses node"

[strk]

I was playing with hexagon grids and knowing irrational numbers approximation would have introduced slight discrepancies I thought to take a look using PostGIS topology.

To be honest I did expect problems

Well, here's the problem: ERROR: Invalid edge (no two distinct vertices exist)

To reproduce:

1. load

​https://github.com/Vizzuality/cartodb/blob/develop/lib/sql/CDB_Hexagon.sql

1. create a grid table: create table hgrid10 as select CDB_HexagonGrid(ST_MakeEnvelope(-180, -90, 180, 90), 10) as geom;
2. create a topology out of it: select createtopology('hexagon_topo'); select ST_CreateTopoGeo('hexagon_topo', st_collect(geom)) from hgrid10;

[strk]

I just verified the error occurs also with toTopoGeom

[strk]

Simpler input:

0106000000020000000103000000010000000700000000000000004065C0031AD965BE5554C00000000000E065C000000000008056C0FFFFFFFFFF1F67C001000000008056C00000000000C067C0051AD965BE5554C001000000002067C00734B2CB7C2B52C00200000000E065C00434B2CB7C2B52C000000000004065C0031AD965BE5554C00103000000010000000700000000000000006063C000000000008056C000000000000064C0FDE5269A41AA58C0FFFFFFFFFF3F65C0FEE5269A41AA58C00000000000E065C002000000008056C001000000004065C0041AD965BE5554C002000000000064C0011AD965BE5554C000000000006063C000000000008056C0

Can't use text or the error goes away

[strk]

Smaller:

01050000000200000001020000000500000000000000004065C0031AD965BE5554C00000000000E065C000000000008056C0FFFFFFFFFF1F67C001000000008056C00200000000E065C00434B2CB7C2B52C000000000004065C0031AD965BE5554C0010200000004000000FFFFFFFFFF3F65C0FEE5269A41AA58C00000000000E065C002000000008056C001000000004065C0041AD965BE5554C002000000000064C0011AD965BE5554C0

[strk]

And smaller:

0107000000020000000105000000010000000102000000030000000000000000E065C002000000008056C001000000004065C0041AD965BE5554C002000000000064C0011AD965BE5554C001050000000100000001020000000200000000000000004065C0031AD965BE5554C00000000000E065C000000000008056C0

Can't get it smaller than this.

[strk]

For the record: the CDB_HexGrid now snaps to a grid to perfectly match hexagons, so can't use the current version for testing anymore. The simplified inputs in previous comments are still good to trigger the issue

[strk]

Well, smaller (not much, just homogenized):

0105000000020000000102000000030000000000000000E065C002000000008056C001000000004065C0041AD965BE5554C002000000000064C0011AD965BE5554C001020000000200000000000000004065C0031AD965BE5554C00000000000E065C000000000008056C0

[strk]

Further analisys points to a singe line with 3 points:

01020000000300000000000000009865C04FE5D4AD958655C001000000004065C0041AD965BE5554C000000000004065C0041AD965BE5554C0

The last 2 points are _very_ close togheter.

[strk]

Ok, problem is ST_Azimuth returns NULL for these input points:

010100000000000000004065C0041AD965BE5554C0
010100000001000000004065C0041AD965BE5554C0

You can see the points are _not_ equal. ST_Distance returns 2.8421709430404e-14.

[strk]

note that azimuth is used to order edge ends around a node. There's probably a more robust way to do that (as done in GEOS/JTS) but before going there it'll be useful to remove some code duplication (the same azimuth-based code is in ST_AddEdgeModFace and ST_AddEdgeNewFace, and a similar one is in GetNodeEdges)

[strk]

Filed #1791 for the ST_Azimuth issue

[strk]

Alright so with #1791 closed ST_CreateTopoGeo succeeds at creating the topology. BUT:

1. Bad news: the topology is invalid
2. Good news: ValidateTopology correctly detects the invalidity

This is with the very first input in comment 6

[strk]

Funny: magnifying the topology (scaling nodes and edges by 1e14) makes ValidateTopology happy. Trying to "manually" validate the topology I didn't find any obvious problem.

The validation with original (non scaled) topology was reporting an edge crossing a node, but I couldn't confirm _visually_. Looking too much visually also hurts (qgis segfaulted…) so will have to recheck in some better way.

[strk]

Oh, btw, my "Bad news" above was based on being cheated by visually looking at things. After enough zoom in I saw no invalidity. Will attach a diagram of the resulting topology.

[strk]

Here's the diagram of the resulting topology, with magnified adjacencies.

In true dimensions edges 1 and 2 are almost cohincident, as are edges 3 and 4. Nodes 1 and 2 are also very close as well as nodes 4 and 5 (making edge 5 very short).

[strk]

As you can see the topology is good. So I guess the pointer goes back to pre-scaled topology and inspecting of whether ValidateTopology is mistaken!

[strk]

So, validity issue report is: ("edge crosses node",2,1) But:

=# select st_relate(e.geom, n.geom), st_crosses(e.geom, n.geom) from t0.edge e, t0.node n where e.edge_id = 2 and n.node_id = 1;
 st_relate | st_crosses 
-----------+------------
 FF1FF00F2 | f
(1 row)

[strk]

a-ah! we use ST_DWithin with tolerance 0 rather than intersection matrix… see #1625

[strk]

Using ST_Within rather than ST_DWithin fixes the false positive (r9683). The bug in the original submission was fixed by r9682.