Split ... Sort .... Slice ... the Top X in Y

Blog Post created by Dan_Patterson Champion on Jan 28, 2019


Identify the 2 highest elevation points in a series of polygons.



To win a bet amongst neighbors... to locate something like a tower... find observation points for visibility analysis... lots of reasons



  • Intersect the points with the polygons, retaining all attributes... output type... point
  • Add an X and Y field to the result so you have the point coordinates for later use (if you don't have them)
  • Delete any 'fluff' fields that you don't need, but keep the ID, KEY, OBJECTID fields from each. They will have a unique name for identification.
  • Split the result into groupings using the key fields associated with the polygons (see code)
  • Sort the groupings on a numeric field, like elevation (ascending/descending)     (code)
  • Slice the number that you need from each grouping.   (code)
  • Send it back out to a featureclass if you need it.    (code)


For Picture People...

A sample table resulting from the intersection of points with polygons.


The points (in red) and the polygons (a grid pattern) with the two points with the highest 'Norm' value in each polygon


An upclose look of the result


For the Pythonistas....

Fire up Spyder or your favorite python IDE


Table to NumPy Array....

>>>  a = arcpy.da.TableToNumPy("path to the featureclass or table", field_names="*")  # you now have an array


with me so far?  now you have an array


Make a script... call it something (splitter.py for example)


The work code... 

import numpy as np
import arcpy

def split_sort_slice(a, split_fld=None, val_fld=None, ascend=True, num=1):
    """Does stuff  Dan_Patterson@carleton.ca 2019-01-28

    def _split_(a, fld):
        """split unsorted array"""
        out = []
        uni, idx = np.unique(a[fld], True)
        for _, j in enumerate(uni):
            key = (a[fld] == j)
        return out
    err_0 = "The split_field {} isn't present in the array"
    if split_fld not in a.dtype.names:
        return a
    subs = _split_(a, split_fld)
    ordered = []
    for i, sub in enumerate(subs):
        r = sub[np.argsort(sub, order=val_fld)]
        if not ascend:
            r = r[::-1]
        num = min(num, r.size)
    out = np.asarray(ordered)
    return out

# ---- Do this ----

in_fc = r"C:\path_to_your\Geodatabase.gdb\intersect_featureclass_name"

a = arcpy.da.TableToNumPyArray(in_fc, "*")

out = split_sort_slice(fn, split_fld='Grid_codes', val_fld='Norm', ascend=False, num=2)

out_fc = r"C:\path_to_your\Geodatabase.gdb\output_featureclass_name"

arcpy.da.NumPyArrayToFeatureClass(out, out_fc, ['Xs', 'Ys'], '2951')


Lines 40 - 42

NumPyArrayToFeatureClass—Data Access module | ArcGIS Desktop 


Open ArcGIS Pro, refresh your database and add the result to your map.


I will put this or a variant into...


Point Tools …

So that you can click away, for those that don't like to type.