Determine overlapping points using a spatial filter in ArcObject

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?

-
Can a point overlap? A simple python script could suffice; just loop through the layer and check all of the other points; around 10 lines of code. Output to screen/list/file. – Hairy Jan 25 '12 at 9:07
but i have to short out this problem by using spatial filter in C# arcobject and i dont know how to use spatial filter plz help me its urgent – Sachin Jan 25 '12 at 10:06
You don't have to do this in ArcObjects. Simply opening the data in Arcpy, you cna loop through, point by point, checking if they are the same, listing the ones which are. Arcpy is really the way to go with this, imo – Hairy Jan 25 '12 at 11:04
what you are saying is fine in arcpy,but i want to solve this problem only in arcobject C#, so please tell me how can i do this in arcobject C# – Sachin Jan 25 '12 at 11:30
"Using a spatial filter" seems like an unnecessary constraint, Sachin. (@Hairy has pointed out that a filter is not necessary.) Is this perhaps a problem for homework or from a textbook? – whuber Jan 25 '12 at 15:21

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 :)

-
I believe you make an excellent point but it seems to get lost in the presentation. (Most people won't be too worried about the difference between a O(n) and a O(n log(n)) algorithm because in practice the implicit constants are more important.) Could you possibly explain how the graph works and what it shows? Without explanatory labels and a legend, it could be displaying anything! – whuber Feb 29 '12 at 16:44
@whuber You are completely correct. A graph without labels simply sucks. I will edit the answer. – Ragi Yaser Burhum Feb 29 '12 at 18:53
Thanks. But what precisely is the difference between the light gray curve and the dark gray one? If I have understood you correctly, the light gray one should be a parabola, but at this scale it looks linear, so asymptotically the two algorithms appear to be the same, just with different constants. That's why such plots are usually made on a log-log-scale. It also seems confusing to refer to your approach as involving "comparisons" when in fact it uses hash table lookups (which require point-point comparisons only when resolving collisions, which should be rare). – whuber Feb 29 '12 at 19:06
@whuber OK, fixed it. Also, if you think about what the "implicit constants" are, you will realize that network I/O plays a huge role. So the time difference becomes even more impressive if you are in a client/server environment since the constant for fetching the features needs to be multiplied by every iteration (1 spatial query = 1 heavy constant cost for I/O through the cursor) – Ragi Yaser Burhum Feb 29 '12 at 19:09
I would love to see you include considerations like network I/O in your answer, because this issue of choosing an algorithm that can scale effectively comes up a lot (especially with ArcGIS). Could I persuade you to consider contributing a blog post about this general subject? – whuber Feb 29 '12 at 19:13

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

-
You should study the reply given by Ragi Yaser Burhum. – whuber Feb 29 '12 at 19:15
@whuber i have read that particular reply..i do not see any looping or use of spatial filter in my code..or is it that IIndexQuery2.NearestFeatures method itself is built on a spatial filter? – vinayan Mar 2 '12 at 7:28
I see at least one loop, vinayan (although I'll grant that every one of the points has to be inspected at least once!). One of the loops is explicit, although hidden by lack of indenting: it is delimited by "Do" and "Loop While Not pFtr is Nothing". The other is implicit in the call to "NearestFeatures". In the worst case this is a O(N) loop over all other features; in the best case, it will require a O(log(N)) search within a spatial data structure. Ragi's main point is that you don't need that search if you hash the point coordinates. – whuber Mar 2 '12 at 14:16
@whuber - now you really made me read the Ragi approach thoroughly..it does look very convincing..thanks for pointing out..let me see if i can post a result with hash approach and a loop approach.. – vinayan Mar 3 '12 at 7:55
vinayan Would share the hash approach code? I am trying to find out the identical features (in point, polyline and polygon)i.e. int hashCode = featureEnvlope.Envelope.GetHashCode(); But, it does not work as it should be. Thanks -jjy – user12665 Nov 10 '12 at 17:18

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);
``````
-
i have already gone through with that but it will be not working – Sachin Jan 25 '12 at 12:20
WHy doesn't it work? – Hairy Jan 25 '12 at 12:58
If I understand the question correctly, I think @Sachin wants to find Spatial Duplicates. – Jakub Jan 25 '12 at 16:17