# Performing a 'far analysis' using ArcGIS for Desktop?

I would like to perform an analysis that is the inverse of `Near_analysis` for ArcGIS. Instead of my result being `NEAR_DIST` and `NEAR_FID`, I'd like each feature to have a `FAR_DIST` and `FAR_FID` populated with the furthest feature information from a 'far feature'. I can do this with some scripting, but it would be slow, and involve iterating through each feature in a feature class and performing a Near_analysis on each, and picking the furthest distance.

Any other ideas? I'm ArcGIS 10.1 and proficient in python/arcpy.

• For 2 point layers (note that you can use copy of layer as 2nd) use point distance. Sort output table in descending order by distance. Remove duplicates in input_fid Dec 23 '14 at 22:08
• Indexing, by its very nature, is optimized to find "close" things. It's unlikely you'll ever find an efficient way to find "far" ones. There is no shortcut to the furthest possible point. Dec 23 '14 at 22:13

## 3 Answers

This code worked for me. I modified a near point distance tool I wrote last year. To find max distance, just change the dist_type parameter from min to max:

``````import arcpy, os, sys, traceback, math
from os import path as p
from datetime import datetime as d
arcpy.env.overwriteOutput = True

def Message(msg):
print str(msg)
arcpy.AddMessage(msg)

def findDistance(a,b):
x = abs(a[0] - b[0])
y = abs(a[1] - b[1])
return math.sqrt((x*x) + (y*y))

def nearPoints(in_points, near_points, out_table, dist_type=min):
"""function to find the minimum or maximimum distance between points

Required:
in_points -- source points for near analysis
near_points -- points to find distance to (can be same as in_points)
out_table -- output near table

Optional:
dist_type -- set to min or max to find the closest or farthest distnace
between points.
"""

# fix args if coming from script tool
if str(dist_type).lower() == 'max':
dist_type = max
else:
dist_type = min

startTime = d.now()
# grab xy coords
with arcpy.da.SearchCursor(in_points, ['OID@','SHAPE@XY']) as rows:
point_dict = dict((r[0],r[1]) for r in rows)

# grab xy coords near points
with arcpy.da.SearchCursor(near_points, ['OID@','SHAPE@XY']) as rows:
npoint_dict = dict((r[0],r[1]) for r in rows)

# create dictionary to find nearest point
same = in_points == near_points
near_dict = {}
for key in point_dict.keys():
this_pt = point_dict[key]
distList = {}
for oid,coords in npoint_dict.iteritems():
distList[oid] = findDistance(this_pt,coords)
if same:
closest = dist_type(filter(None, distList.values()))
else:
closest = dist_type(distList.values())
near_id = [k for k,v in distList.items() if v==closest][0]
near_dict[key] = [near_id,closest]
del point_dict, distList

# create output table
path = os.path.dirname(out_table)
name = os.path.basename(out_table)
if arcpy.Exists(out_table):
arcpy.Delete_management(out_table)
arcpy.TableToTable_conversion(in_points,path,name)
arcpy.AddField_management(out_table,'NEAR_ID','LONG')
arcpy.AddField_management(out_table,'NEAR_DIST','DOUBLE')
fields = ['OID@','NEAR_ID','NEAR_DIST']
with arcpy.da.UpdateCursor(out_table,fields) as rows:
for row in rows:
if row[0] in near_dict:
row[1] = near_dict[row[0]][0]
row[2] = near_dict[row[0]][1]
rows.updateRow(row)
Message('Created: %s' %os.path.basename(out_table))
Message('(Elapsed time: %s)' %(str(d.now() - startTime)[:-3]))
return out_table

if __name__ == '__main__':

##    # stand alone
##    pts = r'C:\TEMP\carlton_cities3.shp'
##    out = r'C:\TEMP\carlton_max_dist.dbf'
##    nearPoints(pts, pts, out, max)

# Get Args
argv = tuple(arcpy.GetParameterAsText(i) for i in range(arcpy.GetArgumentCount()))

# Run it
nearPoints(*argv)
``````

If you have access to an ArcInfo license, you could use the Point Distance tool, which:

Determines the distances from input point features to all points in the near features.... The tool creates a table with distances between two sets of points. if the default search radius is used, distances from all input points to all near points are calculated. The output table can be quite large. For example, if both input and near features have 1,000 points each, then the output table can contain one million records.

Run this with no search radius, then filter to find the pairs with the greatest distance.

• If the asker's license level is an issue then "Point Distance for Basic" code is available from gis.stackexchange.com/questions/113446/…
– PolyGeo
Dec 23 '14 at 22:13
• @PolyGeo handy code. You could post that into a new answer, and do an analysis within the cursor, only writing the distance if it's greater than the current distance Dec 23 '14 at 22:16
• Excellent idea but I don't see me finding time at the moment - you (or anyone else) is welcome to use it as starting code for a "Far Points" tool.
– PolyGeo
Dec 23 '14 at 22:29

I've created a python toolbox which performs a `far analysis` as described above.

Inputs:

• Input Feature Class - The input feature class
• Far Feature Class - The feature class with far feature candidates
• Output Geodatabase - The geodatabase that will house the output table
• Output Table Name - The name of the output table
• Distance Determination Method - Indicates if a far feature's distance will be determined by its nearest point or its farthest point
• True Curve Vertices - The number of (equally spaced) vertices used to represent any true curves in the far features.

The tool creates a geodatabase table:

Advanced license is required if `FARTHEST POINT OF FARTHEST FEATURE` is selected for Distance Determination Method. The code:

``````print "importing"
import arcpy
import os

def FeatureClassToPoints (input_fc, sr, addedPoints):
desc = arcpy.Describe (input_fc)
shape = desc.shapeType
di = {}
with arcpy.da.SearchCursor (input_fc, ["OID@", "SHAPE@"], "", sr) as curs:
for oid, geom in curs:
lines = None
pnts = set ()
if shape == "Polygon":
lines = arcpy.FeatureToLine_management (geom, arcpy.Geometry ())
elif shape == "Polyline":
lines = [geom]
if lines:
for line in lines:
if "curve" in line.JSON:
vPnts = arcpy.FeatureVerticesToPoints_management (line, arcpy.Geometry ())
spLines = arcpy.SplitLineAtPoint_management (line, vPnts, arcpy.Geometry (), "1 feet")
for spLine in spLines:
if "curve" in spLine.JSON:
segLen = spLine.length / addedPoints + 1
for i in range (1, addedPoints + 1):
pnt = spLine.positionAlongLine (segLen * i)
pnts.add (arcpy.Geometry ("point", pnt.firstPoint, sr))
else:
for array in spLine:
for pnt in array:
pnts.add (arcpy.Geometry ("point", pnt, sr))
else:
for array in line:
for pnt in array:
pnts.add (arcpy.Geometry ("point", pnt, sr))
elif shape == "Multipoint":
for pnt in geom:
pnts.add (arcpy.Geometry ("point", pnt, sr))
elif shape == "Point":
pnts.add (geom)
else:
print "Unsupported shape type for feature class {}".format (input_fc)
arcpy.AddMessage ("Unsupported shape type for feature class {}".format (input_fc))
print "{} feature classes are not supported""".format (shape)
arcpy.AddMessage ("{} feature classes are not supported""".format (shape))
return {}
di [oid] = pnts
return di

class Toolbox(object):
def __init__(self):
"""Define the toolbox (the name of the toolbox is the name of the
.pyt file)."""
self.label = "Custom Toolbox"
self.alias = ""

# List of tool classes associated with this toolbox
self.tools = [Tool]

class Tool(object):
def __init__(self):
"""Define the tool (tool name is the name of the class)."""
self.label = "Far (analysis)"
self.description = "Identifies the farthest feature from an input features to another layer or feature class."
self.canRunInBackground = False

def getParameterInfo(self):

inFc = arcpy.Parameter (
displayName = "Input Feature Class",
name = "infc",
datatype = "GPFeatureLayer",
parameterType = "Required",
direction = "Input")

farFc = arcpy.Parameter (
displayName = "Far Feature Class",
name = "farfc",
datatype = "GPFeatureLayer",
parameterType = "Required",
direction = "Input")

outGdb = arcpy.Parameter(
displayName = "Output Geodatabase",
name = "outgdb",
datatype = "Workspace",
parameterType = "Required",
direction = "Input")

outName = arcpy.Parameter(
displayName = "Output Table Name",
name = "outname",
datatype = "GPString",
parameterType = "Required",
direction = "Input")

farPnt = arcpy.Parameter(
displayName = "Distance Determination Method",
name = "farpoint",
datatype = "GPString",
parameterType = "Required",
direction = "Input")

farPnt.filter.list = ["FARTHEST POINT OF FARTHEST FEATURE",
"NEAREST POINT OF FARTHEST FEATURE"]

curvePnts = arcpy.Parameter(
displayName = "True Curve Vertices",
name = "calc",
datatype = "GPLong",
parameterType = "Required",
direction = "Input")

outTab = arcpy.Parameter (
displayName = "Out Table",
name = "outtab",
datatype = "DETable",
parameterType = "Derived",
direction = "Output")

return [inFc, farFc, outGdb, outName, farPnt, curvePnts, outTab]

def isLicensed(self):
"""Set whether tool is licensed to execute."""
return True

def updateParameters(self, parameters):
"""Modify the values and properties of parameters before internal
validation is performed.  This method is called whenever a parameter
has been changed."""
# 0 input feature class
# 1 far feature class
# 2 output gdb
# 3 output file name
# 4 calculation type
# 5 true curve vertices
# 6 output feature class
gdb = arcpy.env.workspace
if not parameters [2].altered:
parameters [2].value = gdb
if not parameters [4].altered:
parameters [4].value = "FARTHEST POINT OF FARTHEST FEATURE"
if not parameters [5].altered:
parameters [5].value = 3
if parameters [5].altered and parameters [5].value < 0:
parameters [5].value = abs (parameters [5].value)
if parameters [3].altered and parameters [3].valueAsText and parameters [2].valueAsText:
parameters [3].value = os.path.basename (arcpy.ValidateTableName (parameters [3].valueAsText.strip (),
parameters [2].valueAsText))
if parameters [0].altered and not parameters [3].valueAsText and parameters [2].valueAsText:
outName = os.path.basename ("{}_Far".format (parameters [0].valueAsText))
outName = os.path.basename (arcpy.CreateUniqueName (outName, parameters [2].valueAsText))
parameters [3].value = os.path.basename (arcpy.ValidateTableName (outName,
parameters [2].valueAsText))

def updateMessages(self, parameters):
"""Modify the messages created by internal validation for each tool
parameter.  This method is called after internal validation."""

try:
if parameters [2].valueAsText:
if arcpy.Describe (parameters [2].valueAsText).workspaceType == "FileSystem":
parameters [2].setErrorMessage ("Output location is not a geodatabase.")
except: parameters [2].setErrorMessage ("Output location is not a geodatabase.")

def execute(self, parameters, messages):
"""The source code of the tool."""

inFc = parameters [0].valueAsText
farFc = parameters [1].valueAsText
outGdb = parameters [2].valueAsText
outName = parameters [3].valueAsText
calc = parameters [4].value
addedPoints =  parameters [5].value
if not addedPoints: addedPoints = 3
outTab = os.path.join (outGdb, outName)

print "determining spatial reference"
arcpy.AddMessage ("determining spatial reference")
sr = None
inSr = arcpy.Describe (inFc).spatialReference
if inSr.linearUnitName:
sr = inSr
else:
farSr = arcpy.Describe (farFc).spatialReference
if farSr.linearUnitName:
sr = farSr
if not sr:
print "one input feature class must have a linear unit spatial reference"
arcpy.AddMessage ("one input feature class must have a linear unit spatial reference")

oidDi = {}
if sr and calc == "FARTHEST POINT OF FARTHEST FEATURE":
print "getting far points"
arcpy.AddMessage ("getting far points")
farPntDi = FeatureClassToPoints (farFc, sr, addedPoints)
distDi = {}
print "finding far features"
arcpy.AddMessage ("finding far features")
total = arcpy.GetCount_management (inFc) [0]
arcpy.SetProgressor ("step", "finding far features...", 0, int (total))
with arcpy.da.SearchCursor (inFc, ["OID@", "SHAPE@"], "", sr) as curs:
for i, (inOid, inGeom) in enumerate (curs):
arcpy.SetProgressorPosition (i)
for farOid, farPnts in farPntDi.items ():
for farPnt in farPnts:
try: lastDist = distDi [inOid]
except KeyError:
distDi [inOid] = inGeom.distanceTo (farPnt)
oidDi [inOid] = farOid
continue
dist = inGeom.distanceTo (farPnt)
if dist > lastDist:
distDi [inOid] = dist
oidDi [inOid] = farOid
arcpy.ResetProgressor()

elif sr:
arcpy.AddMessage ("getting far geometries")
farGeomDi = {}
with arcpy.da.SearchCursor (farFc, ["OID@", "SHAPE@"], "", sr) as curs:
for oid, geom in curs:
farGeomDi [oid] = geom

arcpy.AddMessage ("finding far features")
oidDi = {}
distDi = {}
total = arcpy.GetCount_management (inFc) [0]
arcpy.SetProgressor ("step", "finding far features...", 0, int (total))
with arcpy.da.SearchCursor (inFc, ["OID@", "SHAPE@"], "", sr) as curs:
for i, (inOid, inGeom) in enumerate (curs):
arcpy.SetProgressorPosition (i)
for farOid, geom in farGeomDi.items ():
try: lastDist = distDi [inOid]
except KeyError:
distDi [inOid] = inGeom.distanceTo (geom)
oidDi [inOid] = farOid
continue
dist = inGeom.distanceTo (geom)
if dist > lastDist:
distDi [inOid] = dist
oidDi [inOid] = farOid
arcpy.ResetProgressor()

if oidDi:
print "creating table"
arcpy.AddMessage ("creating table")
outPath, outName = os.path.split (outTab)
outTab = arcpy.CreateUniqueName (outName, outPath)
outPath, outName = os.path.split (outTab)
outTab = arcpy.CreateTable_management (outPath, outName) [0]
for fld in ["INPUT_OID", "FAR_OID"]:
print "adding field", fld
arcpy.AddMessage ("adding field {}".format (fld))
arcpy.AddField_management (outTab, fld, "LONG")

with arcpy.da.InsertCursor (outTab, ["INPUT_OID", "FAR_OID"]) as curs:
for inOid, farOid in oidDi.items ():
row = (inOid, farOid)
curs.insertRow (row)
print
print "created:"
arcpy.AddMessage ("\ncreated:")
print outTab
arcpy.AddMessage (outTab)
parameters [6].value = outTab
``````