Performing a 'far analysis' using ArcGIS for Desktop?

Question

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.

Answer 1

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)

Answer 2

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.

Answer 3

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