Area of overlap between two polygons

[JustinConklin]

I am trying to calculate the area (square feet) of overlap (intersection) of polygons in two shapefiles.  Layer 1 has only one polygon.  Layer two has many and I know that all of them intersect with the polygon in layer 1.  I am trying to calculate the area of the overlap.  My code loops through every feature in layer 2 and intersects it with the polygon in Layer 1 using iTopologicalOperator6.IntersectEx.  The problem is that the area it calculates is zero for most of the features.  It only gets an area for some.  As I mentioned, in reality, all polygons intersect.  My code is attached.  I verified that the loop is getting all features, so it seems the function using the .intersectex method has a problem. 

Any help is appreciated.  Thanks.
Justin

Protected Overrides Sub OnClick()                  Try              Dim pMap As IMap = My.ArcMap.Document.FocusMap             Dim pFLIndex As IFeatureLayer = pMap.Layer(0)             Dim pFLSections As IFeatureLayer = pMap.Layer(1)              'get selected features from Sections             Dim pFeatSel As IFeatureSelection = pFLSections             Dim pSelFeats As ISelectionSet = pFeatSel.SelectionSet              'MsgBox(pSelFeats.Count)              'Get feature from Index Layer             Dim pIndexFeature As IFeature = Nothing             Dim pCursor As ICursor = pFLIndex.Search(Nothing, False)              pIndexFeature = pCursor.NextRow              'MsgBox(pBufFeature.OID & "  " & pIndexFeature.OID)              'Loop through features in Buffered sections and calculate overlap sq ft             Dim pCur As IFeatureCursor = Nothing              Dim pBufFeature As IFeature = Nothing             pSelFeats.Search(Nothing, False, pCur)              pBufFeature = pCur.NextFeature              Do Until pBufFeature Is Nothing                 Dim dblOverLap As Double = GetArea(pBufFeature, pIndexFeature) 'Calls function                  pBufFeature.Value(10) = dblOverLap                 pBufFeature.Store()                  pBufFeature = pCur.NextFeature             Loop              My.ArcMap.Document.ActiveView.Refresh()           Catch ex As Exception             MsgBox(ex.ToString & vbCr & vbCr & ex.Message)         End Try          My.ArcMap.Application.CurrentTool = Nothing      End Sub           Private Function GetArea(ByVal FeatureA As IFeature, ByVal FeatureB As IFeature) As Double          Dim SourceArea As IArea         Dim TargetArea As IArea         Dim TopoOp As ITopologicalOperator6         Dim IntersectArea As IArea          SourceArea = FeatureA.ShapeCopy 'FeatureA and FeatureB are the polygons to be checked for overlap         TargetArea = FeatureB.ShapeCopy          TopoOp = TargetArea          IntersectArea = TopoOp.IntersectEx(SourceArea, False, ESRI.ArcGIS.Geometry.esriGeometryDimension.esriGeometry2Dimension)         GetArea = IntersectArea.Area '/ SourceArea.Area) * 100 'returns the percentage of first polygon which overlaps second polygon          SourceArea = Nothing         TargetArea = Nothing         TopoOp = Nothing         IntersectArea = Nothing      End Function

[NeilClemmons]

Do both geometries have the same spatial reference?  If not, then you should project one geometry into the spatial reference of the other before calling Intersect.

[JustinConklin]

Yes.  Spatial references are identical.

[AlexanderGray]

Any reason for using IntesectEx on ITopologicalOperator6 rather than Intersect of ITopologicalOperator? 
Another thing, in your search cursor for the layer 2, if you create a spatial filter with intersect relationship with the layer 1 geometry (single polygon), are you getting all the features you expect?

[JustinConklin]

Any reason for using IntesectEx on ITopologicalOperator6 rather than Intersect of ITopologicalOperator? 
Another thing, in your search cursor for the layer 2, if you create a spatial filter with intersect relationship with the layer 1 geometry (single polygon), are you getting all the features you expect?

I tried both IntersectEx of ITopologicalOperator6 and Intersect of ITopoligicalOperator and they produced the same results.  I was just getting the feature selection, but I tried your test and yes SpatialFilter intersect relationship will get all the features I expected.  So, the spatial filter recognizes there is an intersection, but topological operator does not.  Is there another interface that will get the geometry so i can calculate area?

Thanks for the help.

[AlexanderGray]

This is a shot in the dark but if you use ITopologicalOperator2 to set isKnownSimple_2 to false and then call simplify on the result geometry of the intersect, do you get an area?

[WeifengHe]

This is getting interesting.  Can you please share your data?  Not necessary all the features, just some that overlap but not generating area after the call.

Thanks.

[JustinConklin]

No.  I tried the simplify work around as well.  If I run the intersect geoprocessing tool in ArcMap same results.  Only gets some of the features.  I don't get why a select by location in ArcMap will get all of the intersecting features and programmatically a spatial filter will get them too, yet I can't get all the intersections in the GUI or through code with the intersect tools.  I have to do this for several datasets each containing several hundred intersecting polygons.  Maybe its just someting with this test dataset I am using?

[JustinConklin]

This is getting interesting.  Can you please share your data?  Not necessary all the features, just some that overlap but not generating area after the call.

Thanks.

Sure.  With this set both my code, and the intersect tool find only one of 11 to be an intersecting polygon.

Latest Code

Protected Overrides Sub OnClick()
        
        Try

            Dim pMap As IMap = My.ArcMap.Document.FocusMap
            Dim pFLIndex As IFeatureLayer = pMap.Layer(0)
            Dim pFLSections As IFeatureLayer = pMap.Layer(1)
            Dim pFCSections As IFeatureClass = pFLSections.FeatureClass

            'Get feature from Index Layer
            Dim pIndexFeature As IFeature = Nothing
            Dim pCursor As ICursor = pFLIndex.Search(Nothing, False)

            pIndexFeature = pCursor.NextRow

            'MsgBox(pBufFeature.OID & "  " & pIndexFeature.OID)

            Dim pSpatialFilter As New SpatialFilter
            pSpatialFilter.Geometry = pIndexFeature.Shape
            pSpatialFilter.SpatialRel = esriSpatialRelEnum.esriSpatialRelIntersects

            Dim pSelSet As ISelectionSet = pFCSections.Select(pSpatialFilter, esriSelectionType.esriSelectionTypeIDSet, esriSelectionOption.esriSelectionOptionNormal, Nothing)

            MsgBox(pSelSet.Count)
            'get selected features from Sections
            'Dim pFeatSel As IFeatureSelection = pFLSections
            'Dim pSelFeats As ISelectionSet = pFeatSel.SelectionSet

            'MsgBox(pSelFeats.Count)


            'Loop through features in Buffered sections and calculate overlap sq ft
            Dim pCur As IFeatureCursor = Nothing

            Dim pBufFeature As IFeature = Nothing
            pSelSet.Search(Nothing, False, pCur)

            pBufFeature = pCur.NextFeature

            Do Until pBufFeature Is Nothing
                Dim dblOverLap As Double = GetArea(pBufFeature, pIndexFeature) 'Calls function 
                pBufFeature.Value(10) = dblOverLap
                pBufFeature.Store()

                pBufFeature = pCur.NextFeature
            Loop

            My.ArcMap.Document.ActiveView.Refresh()


        Catch ex As Exception
            MsgBox(ex.ToString & vbCr & vbCr & ex.Message)
        End Try

        My.ArcMap.Application.CurrentTool = Nothing

    End Sub

    
    Private Function GetArea(ByVal FeatureA As IFeature, ByVal FeatureB As IFeature) As Double

        Dim SourceArea As IArea
        Dim TargetArea As IArea
        Dim TopoOp As ITopologicalOperator
        Dim IntersectArea As IArea

        SourceArea = FeatureA.ShapeCopy 'FeatureA and FeatureB are the polygons to be checked for overlap
        TargetArea = FeatureB.ShapeCopy

        TopoOp = TargetArea


        Dim TopoOp2 As ITopologicalOperator2 = TargetArea

        TopoOp2.IsKnownSimple_2 = False
        TopoOp2.Simplify()


        IntersectArea = TopoOp2.Intersect(SourceArea, ESRI.ArcGIS.Geometry.esriGeometryDimension.esriGeometry2Dimension)


        GetArea = IntersectArea.Area '/ SourceArea.Area) * 100 'returns the percentage of first polygon which overlaps second polygon


        'SourceArea = Nothing
        'TargetArea = Nothing
        'TopoOp = Nothing
        'IntersectArea = Nothing

    End Function

[NeilClemmons · Accepted solution]

All of the polygons in your export_output_3 shapefile have a negative area.  ArcGIS will consider these polygons "holes".  The area is negative because the polygons were created in the wrong direction.  You can fix the problem by looping through all of the features, getting the shape as IPolygon, calling IPolygon.ReverseOrientation, and setting the updated geometries back to the features.  You could also put a check into your current code to look at the area prior to performing the Intersect and reverse the orientation of the polygon if the area is negative.