topic Re: Finding furthest distances between two sets of points in Python Questions

Finding furthest distances between two sets of points

MarleneSanchez — Wed, 18 Mar 2015 17:14:14 GMT

I have two sets of points- list A and list B. I need to efficiently select the points from list B that are furthest from points in list A.

The reason Near Table isn't working out is because of the number of (redundant) records showing the distance of each point on list B to every single point on list A. That's 9000 entries (A) x 500 entries (B) for a total of 4,500,000 records generated by the near table.

I would like to end up with an ordered list of points from list B that are furthest away from any points on list A.

It seems like a simple thing to ask but not so simple to do.

How can I approach this?

Re: Finding furthest distances between two sets of points

ChrisDonohue__GISP — Wed, 18 Mar 2015 17:29:41 GMT

Some questions to help clarify this:

1. Is the issue technical - i.e. the analysis bogs down/software crashes due to the large amount of records?

2. Is the issue that their are many resulting Near records and uncertainty on how to query these results down to get the final answer?

Chris Donohue, GISP

Re: Finding furthest distances between two sets of points

MarleneSanchez — Wed, 18 Mar 2015 17:49:24 GMT

1. No technical issue. I did an experiment with a smaller data set. The application produces the table after some time. I really don't mind letting it run overnight to get 4.5 million records if that's what needs to be done.

2. Yes, there are way too many records produced. It's not the amount of them but their redundancy. Each point in B is matched to all points in A, multiplying the amount of records. I suspect that there is a way to leave just the highest-distance records in the near table for each point from the initial list of points B.

Re: Finding furthest distances between two sets of points

RichardFairhurst — Wed, 18 Mar 2015 18:00:57 GMT

It can be done with Python, but I don't know if I understand the end result you are looking for. Your way of describing what you want is very confusing. Here is what I think you want:

1. For any given point in A (shape) find the distance to the closest point in B along with the point B ObjectID (Attributes).

2. Sort the points in A by the distance to B in descending order so that the furthest distance point in A is at the top of the list.

Is that correct? Or perhaps rule 1 should be : For any given point in A find the distance to the furthest point in B along with the Point B ObjectID.

This assumes you want the point shape from A and just the attributes from B. But your description could mean that you want the to keep the point shape in B and attribute it with the distance and ObjectID of a point in A. Or do you want a table with the coordinates of both points that meet your rule 1 so that you could make an XY point layer based on either the location in A or the location in B?

So basically at the end should there be 9000 records (one for each point in A) or 500 records (one for each point in B).

Re: Finding furthest distances between two sets of points

MarleneSanchez — Wed, 18 Mar 2015 18:24:51 GMT

Richard, I think you got it right in your first description.

I must find the closest distances between points on separate lists, then sort them so the furthest distance is at the top of the list as you've said.

Here is a way to illustrate the problem.

Suppose I have a list of potential construction sites and a list of toxic waste facilities. I want to list the construction sites that are furthest away from all waste facilities. I'm showing each site with a point, not a polygon.

So I want to end up with the theoretical 500 sites sorted by distance to the nearest 9000 hazards, for a total of 500 records.

Re: Finding furthest distances between two sets of points

DarrenWiens2 — Wed, 18 Mar 2015 18:40:36 GMT

Not sure I completely understand, but to get the nearest distance between construction and toxic sites, use Near (outputs only nearest features) rather than Generate Near Table (potentially outputs all combinations). Then, sort the feature class by near distance.

Re: Finding furthest distances between two sets of points

JoshuaBixby — Wed, 18 Mar 2015 19:16:37 GMT

I think part of what is confusing people is semantics, particularly the use of the word "furthest." Assuming no two construction sites are exactly the same distance from a waste facility, you will only have one construction site that is furthest away. If by furthest you mean far, then a definition of what is far and what is close would be helpful. Or, are you looking to simply determine the distance from each construction site to each waste facility?

Re: Finding furthest distances between two sets of points

MarleneSanchez — Wed, 18 Mar 2015 19:23:37 GMT

I think this is it. Please let try this and get back to you.

Re: Finding furthest distances between two sets of points

JoshuaBixby — Wed, 18 Mar 2015 19:38:32 GMT

Are you after nearest or farthest points? You say "furthest" a lot but then you also say "closest" at one point.

Re: Finding furthest distances between two sets of points

DarrenWiens2 — Wed, 18 Mar 2015 19:39:57 GMT

My interpretation is that she's looking for "the construction site farthest from the nearest toxic waste facility," which makes sense.

Re: Finding furthest distances between two sets of points

MarleneSanchez — Wed, 18 Mar 2015 19:42:01 GMT

Thank you everyone for all your help.

Since there is a limited number sites and facilities, we can order all the sites by distance from their closest waste facilities.

In practice you would not really have identical distances between points.

So at the very top of the list you would have one construction site that you could say is further from all environmental hazards than any other site you have on the list, right?

Re: Finding furthest distances between two sets of points

JoshuaBixby — Wed, 18 Mar 2015 19:47:15 GMT

OK, I get it, you want to find the nearest hazard site to each construction site, and then rank/order all the construction sites in descending order of that value. So, the construction site at the top has the largest minimum distance to a hazard site.

Re: Finding furthest distances between two sets of points

MarleneSanchez — Wed, 18 Mar 2015 20:04:53 GMT

That's it.

I realize that my original question may have been a bit off.

I was originally using Near Table and was trying to efficiently get rid of the non-relevant records. (Non-relevant meaning those that do not help identify points with the largest minimum distance to a hazard)

I did't realize I could do this with just the Near tool as Darren has suggested.

Re: Finding furthest distances between two sets of points

BlakeTerhune — Sat, 11 Dec 2021 16:12:21 GMT

I've done something similar and the Near table worked for me. When I was doing this, I also found it frustrating that the Near table would only show the ObjectID of the features, requiring additional joins to actually see what features the OIDs were for. I created this Python script that will create a Near table with a key field of your choice in the output. I haven't tested it much but maybe you can get some use out of it.

import arcpy
import os

try:
    # Set local variables
    sourceGDB = r"C:\temp\mytempsource.gdb"
    outputGDB = r"C:\temp\mytempoutput.gdb"
    ## IN and NEAR features can be the same
    in_fc = "FeatureClassIn_Name"
    in_fc_keyfield_name = "KeyField_Name"
    in_fc_keyfield_type = "Long"
    in_fc_path = os.path.join(sourceGDB, in_fc)
    near_fc = "FeatureClassNear_Name"
    near_fc_keyfield_name = "KeyField_Name"
    near_fc_keyfield_type = "Long"
    near_fc_path = os.path.join(sourceGDB, near_fc)
    search_radius = "25 Feet"

    # Create near table and make as TableView
    neartable_inmem = os.path.join("in_memory", "neartable_inmem")
    arcpy.GenerateNearTable_analysis(
        in_fc_path,  ## in_features
        near_fc_path,  ## near_features
        neartable_inmem,  ## out_table
        search_radius,  ## search_radius
        "NO_LOCATION",  ## location
        "NO_ANGLE",  ## angle
        "ALL",  ## closest
        "0",  ## closest_count
        "PLANAR"  ## method
    )
    arcpy.MakeTableView_management(neartable_inmem, "neartable_view")
    print "Near table view created"

    # Add new IN and NEAR key fields to near table
    in_kf = "IN_{}".format(in_fc_keyfield_name)
    arcpy.AddField_management(
        "neartable_view",  ## in_table
        in_kf,  ## field_name
        in_fc_keyfield_type,  ## field_type
    )
    near_kf = "NEAR_{}".format(near_fc_keyfield_name)
    arcpy.AddField_management(
        "neartable_view",  ## in_table
        near_kf,  ## field_name
        near_fc_keyfield_type,  ## field_type
    )
    print "Key fields added"

    # Join and calc IN key field
    arcpy.AddJoin_management(
        "neartable_view",  ## in_layer_or_view
        "IN_FID",  ## in_field
        in_fc_path,  ## join_table
        "OBJECTID",  ## join_field
        "KEEP_COMMON"  ## join_type
    )
    arcpy.CalculateField_management(
        "neartable_view",  ## in_table
        in_kf,  ## field
        "!{}!".format(in_fc_keyfield_name),  ## expression
        "PYTHON_9.3"  ## expression_type
    )
    arcpy.RemoveJoin_management("neartable_view")
    print "{} field populated".format(in_kf)

    # Join and calc NEAR key field
    arcpy.AddJoin_management(
        "neartable_view",  ## in_layer_or_view
        "NEAR_FID",  ## in_field
        near_fc_path,  ## join_table
        "OBJECTID",  ## join_field
        "KEEP_COMMON"  ## join_type
    )
    arcpy.CalculateField_management(
        "neartable_view",  ## in_table
        near_kf,  ## field
        "!{}!".format(near_fc_keyfield_name),  ## expression
        "PYTHON_9.3"  ## expression_type
    )
    arcpy.RemoveJoin_management("neartable_view")
    print "{} field populated".format(near_kf)

    # Export final near table
    kf_near_table = "{}_NEAR".format(in_fc.replace('.', '_'))  ## Replace any dot with underscore and add _NEAR suffix
    kf_near_table_path = os.path.join(outputGDB, kf_near_table)
    arcpy.CopyRows_management(
        "neartable_view",  ## in_rows
        kf_near_table_path  ## out_table
    )
    print "Final near table output at {}".format(kf_near_table_path)

except Exception as err:
    print err
    ## If there was an ArcPy error, write all of the messages returned by the last tool
    if arcpy.GetMessages(2):
        print arcpy.GetMessages()

finally:
    # Cleanup
    arcpy.Delete_management("in_memory")
    arcpy.ClearWorkspaceCache_management()

Re: Finding furthest distances between two sets of points

MarleneSanchez — Wed, 18 Mar 2015 22:45:52 GMT

Good god. I have no experience with Python yet, so I'm going to take this down and chew on it later ^^

Re: Finding furthest distances between two sets of points

BlakeTerhune — Wed, 18 Mar 2015 22:51:10 GMT

Haha! Well, you did post in the Python section so I thought I would go for it. It's not as complicated as it looks.

Create Near table
Add two new fields for in and near key fields
Join the original feature classes to field calc the in and near key field values
Export final near table

Re: Finding furthest distances between two sets of points

JoshuaBixby — Sat, 11 Dec 2021 16:12:24 GMT

If you were interested in the farthest instead of nearest, then you would have to use Point Distance or Generate Near Table in your workflow, unless you wanted to role your own NumPy-based function. When working with data sets in the tens of thousands, which much from a practical stand point, I have found the Point Distance tool to be very slow and sometimes Generate Near Table fails with unspecified errors.

For situations where I want farthest points or quicker outputs similar to Point Distance, I have roled my own NumPy-based functions. For example:

def stat_point(in_features, other_features, stat='MINIMUM'):
    import arcpy
    import numpy
   
    stats = {
        'MINIMUM': {'FID': 'MIN_FID',
                    'DIST': 'MIN_DIST',
                    'INDEX': lambda x: 0},
        'MAXIMUM': {'FID': 'MAX_FID',
                    'DIST': 'MAX_DIST',
                    'INDEX': lambda x: x - 1},
        'MEDIAN_HIGH': {'FID': 'MEDH_FID',
                        'DIST': 'MEDH_DIST',
                        'INDEX': lambda x: x / 2},
        'MEDIAN_LOW': {'FID': 'MEDL_FID',
                    'DIST': 'MEDL_DIST',
                    'INDEX': lambda x: x / 2 - 1 if x % 2 == 0 else x / 2 }
    }
   
    desc = arcpy.Describe(in_features)
    SR = desc.spatialReference
   
    desc = arcpy.Describe(other_features)
    OID_name_other = desc.OIDFieldName
    shape_name_other = desc.ShapeFieldName
    narr_other = arcpy.da.FeatureClassToNumPyArray(
        other_features,
        [OID_name_other, desc.ShapeFieldName],
        spatial_reference = SR
    )
    xy_other = narr_other[shape_name_other]
    idx = stats[stat]['INDEX'](numpy.shape(xy_other)[0])
   
    arcpy.AddField_management(in_features, stats[stat]['FID'], 'LONG')
    arcpy.AddField_management(in_features, stats[stat]['DIST'], 'DOUBLE')
   
    with arcpy.da.UpdateCursor(
        in_features,
        ["SHAPE@XY", stats[stat]['FID'], stats[stat]['DIST']]
    ) as cur:
        for xy, FID, DIST in cur:
            xy_in = numpy.array(xy)
            d0 = numpy.subtract.outer(xy_in[0], xy_other[:,0])
            d1 = numpy.subtract.outer(xy_in[1], xy_other[:,1])
            dist = numpy.hypot(d0, d1)
            i = dist.argsort()[idx]
            FID = narr_other[OID_name_other]
            DIST = dist
            cur.updateRow([xy, FID, DIST])

When used with "MINIMUM", the stat_point function mimics the Near tool. The Near tool is more performant than using stat_point to find the minimum distance point, but stat_point is orders of magnitude faster at finding non-nearest points compared to Point Distance or Generate Near Table. Of course, stat_point is doing planar measurements so that is something to be aware of. I threw MEDIAN parameters in as an example of how simple it is to extend this methodology beyond minimum or maxiumum points.

The NumPy heavy lifting, lines 43-45, comes from Alex Martelli in response to a Stack Overflow post: Euclidean distance between points in two different Numpy arrays, not within...

Re: Finding furthest distances between two sets of points

DarrenWiens2 — Tue, 24 Mar 2015 21:27:04 GMT

I should probably read this more closely, but I've used scipy.spatial.distance.pdist for approximately what I think may be happening in lines 43-45 (pairwise distances?).

Re: Finding furthest distances between two sets of points

JoshuaBixby — Tue, 24 Mar 2015 21:41:11 GMT

I stuck with NumPy because SciPy isn't shipped with ArcGIS Desktop, yet, and SciPy isn't part of our standard data center Python deployment either. From a performance perspective, it would be interesting to compare the two approaches.

Re: Finding furthest distances between two sets of points

DarrenWiens2 — Sat, 11 Dec 2021 16:12:27 GMT

Just ran a quick test between stat_point() and the following function for 1,000 points and the results are: stat_point() - 29s, numpyNear() - 0.5s. Almost of the stat_point() time is due to calls to arcpy.Describe (14.1s) and arcpy.AddField (14.7s), so the numpy stuff isn't the problem.

... def numpyNear(fc):
...     npArray = arcpy.da.FeatureClassToNumPyArray(fc,["OID@","SHAPE@XY"])
...     distArray = scipy.spatial.distance.pdist(npArray['SHAPE@XY'],'euclidean')
...     sqDistArray = scipy.spatial.distance.squareform(distArray)
...     nearFID = npArray['OID@'][numpy.argsort(sqDistArray)[:,1]].transpose()
...     npAppend = numpy.lib.recfunctions.append_fields(npArray,'NearFID',nearFID)
...     nearDist = numpy.sort(sqDistArray)[:,1].transpose()
...     npAppend = numpy.lib.recfunctions.append_fields(npAppend,'NearDist',nearDist)
...     outFC = r'in_memory\outPts'+str(random.random()*1000)
...     outFeat = arcpy.da.NumPyArrayToFeatureClass(npAppend,outFC,"SHAPE@XY")
...     arcpy.CopyFeatures_management(outFC,r'in_memory\numpy_Pts')