Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Polygon disappears in unaryUnion #965

Open
mortewle opened this issue Sep 27, 2023 · 22 comments
Open

Polygon disappears in unaryUnion #965

mortewle opened this issue Sep 27, 2023 · 22 comments
Labels

Comments

@mortewle
Copy link

I have a list of polygons where one polygon dissappears after unaryUnion. Tested in geosop with GEOS 3.13.0dev (also tested in shapely with 3.11 and QGIS with 3.12).

The disappearing polygon is this one:

POLYGON ((33668.77730000019 6953106.171599999, 33649.936599999666 6953078.862100001, 33650.19849999994 6953082.8554, 33652.10510000028 6953089.261399999, 33654.99129999988 6953094.7304, 33659.10290000029 6953099.521200001, 33664.003800000064 6953103.493299998, 33668.77730000019 6953106.171599999))

dissolve_error.txt

@dr-jts
Copy link
Contributor

dr-jts commented Sep 27, 2023

Weird, I'll look into this.

It doesn't happen in JTS.

@dr-jts dr-jts added the Bug label Sep 27, 2023
@dbaston
Copy link
Member

dbaston commented Sep 27, 2023

Looks like the area of the result changed from 517747 in 3.9 to 517642 in 3.10. I can try bisecting on that.

@dr-jts
Copy link
Contributor

dr-jts commented Sep 27, 2023

The missing polygons are on the extreme right of the input geometry:

image

@dr-jts
Copy link
Contributor

dr-jts commented Sep 27, 2023

While it looks like a single polygon, it's actually 3, one of which is extremely narrow.

Perhaps there's a robustness problem causing an exception to be thrown, and somehow the polygons involved are being dropped from the processing as a result of that?

However, unioning ONLY those polygons does succeed.

image image

@dbaston
Copy link
Member

dbaston commented Sep 27, 2023

Breaking commit is 0b2c29a

@dbaston
Copy link
Member

dbaston commented Sep 27, 2023

It seems that the commit produces a different sequence of pairwise unions, one of which produces an incorrect result both in JTS and GEOS:

A

MULTIPOLYGON (((32716.64690000005 6953061.072500002, 32728.48959999997 6953063.269299999, 32734.278900000267 6953064.8935, 32739.656799999997 6953067.180199999, 32745.473699999973 6953068.993099999, 32751.428299999796 6953069.585000001, 32762.094200000167 6953067.693599999, 32767.982300000265 6953069.298500001, 32774.26150000002 6953068.6219, 32780.27089999989 6953066.671799999, 32786.03369999956 6953064.4527, 32792.24100000039 6953062.574700002, 32798.685700000264 6953060.865899999, 32818.3320000004 6953056.606199998, 32824.6584999999 6953054.8178, 32830.520100000314 6953052.579399999, 32835.6333999997 6953049.333799999, 32838.19000000041 6953043.719300002, 32837.36629999988 6953037.051100001, 32835.28870000038 6953030.862300001, 32832.65309999976 6953025.0579, 32829.450000000186 6953019.537999999, 32825.71710000001 6953014.4903, 32822.153599999845 6953009.316100001, 32819.80709999986 6953003.3739, 32820.83899999969 6952996.955400001, 32824.83980000019 6952992.2531, 32829.79609999992 6952988.619399998, 32834.11930000037 6952984.783, 32837.27190000005 6952979.082699999, 32838.52109999955 6952972.835200001, 32837.626299999654 6952966.375, 32831.765200000256 6952963.549400002, 32827.18130000029 6952959.738899998, 32820.92449999973 6952958.389899999, 32815.08430000022 6952958.592700001, 32808.56010000035 6952960.4197, 32802.83899999969 6952962.222100001, 32799.65859999973 6952967.7337, 32794.04260000028 6952970.230999999, 32788.35570000019 6952973.046999998, 32781.352099999785 6952973.620099999, 32775.8236999996 6952978.2434, 32769.22369999997 6952980.983600002, 32762.215300000273 6952982.261799999, 32755.95079999976 6952983.742400002, 32749.51080000028 6952984.746100001, 32743.37000000011 6952987.121199999, 32738.137600000016 6952990.277399998, 32731.69830000028 6952992.690400001, 32726.76379999984 6953000.872499999, 32721.422299999744 6953003.938099999, 32714.61730000004 6953004.4826, 32704.637900000438 6953006.179000001, 32699.86110000033 6953009.6853, 32693.908900000155 6953012.0326999985, 32691.319099999964 6953016.643599998, 32692.87770000007 6953019.860599998, 32692.272800000384 6953032.9037999995, 32690.786299999803 6953038.9822, 32690.82899999991 6953045.048799999, 32695.238099999726 6953049.791700002, 32700.498499999754 6953053.4081999995, 32705.67399999965 6953056.438700002, 32710.909500000067 6953059.141399998, 32716.64690000005 6953061.072500002)), 
  ((33652.10162677833 6953082, 33651 6953080.403199323, 33651.95034120033 6953081.780712563, 33652.10162677833 6953082)))

B

POLYGON ((33668.7616367347 6953106.1488961745, 33668.762000000104 6953106.149099998, 33649.952600000426 6953078.885000002, 33649.952621165845 6953078.885322602, 33649.936599999666 6953078.862100001, 33650.19849999994 6953082.8554, 33652.10510000028 6953089.261399999, 33654.99129999988 6953094.7304, 33659.10290000029 6953099.521200001, 33664.003800000064 6953103.493299998, 33668.77730000019 6953106.171599999, 33668.7616367347 6953106.1488961745))

image

Result

POLYGON ((32728.48959999997 6953063.269299999, 32734.278900000267 6953064.8935, 32739.656799999997 6953067.180199999, 32745.473699999973 6953068.993099999, 32751.428299999796 6953069.585000001, 32762.094200000167 6953067.693599999, 32767.982300000265 6953069.298500001, 32774.26150000002 6953068.6219, 32780.27089999989 6953066.671799999, 32786.03369999956 6953064.4527, 32792.24100000039 6953062.574700002, 32798.685700000264 6953060.865899999, 32818.3320000004 6953056.606199998, 32824.6584999999 6953054.8178, 32830.520100000314 6953052.579399999, 32835.6333999997 6953049.333799999, 32838.19000000041 6953043.719300002, 32837.36629999988 6953037.051100001, 32835.28870000038 6953030.862300001, 32832.65309999976 6953025.0579, 32829.450000000186 6953019.537999999, 32825.71710000001 6953014.4903, 32822.153599999845 6953009.316100001, 32819.80709999986 6953003.3739, 32820.83899999969 6952996.955400001, 32824.83980000019 6952992.2531, 32829.79609999992 6952988.619399998, 32834.11930000037 6952984.783, 32837.27190000005 6952979.082699999, 32838.52109999955 6952972.835200001, 32837.626299999654 6952966.375, 32831.765200000256 6952963.549400002, 32827.18130000029 6952959.738899998, 32820.92449999973 6952958.389899999, 32815.08430000022 6952958.592700001, 32808.56010000035 6952960.4197, 32802.83899999969 6952962.222100001, 32799.65859999973 6952967.7337, 32794.04260000028 6952970.230999999, 32788.35570000019 6952973.046999998, 32781.352099999785 6952973.620099999, 32775.8236999996 6952978.2434, 32769.22369999997 6952980.983600002, 32762.215300000273 6952982.261799999, 32755.95079999976 6952983.742400002, 32749.51080000028 6952984.746100001, 32743.37000000011 6952987.121199999, 32738.137600000016 6952990.277399998, 32731.69830000028 6952992.690400001, 32726.76379999984 6953000.872499999, 32721.422299999744 6953003.938099999, 32714.61730000004 6953004.4826, 32704.637900000438 6953006.179000001, 32699.86110000033 6953009.6853, 32693.908900000155 6953012.0326999985, 32691.319099999964 6953016.643599998, 32692.87770000007 6953019.860599998, 32692.272800000384 6953032.9037999995, 32690.786299999803 6953038.9822, 32690.82899999991 6953045.048799999, 32695.238099999726 6953049.791700002, 32700.498499999754 6953053.4081999995, 32705.67399999965 6953056.438700002, 32710.909500000067 6953059.141399998, 32716.64690000005 6953061.072500002, 32728.48959999997 6953063.269299999))

image

@dbaston
Copy link
Member

dbaston commented Sep 27, 2023

The second component of A is a tiny sliver that intersects B:

image

@dr-jts
Copy link
Contributor

dr-jts commented Sep 28, 2023

Thanks for the analysis, @dbaston . It's possible this is another instance of the OverlayNG robustness issue in locationtech/jts#1000

This issue occurs in part because of situations which pass the (fairly coarse) heuristics in place to detect the failure. I'm thinking about replacing them with an area-based heuristic. That might resolve this case.

@mortewle
Copy link
Author

@dr-jts I would very much favor the expensive and robust area check you mention in locationtech/jts#1000.

These bugs seem to happen quite a lot in the land use dataset I am trying to produce, earlier made with ArcMap. Especially when doing intersection and difference (#942). Rounding the coordinates is not an appealing option in this case, since we can't have gaps between the polygons.

I can try to track down some more examples, if this is of interest. I came over this one, which results in two holes in the middle when dissolving in QGIS (GEOS 3.12) and shapely (GEOS 3.11.1), but not in GEOSOP (3.13.0dev). Don't know if this is related.

dissolve_error2.txt

@dr-jts
Copy link
Contributor

dr-jts commented Sep 28, 2023

I came over this one, which results in two holes in the middle when dissolving in QGIS (GEOS 3.12) and shapely (GEOS 3.11.1), but not in GEOSOP (3.13.0dev). Don't know if this is related.

The holes in dissolve_error2.txt are legitimate - they are gaps between the polygons being unioned. Here's a picture of one of them:
image

My testing with geosop in 3.13dev produces a result with the holes, so not sure why you're not seeing that.

@dr-jts
Copy link
Contributor

dr-jts commented Sep 28, 2023

@dr-jts I would very much favor the expensive and robust area check you mention in locationtech/jts#1000.

That's a bit of a research project, so not sure when (or if) it will be added. It is also possible to add an envelope-based check, but that probably will only detect a subset of errors.

One thing that would help is to remove extremely small-area polygons (slivers) before processing. This won't solve all problems, but might help with some. As per examples in #942 removing spikes would also help - but that will require custom code.

These bugs seem to happen quite a lot in the land use dataset I am trying to produce, earlier made with ArcMap. Especially when doing intersection and difference (#942). Rounding the coordinates is not an appealing option in this case, since we can't have gaps between the polygons.

At some point (soon hopefully) I'll be working on a Coverage Cleaning operation. Perhaps that will also help solving this issues.

@mortewle
Copy link
Author

@dr-jts Yes, the small holes returned from GEOSOP are legitimate, but not the larger ones returned from QGIS/shapely with different GEOS versions (so probably not relevant here).

Removing slivers would solve some issues, but it will also leave gaps which would need dissolving later. Maybe an option.

Removing spikes is tricky. Are there any good ways to do this other than buffering (which slows things down way too much)?

Currently, I am trying to bypass the opposite overlay results by doing both intersection and difference for each overlay operation (then picking the relevant polygons based on intersects of their representative point). This is slow, but seemingly effective.

Dissappearing or closed polygons after unaryUnion is harder to guard against.

@dr-jts
Copy link
Contributor

dr-jts commented Sep 28, 2023

Removing slivers would solve some issues, but it will also leave gaps which would need dissolving later. Maybe an option.

Right, I see that in the dissolve_error2 dataset. There is a range of narrow polygons of varying areas, so would be hard deciding what threshold to remove at.

You seem to have the perfect dataset to expose these errors - unluckily for you! Too bad the original data isn't of higher quality.

@dr-jts
Copy link
Contributor

dr-jts commented Sep 28, 2023

Removing spikes is tricky. Are there any good ways to do this other than buffering (which slows things down way too much)?

Not really in GEOS at the moment. There are various algorithms around for removing spikes, which could be implemented.

In JTS the PolygonHullSimplifier can be used to remove spikes by computing an inner hull using a very small "area delta ratio" as the simplification parameter. (It would be even better if PolygonHullSimplifier supported a "minimum spike width" parameter to specifically remove spikes, without the risk of removing "nearly flat" vertices as well).

image

Unfortunately GEOSPolygonHullSimplify() doesn't expose this option.

@mortewle
Copy link
Author

The original data is of high quality, but we are combining a bunch of datasets, which results in a whole lot of edge cases. Lucky indeed!

I cannot really use JTS, since this is all done in Python. But intersesting to hear.

Off-topic, out of curiosity: why are JTS and GEOS separate? Do they have different use cases?

@dr-jts
Copy link
Contributor

dr-jts commented Sep 29, 2023

The original data is of high quality, but we are combining a bunch of datasets, which results in a whole lot of edge cases. Lucky indeed!

Ah, so the slivers are the result of an overlay?

Off-topic, out of curiosity: why are JTS and GEOS separate? Do they have different use cases?

No, same use case. As you have noticed, the main reason GEOS exists is that it is usable in a lot of places where JTS isn't. It is a separate project for many reasons: different technology, different set of contributors, different governance model. But the projects are pretty tightly linked. Some improvements start in JTS and then are ported into GEOS, and there's also things that go the other way.

@mortewle
Copy link
Author

Ah, so the slivers are the result of an overlay?

Yes, or about 10-15 overlays.

I came across another example where poygons dissappear in shapely and QGIS but not GEOSOP. Are there any changes from GEOS 3.12 to 3.13.0dev that solve some, but not all, of these issues, or is this a QGIS and shapely issue?

dissolve_error3.txt

@dr-jts
Copy link
Contributor

dr-jts commented Sep 29, 2023

I came across another example where polygons disappear in shapely and QGIS but not GEOSOP. Are there any changes from GEOS 3.12 to 3.13.0dev that solve some, but not all, of these issues, or is this a QGIS and shapely issue?

Right, I see that behaviour. I'm not sure what might have changed in 3.13dev to improve this. Possibly #937? Would have to bisect to see if that's the cause of the improvement. If so, that is an excellent result, for the cases it fixes.

@dr-jts
Copy link
Contributor

dr-jts commented Sep 29, 2023

Ah, so the slivers are the result of an overlay?

Yes, or about 10-15 overlays.

The dissolve_error3.txt dataset is FULL of very narrow sliver polygons. The red below shows all polygons with area < 0.0001:

image

It would be good to figure out a way to remove these prior to unioning. That sounds like another coverage cleaning operation (perhaps more tractable than full coverage cleaning, too).

@mortewle
Copy link
Author

It would be good to figure out a way to remove these prior to unioning. That sounds like another coverage cleaning operation (perhaps more tractable than full coverage cleaning, too).

I have made some functions for cleaning (eliminating slivers, closing gaps, closing small holes), but unary_union is often involved at some point, so it doesn't always fix this problem. Are there other ways to coverage clean?

@mortewle
Copy link
Author

mortewle commented Oct 11, 2023

@dr-jts Not to nag, but I'm just wondering if there will (or might, or will not) be a fix for this of any sort in the next few months?

I'm asking since we're publishing this dataset yearly in the spring, hopefully based on GEOS next year.

Meanwhile, I'm working on a coverage cleaning/snapping function that almost works, but it's not worth getting it perfect if this is added to JTS then GEOS shortly.

@dr-jts
Copy link
Contributor

dr-jts commented Oct 11, 2023

I'm just wondering if there will (or might, or will not) be a fix for this of any sort in the next few months?

I think it's unlikely there will be a fix for this in the next few months. It's pretty high priority, but it's also a tough problem to solve. If you have a potential solution in the works, then it's worth pursuing that.

@dr-jts dr-jts changed the title Polygon dissappears in unaryUnion Polygon disappears in unaryUnion Dec 19, 2023
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

No branches or pull requests

3 participants