Intersection-geometry has negative area

[CarstenSchumann]

I have a code that looks for parcels that touch another one:

IGeometry parcelGeometry = parcel.GetEsriGeometry();
ISpatialFilter filter = new SpatialFilter
{
    SpatialRel = esriSpatialRelEnum.esriSpatialRelIntersects, 
    Geometry = parcelGeometry
};

IFeatureCursor cursor = parcelFeatureClass.Search(filter, true);
IFeature feature;


while ((feature = cursor.NextFeature()) != null)
{
    var geom = (ITopologicalOperator) feature.Shape;
    var intersectGeom = geom.Intersect(parcelGeometry, esriGeometryDimension.esriGeometry2Dimension) as IArea;
}



(geom as ESRI.ArcGIS.Geometry.IRelationalOperator).Overlaps(parcel.GetEsriGeometry())

When I execute this chunk of code I get an area of intersection of -0.0000000014294918595190997 which is strange for two reasons:

1. it is negative
2. it is far beyond the actual cluster-tolerance (set to 0.1mm)

So I assume this is due to the inaccuarcy of type double.

As my actual purpose is on identifying touching geometries I initially used the esriSpatialRelTouches-operator instead of an intersect for the spatial filter which did not return any feature. esriSpatialRelOverlaps however returns the feature. Quite more weird is that when I filter the geometries in ArcMap using select features from ... that touch the boundary of the source-layer I also get the feature of which I thought it touches my source-feature.

When this problem comes due to the inaccuracy of double I wonder if there are any opportunities to trust the touch-operator anyway as it is a quite strict operator asking for an exact solution which double can´t (obviously) handle. However I´d expected that if an intersection is detected between two features that is far beyond the cluster-tolerance (and moreover is negative) the functions of ITopolgicalOperator are smart enough to filter those geometries out.

I provide a shape-file containing the two geometries but I doubt the problem is reproducible as I exported the geometries from an SDE.

EDIT: I also wrote this arcpy-snipped which gave me expected result:

geometries = []
i = 0
for row in arcpy.SearchCursor('LayerWithSelectedFeatures'):
    geometries.append(row.shp)
    i = i + 1
area = geometries[0].intersect(geometries[1], 2)
print area.area

This prints out 0.0 as I expected.

[DanPatterson_Retired]

negative geometry has been found to be associated with rings that go counter-clockwise instead of clockwise or for geometry with geometry errors like self-intersections and the like.  check your geometry

[CarstenSchumann]

However as the geometry with negative area comes directly from an ITopologicalOperator-method which is guaranteed to return only simple geometries (as by the docs, have to look further to get the source of this statement) the area must not be negative, mustn´t it?

However I checked both input-geometries if they are simple. To be safe I also set the IsKnownSimple-Property for both to false to enforce a Simplify which didn´t solve the problem. Anyway when I do the same with the result-geometry its area is zero as expected. This doesn´t neccessarily mean that the result-geometry per se is topologically incorrect as it is simply smaller than the CRS´ grid. However this also applies for the source-geometries which seem to be snapped to different coordinates.

[DanPatterson_Retired]

Robert Scheitlin, GISP​ may recollect the thread where ring order was important, I can only demonstrate this in the standard method of calculating area which is based, in simple form on Green's theorem (see the Shoelace formula)

Python as an example

cw clockwise, ccw counter...

def ring_area(x, y):
    """x and y as separate lists or arrays from references above
    - assumes first and last points are the same so uses slices
    """
    a0 = np.dot(x[1:], y[:-1])
    a1 = np.dot(x[:-1], y[1:])
    area = (a0-a1) * 0.5 # np.abs(a0-a1) * 0.5 for unsigned area
    print("\nx: {}\ny: {}\na0: {} a1: {}  area: {}".format(x,y,a0,a1,area))
    return area

>>> cw = np.array([[0,0], [0,10], [10,10], [10,0],[0,0]])
>>> xs = cw[...,0]
>>> ys = cw[...,1]
>>> ring_area(xs,ys)
x: [ 0  0 10 10  0]
y: [ 0 10 10  0  0]
a0: 200 a1: 0  area: 100.0
100.0
>>> ccw = np.array([[0,0],[10,0], [10,10], [0,10],  [0,0]])
>>> xs = ccw[...,0]
>>> ys = ccw[...,1]
>>> ring_area(xs,ys)
x: [ 0 10 10  0  0]
y: [ 0  0 10 10  0]
a0: 0 a1: 200  area: -100.0
-100.0
>>>

[CarstenSchumann]

I´ve updated my question. The arcpy-code that I´ve attached does not return any negative value, it simply retuns 0.0. So we´re back at ITopologicalOperator, I guess.

[RobertScheitlin__GISP]

Carsten,

Did you say that you explicitly called geom.Simplify() in your code at some point in your testing?

[CarstenSchumann]

Yeap, within the loop I call Simplify for every feature. Something like this:

while ((feature = cursor.NextFeature()) != null)  
{  
     var topo = (ITopologicalOperator2) feature.Shape;
     topo.IsKnownSimple_2 = false;
     topo.Simplify();
     var intersectGeom = topo.Intersect(parcelGeometry, esriGeometryDimension.esriGeometry2Dimension) as IArea;  
}  

Just for testing I did the same with my parcelGeometry and the actual intersectionGeometry. When simplifying the latter its area surely is 0.0 which makes me think that Intersect results in a non-simple geometry.

[RobertScheitlin__GISP]

What about casting the topo.Intersect result to ITopologicalOperator and Simplify and then get the area?

[CarstenSchumann]

I already checked this as mentioned in the commment from 23.05.2016 23:47. ​After simplifying the result-geometry its area is surely 0.0 which is strange a bit as it might mean that ITopologicalOperator produced some non-simply geometries. But I rather assume this is because the geometry being so small (remember: smaller than the grid itself) making it impossible for ITopo to process appropriately.