BUG: valid geometry fails with topology error in GEOSGeom_setPrecision

[brendan-ward]

Using the default flag to produce valid outputs, an input polygon that is valid, and a precision of 0.00001, I'm getting

TopologyException: Ring edge missing at 127.11725000000001 34.562390000000001

(via shapely 2.0 & latest build from GEOS main)

for this geometry:

POLYGON ((127.11125000000001 34.55639, 127.11125000000001 34.56839, 127.12325000000001 34.56839, 127.12325000000001 34.559037113624825, 127.123156453 34.5589634680001, 127.122811919 34.5587243150001, 127.122414652 34.5584866740001, 127.122002904 34.5582775410001, 127.121541652 34.5580788670001, 127.121044598 34.5579037090001, 127.12055893 34.5577654690001, 127.120676924 34.557430595, 127.120795189 34.5569197390001, 127.12083733777376 34.55639, 127.11125000000001 34.55639), (127.117461568 34.562519572, 127.117607152 34.562312309, 127.117578583 34.5623224060001, 127.117254733 34.5621607510001, 127.117314277 34.5620469410001, 127.118203396 34.561428344, 127.118190156 34.5613418900001, 127.118384564 34.5613216470001, 127.118737508 34.561522295, 127.119362293 34.561425193, 127.119599297 34.5614524040001, 127.119733456 34.561485114, 127.11992963 34.561673145, 127.120093027 34.5620903450001, 127.120079113 34.562106242, 127.120189288 34.562190672, 127.120516852 34.5621474980001, 127.120628803 34.5621547800001, 127.120618491 34.5622284200001, 127.120774952 34.562272916, 127.12096395 34.5623545190001, 127.121092847 34.5622872980001, 127.121047121 34.562239017, 127.120954155 34.5622026450001, 127.120956522 34.562114165, 127.121016436 34.562116614, 127.121027842 34.5620463710001, 127.121154732 34.561871357, 127.121171039 34.561875641, 127.121187768 34.561831698, 127.121292353 34.5617865670001, 127.12155222 34.5617921530001, 127.121600392 34.561753272, 127.121769634 34.5612906470001, 127.121871731 34.5613190400001, 127.121888898 34.5613071240001, 127.12189454 34.5613286010001, 127.122906822 34.5618186950001, 127.122965601 34.5620002050001, 127.122436258 34.5623006220001, 127.122388016 34.562411323, 127.122163789 34.562609222, 127.122195096 34.5626851820001, 127.122188291 34.5627269080001, 127.122039124 34.562739865, 127.121948872 34.562636766, 127.121954912 34.562624176, 127.12164618 34.56258689, 127.121555526 34.562560648, 127.121529786 34.563089668, 127.12132344 34.5634194470001, 127.121323793 34.5635131510001, 127.12137824 34.563570618, 127.121008002 34.5639598650001, 127.120796029 34.5640138700001, 127.120592147 34.5636423400001, 127.120479325 34.5633902170001, 127.119916637 34.5626650310001, 127.119816205 34.562558858, 127.119573987 34.562359648, 127.119427371 34.5622498590001, 127.119269955 34.562177192, 127.118943752 34.5620657370001, 127.118614649 34.5620455640001, 127.11852493 34.5620550880001, 127.117776763 34.562382928, 127.11769759 34.5623760980001, 127.117683778 34.5624113220001, 127.117601744 34.562384149, 127.117926302 34.5618559370001, 127.117923867 34.5617969850001, 127.117827728 34.5619306070001, 127.117603496 34.5623196400001, 127.117462315 34.562520066, 127.117245225 34.562385186, 127.117461568 34.562519572))

The part that is running into an issue has a very narrow extension of the hole in the polygon:

In this case, the desired behavior is to eliminate that very narrow cut.

It seems that it is sensitive to the order of magnitude of the precision parameter: 0.0001 and 0.000001 both work fine, but 0.00001 - 0.00008 fails

---

I'm getting entirely different results in geosop, so I'm either not using that correctly or there is a different error there:

❯ geosop reducePrecision 0.00001 -a /tmp/test.wkt
POLYGON EMPTY

and likewise the pointwise operation returns invalid looking data:

❯ geosop reducePrecisionPointwise 0.00001 -a /tmp/test.wkt                                                                                             10ms
POLYGON ((0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0), (0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 0))

[brendan-ward]

This produces the same error if you convert the hole to a polygon, so it isn't specific to the ring being part of the hole rather than exterior:

POLYGON ((127.117461568 34.562519572, 127.117483252 34.5624884690001, 127.117603304 34.562319127, 127.117607152 34.562312309, 127.117607012 34.562312359, 127.117594235 34.562316874, 127.117578583 34.5623224060001, 127.117254733 34.5621607510001, 127.117314277 34.5620469410001, 127.118203396 34.561428344, 127.118190156 34.5613418900001, 127.118384564 34.5613216470001, 127.118737508 34.561522295, 127.119362293 34.561425193, 127.119599297 34.5614524040001, 127.119733456 34.561485114, 127.119896569 34.561641464, 127.119898945 34.5616420740001, 127.119898457 34.5616432640001, 127.11992963 34.561673145, 127.120093027 34.5620903450001, 127.120079113 34.562106242, 127.120091809 34.5621152430001, 127.120130491 34.5621450080001, 127.120189288 34.562190672, 127.120516852 34.5621474980001, 127.120628803 34.5621547800001, 127.120618491 34.5622284200001, 127.120650265 34.562236179, 127.120664228 34.562239409, 127.120738619 34.5622620890001, 127.120774952 34.562272916, 127.120824466 34.5622955210001, 127.120870952 34.5623195540001, 127.12096395 34.5623545190001, 127.121019401 34.5623309100001, 127.121092847 34.5622872980001, 127.121047121 34.562239017, 127.121011428 34.5622189860001, 127.120954155 34.5622026450001, 127.120956522 34.562114165, 127.121016436 34.562116614, 127.121027842 34.5620463710001, 127.121154732 34.561871357, 127.121171039 34.561875641, 127.121187768 34.561831698, 127.121292353 34.5617865670001, 127.121329769 34.5617935450001, 127.121335827 34.561793213, 127.121484042 34.5617850630001, 127.12155222 34.5617921530001, 127.121600392 34.561753272, 127.121600792 34.561735577, 127.121636655 34.5616471810001, 127.121712511 34.5614312930001, 127.121769634 34.5612906470001, 127.121871731 34.5613190400001, 127.121888898 34.5613071240001, 127.121896931 34.5613243710001, 127.12189454 34.5613286010001, 127.122906822 34.5618186950001, 127.122965601 34.5620002050001, 127.122436258 34.5623006220001, 127.122388016 34.562411323, 127.122163789 34.562609222, 127.122182142 34.5626415020001, 127.122195096 34.5626851820001, 127.122188291 34.5627269080001, 127.122122085 34.5627277030001, 127.122039124 34.562739865, 127.121996808 34.5626985660001, 127.121948872 34.562636766, 127.121954912 34.562624176, 127.12164618 34.56258689, 127.121555526 34.562560648, 127.121529786 34.563089668, 127.12132344 34.5634194470001, 127.121323793 34.5635131510001, 127.12137824 34.563570618, 127.121008002 34.5639598650001, 127.120796029 34.5640138700001, 127.120777872 34.563980777, 127.120777034 34.563981815, 127.120775687 34.5639768130001, 127.12073337 34.563899699, 127.120592147 34.5636423400001, 127.120588848 34.5636349060001, 127.120586958 34.563630684, 127.120573063 34.5635996390001, 127.120521851 34.563487009, 127.120479325 34.5633902170001, 127.120441651 34.563333986, 127.120386654 34.5632641780001, 127.120294068 34.5631432760001, 127.120177873 34.5629953550001, 127.120049986 34.5628308930001, 127.119916637 34.5626650310001, 127.119816205 34.562558858, 127.119663441 34.562432173, 127.119573987 34.562359648, 127.119427371 34.5622498590001, 127.119269955 34.562177192, 127.119169572 34.562135174, 127.118943752 34.5620657370001, 127.118848398 34.562058581, 127.118754535 34.562050991, 127.118614649 34.5620455640001, 127.11852493 34.5620550880001, 127.118137453 34.562221234, 127.117941432 34.5623049490001, 127.117776763 34.562382928, 127.11769759 34.5623760980001, 127.117684764 34.5624088070001, 127.11768425 34.5624101180001, 127.117683778 34.5624113220001, 127.117637142 34.5623963030001, 127.117601744 34.562384149, 127.117602256 34.5623833220001, 127.117630384 34.5623379620001, 127.117651691 34.562303603, 127.117651993 34.5623017720001, 127.11771705 34.5621953670001, 127.117849893 34.5619790390001, 127.117892326 34.5619096950001, 127.117926302 34.5618559370001, 127.117926295 34.5618557600001, 127.117925874 34.5618455640001, 127.117923867 34.5617969850001, 127.117827728 34.5619306070001, 127.117746661 34.5620659730001, 127.117603496 34.5623196400001, 127.117484065 34.562488982, 127.117462315 34.562520066, 127.117245225 34.562385186, 127.117461568 34.562519572))

One workaround I discovered for this polygon was to use GEOSDensify to add more points before setting the precision.

[pramsey]

I've confirmed this is reproducible in C++ land, with a test 24 in GeometryPrecisionReducerTest,

template<>
template<>
void object::test<24> ()
{
    std::string wkt("POLYGON ((127.11125000000001 34.55639, 127.11125000000001 34.56839, 127.12325000000001 34.56839, 127.12325000000001 34.559037113624825, 127.123156453 34.5589634680001, 127.122811919 34.5587243150001, 127.122414652 34.5584866740001, 127.122002904 34.5582775410001, 127.121541652 34.5580788670001, 127.121044598 34.5579037090001, 127.12055893 34.5577654690001, 127.120676924 34.557430595, 127.120795189 34.5569197390001, 127.12083733777376 34.55639, 127.11125000000001 34.55639), (127.117461568 34.562519572, 127.117607152 34.562312309, 127.117578583 34.5623224060001, 127.117254733 34.5621607510001, 127.117314277 34.5620469410001, 127.118203396 34.561428344, 127.118190156 34.5613418900001, 127.118384564 34.5613216470001, 127.118737508 34.561522295, 127.119362293 34.561425193, 127.119599297 34.5614524040001, 127.119733456 34.561485114, 127.11992963 34.561673145, 127.120093027 34.5620903450001, 127.120079113 34.562106242, 127.120189288 34.562190672, 127.120516852 34.5621474980001, 127.120628803 34.5621547800001, 127.120618491 34.5622284200001, 127.120774952 34.562272916, 127.12096395 34.5623545190001, 127.121092847 34.5622872980001, 127.121047121 34.562239017, 127.120954155 34.5622026450001, 127.120956522 34.562114165, 127.121016436 34.562116614, 127.121027842 34.5620463710001, 127.121154732 34.561871357, 127.121171039 34.561875641, 127.121187768 34.561831698, 127.121292353 34.5617865670001, 127.12155222 34.5617921530001, 127.121600392 34.561753272, 127.121769634 34.5612906470001, 127.121871731 34.5613190400001, 127.121888898 34.5613071240001, 127.12189454 34.5613286010001, 127.122906822 34.5618186950001, 127.122965601 34.5620002050001, 127.122436258 34.5623006220001, 127.122388016 34.562411323, 127.122163789 34.562609222, 127.122195096 34.5626851820001, 127.122188291 34.5627269080001, 127.122039124 34.562739865, 127.121948872 34.562636766, 127.121954912 34.562624176, 127.12164618 34.56258689, 127.121555526 34.562560648, 127.121529786 34.563089668, 127.12132344 34.5634194470001, 127.121323793 34.5635131510001, 127.12137824 34.563570618, 127.121008002 34.5639598650001, 127.120796029 34.5640138700001, 127.120592147 34.5636423400001, 127.120479325 34.5633902170001, 127.119916637 34.5626650310001, 127.119816205 34.562558858, 127.119573987 34.562359648, 127.119427371 34.5622498590001, 127.119269955 34.562177192, 127.118943752 34.5620657370001, 127.118614649 34.5620455640001, 127.11852493 34.5620550880001, 127.117776763 34.562382928, 127.11769759 34.5623760980001, 127.117683778 34.5624113220001, 127.117601744 34.562384149, 127.117926302 34.5618559370001, 127.117923867 34.5617969850001, 127.117827728 34.5619306070001, 127.117603496 34.5623196400001, 127.117462315 34.562520066, 127.117245225 34.562385186, 127.117461568 34.562519572))");
    checkReduce(100000, wkt, "POLYGON EMPTY");
}

It does seem very sensitive to the particular scale factor used, but the one in the test does generate an error

---> group: geos::precision::GeometryPrecisionReducer, test: test<24>
     problem: unexpected exception
     exception typeid: N4geos4util17TopologyExceptionE
     message: `TopologyException: Ring edge missing at 127.11725 34.562390000000001`

[dr-jts]

I'm getting entirely different results in geosop, so I'm either not using that correctly or there is a different error there

geosop reducePrecision uses a scale factor, not a grid size. So the command should be:

geosop  -a /tmp/test.wkt reducePrecision 100000

This produces the same error as the C API.

[dr-jts]

Here's a reduced polygon which produces the error:

POLYGON ((127.117461568 34.562519572, 127.117483252 34.5624884690001, 127.117603304 34.562319127, 127.117607152 34.562312309, 127.117607012 34.562312359, 127.117254733 34.5621607510001, 127.117746661 34.5620659730001, 127.117603496 34.5623196400001, 127.117484065 34.562488982, 127.117462315 34.562520066, 127.117245225 34.562385186, 127.117461568 34.562519572))

bin/geosop -a "POLYGON ((127.117461568 34.562519572, 127.117483252 34.5624884690001, 127.117603304 34.562319127, 127.117607152 34.562312309, 127.117607012 34.562312359, 127.117254733 34.5621607510001, 127.117746661 34.5620659730001, 127.117603496 34.5623196400001, 127.117484065 34.562488982, 127.117462315 34.562520066, 127.117245225 34.562385186, 127.117461568 34.562519572))" reducePrecision 100000

[dr-jts]

This is due to a bug in a CoordinateSequence::add method, which is called from SegmentNodeList.addEdgeCoordinates during snap-rounding noding. I'll submit a PR soon.