Problem doing Intersection and Difference on small polygons

[Keld Laursen]

Hello.

First, thank you for your effort in creating NetTopologySuite. We have used it for some time, and it looks very sturdy and fits the bill very much for us.

I have run into somewhat of a snag, though. We are looking at farmers field polygons, and we have some polygons representing various legislation based no-go zones.

In order to correctly calculate areas for the field polygons, we subtract the areas of the no-go zones from the field polygons. Locally, the no-go zones are represented as only the zone area that overlaps the field, which some times can be very small areas.

When we calculate the area to subtract, we have to remove areas of a specific no-go zone that has already been covered by another no-go zone. From time to time, this work ends us with very small multipolygons as a remainder. Some times these remains wreak havoc with the calculations, thus ending with exceptions being thrown.

I have created an example (Accompanying .zip file), built on version 1.14 of NetTopologySuite. My problem were found in a larger application, using NTS version 1.7

The basic code is like this:

1 static void Main(string[] args)

2 {

3 var pol1String= "MULTIPOLYGON (((636192.59787309519 6154321.699591

4 var pol2String = "MULTIPOLYGON (((636167.79101204267 6154298.27038

5 var polygon1 = (MultiPolygon) new WKTReader().Read(pol1String);

6 var polygon2 = (MultiPolygon) new WKTReader().Read(pol2String);

7 var result = polygon1.Difference(polygon2);

8 Console.WriteLine(result.Area);

9 Console.ReadKey();

10 }

11 }

When I go and find the difference in line 7, the program throws an error. The error stems from one of the polygons in multipolygon polygon1. Asking the polygon if it is valid yields true.

I have looked at the offending polygon, and the sides are VERY close to each others, thus yielding very small areas to consider.

Can you confirm whether or not the problem is my polygon (and suggest ways to better the polygon without hurting the area of it)?

Can you suggest other ways to allow the program to do the correct calculation?

If none of the above: Is this a bug in NTS?

MS SQL server geometry can do the above calculations with no errors.

Best regards,

Keld Laursen

[FObermaier]

This is an issue in both NTS and JTS. Most likely a precision thing.

What you could try is whenever you get this error, insert the topology exception location in the geometries by adding it at the apropriate place or replace it when necessary.

If I do that on your geometries, I get POLYGON EMPTY as result.

What does SqlServer return?

[FObermaier]

I'm sorry, I overlooked a small glich in the code

I replaced polygon2 with polygon1, sorry. Once you fix that, you get the following results:

WKT : MULTIPOLYGON (((636184.8487384459 6154321.8321918845, 636192.59787309519 6154321.6995911133, 636185.751571696 6154321.8162046, 636184.84907004656 6154321.8315769676, 636182.89523056 6154321.86485684, 636184.8490700468 6154321.8315769676, 636184.8487384459 6154321.8321918845)), ((636211.11390913755 6154314.0155097777, 636207.43681213306 6154321.4456718238, 636207.437031418 6154321.4456680715, 636211.114677974 6154314.0171801187, 636211.11390913755 6154314.0155097777)), ((636203.68238669413 6154308.2698675757, 636208.89554388 6154309.19807967, 636209.84867722925 6154311.2667209692, 636208.89654585847 6154309.1981610162, 636203.68238669413 6154308.2698675757)), ((636152.56918352435 6154316.2767154472, 636152.56918352365 6154316.2767155087, 636152.581231977 6154315.22303587, 636152.56918352435 6154316.2767154472)))

Area: 0,0103460097064597

Am Donnerstag, 17. März 2016 14:48:56 UTC+1 schrieb Keld Laursen:

[Keld Laursen]

Thank you for this. I would never have come up with that set of InsertExceptionTopologyPoint functions on my own. I couldn't recognise the coordinate, so I couldn't determine where to stick the coordinate.

I will try and use your code in my program and see what happens.

For reference, the result from SQL server is:

MULTIPOLYGON (((636192.59787309519 6154321.6995911133, 636184.84873844613 6154321.8321918845, 636184.8490700468 6154321.8315769676, 636185.751571696 6154321.8162046, 636192.59787309519 6154321.6995911133)), ((636211.11390913755 6154314.0155097777, 636211.114677974 6154314.0171801187, 636207.437031418 6154321.4456680715, 636207.43681213306 6154321.4456718238, 636211.11390913755 6154314.0155097777)), ((636203.68238669413 6154308.2698675757, 636208.89654585847 6154309.1981610162, 636209.84867722925 6154311.2667209692, 636208.89554388 6154309.19807967, 636203.68238669413 6154308.2698675757)))

Area: 0,010498046875

Best regards,

Keld List Laursen

--
You received this message because you are subscribed to a topic in the Google Groups "NetTopologySuite" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/nettopologysuite/6YTDe0pMMgI/unsubscribe.
To unsubscribe from this group and all its topics, send an email to nettopologysui...@googlegroups.com.
To post to this group, send email to nettopol...@googlegroups.com.
Visit this group at https://groups.google.com/group/nettopologysuite.
For more options, visit https://groups.google.com/d/optout.

[FObermaier]

Another option is to reduce the precision of the geometries:

//2nd Attempt

var gpr = new GeometryPrecisionReducer(new PrecisionModel(10000000000));

var p1 = gpr.Reduce(polygon1);

var p2 = gpr.Reduce(polygon2);

result = null;

Assert.DoesNotThrow(() => result = p1.Difference(p2));

area = 0;

Assert.DoesNotThrow(() => area = result.Area);

Assert.AreEqual(0.01025390625, area, 0.01);

Console.WriteLine("WKT : {0}", result.AsText());

Console.WriteLine("Area: {0}", area);

This works well, too.

Am Donnerstag, 17. März 2016 14:48:56 UTC+1 schrieb Keld Laursen: