rings not closed after update

[larsop]

We started out with this map.

Then we added a line and we get this map

Here we se that line 219 is "new" data and we got from the web client.

Line 221, 220, 218, 222 are data that have change/added that we also have to send back to NGIS OpenAPi.

And this lines are causing us the problem described in this issues. In this issue you also find files that we get from NGIS-OpenAPI.

Then problem is where line 31, 117 and 211 meets in single point.

Line 31 and 211 are part off ring that that should be closed which it not is because of this 2 points (654032.5399999999, 7812386.27) and (654032.54,7812386.27) that we get back from NGIS-OpenApi. This is a difference that is not possible to se on a any map but it forces our server side to do some kind off snapping.

I had expected line 221 had snapped to the common point for line 31 and 117 that I assume already exists in the database.

[meastp]

We have removed snapping in QMS 13. What tolerances are you using? Are you expecting floating point numbers to be exactly equal at the binary level?

[larsop]

We have removed snapping in QMS 13. What tolerances are you using? Are you expecting floating point numbers to be exactly equal at the binary level?

We run with tolerance 1e-09 in Postgis Topology.

In QMS 12 when we got one single payload where all points matched I mean and then did we not need to do any snapto, when doing pre checks of the input lines, as we have to do now.

Yes we are able to handle this, but it makes things more complicated.

[meastp]

This is an intended change, unfortunately.

[larsop]

Ok, I think we need discuss this more closely. I will try to make add more comments on Monday.

[larsop]

The main reason for not snapping is related to performance when adding data as I understand it.

I think there are some major problems related to not doing snapping when updating/adding new data in QMS.

This means that in the QMS database there will raw data with invalid rings (non-closed) that are to be used in polygons just like we got through NGIS-OpenApi when we tested.

If the QMS database then are going to do inconsistency check it has to do some kind of snapping before doing that check. The same when QMS needs to deliver valid simple feature polygons it also also have to some kind of snapping.

That can be solved in different ways

• Before running ST_MakePolygon check that all rings are closed if not try to run for instance ST_SnapToGrid and check if it closes (like we do here https://gitlab.com/nibioopensource/pgtopo_update_sql/-/commit/71b6166967b4c4db4126be22610e7e099d60c132)
• Run a try catch in PostgresSQL do some snapping if you get an exception when calling. ST_MakePolygon
• Always do some kind of snapping for all rings to avoid exceptions and testing on ST_IsClosed.

All this issues will have a performance impact and that will happen every time the server has to deliver data. (You may make cache on top to help out this but that adds even more complexity to the QMS database)

Another problem is what should you do if you are not able to make closed ring, when building a polygon. For instance should you just continue increase tolerance in a loop and break when the ring is closed, but then you have no control of what the resulting ring will lock like. Since this may happens a long time after data for instance related to a new software version that affects snapping.

Another problem is that when doing snap off a single ring that have edges that have a relation 2 different polygons, there are no guarantee that will snap to the same coordinates. You will end with tiny spikes, if the receiver are using plain simple feature database, with no global snapping.

And it will not make thing any better when this snapping is done different system and different coordinate systems around in Norway.

So I think it’s important that this works the same way as it does in production today (with snapping) to avoid a lot extra work and to secure data quality.

[larsop]

As understand this, the new rule in QMS 13 now, is that if two points have as distance less than 0.1 millimeter they are considered to be equal, but they are not snapped in to a single point.

Lets say we 2 points A and B and the distance between A and B are less than 0.1 millimeter. This is then valid data in QMS 13 and when doing checks this will be considered as one single point.

A   B
.   .

Then we add point C, which have distance less than 0.1 millimeter to B so this is also as valid point. Then we get this dataset. The distance between A and C is more than 0.1 millimeter.

A   B   C
.   .   .

This 3 points are endpoints for 3 different edges that logically should meet in one single point. That may happen if A snaps B and C snaps B, like this

    (A,B,C  )
    .   

But we might as well end up with case like this, where B snaps to A and then C can not snap to AB because if thats out of tolerance. We then end up with tiny gaps.

(A,B)       C
.       .

With 4 points this gets even more complicated and we may end having to lower the tolerance to get all points into a single point.

Here is some sample code trying to show the problem that may arise when based on a very simple distance check and not doing any snap to

do $$
declare
point_num integer := 1;
x real := 631451; 
y real := 6604232;      
max_snap_to real := 0.0001;

master_point geometry;
before_point geometry;
new_point geometry;
point_collection geometry;
srid int = 25832;
v_point_array geometry[];
inv_point_array geometry[];
status varchar;
begin

new_point = ST_SetSrid(ST_MakePoint(x,y+((max_snap_to/3)*point_num)),srid);
master_point = new_point;
before_point = new_point;

v_point_array = ARRAY_APPEND(v_point_array,new_point);

while point_num <10 loop
    point_num := point_num+1;
    new_point = ST_SetSrid(ST_MakePoint(x,y+((max_snap_to/3)*point_num)),srid);
    IF ST_Distance(master_point,new_point) < max_snap_to THEN 
    status = 'VALID POINT ADDED  ';
    v_point_array = ARRAY_APPEND(v_point_array,new_point);
    ELSE 
    status = 'INVALID POINT ADDED';
    inv_point_array = ARRAY_APPEND(inv_point_array,new_point);
    END IF;

    raise notice 'Addding point num % % , distance to master % , distance to last point %', 
    point_num,status,ST_Distance(master_point,new_point), ST_Distance(before_point,new_point);  

    before_point = new_point;
    before_point = new_point;
end loop;

raise notice 'Valid Points  % ', ST_AsText(ST_Collect(v_point_array));
raise notice 'Invalid Points  % ', ST_AsText(ST_Collect(inv_point_array));

end$$;

Here we get a list of valid and invalid points but all of then are valid based on the distance check if I select the last point added when comparing

If check this at scale 1:1 it locks OK.

but we know that this is wrong and with this points we will get overlap and gaps in AR5.

When zoom into a scale at 1000:1 we se problem visualy.

[larsop]

To compute distance has differences when working in degrees and UTM zones

This is sample data we get from NGIS Open API

{ "navnerom": "http:\/\/data.geonorge.no\/SFKB\/FKB-AR5\/so", "lokalId": "0d767a8e-a8a4-4aad-916f-707fcb7cebca", "versjonId": "2023-02-17 03:15:22.979076000" }

 "coordinates": [
          [
            654033.2899999999,
            7812399.39
          ],
          [
            654032.5399999999,
            7812386.27
          ]
        ]

{ "lokalId": "6f716f91-09eb-4165-b89f-bfaf73f177f9", "navnerom": "http:\/\/data.geonorge.no\/SFKB\/FKB-AR5\/so", "versjonId": "2023-04-20 16:39:24.853336000" }

"coordinates": [
          [
            654032.54,
            7812386.27
          ],
          [
            654033.34,
            7812385.42
          ],
          [
            654033.77,
            7812384.27
          ],
          [
            654034.02,
            7812382.83
          ],
          [
            654034.2,
            7812381.15
          ],
          [
            654034.51,
            7812378.45
          ],
          [
            654035.39,
            7812370.56
          ],
          [
            654035.975,
            7812364.747
          ]
        ]

Below we the complete line string where this two lines are used after converted to degrees

SELECT ST_IsClosed('SRID=4258;LINESTRING(31.11102944213983 70.37038782385574,31.11104919371457 70.37037973728052,31.11105856188665 70.37036919159968,31.111062619238577 70.37035616121248,31.111064382149294 70.3703410267053,31.111067766453775 70.37031669084811,31.11107696738457 70.37024559174006,31.111082061039635 70.37019324744253,31.111082882793998 70.37018472017874,31.111096875182398 70.37010291000504,31.11110457240721 70.37003894307016,31.111104315121306 70.37002439362642,31.111103859731852 70.37000886036874,31.11109776636949 70.36999983460704,31.111102438492438 70.36990359141274,31.11111140867866 70.3698963791409,31.11112193486194 70.36988760404529,31.111134919444464 70.36981273505832,31.111147533399393 70.3697094826833,31.111153649647783 70.36963899285314,31.11115521846009 70.36962422214988,31.111153669659767 70.36961653051394,31.11114873876065 70.36961062296068,31.111153011336103 70.36952301421394,31.111160173238922 70.36951611263676,31.111165270786096 70.36950692201015,31.111175228540297 70.36945152897289,31.111189074982942 70.3693570535903,31.11120768115925 70.3692357851451,31.111224880484013 70.36913009234024,31.111234848095915 70.36905484275994,31.11123713869192 70.3690383485209,31.1112359999956 70.3690313663313,31.111226236735483 70.36902269365459,31.11121399757623 70.36901632352206,31.11119650045556 70.36901303803796,31.111449282324152 70.36901294450949,31.111435890219145 70.36901810112336,31.111428722980037 70.36902895613116,31.111422605394065 70.36907554657469,31.111407855287496 70.36918145321174,31.11139314372625 70.36929418738093,31.111379827349037 70.36939925274511,31.111359754685388 70.36952783225651,31.111346549036984 70.36962283216427,31.111332424757833 70.36971327075271,31.11131832101824 70.36981044743575,31.1113111355902 70.36988491490932,31.111311343118818 70.36990054365889,31.11131705832819 70.36990769122355,31.11130076276515 70.3699900030964,31.11129544904048 70.36999812047536,31.111289091509597 70.37003025090684,31.11127710057537 70.37011668762902,31.11126517746239 70.37020213448135,31.11125539074223 70.37029022813076,31.111253417105978 70.37031758670452,31.111245216729674 70.37034454791782,31.111254916875275 70.37041925998398,31.11107302313176 70.37050471192268,31.111029442139824 70.37038782385574)');
 st_isclosed 
-------------
 f
(1 row)

If we lock at the start and end point this line in the original data from NGIS-OpenAPI we find this 2 points and then we check the distance we see that we have a distance and that is less than 1 millimeter that make sense.

select ST_Distance(
'SRID=25835;POINT (654032.5399999999 7812386.27)',
'SRID=25835;POINT (654032.54 7812386.27)'
);

      st_distance       
------------------------
 1.1641532182693481e-10
(1 row)

In QMS they have 3 different archives one for each UTM for AR5 in Norway. This split splitting affects data consistency checks and complicates client code for updating in the border areas between the zones.

So in Postgis Topology we have one single archive to avoid this problems and then we have to use degrees to get correct transformations across all of Norway.

Lets see what to distance when working in degrees with the same two points.

select ST_Distance(
ST_transform('SRID=25835;POINT (654032.5399999999 7812386.27)',4258),
ST_transform('SRID=25835;POINT (654032.54 7812386.27)',4258),
true
);

st_distance 
-------------
           0
(1 row)

In degrees we then have a start and end point with zero distance but the linestring is still not closed.

[larsop]

This evening Karverket/NorKart added this code and then the non-closed rings seems to be gone for AR5.

UPDATE q_geometry SET geo = ST_SnapToGrid(geo, ST_MakePoint(0.0, 0.0), resolution, resolution, 0, 0)

[larsop]

This fix was added to the archive for EPSG:25835 in QMS 13 and after that I was not able reproduce the error.

For the archive EPSG:25832 the fix was not applied and there we still had the the error.

[larsop]

There is internal discussion around ST_SnapTo that a point may snap to "wrong?" a grid point and then cause a gap.

This can be checked easy with a very a low cost in the pre-commit code when data are changed, since QMS logically knows the line strings before and after the line's added.

It can just use ST_Touches on the actual lines confirm that the new point's are OK related to existing lines in the database.

select  ST_Touches(
'SRID=25835;POINT (654032.5399999999 7812386.27)',
'SRID=25835;POINT (654032.54 7812386.27)'
);

 st_touches 
------------
 f

select ST_Touches(
ST_transform('SRID=25835;POINT (654032.5399999999 7812386.27)',4258),
ST_transform('SRID=25835;POINT (654032.54 7812386.27)',4258)
);

 st_touches 
------------
 f
(1 row)

[larsop]

Here is some simple code to test performance related to secure that a point in new line is connected to the next edge . This seems to take less that 100 ms to check that 100000 new lines with two points match the start point of a master line. We also have to check the last point the new line let's double the time 200 ms for 100000.

The code is here

do $$
declare
point_num integer := 1;
x real := 631451; 
y real := 6604232;      
max_snap_to real := 0.00000001;

master_point geometry;
new_start_point_from_added_line geometry;
new_end_point_from_added_line geometry;

master_line geometry;
new_line geometry;

before_point geometry;
new_point geometry;
point_collection geometry;
srid int = 25832;

v_point_dist_array geometry[];
inv_point_dist_array geometry[];

v_point_touch_array geometry[];
inv_point_touch_array geometry[];

status varchar;
origin_point geometry;
resolution numeric = 0.0001;
begin

origin_point = ST_SetSrid(ST_MakePoint(0.0, 0.0),srid);

new_point = ST_SetSrid(ST_MakePoint(x,y+((max_snap_to/3)*point_num)),srid);
new_point = ST_SetSrid(ST_SnapToGrid(new_point, origin_point, resolution, resolution, 0, 0),srid);
master_point = new_point;
before_point = new_point;
v_point_dist_array = ARRAY_APPEND(v_point_dist_array,new_point);
v_point_touch_array = ARRAY_APPEND(v_point_touch_array,new_point);

master_line = ST_SetSrid(ST_MakeLine(new_point, ST_SetSrid(ST_MakePoint(x,y-((max_snap_to/3)*10)),srid)),srid);

while point_num <100000 loop
    point_num := point_num+1;
    new_point = ST_SetSrid(ST_MakePoint(x,y+((max_snap_to/3)*point_num)),srid);

    new_point = ST_SetSrid(ST_SnapToGrid(new_point, origin_point, resolution, resolution, 0, 0),srid);

    IF ST_Distance(master_point,new_point) < max_snap_to THEN 
        status = 'VALID POINT ADDED  ';
        v_point_dist_array = ARRAY_APPEND(v_point_dist_array,new_point);
    ELSE 
        status = 'INVALID POINT ADDED';
        inv_point_dist_array = ARRAY_APPEND(inv_point_dist_array,new_point);
    END IF;

    new_line = ST_SetSrid(ST_MakeLine(new_point, ST_SetSrid(ST_MakePoint(x,y-((max_snap_to/3)*10)),srid)),srid);
    new_start_point_from_added_line = ST_StartPoint(new_line);

    IF ST_Equals(ST_StartPoint(master_line),new_start_point_from_added_line) THEN -- real   0m0.372s
    -- IF TRUE THEN -- real 0m0.308s
        status = 'VALID POINT TOUCH ADDED  ';
        v_point_touch_array = ARRAY_APPEND(v_point_touch_array,new_point);
    ELSE 
        status = 'INVALID POINT TOUCH ADDED '|| ST_ASText(master_point) || ' ' || ST_ASText(new_point) ;
        inv_point_touch_array = ARRAY_APPEND(inv_point_touch_array,new_point);
    END IF;

--  raise notice 'Addding point num % % , distance to master % , distance to last point %', 
--  point_num,status,ST_Distance(master_point,new_point), ST_Distance(before_point,new_point);  

    before_point = new_point;
    before_point = new_point;
end loop;

--raise notice 'Valid Points  % ', ST_AsText(ST_Collect(v_point_dist_array));
--raise notice 'Invalid Points  % ', ST_AsText(ST_Collect(inv_point_dist_array));

raise notice 'Valid length points  % ', ARRAY_LENGTH(v_point_dist_array,1);
raise notice 'Invalid length points  % ', ARRAY_LENGTH(inv_point_dist_array,1);

raise notice 'Valid touch Points  % ', ARRAY_LENGTH(v_point_touch_array,1);
raise notice 'Invalid touch Points  % ', ARRAY_LENGTH(inv_point_touch_array,1);

end$$;

[larsop]

Snapping code was added to production in late May