Linearly comparing one point to all points within 20' foot area, for every point in dataset using ArcPy?

Question

I have a very large dataset of GPS point that I am trying to write a python script for that will linearly compare one points to the ones around it. Original I was using the ID to determine which point to compare to the next one but do to how the points was collected the ID are not always in order, results in a lot of extra points being selected.

I believe the best solution would be to compare one point to all points within a 20’ foot buffer area of it. I am not sure how to best write a code that does this, I was thinking of creating an in-memory buffer of all points and then using that to select the points that I want to compare to each other.

Here is part of my original script that I am using now (this not using a buffer but an ID# to compare points):

JONum = arcpy.GetParameterAsText(0)
JO = "{}" "{}" "{}" "{}".format('JONumber = ',"'",JONum,"'")
OutPut = arcpy.GetParameterAsText(1)
arcpy.MakeFeatureLayer_management('GPSPipePoint', 'Temp', JO, OutPut)
arcpy.AddMessage('MN Feature Layer')

layer = 'GPS_GasMain_Temp'
fields =  ['ID', 'NAME', 'SIGNAL', 'DEPTH', 'Z_COORD','OPERATORNAME','COLLECTIONDATE']
all_records = [i for i in arcpy.da.SearchCursor(layer,fields)] 
all_records = sorted(all_records, key=lambda x: (x[6],x[7],x[1]))

ids_to_select2 = []
for index in range(1,len(all_records)-1): #Skip first and last record
    if (abs(all_records[index][2] >= (1.2 * all_records[index-1][2]))):
        ids_to_select2.append(all_records[index][0]) #If record >+-20% 

sql2 = """{0} IN({1})""".format(arcpy.AddFieldDelimiters(datasource=layer, field=fields[0]), ','.join([str(i) for i in ids_to_select2])) #Build a query using ids_to_select list. For example 'GISID' IN(2023,2024)

Answer 1

You can Spatial Join the Points to themselves using WITHIN_A_DISTANCE and a search radius. This method should be fast since spatial join usually is. Multiple select by location is often slow.

Then use Collections.defaultdict(list) to store each Point id as key and all Closest Point ids as values. It is unclear what comparison you want to do but you could select by attributes or use dictionaries of stored attributes etc. to compare the attributes.

import arcpy
from collections import defaultdict
arcpy.env.overwriteOutput = True

fc = r"X:\sve_1milj_Sweref_99_TM_shape\svk\riks\bs_riks.shp" #Change
search_distance = '10 Kilometers' #Change

memjoin = r'in_memory\spatjoin'
arcpy.SpatialJoin_analysis(target_features=fc, join_features=fc,
                           out_feature_class=memjoin, join_operation='JOIN_ONE_TO_MANY',
                           match_option='WITHIN_A_DISTANCE', search_radius=search_distance)
closestdict = defaultdict(list)
with arcpy.da.SearchCursor(memjoin,['TARGET_FID','JOIN_FID']) as cursor:
    for target, join in cursor:
        if target!=join:
            closestdict[target].append(join)

For example, Point with objectid 0 has these 6 Points within the distance:

closestdict[0]
[5, 17, 18, 363, 562, 2223]

To calculate mean of the closest points you can do:

all_values = {k:v for k,v in arcpy.da.SearchCursor(fc,['OID@','somevaluefield'])}

for key,values in closestdict.iteritems(): #.items() in python 3
    mean_value = sum([all_values[value] for value in values])/len(values)
    print('Mean value of the points closest to oid {0} is {1}'.format(key, mean_value))

Answer 2

In this type of analysis it can be nice to jump out of nested cursors and buffers and work directly on the data with a spatial index. In doing similar work I've found that using the rtree library cuts the processing time in half. The only catch is that you have to be willing and able to install the rtree package.

import math
from rtree import index

# Grab a list of all the points.
# Each record should be [[x,y], id]

all_records = [i for i in arcpy.da.SearchCursor(layer, ["SHAPE@XY", "ID"])]

# Lets make an output that's a dictionary of each point id,
# storing a list of nearby point ids.
# Each record would be something like 15: [23, 47, 124]

output = {}

# rtree indexes work on square bounding boxes, 
# so to filter by the cartesian distance
# we'll set up a distance function to use
def distance(p0, p1):
    return math.sqrt((p0[0] - p1[0])**2 + (p0[1] - p1[1])**2)

# Create an index and fill with the points, storing the object record
idx = index.Index()
for i, record in enumerate(all_records):
    x,y = record[0]
    idx.Insert(i, [x, y, x, y], record)

# Search the index for every record
for record in all_records:
    source_point, source_id = record
    # find everything in a 20 foot square
    x,y = source_point
    results = idx.intersection([x-10.0, y-10.0, x+10.0, y+10.0], objects=True)
    # with objects=True, the results are the records we stored
    for result in results:
        result_point, result_id = result
        # skip the point we're searching
        if source_id == result_id:
            continue
        # skip points that aren't within 20 feet 
        if distance(source_point, result_point) > 20.0:
            continue
        # populate the output dictionary
        if source_id not in output:
            output[source_id] = [result_id]
        else:
            output[source_id].append(result_id)

Answer 3

Update: I wasn't aware that the tool was under the advanced license. This can still be done easily with a loop and a select by location. The near tool is the same tool for the standard license as well.

Depending on what you're trying to do you either just need select by location with WITHIN_A_DISTANCE or you need to just put that inside of a loop.

1. If you just need to find ALL points in a layer that are within 20m of any point in another layer just use a single select by location.

2. If you truly need distances of specific IDs from each specific point you just need to loop it.

Make a layer of your points you need distances for. Make another layer of points you're wondering if they are within 20m of. Selections are feature layer specific so you'll need to have two even if it's conceptually the same feature class.

SearchCursor your first layer selecting each point by uid then select by location on your other layer. Then extract UIDs of that selection to store IDS that are within a distance of your point. You could do this with a dictionary easily {ID:[list of ids within 20m]}

3. If you need to know all the points that are within 20m of all other points progressively you could just make a single selection of a point, and add to the selection where it is within 20m of selection. If the selection count goes up, do that again, and select all points that are within 20m of your selection. Repeat until that count doesn't increase.

It should be noted that this will not work well with a random seed point, as it may just be an outlier. This logic can be used to cluster points though if you repeat the logic selecting points that haven't already been clustered.

[ORIGINAL ADVICE/Advanced license] Arc makes a tool for this: http://desktop.arcgis.com/en/arcmap/10.3/tools/analysis-toolbox/point-distance.htm

Just run that on a selection of points to get other points within 20 meters. Then just iterate through that selection to do whatever it is you're trying to do.

This tool is actually deprecated sorry. The tool you want is Near and you should be able to use it with a standard license.