Dan_Patterson

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

Blog Post created by Dan_Patterson Champion on Jan 28, 2019

Task...

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

 

Purpose...

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

 

Conceptualization...

  • 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)
            out.append(a[key])
        return out
    #
    err_0 = "The split_field {} isn't present in the array"
    if split_fld not in a.dtype.names:
        print(err_0.format(split_fld))
        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)
        ordered.extend(r[:num])
    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.

Outcomes