Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I have a script that calculates the average distance between points. What I want is to rewrite it so I can input the point layer, input a distance, and the output field will tell me a count of the number of points at that distance. Can anyone help? so basically I just want a script that will tell me a count of the number of points from a point at a certain distance.

#Import standard library modules
import sys, string, arcpy

    inputFC = arcpy.GetParameterAsText(0)# Input Point Layer
    inputdistance = arcpy.GetParameterAsText(1) #input distance
    outputfield = arcpy.GetParameterAsText(1) #output field

# Read the point data into list.  
#
    xs = []
    ys = []
 rowCursor = arcpy.SearchCursor(inputFC)
 for row in rowCursor:
    geom = row.getValue (properties.shapeFieldName)
    cen = geom.centroid
    xs.append(cen.X)
    ys.append(cen.Y)
 del rowCursor

 # compute distance 
 sumdist = 0.0
 count = 0.0
 for i in range(0,len(xs)):
    for j in range(i,len(xs)):
        if (i <> j):
           dist = ((xs[i] -xs[j])**2 + (ys[i] - ys[j])**2)**0.5
           sumdist = sumdist + dist
           count = count + 1
  print "The average inter-point distance is " + str(sumdist/count)
  arcpy.AddMessage("The average inter-point distance is " + str(sumdist/count))

With the following script. Ive tried to run it as Distance.add_nearby_points_count_field('gauges.shp', 100)

 import sys
 import arcpy
 import math

def add_nearby_points_count_field(inputFC, inputdistance):

# Read the point data into list.
coord_pairs = []
desc = arcpy.Describe(inputFC)
rowCursor = arcpy.SearchCursor(inputFC)
for row in rowCursor:
    geom = row.getValue(desc.shapeFieldName)
    coord_pairs.append(geom.centroid)

# Add a field called 'cnt' and calculate nearby points
arcpy.AddField_management(inputFC, 'cnt', 'LONG') #i also tried cutting   out this part and creating the "cnt" field manually
rowCursor = arcpy.UpdateCursor(inputFC)
for row in rowCursor:
    geom = row.getValue(desc.shapeFieldName)
    from_point = geom.centroid
    near_pts = 0
    for to_point in coord_pairs:
        distance = math.sqrt(pow((to_point.X - from_point.X), 2) + pow((to_point.Y - from_point.Y), 2))
        if distance <= inputdistance:
            near_pts += 1
    row.cnt = near_pts - 1 # Subtract 1 to remove the measurement to itself
    rowCursor.updateRow(row)
del rowCursor

if __name__ == '__main__':
    inputFC ='guages.shp' #originally its a featureclass but I converted it to a shapefile to try both methods
    inputdistance = 100
    add_nearby_points_count_field(inputFC, inputdistance)
share|improve this question
    
Are you looking for points that are exactly the specified distance from the input point, near that specified distance (i.e. with tolerance), or all points within or beyond that specified distance (buffer/intersect operation)? –  Allan Adair Sep 29 '11 at 8:51
    
what I have here is "points within the vicinity of other points". So I'm assuming if you have an input distance of say 10km, it will give you all the points within that distance (10km and less) –  Nevu Sep 29 '11 at 9:06
    
I know buffer seems like an obvious choice but they don't want it done that way –  Nevu Sep 29 '11 at 9:23
    
So for each point feature of a feature class, you will want a list of all other point features within a certain distance from the same originating feature class? –  Allan Adair Sep 29 '11 at 9:31
    
thats close, but rather an output field in the layer that will show a count of how many points are in the vicinity of that point. –  Nevu Sep 29 '11 at 10:47
show 5 more comments

1 Answer

up vote 0 down vote accepted

EDIT based on comment...

Kinn, you can use this in a number of ways. One way is to save the script as "count_nearby_points.py", and as long as it is in your python path, you can use it in any other script by importing it.

import count_nearby_points

Then you should be able to call the count_nearby_points.add_nearby_points_count_field(inputFC, inputdistance) function from within your script.

Or, if you run this script directly, you can modify the parameters under if __name__ == '__main__':

You can also copy the part(s) that you want into your own script.

Feel free to ask if you have further questions.

#-------------------------------------------------------------------------------
# Name:        count_nearby_points.py
#
# Purpose:     Return a count of nearby point features
#
# Created:     29/09/2011
#-------------------------------------------------------------------------------
#!/usr/bin/env python

import sys
import arcpy
import math

def add_nearby_points_count_field(inputFC, inputdistance):

    # Read the point data into list.
    coord_pairs = []
    desc = arcpy.Describe(inputFC)
    rowCursor = arcpy.SearchCursor(inputFC)
    for row in rowCursor:
        geom = row.getValue(desc.shapeFieldName)
        coord_pairs.append(geom.centroid)

    # Add a field called 'cnt' and calculate nearby points
    arcpy.AddField_management(inputFC, 'cnt', 'LONG')
    rowCursor = arcpy.UpdateCursor(inputFC)
    for row in rowCursor:
        geom = row.getValue(desc.shapeFieldName)
        from_point = geom.centroid
        near_pts = 0
        for to_point in coord_pairs:
            distance = math.sqrt(pow((to_point.X - from_point.X), 2) + pow((to_point.Y - from_point.Y), 2))
            if distance <= inputdistance:
                near_pts += 1
        row.cnt = near_pts - 1 # Subtract 1 to remove the measurement to itself
        rowCursor.updateRow(row)
    del rowCursor

if __name__ == '__main__':
    # Test data for when running as __main__, make sure that 'points.shp' exists.
    inputFC ='points.shp'
    inputdistance = 100
    add_nearby_points_count_field(inputFC, inputdistance)
share|improve this answer
    
thanks! that looks great. Only problem is when I run it I get an error "<type 'exceptions.TypeError'>: add_nearby_points_count_field() takes exactly 2 arguments (4 given)". But I remove the if_name part the script completes but nothing happens. –  Nevu Sep 29 '11 at 15:32
    
@Kinn, I made a small change to the code and edited my response with some supporting information. Please mark my answer as the solution if it fits your needs. Thanks. –  Allan Adair Sep 29 '11 at 16:05
    
hmm I've tried to run it as a script and from the command line and I get a 99999 Error executing function. –  Nevu Sep 30 '11 at 3:31
    
@Kinn, make sure there isn't a schema lock on the inputFC (i.e. open in ArcMap or selected in ArcCatalog), and that there isn't already a field called "cnt" in inputFC. –  Allan Adair Sep 30 '11 at 5:08
    
thanks, yes I checked for both of those before. –  Nevu Sep 30 '11 at 8:12
show 1 more comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.