GEOSGeom_setPrecision returns error for valid geometry

[kalenik]

GEOSGeom_setPrecision return error on valid polygon geometry. For geometry provided in the example only 612 gridSize will cause an error (maybe 612 is not the only value that allows to reproduce the problem but I didn't manage to find another example).

00:03:04 [gis] > select ST_ReducePrecision(ST_GeomFromEWKT('SRID=3857;POLYGON((3668147.751966311 3401469.706472568,3668165.0621471293 3401466.722306755,3668185.8121002135 3401460.6009431183,3668207.4637411726 3401441.930801724,3668265.383272232 3401363.1186134624,3668300.771738355 3401321.901813607,3668313.462160306 3401306.381783725,3668333.911550765 3401276.82108756,3668366.7619324974 3401213.8231056267,3668432.740994691 3401094.2425486096,3668541.4222135525 3400913.2988541853,3668667.1353145055 3400718.0242803628,3668737.6005521775 3400589.0372875025,3668754.020177069 3400537.3790869503,3668779.6904516467 3400497.3763204413,3668794.317832737 3400456.0092274295,3668806.718824011 3400412.4234491913,3668815.524195733 3400388.9856331884,3668828.9604582717 3400365.9941711505,3668877.6293396465 3400274.6280562477,3668923.3148586676 3400208.0516982456,3669016.4670085637 3400094.881593592,3669100.0790780988 3399987.9861564394,3669169.6426278944 3399902.5902982233,3669257.1286157095 3399752.1128553743,3669321.026003425 3399626.7056830493,3669355.3569343863 3399556.8428698564,3669396.266847252 3399492.6418365375,3669413.888722645 3399460.700824149,3669430.608910162 3399414.7084535174,3669444.0451727 3399386.94989833,3669480.4689100883 3399344.2474863273,3669535.7056414196 3399255.8721699594,3669636.627891773 3399058.7990224855,3669663.5004168507 3399001.179658935,3669704.1097670915 3398924.7409967543,3669737.249579501 3398848.302781244,3669828.0194922932 3398690.3531348244,3669846.2313609873 3398634.520405572,3669862.3504232545 3398538.9716382483,3669865.344917557 3398483.738802024,3669868.027717285 3398446.407044794,3669885.939023354 3398401.9227231788,3669908.637067526 3398365.2031907514,3669933.1162235523 3398345.491984257,3669962.682680307 3398328.7642566375,3670002.6909052976 3398323.3965985705,3670037.6229615086 3398318.9214269016,3670071.363899168 3398301.0080057825,3670101.820911849 3398281.5901399856,3670140.036893038 3398247.5611327984,3670202.142036952 3398154.093378526,3670236.328252575 3398102.1515541063,3670255.141246519 3398077.506671257,3670273.4978305507 3398063.1761968303,3670307.5393308355 3398046.6145737856,3670343.963068223 3398015.5568060298,3670386.364662266 3397951.3633366358,3670424.1687613395 3397873.6558740186,3670455.638781387 3397794.0109847677,3670501.023737783 3397683.2331062667,3670527.29513761 3397626.2062593116,3670539.5402815975 3397571.563732348,3670542.5236439505 3397542.903907426,3670547.2992501054 3397520.503901133,3670564.921125498 3397500.806716133,3670598.361500533 3397464.3701892784,3670660.767207071 3397372.4123940207,3670716.0039384025 3397246.1101231193,3670738.991413251 3397200.724725436,3670754.5204822174 3397171.7597449264,3670827.078526316 3397066.966204879,3670902.619932769 3396967.2346406784,3670913.9633888807 3396951.1078516874,3670939.6336634574 3396937.3777869204,3670995.4715200397 3396926.0316904164,3671077.280213823 3396905.4302639295,3671203.8838707027 3396908.120176068,3671334.962571111 3396904.8310892633,3668137.299066126 3396904.8310892633,3668137.299066126 3401475.9808747065,3668147.751966311 3401469.706472568))'), 612);
ERROR:  XX000: lwgeom_reduceprecision: GEOS Error: TopologyException: side location conflict at 3670776 3397212. This can occur if the input geometry is invalid.
LOCATION:  pg_error, lwgeom_pg.c:342

GEOS Version: 3.10.0dev

[mdavis]

Fails in JTS as well, using PrecisionModel scale factor of 1/612 = 0.001633986928105. (Note that the GEOS grid size is the reciprocal of the JTS scale factor).

It works using a scale factor of 0.00163398.

Precision Reduction works better with scale factors/grid sizes with few digits of precision (e.g. 600 or 610 rather than 612). What is the reason for using such a specific grid size? Is it possible to use a lower-precision one?

It's going to be tricky to identify why this is failing, and also to fix it, if that's even possible. So it will probably be necessary to explore perturbing the grid size slightly.

[komzpa]

What does "works better" even mean?

The C api promise is that the output will be a valid geometry: ​https://github.com/libgeos/geos/blob/9e3aa68491ca56a93181e6a11272bb1feef06161/capi/geos_c.h.in#L3709

612 is pre-rounded mercator pixel size at resolution 8. Using non-rounded value failed in a different corner in a different fashion:

 error: ERROR: lwgeom_reduceprecision: GEOS Error: TopologyException: Ring edge missing at 2027415.73823134 8183195.3741355296 (SQLSTATE XX000)

[kalenik]

Reproducer for the "ring edge missing" exception above:

21:46:32 [gis] > select ST_ReducePrecision('SRID=3857;POLYGON((2028943.4369208599 8183234.881054974,2027337.0600453038 8183234.881054974,2027336.6513907493 8183236.913533773,2027332.8442641646 8183237.367598263,2027322.7475863493 8183235.7459395,2027317.6268897727 8183238.232483071,2027316.7474657956 8183240.005497266,2027313.0850545487 8183240.437939815,2027306.9402186568 8183236.200003917,2027306.1738633113 8183234.881054974,2027112.7433669677 8183234.881054974,2027112.7433669677 8183529.661146606,2028943.4369208599 8183529.661146606,2028943.4369208599 8183234.881054974),(2027367.9766954589 8183235.7459395,2027369.880258751 8183237.367598263,2027375.145670666 8183245.411030931,2027376.0250946432 8183248.935440996,2027375.7356639667 8183252.740940863,2027374.5668093138 8183255.508578343,2027370.3144047656 8183262.103343689,2027369.5908280753 8183264.5898960745,2027369.5908280753 8183268.395403939,2027370.6149673907 8183270.449514109,2027375.0009553276 8183272.482002624,2027382.169930535 8183272.78471373,2027385.0976331427 8183271.314402774,2027387.1459117732 8183266.773738418,2027388.4706137138 8183266.0602057045,2027390.0736143813 8183266.773738418,2027390.8083230206 8183268.546759413,2027390.5188923448 8183270.730602918,2027388.1700510888 8183273.952313813,2027388.0253357503 8183275.855069901,2027391.398316322 8183279.37949439,2027391.5430316597 8183281.130896133,2027389.049475066 8183284.633700866,2027384.6634871287 8183286.536459686,2027381.7357845209 8183289.325732009,2027378.2180886115 8183296.050261064,2027379.9769365662 8183298.839536964,2027380.266367242 8183302.645062346,2027379.6875058904 8183304.396469572,2027376.9045186206 8183307.618195051,2027373.8321006747 8183309.088513183,2027371.483259419 8183309.088513183,2027366.0731321662 8183307.034392344,2027365.4831388649 8183301.477457763,2027364.3142842117 8183300.158497241,2027360.9524355896 8183298.536824796,2027355.675891726 8183301.02338936,2027349.9763337977 8183302.493706188,2027348.7963471948 8183294.882657515,2027346.7480685643 8183288.0067736525,2027343.8203659565 8183290.493334692,2027342.5067959654 8183289.325732009,2027338.5549540422 8183295.466459267,2027335.6272514344 8183298.104378869,2027333.7236881417 8183298.536824796,2027332.6995488266 8183297.217864807,2027326.6994282724 8183296.201617097,2027336.5066754117 8183274.103669395,2027338.5549540422 8183271.617113829,2027340.0243713208 8183271.314402774,2027339.723808695 8183273.065802623,2027340.8926633487 8183267.6602488635,2027346.8927839026 8183265.173695451,2027351.7240498029 8183260.351946419,2027350.6999104875 8183252.135520296,2027353.9281757206 8183247.897577736,2027355.3864610502 8183244.6758781215,2027357.1453090047 8183244.524523144,2027358.603594334 8183245.411030931,2027359.4830183114 8183246.729981743,2027361.0971509276 8183246.881336766,2027361.0971509276 8183245.129943086,2027359.3383029734 8183241.756889739,2027359.7835809365 8183238.535192785,2027361.9765749052 8183236.329736616,2027365.193708189 8183234.881054974,2027367.9766954589 8183235.7459395))', 152.87406);
ERROR:  XX000: lwgeom_reduceprecision: GEOS Error: TopologyException: Ring edge missing at 2027415.7837199997 8183195.5577399991
LOCATION:  pg_error, lwgeom_pg.c:342

Thanks for the answer: limiting scale precision seems to help for some cases but for my use case, large variety of polygons need to be processed with ST_ReducePrecision and choosing "right" scale factor feels like guessing game -- there are always polygons that causes an error for chosen scale factor.

[mdavis]

"Works better" means "will produce valid output".

I know the C API promises to produce valid output. That is the intent of this function. But the implementation of this kind of topological process is subject to robustness errors caused by numerical precision issues. Not an easy problem to solve.

[mdavis]

Here's something to try: for grid sizes which are greater than 1, first reduce precision to a gridsize of 1, and then to the desired grid size. That eliminates the decimal places in the input first, which provides more precision "headroom" to correctly compute the reduction to the desired gridsize.

I am keen to hear if this works. If it does, then this could be built into the API as a fix for this problem.

[kalenik]

Example of query that fails when gridSize equal to 1 is specified:

select ST_ReducePrecision(ST_GeomFromEWKT('SRID=3857;POLYGON((-6525887.72688684 -3639625.53883344,-6523441.74198171 -3637179.5539283105,-6523441.74198171 -3637179.55392831,-6520995.75707658 -3634733.56902318,-6516103.78726632 -3634733.56902318,-6514880.794813755 -3633510.5765706147,-6516103.78726632 -3632287.58411805,-6516103.78726632 -3629841.59921292,-6513657.80236119 -3627395.6143077905,-6513657.80236119 -3624949.62940266,-6511211.81745606 -3622503.6444975296,-6511211.81745606 -3620057.6595924,-6508765.83255093 -3617611.6746872696,-6508765.83255093 -3615165.68978214,-6506319.8476458 -3612719.7048770096,-6506319.8476458 -3610273.71997188,-6503873.862740669 -3607827.7350667496,-6496535.90802528 -3607827.7350667496,-6494089.923120149 -3610273.71997188,-6484305.983499629 -3610273.71997188,-6481859.998594499 -3612719.70487701,-6481859.998594499 -3615165.68978214,-6484305.983499629 -3617611.67468727,-6484305.983499629 -3620057.6595924,-6481859.998594499 -3622503.64449753,-6481859.998594499 -3629841.5992129208,-6479414.01368937 -3632287.58411805,-6479414.01368937 -3637179.5539283096,-6476968.02878424 -3637179.5539283096,-6474522.043879109 -3639625.53883344,-6472076.05897398 -3639625.53883344,-6469630.07406885 -3642071.5237385696,-6464738.104258589 -3642071.5237385696,-6462292.119353459 -3644517.5086437,-6459846.13444833 -3644517.5086437,-6457400.1495432 -3646963.4935488296,-6454954.164638069 -3649409.47845396,-6454954.164638069 -3650842.303238026,-6453224.5921249725 -3652571.8757511224,-6452508.179732939 -3654301.44826422,-6452508.179732939 -3654301.4482642207,-6450062.19482781 -3656747.43316935,-6450062.19482781 -3659193.41807448,-6447616.209922681 -3661639.4029796096,-6447616.20992268 -3661639.4029796096,-6445170.2250175495 -3664085.38788474,-6445170.2250175495 -3666531.37278987,-6447616.20992268 -3668977.3576950002,-6447616.20992268 -3671423.34260013,-6450062.19482781 -3673869.3275052602,-6452508.179732939 -3676315.312410389,-6452508.179732939 -3678761.29731552,-6453224.5921249725 -3679477.709707553,-6452508.179732939 -3681207.28222065,-6452508.179732939 -3683653.2671257798,-6453224.5921249725 -3685382.839638877,-6455670.577030103 -3687828.8245440074,-6457400.1495432 -3688545.23693604,-6459846.13444833 -3690991.22184117,-6462292.119353459 -3690991.22184117,-6462292.119353459 -3695883.1916514295,-6459846.13444833 -3695883.1916514295,-6457400.1495432 -3698329.1765565597,-6454954.164638069 -3700775.16146169,-6454954.164638069 -3703221.1463668197,-6452508.179732939 -3705667.13127195,-6452508.179732939 -3715451.07089247,-6454954.164638069 -3717897.0557976,-6454954.164638069 -3722789.0256078597,-6457400.1495432 -3725235.01051299,-6457400.1495432 -3727680.9954181197,-6459846.13444833 -3730126.98032325,-6462292.11935346 -3732572.96522838,-6464738.104258589 -3732572.96522838,-6467184.08916372 -3730126.98032325,-6469630.07406885 -3730126.98032325,-6472076.05897398 -3732572.96522838,-6479414.01368937 -3732572.96522838,-6481859.9985945 -3730126.98032325,-6484305.983499629 -3730126.98032325,-6486751.968404759 -3727680.99541812,-6486751.96840476 -3727680.99541812,-6489197.95330989 -3725235.01051299,-6491643.93821502 -3722789.0256078597,-6491643.93821502 -3713005.0859873397,-6490927.525822987 -3711275.5134742428,-6491643.93821502 -3710559.10108221,-6491643.93821502 -3705667.13127195,-6494089.923120149 -3705667.13127195,-6496535.90802528 -3703221.1463668197,-6496535.90802528 -3698329.17655656,-6498265.480538377 -3697612.7641645274,-6500540.204014536 -3695338.040688368,-6507676.635952639 -3690580.419396299,-6508765.83255093 -3688545.2369360398,-6508765.83255093 -3687112.412151974,-6510495.405064027 -3685382.839638877,-6510495.405064027 -3681923.6946126823,-6508049.420158897 -3679477.7097075526,-6506319.8476458 -3678761.29731552,-6506319.8476458 -3674882.487626324,-6508049.420158897 -3673152.9151132274,-6508765.83255093 -3671423.34260013,-6511211.81745606 -3668977.357695,-6511211.81745606 -3666531.37278987,-6508765.83255093 -3664085.38788474,-6506319.8476458 -3661639.4029796096,-6503873.862740669 -3661639.4029796096,-6501427.87783554 -3659193.41807448,-6501427.87783554 -3661639.4029796096,-6501427.87783554 -3661639.40297961,-6501427.87783554 -3662862.3954321747,-6502650.870288106 -3664085.38788474,-6503873.86274067 -3664085.38788474,-6503873.86274067 -3661639.40297961,-6503873.86274067 -3661639.4029796096,-6506319.847645799 -3659193.4180744803,-6506319.8476458 -3659193.4180744803,-6508765.83255093 -3656747.43316935,-6513657.802361189 -3656747.43316935,-6516103.78726632 -3654301.4482642203,-6518549.77217145 -3651855.46335909,-6520995.75707658 -3649409.4784539603,-6523441.74198171 -3646963.49354883,-6523441.74198171 -3644517.5086436993,-6525887.72688684 -3642071.52373857,-6525887.72688684 -3639625.53883344))'), 1);

[mdavis]

Thanks. This is a simpler example, so will be easier to investigate if this is a topology bug or a robustness problem.

[mdavis]

The above example fails because the input polygon is invalid. It has a self-crossing loop.

You can fix it using MakeValid. You may wish to use the new GEOSMakeValidWithParams_r, to eliminate the self-crossing "hole".

[kalenik]

right, sorry, it is not related to GEOSGeom_setPrecision but this geometry was produced by ST_Buffer so I reported another ticket #1131

[mdavis]

Have you found any issues using the hack of reducing to gridSize=1 and then the desired gridSize?

[kalenik]

Yes, unfortunately I still having issues. Example of geometry with precision reduced to gridSize = 1 that is failing on desired gridSize is here ​https://gist.github.com/kalenikaliaksandr/e5c1b959335ac8dbaacac6b11fe16317.

Another thing is that I would like to use ST_Union in PostGIS specifying gridSize but seems like it is suffering from the same sort of problem with numerical precision and throws error on valid input geometries. I can provide you with a reproducer.

[mdavis]

Replying to kalenik:

Yes, unfortunately I still having issues. Example of geometry with precision reduced to gridSize = 1 that is failing on desired gridSize is here ​https://gist.github.com/kalenikaliaksandr/e5c1b959335ac8dbaacac6b11fe16317.

Thanks for the failure case. I've confirmed that this fails in JTS as well.

One thing that did work in JTS is to use a bigger gridsize for the initial reduction, and compensate the secondary reduction gridsize accordingly. So in this case the reductions used 100 and then 12.2299245256.

Obvious then the question is which strategy to use. This could be done by iterating using larger primary reduction sizes until a working one is found.

Not great, but the best solution I have at the moment.

[mdavis]

Replying to kalenik:

Another thing is that I would like to use ST_Union in PostGIS specifying gridSize but seems like it is suffering from the same sort of problem with numerical precision and throws error on valid input geometries. I can provide you with a reproducer.

It's not surprising that ST_Union suffers from the same issue, since it uses the same codebase.

One workaround is to precision reduce both geometries first (using the strategies discussed here). That should allow union to work correctly.

[robe]

Retargeting in prep for GEOS 3.10.0 release

[kalenik]

broken noded segments

[kalenik]

normal noded segments

[kalenik]

I tried to investigate this problem a bit myself. I got visualisations of how noded segments look for "broken" scaleFactor vs how noded segments look for "normal" scaleFactor. Geometry from ticket description is used. Not sure if that is helpful at all.

[mdavis]

I can reproduce this in JTS, and the noded segment situation is the same. Clearly something is wrong. But it could be a while (if ever) until a fix to the existing algorithm is identified.

As mentioned before, the precision reduction code is tailored towards situations where the actual numerical precision of the input numbers is reduced. Also, the gridsize argument is actually inverted for use by the actual code - so that 612 becomes 0.001633986928105. This is not reducing the precision of the input numbers, only changing them.

Instead, try using a slightly larger grid size which is a "simple" multiple of 5, 10, and 2. I have confirmed this works for all 3 cases provided so far on this ticket.

• 612 -> 800
• 152.87406 -> 200
• 1222.99245256 -> 1250

[mdavis]

Further analysis indicates this problem is caused by a ​heuristic (really a hack) in the SnapRoundingNoder logic. This was intended to increase robustness, but it appears to have some (now) obvious failure modes.

The fix is not clear right now.

[kalenik]

Thanks for the progress on this problem! I am wondering if NEARNESS_FACTOR can be computed based on input geometry somehow.

[mdavis]

Replying to kalenik:

I am wondering if NEARNESS_FACTOR can be computed based on input geometry somehow.

It's not clear to me how. Input geometry can be arbitrarily complex, so attempting to derive some sort of tolerance from it has always seemed to be (a) difficult and (b) slow. The NEARNESS_FACTOR is already used in conjunction with the requested grid size (scale factor).

I'm a bit more hopeful that the tests for when to invoke the "nearby snapping" logic can be improved.

[sergeish]

simplified example

[sergeish]

Simplified polygon causing same exception in GEOSGeom_setPrecision(polygonGeometry, 612, 0):

POLYGON((3670939.6336634574 3396937.3777869204, 3670995.4715200397 3396926.0316904164, 3671077.280213823 3396905.4302639295, 3671203.8838707027 3396908.120176068, 3671334.962571111 3396904.8310892633, 3670037.299066126 3396904.8310892633, 3670037.299066126 3398075.9808747065, 3670939.6336634574 3396937.3777869204))

TopologyException: side location conflict at 3670776 3397212. This can occur if the input geometry is invalid.

As explained by mdavis, an intersection is detected by noder at point 4 (see image) by a heuristic because point 4 happens to be close enough to segment (5, 6), within distance of 612/100.

[sergeish]

Partially ​fixed, the case in the description and simplified one no longer cause an exception.
Cases that still break GEOSGeom_setPrecision are in comment:4 and comment:7

[komzpa]

For the geometry in comment:7 it is okay to fail, as it fails IsValid test - it should have never existed and the bug that caused its generation is in #1131.

So, only comment:4 left?

[komzpa]

A piece of discussion of this issue is in now closed ​https://github.com/libgeos/geos/issues/511

[sergeish]

Issue ​https://github.com/libgeos/geos/issues/511 now reopened, simplified case for comment:4 added: ​https://github.com/libgeos/geos/issues/511#issuecomment-984494242.