I have a point layer with thousands of points. Using a spatial filter, I have to find out those points which are intersecting each other. Then the objectid's of the intersecting points should be displayed in a list box.
I'm using ArcObjects and C# and need a solution in this environment.
How should I approach this problem?
Spatial Filter is not the way to go. Think about this.
Since in ArcObjects (unlike in Oracle Spatial or PostGIS) there is no way to do a real self join, you are left with having to create a Spatial Query for each point feature and check for intersection. n spatial queries where n is the number of features!
If you did this your time complexity would be O(n) (the cost for looping through every feature) * O(log n) (the cost for searching through the index - and at "log n", I am actually being too kind because the ESRI index is a multi-level grid - not an Rtree).
To get an idea how slow it degrades over time, take a look at this graph.
X is the number of points. Y is the number of comparisons (worst case) using the SpatialFilter approach. 250 points, with the approach I describe below, will be done in 250 comparisons. With the spatial filter approach, the number of comparisons shoots out of the chart (i.e, as the problem size increases, it gets slower really really fast).
If you are not shocked yet, you should be.
Another example is that with 1000 point features, the approach below will be done in 1000 comparisons. With the SpatialFilter approach, you would need ~9965 comparisons. That is 9x slower.
Don't use the spatial filter approach for finding duplicates.
Add a 200,000 features and prepare to wait for several hours (or even days) for something that should take minutes instead.
Better approach is to:
Easy approach, much better than what you would do without having the luxury of joins.
Yes, the space complexity is O(n) as opposed to O(1) for your first approach, but in this case I would not be concerned and we could get into all kinds of other discussions that would last forever :)
Another approach is to to create a cursor over your point data so you pass through your point dataset only once.
For each point get the X and Y coordinates and turn these into a string and concatenate into a single string thus a point whose x is 11.6 and y is 12.8 would become a string of "11.6,12.8" .
Then add this to a dictionary, the key is the XY string and the item is a collection object with ObjectID's. But when you add it to a dictionary you check if it exists already, if if does then return the item (which is the collection object) and add the objectID to the collection and reinsert back into the dictionary.
At the end of the run you will have a dictionary where the key is the XY string and the item is a collection of 1 or more objectID's. You can then step through the dictionary and display only the collections that have 2 or more entries.
No spatial filter used and it can be implemented in C# and ArcObjects.
i did this using IIndexQuery2 Interface in ArcObjects. Code is in VBA. I think you could easily convert this to C#. I find this quite fast. Please correct me if i am doing this wrong..
Option Explicit
Private Sub DuplicatePoints()
On Error GoTo EH
Dim m_IQ2 As IIndexQuery2
Dim pFI2 As IFeatureIndex2
Set pFI2 = New FeatureIndex
Dim pMxDoc As IMxDocument
Set pMxDoc = ThisDocument
Dim pFLayer As IFeatureLayer
Set pFLayer = pMxDoc.FocusMap.Layer(0)
Set pFI2.FeatureClass = pFLayer.FeatureClass
Set pFI2.OutputSpatialReference("Shape") = pMxDoc.FocusMap.SpatialReference
pFI2.Index Nothing, pFLayer.AreaOfInterest
Set m_IQ2 = pFI2
Dim pAV As IActiveView
Set pAV = pMxDoc.FocusMap
Dim pPoint As IPoint
Dim pFtrCursor As IFeatureCursor
Dim pFC As IFeatureClass
Set pFC = pFLayer.FeatureClass
Set pFtrCursor = pFC.Search(Nothing, False)
Dim pFtr As IFeature
Set pFtr = pFtrCursor.NextFeature
Do
Set pPoint = pFtr.Shape
Dim lOID As Long, dDist As Double
Dim objIDs As Variant
Dim distNearest As Double
m_IQ2.NearestFeatures pPoint, objIDs, distNearest
If UBound(objIDs) > 0 Then
Debug.Print "ObjectID:" & pFtr.OID & " has duplicates!"
End If
Set pFtr = pFtrCursor.NextFeature
Loop While Not pFtr Is Nothing
Exit Sub
EH:
MsgBox Err.Description
End Sub
Here's an example of a C# spatialFilter with an intersect. Supplied from the excellent Esri docs. You should be able to adjust it for your use..
[C#]
// Get the feature and its geometry given an Object ID.
IFeature stateFeature = stateFeatureClass.GetFeature(14);
IGeometry queryGeometry = stateFeature.Shape;
// Create the spatial filter. "highwayFeatureClass" is the feature class containing
// the highway data. Set the SubFields property to "FULL_NAME", as only that field is
// going to be displayed.
ISpatialFilter spatialFilter = new SpatialFilterClass();
spatialFilter.Geometry = queryGeometry;
spatialFilter.GeometryField = highwayFeatureClass.ShapeFieldName;
spatialFilter.SpatialRel = esriSpatialRelEnum.esriSpatialRelIntersects;
spatialFilter.SubFields = "FULL_NAME";
// Find the position of the "FULL_NAME" field in the highway feature class.
int nameFieldPosition = highwayFeatureClass.FindField("FULL_NAME");
// Execute the query and iterate through the cursor's results.
IFeatureCursor highwayCursor = highwayFeatureClass.Search(spatialFilter, false);
IFeature highwayFeature = null;
while ((highwayFeature = highwayCursor.NextFeature()) != null)
{
String name = Convert.ToString(highwayFeature.get_Value(nameFieldPosition));
Console.WriteLine("Highway found: {0}", name);
}
// The cursors is no longer needed, so dispose of it.
Marshal.ReleaseComObject(highwayCursor);
