nav[aria-label="Primary Navigation"] { padding: 0; & ul { list-style: none; width: 100%; display: flex; flex-direction: row; justify-content: start; align-items: start; gap: 30px; padding: 0; & li { margin: 0; } & ul li { list-style: none; } } }

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:

it is negative
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.

negative_intersection_area.zip

Vind meer berichten getagd met

Accepted answers

All comments

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.

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

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

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

Doe mee aan het gesprek

Aangemelde leden kunnen berichten plaatsen, updates volgen en meer. Nieuw hier? Registreer een gratis account.

Inloggen
Registreren

Snelle Links

Community Help
Find useful guides, FAQs, and documents to help you navigate and make the most of Esri Community.