1

Given a set of points representing villages stored in a .shp file, is there a way to compute the pair-wise distance (Euclidean) between them? What I would like to come up with is a sort of symmetric matrix to be further analysed statistically.

I am using ArcGIS 10.1 with an Advanced license.

3 Answers 3

3

Alternatively you can use the Generate Near Table tool to do the same thing as GISGe's suggestion, but with some additional attributes (such as xy of the near feature and the angle between source and near features) if those would be useful in your study. It also allows you to specify a number of nearest matches so it isn't all or just nearest.

2

Use the Point Distance tool to do this, with your point shapefile as input features and near features parameters.

enter image description here

10
  • Thank. Do you mean using the same shape file as input in the two fields? I tried this solution, but arcGIS does not complete the task, reporting "999999 : Error executing function". Any idea?
    – NewAtGis
    Jun 2, 2015 at 6:33
  • Yes that's what I mean, how many points do you have in your shapefile? And what is its coordinate system?
    – GISGe
    Jun 2, 2015 at 6:41
  • The number of points is 15. As for the second question, here are the details: Projected Coordinate System: ED_1950_UTM_Zone_33N; Projection: Transverse_Mercator
    – NewAtGis
    Jun 2, 2015 at 8:18
  • Weird, can you upload the shapefile somehow so we can make a test? Do you have the same issue with the Generate Near Table tool mentioned by @ChrisW?
    – GISGe
    Jun 2, 2015 at 8:32
  • Hello. I am wondering how can I upload something here....
    – NewAtGis
    Jun 3, 2015 at 7:46
1

I'm sorry, I would reply to previous comments but I don't have 50 points on my profile (new member) and I'm getting blocked. Since your "Input" and "Near" feature classes are both the same, why not make a copy of your original FC and use the original as the "Input" and the copy as the "Near" for the tool? This might circumvent the issues you're having with using the same feature class as Input/Near.

The Near tools and Point Distance tools require an Advanced license, so it seems you have the necessary license level. For those who have Basic/Standard, ESRI's Near algorithm isn't too hard to implement yourself in Python. This answer is similar to one that I posted in Calculate nearest distance in ArcView, XY coordinates. But since you want all pairwise comparisons (not just the nearest feature), the tool will create a table with all pairwise comparisons. It will also create a "plumb lines" feature class that connects all pairs and has the same info as the table (this output is more useful IMO). For the Target and Source ID Fields, you can use "OBJECTID" or any other unique ID of your choice. It will work on any license level. If you don't have luck with the above suggestions, try this script.

__author__ = "John K. Tran"
__contact__ = "[email protected]"
__version__ = "1.0"
__created__ = "7/1/15"

"""Compares all target point features to all source point features and populates
the ID and distance of all pairwise comparisons between each target to each source.
Can optionally create plumb lines to visualize the links between source and target featues."""

import arcpy
import os
import math
import random

arcpy.SetProgressor("Default", "Firing up script...")

# Set up initial parameters for script.
source = arcpy.GetParameterAsText(0) # Point feature class.
target = arcpy.GetParameterAsText(1) # Point feature class.
outputtable = arcpy.GetParameterAsText(2) # Table feature class that stores distance and ID info.
sourceIDfield = arcpy.GetParameterAsText(3) # Field to identify unique ID of Source values.
targetIDfield = arcpy.GetParameterAsText(4) # Field to identify unique ID of Target values.
createplumblines = arcpy.GetParameterAsText(5) # Optional output location for plumb lines. HIGHLY suggest you specify a path for this!
arcpy.env.parallelProcessingFactor = "100%"

#Convert arcpy field types for "Field Type" argument in AddField_management.
fieldtypedict = {"Blob":"BLOB", "Date":"DATE", "Double":"DOUBLE", "Guid":"GUID", "Integer":"LONG", "Raster":"RASTER", "Single":"FLOAT", "SmallInteger":"SHORT", "String":"TEXT", "OID":"LONG"}

# Get the field type for the source ID field, so we can make analogous field for NearID in target.
sourcefields = arcpy.ListFields(source)
targetfields = arcpy.ListFields(target)
sourceIDfieldtype = fieldtypedict[[field.type for field in sourcefields if field.name == sourceIDfield].pop()]
targetIDfieldtype = fieldtypedict[[field.type for field in targetfields if field.name == targetIDfield].pop()]

# Get the spatial references for source and target and make sure they're the same.
sourcespatialref = arcpy.Describe(source).spatialReference
targetspatialref = arcpy.Describe(target).spatialReference
assert sourcespatialref.name == targetspatialref.name, "The source and target must have the same coordinate system. Ensure both data are in an identical projected coordinate system."

# Create a Point class with a distance method. More efficient than native arcpy.Geometry class.
class Point(object):
    def __init__(self, x, y, ID):
        """Create the constructor to initialize the point."""
        self.x = x
        self.y = y
        self.ID = ID
    def GetDist(self, OtherPoint):
        """Determines the distance between two points using Pythagorean Theorem."""
        return math.sqrt((self.x-OtherPoint.x)**2+(self.y-OtherPoint.y)**2)

# Create a list of source points to compare with each target point.
scount = 0
sourcepoints = list()
with arcpy.da.SearchCursor(source, ["SHAPE@XY", sourceIDfield]) as sourcecursor:
    for row in sourcecursor:
        arcpy.SetProgressorLabel("Gathering Source Points: {0} features complete".format(str(scount)))
        sourcepoint = Point(row[0][0], row[0][1], row[1])
        sourcepoints.append(sourcepoint)
        scount += 1
del sourcecursor

# Create several startup variables if creating plumb lines was enabled.
if createplumblines:
    arcpy.SetProgressorLabel("Making Plumb Line Feature Class")
    plumblines = list()
    arcpy.CreateFeatureclass_management(os.path.dirname(createplumblines), os.path.basename(createplumblines), "POLYLINE", None, "DISABLED", "DISABLED", targetspatialref)
    arcpy.AddField_management(createplumblines, "SourceID", sourceIDfieldtype)
    arcpy.AddField_management(createplumblines, "TargetID", targetIDfieldtype)

# Create a table to store distance info.
arcpy.SetProgressorLabel("Making Distance Table")
arcpy.CreateTable_management(os.path.dirname(outputtable), os.path.basename(outputtable))
arcpy.AddField_management(outputtable, "SourceID", sourceIDfieldtype)
arcpy.AddField_management(outputtable, "TargetID", targetIDfieldtype)
arcpy.AddField_management(outputtable, "Distance", "DOUBLE")

# Create a list to store each comparison, so we can add it into the table afterwards.
distlist = []

# Loop through each target and compare with each source and populate ID and Distance fields in table.
count = 0
with arcpy.da.SearchCursor(target, ["SHAPE@XY", targetIDfield]) as targetcursor:
    for row in targetcursor:
        arcpy.SetProgressorLabel("Computing Distance: {0} features complete".format(str(count)))
        targetpoint = Point(row[0][0], row[0][1], row[1])
        for sourcepoint in sourcepoints:
            currentdist = targetpoint.GetDist(sourcepoint)
            distlist.append((sourcepoint.ID, targetpoint.ID, currentdist))
            if createplumblines:
                plumblines.append((sourcepoint, targetpoint))
        count = count + 1
del targetcursor

# Update the table with the info from the distlist.
tnum = 0
with arcpy.da.InsertCursor(outputtable, ["SourceID", "TargetID", "Distance"]) as tablecursor:
    for row in distlist:
        arcpy.SetProgressorLabel("Adding distances to table: {0} features complete".format(str(tnum)))
        tablecursor.insertRow(row)
del tablecursor

# Create the plumb lines if the option was enabled.
if createplumblines:
    plnum = 0
    plumbcursor = arcpy.da.InsertCursor(createplumblines, ["SHAPE@", "SourceID", "TargetID"])
    for sourcepnt, targetpnt in plumblines:
        arcpy.SetProgressorLabel("Generating plumb lines: {0} features complete".format(str(plnum)))
        plumbarray = arcpy.Array([arcpy.Point(sourcepnt.x, sourcepnt.y), arcpy.Point(targetpnt.x, targetpnt.y)])
        plumbline = arcpy.Polyline(plumbarray)
        plumbcursor.insertRow((plumbline, sourcepnt.ID, targetpnt.ID))
        plnum += 1
    del plumbcursor

# Done!
arcpy.ResetProgressor()

Let me know if you have any issues running the script. It works for point-to-point geometry comparisons. A more challenging exercise is writing an algorithm to compare the higher geom types to each other.

Best of luck!

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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