TopoGeo_addLinestring corrups topology: adjacent edges 33 and -32 bind different face (1 and 0)

[Lars Aksel Opsahl]

When running this code on

select topology.CreateTopology ('rmp_temp_data_01_5526_2', 4258, 1e-06);

SELECT abs(topogeo_addlinestring) FROM topology.TopoGeo_addLinestring('rmp_temp_data_01_5526_2','0102000020A21000001D000000C0B380A8B5ED31406FD0A859974251404D4BDD9059EE314001698EEF8542514025A7A28A1CEF31404BEE5D0089425140B647DCB92AEF3140F91A58BB88425140878B46E737EF31400CF351CB8642514026A9D32B45EF3140173B6ED784425140FCFAF3556FEF31408A4851A37E425140D53A5284A9F03140005B796950425140C84AD83876F13140B4626EEB31425140F56FE7A72DF3314095AF7E06F14151409BD40C6BBEF33140A6A28A79DB4151403D887331EDF33140B77DA182D44151401FFFF38F55F4314018976FF7C441514091407B9F76F63140380F100EA641514078321D710AF7314061F3C7AB9D415140801488974BF731408CD435C49B4151405AEB4E273BF831405F985B07A14151402C2C4D6B5AF9314044AB3CC89A415140D7F2B7D786FA3140100A8988914151405500D5FADAFA314050CB49418F41514086878F6946FB3140C630D96E87415140577F478F94FB3140B203F92D804151407A9B1193FAFB3140CCA9897F774151401CB8BD1C7DFC31402C9556ED6C415140EFA8F8EC66FD314053AF581459415140A101AD3649FE31404BEC0C7146415140BD0FE29813FF314000285E3A354151400A5EA0548AFF3140ED0E8EB725415140FD6D3FA5D70532402C9DFD7B2D405140',0);
SELECT abs(topogeo_addlinestring) FROM topology.TopoGeo_addLinestring('rmp_temp_data_01_5526_2','0102000020A2100000030000007E32E78A3EEF314049272DD185425140D469DD2B45EF3140424E6ED7844251401AA1CDBA52EF31407F21C1D882425140',0);
SELECT abs(topogeo_addlinestring) FROM topology.TopoGeo_addLinestring('rmp_temp_data_01_5526_2','0102000020A2100000030000001AA1CDBA52EF31407F21C1D882425140D469DD2B45EF3140424E6ED7844251407E32E78A3EEF314049272DD185425140',0);
SELECT abs(topogeo_addlinestring) FROM topology.TopoGeo_addLinestring('rmp_temp_data_01_5526_2','0102000020A21000000300000038955FC948F3314046DA87FCEC4151408D3FDDA72DF331405E0B8106F14151400AC4A273D2F23140B16F7F7EFE415140',0);
SELECT abs(topogeo_addlinestring) FROM topology.TopoGeo_addLinestring('rmp_temp_data_01_5526_2','0102000020A2100000030000000AC4A273D2F23140B16F7F7EFE4151408D3FDDA72DF331405E0B8106F141514038955FC948F3314046DA87FCEC415140',0);
SELECT abs(topogeo_addlinestring) FROM topology.TopoGeo_addLinestring('rmp_temp_data_01_5526_2','0102000020A210000003000000829FE21B0BF83140ED302FF99F415140015C87974BF73140FE133AC49B415140D6F13F5D1FF73140FCA4330F9D415140',0);
SELECT abs(topogeo_addlinestring) FROM topology.TopoGeo_addLinestring('rmp_temp_data_01_5526_2','0102000020A210000003000000D6F13F5D1FF73140FCA4330F9D415140015C87974BF73140FE133AC49B415140829FE21B0BF83140ED302FF99F415140',0);
SELECT abs(topogeo_addlinestring) FROM topology.TopoGeo_addLinestring('rmp_temp_data_01_5526_2','0102000020A2100000030000000A1BAAA8B5ED31406FD0A85997425140653CEF9059EE3140AA8591EF85425140CBA0469CF6EE3140F40FB06788425140',0);
SELECT abs(topogeo_addlinestring) FROM topology.TopoGeo_addLinestring('rmp_temp_data_01_5526_2','0102000020A210000003000000CBA0469CF6EE3140F40FB06788425140653CEF9059EE3140AA8591EF854251400A1BAAA8B5ED31406FD0A85997425140',0);
SELECT abs(topogeo_addlinestring) FROM topology.TopoGeo_addLinestring('rmp_temp_data_01_5526_2','0102000020A2100000030000006406657B10F13140AA28181341425140B3964C84A9F03140EC847D69504251406A5784D908F03140E6C41C0D68425140',0);
SELECT abs(topogeo_addlinestring) FROM topology.TopoGeo_addLinestring('rmp_temp_data_01_5526_2','0102000020A2100000030000006A5784D908F03140E6C41C0D68425140B3964C84A9F03140EC847D69504251406406657B10F13140AA28181341425140',0);
SELECT abs(topogeo_addlinestring) FROM topology.TopoGeo_addLinestring('rmp_temp_data_01_5526_2','0102000020A2100000040000002279EF8639F431409BB84F24C94151400DAF039055F43140B0F46EF7C4415140CDE2729F76F6314070FE140EA6415140AB506A8DDEF63140AF9D0A29A0415140',0);
SELECT abs(topogeo_addlinestring) FROM topology.TopoGeo_addLinestring('rmp_temp_data_01_5526_2','0102000020A210000004000000AB506A8DDEF63140AF9D0A29A0415140CDE2729F76F6314070FE140EA64151400DAF039055F43140B0F46EF7C44151402279EF8639F431409BB84F24C9415140',0);
SELECT abs(topogeo_addlinestring) FROM topology.TopoGeo_addLinestring('rmp_temp_data_01_5526_2','0102000020A21000000900000048A235EF1FFE314050AB6CD749415140B751F8EC66FD3140416B5C14594151402F8CBA1C7DFC3140BA4357ED6C41514066F00793FAFB3140103A8C7F774151405EF2478F94FB31409ECEFB2D804151400F18A06946FB3140CD0CD86E87415140C327E7FADAFA314022DF49418F4151408F57A4D786FA3140DF688C8891415140D45198BF82F93140F4B5648A99415140',0);
SELECT abs(topogeo_addlinestring) FROM topology.TopoGeo_addLinestring('rmp_temp_data_01_5526_2','0102000020A210000009000000D45198BF82F93140F4B5648A994151408F57A4D786FA3140DF688C8891415140C327E7FADAFA314022DF49418F4151400F18A06946FB3140CD0CD86E874151405EF2478F94FB31409ECEFB2D8041514066F00793FAFB3140103A8C7F774151402F8CBA1C7DFC3140BA4357ED6C415140B751F8EC66FD3140416B5C145941514048A235EF1FFE314050AB6CD749415140',0);
SELECT abs(topogeo_addlinestring) FROM topology.TopoGeo_addLinestring('rmp_temp_data_01_5526_2','0102000020A21000000400000063ABC5BD19023240DDE24BDFC04051408173B1548AFF314037B38DB7254151401D5AE89813FF3140478F623A3541514025ED461466FE3140EA6988FC43415140',0);
SELECT abs(topogeo_addlinestring) FROM topology.TopoGeo_addLinestring('rmp_temp_data_01_5526_2','0102000020A21000000400000025ED461466FE3140EA6988FC434151401D5AE89813FF3140478F623A354151408173B1548AFF314037B38DB72541514063ABC5BD19023240DDE24BDFC0405140',0);
SELECT abs(topogeo_addlinestring) FROM topology.TopoGeo_addLinestring('rmp_temp_data_01_5526_2','0102000020A21000000600000030A2B312C1013240A5B7CAE695425140A349BD27BC0132407B0E54EC954251406EC47A65BC0132400FC3FC5C9642514002BAF5E0BC013240C8A5633E974251408E12ECCBC1013240F14EDA38974251409437EBDDC10132406FD0A85997425140',0);
SELECT abs(topogeo_addlinestring) FROM topology.TopoGeo_addLinestring('rmp_temp_data_01_5526_2','0102000020A2100000090000009437EBDDC10132406FD0A859974251408E12ECCBC1013240F14EDA389742514002BAF5E0BC013240C8A5633E97425140383F38A3BC0132406C34B0CD96425140A349BD27BC0132407B0E54EC9542514030A2B312C1013240A5B7CAE6954251406527F6D4C00132404846177695425140B35F5E42E3013240FB5F6B4F954251407D7EF737E30132402E46963C95425140',0);
SELECT abs(topogeo_addlinestring) FROM topology.TopoGeo_addLinestring('rmp_temp_data_01_5526_2','0102000020A21000000A0000007D7EF737E30132402E46963C95425140A6D2F4A3D80132409C40294E95425140E7C82D7FCE013240A5E8784495425140C97C40C5C3013240A5D5B85D954251405331E7A1B501324046381B7295425140AB4F492DAC013240728C84C3954251403C56C637A8013240109CCF9D95425140CE87F3D7AE013240E001D78596425140AC72676CB601324050C6F756974251400F7E7782B60132406FD0A85997425140',0);
SELECT abs(topogeo_addlinestring) FROM topology.TopoGeo_addLinestring('rmp_temp_data_01_5526_2','0102000020A21000002A000000855F31A8D00132406FD0A859974251406F02258DD0013240360749289742514050F25D4EDF013240B302AD179742514013EF6B0FEE013240F7BA1B0797425140499383E5F70132404B0D09FC964251409D2AC6BB01023240661C01F196425140B5BCDD7D0102324009AB4D8096425140EC60F5530B0232405DFD3A759642514022E637160B023240018C870496425140CC31590110023240F2F108FF95425140588A4FEC140232401C9B7FF9954251408E0F92AE14023240C029CC88954251401A68889919023240E9D242839542514050EDCA5B190232408D618F1295425140FA38EC461E023240B70A060D95425140A15B221B1E02324025BC34BD9442514012E9F17F1D0232400D3994B694425140E35431DB1C0232404932AA9D94425140D9527AAA19023240D22F5A2294425140B1E9BD2D15023240AC807253934251400566741B15023240203F2250934251400B671AE61402324020AABD5393425140F909F38513023240F9BAB06B934251409D0F19C50E023240BD7151BE93425140F6AA63F40D02324011B17ECC934251405E086A040A02324049E5FC1094425140E2EEAC0207023240C9744A45944251404E9D0D3F06023240C0379652944251408490873E0502324019E74D5A944251400FFC0365000232409627CD7F94425140EE9736DEFB0132405529C8A294425140CCC97A45F80132408D9ACFC69442514073F631B8F6013240962F5FD6944251402E92B8E4F1013240DB6DBB06954251408993A8FAEF013240AB2AF11995425140B6D8921FE8013240430FCD3095425140D0D6873BE5013240DA473B39954251407D7EF737E30132402E46963C95425140B35F5E42E3013240FB5F6B4F954251407CBB466CD9013240E050735A954251406527F6D4C0013240484617769542514030A2B312C1013240A5B7CAE695425140',0);

is tested in this this versions

POSTGIS="3.5.0dev 3.4.0rc1-1272-g80806376a" [EXTENSION] PGSQL="160" GEOS="3.12.1-CAPI-1.18.1" PROJ="9.4.0 NETWORK_ENABLED=OFF URL_ENDPOINT=https://cdn.proj.org USER_WRITABLE_DIRECTORY=/Users/lop/Library/Application Support/proj DATABASE_PATH=/opt/homebrew/Cellar/proj/9.4.0/share/proj/proj.db" (compiled against PROJ 9.12.1) LIBXML="2.13.0" LIBJSON="0.17" TOPOLOGY

and also on

POSTGIS="3.5.0dev 3.4.0rc1-1310-g8349e1670" [EXTENSION] PGSQL="160" GEOS="3.13.0-CAPI-1.19.0" PROJ="9.3.0 NETWORK_ENABLED=OFF URL_ENDPOINT=https://cdn.proj.org USER_WRITABLE_DIRECTORY=/tmp/proj DATABASE_PATH=/usr/local/share/proj/proj.db" (compiled against PROJ 9.13.0) LIBXML="2.9.13" LIBJSON="0.15" LIBPROTOBUF="1.3.3" WAGYU="0.5.0 (Internal)" TOPOLOGY

and

POSTGIS="3.4.2 c19ce56" [EXTENSION] PGSQL="160" GEOS="3.12.1-CAPI-1.18.1" (compiled against GEOS 3.10.2) SFCGAL="SFCGAL 1.4.1, CGAL 5.3.1, BOOST 1.74.0" PROJ="8.2.1 NETWORK_ENABLED=OFF URL_ENDPOINT=https://cdn.proj.org USER_WRITABLE_DIRECTORY=/tmp/proj DATABASE_PATH=/usr/share/proj/proj.db" GDAL="GDAL 3.4.1, released 2021/12/27" LIBXML="2.9.13" LIBJSON="0.15" LIBPROTOBUF="1.3.3" WAGYU="0.5.0 (Internal)" TOPOLOGY RASTER

it fails on all

[Lars Aksel Opsahl]

With higher tolerance like this it works for this case.

select topology.CreateTopology ('rmp_temp_data_01_5526_2', 4258, 1e-05);

[Lars Aksel Opsahl]

When I test with 1e-05 I now see this this errors, there error below are also present with tolerance 1e-05

XX000 message: SQL/MM Spatial exception - geometry crosses edge 
XX000 message: SQL/MM Spatial exception - end node not geometry end point.

I try create testcase for this errors.

[Lars Aksel Opsahl]

It also fails on this

PostgreSQL 16.4 (Debian 16.4-1.pgdg110+1) on x86_64-pc-linux-gnu, compiled by gcc (Debian 10.2.1-6) 10.2.1 20210110, 64-bit
  Postgis 3.5.0dev - (3.5.0alpha2-55-g8349e1670) - 2024-09-23 05:55:28
  scripts 3.5.0dev 3.5.0alpha2-55-g8349e1670
  GEOS: 3.14.0dev-CAPI-1.20.0
  PROJ: 9.5.0 NETWORK_ENABLED=OFF URL_ENDPOINT=https://cdn.proj.org USER_WRITABLE_DIRECTORY=/var/lib/postgresql/.local/share/proj DATABASE_PATH=/usr/local/share/proj/proj.db

pipeline here

​https://gitlab.com/nibioopensource/resolve-overlap-and-gap/-/jobs/7900975113

to run the test

time sh src/test/sql/regress/run_resolve_test.sh --nodrop  ./src/test/sql/regress/topology_test_01.sql;  

this test does not involve any code from ​https://gitlab.com/nibioopensource/resolve-overlap-and-gap

[strk]

I can reproduce locally with POSTGIS="3.5.0dev 3.4.0rc1-1231-g932e4936c" GEOS="3.13.0beta2-CAPI-1.19.0".

What happens is that one of the TopoGeo_addLinestring calls (the one before the last in the testcase) makes the topology invalid. I'm working on reducing the testcase.

[strk]

I can also reproduce with POSTGIS="3.4.3dev 3.4.2-39-gd9f3bf593" and POSTGIS="3.3.8dev 3.3.7-1-ge035797e5" (same GEOS version)

[strk]

Simplified testcases:

select topology.DropTopology ('t5782');                                                                                                                                                        
                                                                                                                                                                                               
select topology.CreateTopology ('t5782');                                                                                                                                                      
                                                                                                                                                                                               
SELECT NULL FROM topology.TopoGeo_addLinestring('t5782',                                                                                                                                       
'LINESTRING(                                                                                                                                                                                   
18.006852310996862 69.0403992633497,                                                                                                                                                           
18.00677727099686 69.0404005833497,                                                                                                                                                            
18.006780950996863 69.04042744334969,                                                                                                                                                          
18.00678831099686 69.0404811833497,                                                                                                                                                            
18.00686335099686 69.04047986334969,                                                                                                                                                           
18.006864423681307 69.04048768506776)'                                                                                                                                                         
);                                                                                                                                                                                             
                                                                                                                                                                                               
SELECT NULL FROM topology.TopoGeo_addLinestring('t5782',                                                                                                                                       
'LINESTRING(18.006864423681307 69.04048768506776,                                                                                                                                              
18.00686335099686 69.04047986334969,                                                                                                                                                           
18.00678831099686 69.0404811833497,                                                                                                                                                            
18.00678463099686 69.0404543133497,                                                                                                                                                            
18.00677727099686 69.0404005833497,                                                                                                                                                            
18.006852310996862 69.0403992633497,                                                                                                                                                           
18.00684863099686 69.04037239334968,                                                                                                                                                           
18.00737395099686 69.04036317334969,                                                                                                                                                           
18.00737333099686 69.0403586833497)'                                                                                                                                                           
);                                                                                                                                                                                             
                                                                                                                                                                                               
SELECT * FROM topology.ValidateTopology('t5782');                                                                                                                                              
                                                                                                                                                                                               
BEGIN;                                                                                                                                                                                         
                                                                                                                                                                                               
-- Breaks validity                                                                                                                                                                             
SELECT NULL FROM topology.TopoGeo_addLinestring('t5782',                                                                                                                                       
'LINESTRING(18.00737333099686 69.0403586833497,                                                                                                                                                
18.00721192099686 69.0403628733497,                                                                                                                                                            
18.00705714099686 69.0403605633497,                                                                                                                                                            
18.00689347099686 69.0403665833497,                                                                                                                                                            
18.00667774099686 69.0403714433497,                                                                                                                                                            
18.00653346099686 69.04039085334969,                                                                                                                                                           
18.00647305099686 69.0403818633497,                                                                                                                                                            
18.00657415099686 69.0404371833497,                                                                                                                                                            
18.00668981099686 69.04048704334969,                                                                                                                                                           
18.006691126034692 69.04048768506776)'                                                                                                                                                         
);                                                                                                                                                                                             
                                                                                                                                                                                               
SELECT * FROM topology.ValidateTopology('t5782');

[strk]

The starting topology has two edges that are almost coincident (3 and 4, in the picture):

It looks like those two edges influence the code as if I change them to be further apart, as in the next picture, the subsequent addition of a line does not result in a broken topology:

The last incoming line (in yellow in the picture) only touches edge 5 and is otherwise open, not closing any ring. For some reason the code thinks one of the side faces of it is face 1 which is the tiny face between edges 3 and 4. It must be the code that determines which face contains a new point (endpoinnt of incoming line) by lookinng at the closest edge: in this case the "closest edge" is probably 3 or 4 and we get the side wrong. I'll verify by testing addition of a single point.

[strk]

Confirmed, adding the last point of incoming line finds edge 3 to be the closest to it while instead edge 4 would be closest:

_lwt_AddPoint:6559] Adding point: POINT(18.006691126034692 69.04048768506776)
lwt_GetFaceContainingPoint:7502] Closest segment on edge 3 is 1 (dist 9.71684e-05)
lwt_GetFaceContainingPoint:7505] Closest segment on edge 3 is LINESTRING(18.0068 69.0404, 18.0068 69.0405)
lwt_GetFaceContainingPoint:7642] Closest point is NOT a node
lwt_GetFaceContainingPoint:7760] Closest point is internal to closest segment, calling lw_segment_side((18.0068,69.0404),(18.0068,69.0405),(18.0067,69.0405)
lwt_GetFaceContainingPoint:7770] Side of closest segment query point falls on: -1 (left)                

I guess the code could be made more robust to check that the first and last node of a line being added are found to be on the same face or throw an exception if that's not the case

[Lars Aksel Opsahl]

Nice, I have added new test in same branch 84-error-xx000-corrupted-topology-adjacent-edges-33-and-32-bind-different-face-1-and- on gitlab

ERROR:  Edge changed disposition around start node 15

It'stest src/test/sql/regress/topology_test_02.sql

[strk]

This query confirms edge 3 is found closer than edge 4:

=# select edge_id, st_distance(geom, 'POINT(18.006691126034692 69.04048768506776)') from t5782.edge order by geom <-> 'POINT(18.006691126034692 69.04048768506776)' limit 3;
 edge_id |      st_distance      
---------+-----------------------
       3 | 9.716835623092273e-05
       4 | 9.716835623285099e-05
       2 | 9.740220330897892e-05
(3 rows)

We don't even get an equidistance. On the other hand order by angle around nodes confirms edge 4 is on the left of edge 3. Puzzling.

[strk]

For the record: adding ONLY the point:

select topogeo_addPoint('t5782', 'POINT(18.006691126034692 69.04048768506776)');

Sets containing_face of that point to 1, which is wrong, but ValidateTopology will not notice because it will use the same algorithm used during insertion to verify the face containing the point, so will find again edge 3 being the closest edge it and thus consider the containment valid. Maybe ValidateTopology should use a different method to verify containment…

[strk]

Of course the same problem occurs with GetFaceContainingPoint:

=# select GetFaceContainingPoint('t5782', 'POINT(18.006691126034692 69.04048768506776)');
 getfacecontainingpoint 
------------------------
                      1
(1 row)

Both edges are found to have no interior intersection with the shortest line between the query point and edge 3:

=# with candidate as ( select ST_ShortestLine('POINT(18.006691126034692 69.04048768506776)', geom) sl from t5782.edge where edge_id = 3) select edge_id, ST_Relate(geom, sl) from t5782.edge, candidate where ST_Intersects(geom, sl)
;
 edge_id | st_relate 
---------+-----------
       4 | F01FF0102
       3 | F01FF0102
(2 rows)

And both edges are found to have an interior intersection with the shortest line between query point and edge 4:

=# with candidate as ( select ST_ShortestLine('POINT(18.006691126034692 69.04048768506776)', geom) sl from t5782.edge where edge_id = 4) select edge_id, ST_Relate(geom, sl) from t5782.edge, candidate where ST_Intersects(geom, sl)
;
 edge_id | st_relate 
---------+-----------
       4 | 0F1FF0102
       3 | 0F1FF0102
(2 rows)

I don't feel I can trust any of the above

Here's a picture of the shortest line between query point and any of the two edges:

[Lars Aksel Opsahl]

Replying to strk:

Of course the same problem occurs with GetFaceContainingPoint:

=# select GetFaceContainingPoint('t5782', 'POINT(18.006691126034692 69.04048768506776)');
 getfacecontainingpoint 
------------------------
                      1
(1 row)

Both edges are found to have no interior intersection with the shortest line between the query point and edge 3:

=# with candidate as ( select ST_ShortestLine('POINT(18.006691126034692 69.04048768506776)', geom) sl from t5782.edge where edge_id = 3) select edge_id, ST_Relate(geom, sl) from t5782.edge, candidate where ST_Intersects(geom, sl)
;
 edge_id | st_relate 
---------+-----------
       4 | F01FF0102
       3 | F01FF0102
(2 rows)

And both edges are found to have an interior intersection with the shortest line between query point and edge 4:

=# with candidate as ( select ST_ShortestLine('POINT(18.006691126034692 69.04048768506776)', geom) sl from t5782.edge where edge_id = 4) select edge_id, ST_Relate(geom, sl) from t5782.edge, candidate where ST_Intersects(geom, sl)
;
 edge_id | st_relate 
---------+-----------
       4 | 0F1FF0102
       3 | 0F1FF0102
(2 rows)

I don't feel I can trust any of the above

Can it be related to performance or other fixes in GEOS ?

Does it for instance fail on GEOS from 2023 ?

[strk]

Another small testcase focusing on how topology gets the containing face of a newly added point wrong:

SELECT topology.DropTopology ('t5782');                                                                                                                                                        
SELECT topology.CreateTopology ('t5782');                                                                                                                                                      -- Edge 1 (was: 3)                                                                                                                                                                             
SELECT NULL FROM topology.TopoGeo_addLinestring('t5782', 'LINESTRING(18.00677727099686 69.0404005833497,18.006780950996863 69.04042744334969,18.00678831099686 69.0404811833497)');            
-- Edge 2 (was: 4)                                                                                                                                                                             
SELECT NULL FROM topology.TopoGeo_addLinestring('t5782', 'LINESTRING(18.00678831099686 69.0404811833497,18.006784630996860 69.04045431334970,18.00677727099686 69.0404005833497)');            
-- Finds point to be inside the tiny face (wrong!)                                                                                                                                                     
SELECT NULL FROM topology.TopoGeo_addPoint('t5782', 'POINT(18.006691126034692 69.04048768506776)');                                                                                            
SELECT node_id, containing_face FROM t5782.node ORDER BY node_id;  

Note that ValidateTopology won't catch the problem, but the last added isolated node gets containing_face to be 1 instead of 0.

[strk]

What I can tell is that GEOS finds the correct edge being closer. Here's a simple test:

A='LINESTRING(18.00678831099686 69.0404811833497,18.006784630996860 69.04045431334970,18.00677727099686 69.0404005833497)'
B='LINESTRING(18.00677727099686 69.0404005833497,18.006780950996863 69.04042744334969,18.00678831099686 69.0404811833497)'
Q='POINT(18.006691126034692 69.04048768506776)'

AQ_geos=$(geosop -a "$A" -b "$Q" distance)
AQ_pgis=$(psql -XtAc "SELECT ST_Distance('$A', '$Q')")
BQ_geos=$(geosop -a "$B" -b "$Q" distance)
BQ_pgis=$(psql -XtAc "SELECT ST_Distance('$B', '$Q')")

echo "PostGIS:"
echo " A-Q: $AQ_pgis"
echo " B-Q: $BQ_pgis"

echo "GEOS:"
echo " A-Q: $AQ_geos"
echo " B-Q: $BQ_geos"

The results I get:

PostGIS:
 A-Q: 9.716835623285099e-05
 B-Q: 9.716835623092273e-05
GEOS:
 A-Q: 9.7168356230588148e-05
 B-Q: 9.7168356230827946e-05

Basically PostGIS finds B being closer to Q while GEOS finds A being closer to Q. In our case the GEOS answer would be more appropriate.

What we are using to find the closest edge is the <-> operator.

[Lars Aksel Opsahl]

So if ST_Distance in postgis had used the same GEOS call that you do it would have worked ? I assume there are some reasons why ST_Distance is not doing this today.

[mdavis]

It seems problematic to rely on distance calculation to determine topology, since it's inherently less precise than the actual geometric situation.

Perhaps the fundamental problem is treating nearly-identical linestrings as distinct edges. Would it be possible to specify a distance tolerance, and raise an error if lines contain distinct sections which are closer than the tolerance?

[Lars Aksel Opsahl]

Replying to mdavis:

It seems problematic to rely on distance calculation to determine topology, since it's inherently less precise than the actual geometric situation.

Perhaps the fundamental problem is treating nearly-identical linestrings as distinct edges. Would it be possible to specify a distance tolerance, and raise an error if lines contain distinct sections which are closer than the tolerance?

What you say about tolerance think makes since, because I had big area missing and it did help some to make cell, but al it did not fix all problems and I did not expect that .

But then I did only one change switching from

select topology.CreateTopology ('toplogy_name', 4258, 1e-06);

to

select topology.CreateTopology ('toplogy_name', 4258, 1e-05);

and all area seems to recovered for that part.

[Lars Aksel Opsahl]

Replying to Lars Aksel Opsahl:

Replying to mdavis:

It seems problematic to rely on distance calculation to determine topology, since it's inherently less precise than the actual geometric situation.

Perhaps the fundamental problem is treating nearly-identical linestrings as distinct edges. Would it be possible to specify a distance tolerance, and raise an error if lines contain distinct sections which are closer than the tolerance?

What you say about tolerance think makes since, because I had big area missing and it did help some to make cell, but al it did not fix all problems and I did not expect that .

But then I did only one change switching from

select topology.CreateTopology ('toplogy_name', 4258, 1e-06);

to

select topology.CreateTopology ('toplogy_name', 4258, 1e-05);

and all area seems to recovered for that part.

But for resolve-overlap-and-gap/src/test/sql/regress/topology_test_02.sql i god up to 1e-04

[strk]

Replying to mdavis:

It seems problematic to rely on distance calculation to determine topology, since it's inherently less precise than the actual geometric situation.

The computation that goes wrong is finding which face contains a point. This used to be implemented by constructing the geometry of each face whose MBR covered the point and then computing PIP, but was later converted to finding the closest edge and finding out on which side of the edge the point was. I'm not sure which method would be more robust.

Perhaps the fundamental problem is treating nearly-identical linestrings as distinct edges. Would it be possible to specify a distance tolerance, and raise an error if lines contain distinct sections which are closer than the tolerance?

Tolerance can be already specified, and when it's not specified it is computed based on the extent of the input geometry.

We currently only use the tolerance to snap, which doesn't help in this case (it helps if you swap the first to insertions).

Distinct sections closer than tolerance would include incident edge ends which will always be problematic if the angle between them is too tiny.

[strk]

Replying to Lars Aksel Opsahl:

So if ST_Distance in postgis had used the same GEOS call that you do it would have worked ?

In this specific case YES, I've tried it and if ST_Distance used GEOSDistance the test would pass. I guess the reason to use liblwgeom would be avoidinng the conversion to GEOS but if it gives wrong answers we could switch over.

What Martin says about distance being not robust enough remains as a problem though…

[strk]

A pull request that implements ST_Distance via GEOS is here: ​https://git.osgeo.org/gitea/postgis/postgis/pulls/218

[Lars Aksel Opsahl]

Replying to strk:

A pull request that implements ST_Distance via GEOS is here: ​https://git.osgeo.org/gitea/postgis/postgis/pulls/218

Thanks then I will take this branch and use that.

[Lars Aksel Opsahl]

Replying to Lars Aksel Opsahl:

Replying to strk:

A pull request that implements ST_Distance via GEOS is here: ​https://git.osgeo.org/gitea/postgis/postgis/pulls/218

Thanks then I will take this branch and use that.

I still see a lot errors , but confirms that topology_test_01.sql ​https://gitlab.com/nibioopensource/resolve-overlap-and-gap/-/blob/84-error-xx000-corrupted-topology-adjacent-edges-33-and-32-bind-different-face-1-and-0/src/test/sql/regress/topology_test_01.sql?ref_type=heads is running ok.

and that topology_test_02.sql ​https://gitlab.com/nibioopensource/resolve-overlap-and-gap/-/blob/84-error-xx000-corrupted-topology-adjacent-edges-33-and-32-bind-different-face-1-and-0/src/test/sql/regress/topology_test_02.sql?ref_type=heads runs ok.

But we still have lot errors in the log and area lost

POSTGIS="3.5.0dev 3.5.0rc1-5-gb9d4311f9" [EXTENSION] PGSQL="160" GEOS="3.13.0-CAPI-1.19.0" PROJ="9.3.0 NETWORK_ENABLED=OFF URL_ENDPOINT=https://cdn.proj.org USER_WRITABLE_DIRECTORY=/tmp/proj DATABASE_PATH=/usr/local/share/proj/proj.db" (compiled against PROJ 9.13.0) LIBXML="2.9.13" LIBJSON="0.15" LIBPROTOBUF="1.3.3" WAGYU="0.5.0 (Internal)" TOPOLOGY

[Lars Aksel Opsahl]

POSTGIS="3.5.0dev 3.5.0rc1-5-gb9d4311f9" [EXTENSION] PGSQL="160" GEOS="3.13.0-CAPI-1.19.0" PROJ="9.3.0 NETWORK_ENABLED=OFF URL_ENDPOINT=https://cdn.proj.org USER_WRITABLE_DIRECTORY=/tmp/proj DATABASE_PATH=/usr/local/share/proj/proj.db" (compiled against PROJ 9.13.0) LIBXML="2.9.13" LIBJSON="0.15" LIBPROTOBUF="1.3.3" WAGYU="0.5.0 (Internal)" TOPOLOGY

With the code abouve this is the first error I always see

XX000 message: Programmatic error ? Geometry of face 38716

Then I end with errorss like this

  Corrupted topology: face 559 could not be constructed only 

The problem seems to frequent running higher tolerance, I am now using 1e-04.

[strk]

PR going back to build full face geometry in order to determine containment is here: ​https://git.osgeo.org/gitea/postgis/postgis/pulls/219

[strk]

Another pull to make ValidateTopology catch the wrong containing_face: ​https://git.osgeo.org/gitea/postgis/postgis/pulls/220

[strk]

For the record, SFCGAL agrees with GEOS about which line is closer, although gives different numbers (using CG_Distance):

PostGIS:
 A-Q: 9.716835623285099e-05
 B-Q: 9.716835623092273e-05
SFCGAL:
 A-Q: 9.716835622921602e-05
 B-Q: 9.716835623197789e-05
GEOS:
 A-Q: 9.7168356230588148e-05
 B-Q: 9.7168356230827946e-05

[strk]

I confirm PostGIS minimum distance calculation is currently constructing the closest point first and then computing the distance from that point, which is what is causing this robustness issue. PR ​https://git.osgeo.org/gitea/postgis/postgis/pulls/222 is changing that and fixes the issue at hand, giving the same answers as GEOS:

PostGIS:
 A-Q: 9.716835623058815e-05
 B-Q: 9.716835623082795e-05
SFCGAL:
 A-Q: 9.716835622921602e-05
 B-Q: 9.716835623197789e-05
GEOS:
 A-Q: 9.7168356230588148e-05
 B-Q: 9.7168356230827946e-05

[Lars Aksel Opsahl]

Replying to strk:

I confirm PostGIS minimum distance calculation is currently constructing the closest point first and then computing the distance from that point, which is what is causing this robustness issue. PR ​https://git.osgeo.org/gitea/postgis/postgis/pulls/222 is changing that and fixes the issue at hand, giving the same answers as GEOS:

Nice, tested now your branch mindistance-no-construct-point and confirms that test /topology_test_01.sql works with no errors running on

POSTGIS="3.6.0dev 3.5.0-9-g181859475" [EXTENSION] PGSQL="160" GEOS="3.13.0-CAPI-1.19.0" PROJ="8.2.1 NETWORK_ENABLED=OFF URL_ENDPOINT=https://cdn.proj.org USER_WRITABLE_DIRECTORY=/tmp/proj DATABASE_PATH=/usr/share/proj/proj.db" (compiled against PROJ 8.13.0) LIBXML="2.9.13" LIBJSON="0.13.1" LIBPROTOBUF="1.3.3" WAGYU="0.5.0 (Internal)" TOPOLOGY

For test topology_test_02.sql ​​https://gitlab.com/nibioopensource/resolve-overlap-and-gap/-/blob/84-error-xx000-corrupted-topology-adjacent-edges-33-and-32-bind-different-face-1-and-0/src/test/sql/regress/topology_test_02.sql?ref_type=heads I get this error in the same database as above .

The error I get is

ERROR:  XX000: Side-location conflict: new edge starts in face 37 and ends in face 10

Is that expected ?

[strk]

The mindistance-no-construct-point branch was now merged into the master branch with commit [290fffd11dac294431056e6d152098ac0afe60de/git] with regress tests.

I did not analyze the topology_test_02.sql test yet but if as you report it was succeeding with the GEOS-implemented min distance and it is failing with the PostGIS version of it then it is not expected I guess.

[strk]

I would file a separate ticket for your test02

[Lars Aksel Opsahl]

Replying to strk:

I would file a separate ticket for your test02

Ok I will do.

[strk]

Postgis 3.2.8dev is also affected by this bug

[Sandro Santilli <strk@…>]

In 0758764/git:

Improve robustness of min distance calculation

References #5782 in 3.4 branch (3.4.4dev)
Includes regress tests

[Sandro Santilli <strk@…>]

In bb7665db/git:

Improve robustness of min distance calculation

References #5782 in 3.3 branch (3.3.8dev)
Includes regress tests

[strk]

I'll call this closed in 3.3 and avoid backporting to 3.2

Now to look at #5786

[Sandro Santilli <strk@…>]

In bdecd88/git:

Improve robustness of min distance calculation

References #5782 in 3.5 branch (3.5.1dev)
Includes regress tests

[lnicola]

Drive-by comment: hypot() might be more robust than sqrt(x * x + y * y).