Point Clustering w/ ArcGIS + Python

Here is a quick example of clustering an arcgis feature class of points.  WARNING: This is more of a "gist" than a production ready version.  I am releasing it to support my presentation at the Esri Developer Summit.  Come back soon and it should be updated with additional error handling, much much more efficient use of the tools.

import json
import multiprocessing
import os
import unittest

from collections import defaultdict

import arcpy

def cluster_worker(sourcePointFeatureClass, outputDirectory, levelInfos=None, whereClause=None, workerName=None, includeLevels=None):
    print '.rectify_worker({})'.format(workerName)
    arcpy.env.overwriteOutput = True
    outputGeodatabaseName = "temp{}.gdb".format(workerName)
    outputGeodatabase = os.path.join(outputDirectory, outputGeodatabaseName)
    if arcpy.Exists(outputGeodatabase):
        arcpy.Delete_management(outputGeodatabase)
    arcpy.CreateFileGDB_management(outputDirectory, outputGeodatabaseName)
    arcpy.env.workspace = outputGeodatabase

    points_feature_layer = 'temp_points'
    points_feature_class = 'source_points'
    
    arcpy.MakeFeatureLayer_management(sourcePointFeatureClass, points_feature_layer, whereClause)
    arcpy.CopyFeatures_management(points_feature_layer, points_feature_class)

    level_infos_json = """{"origin":{"x":-20037508.342787,"y":20037508.342787},"spatialReference":{"wkid":102100},"lods":[{"level":0,"resolution":156543.033928,"scale":591657527.591555},{"level":1,"resolution":78271.5169639999,"scale":295828763.795777},{"level":2,"resolution":39135.7584820001,"scale":147914381.897889},{"level":3,"resolution":19567.8792409999,"scale":73957190.948944},{"level":4,"resolution":9783.93962049996,"scale":36978595.474472},{"level":5,"resolution":4891.96981024998,"scale":18489297.737236},{"level":6,"resolution":2445.98490512499,"scale":9244648.868618},{"level":7,"resolution":1222.99245256249,"scale":4622324.434309},{"level":8,"resolution":611.49622628138,"scale":2311162.217155},{"level":9,"resolution":305.748113140558,"scale":1155581.108577},{"level":10,"resolution":152.874056570411,"scale":577790.554289},{"level":11,"resolution":76.4370282850732,"scale":288895.277144},{"level":12,"resolution":38.2185141425366,"scale":144447.638572},{"level":13,"resolution":19.1092570712683,"scale":72223.819286},{"level":14,"resolution":9.55462853563415,"scale":36111.909643},{"level":15,"resolution":4.77731426794937,"scale":18055.954822},{"level":16,"resolution":2.38865713397468,"scale":9027.977411},{"level":17,"resolution":1.19432856685505,"scale":4513.988705},{"level":18,"resolution":0.597164283559817,"scale":2256.994353},{"level":19,"resolution":0.298582141647617,"scale":1128.497176}],"units":"esriMeters"}"""
    level_infos = json.loads(unicode(level_infos_json, 'latin-1'))
    clustered_feature_classes = []
    for level_info in level_infos['lods']:
        if not includeLevels or int(level_info['level']) in includeLevels:
            cluster_feature_class = create_cluster_layer(points_feature_class, level_info)
            if cluster_feature_class:
                clustered_feature_classes.append(os.path.join(outputGeodatabase, cluster_feature_class))
    return (clustered_feature_classes, workerName, outputGeodatabase)

def create_cluster_layer(sourcePointFeatureClass, levelInfo):
    level, res = (levelInfo['level'], levelInfo['resolution'])

    integrate_name = 'integrate_{}'.format(level)
    cluster_name = 'clusters_{}'.format(level)
    output_name = 'clusters_{}_joined'.format(level)

    if not arcpy.Exists(output_name):
        arcpy.CopyFeatures_management(sourcePointFeatureClass, integrate_name)
        distance = "{} Meters".format(res * 10)
        try:
            arcpy.Integrate_management(integrate_name, distance)
        except:
            message = arcpy.GetMessage(0)
            if 'maximum tolerance exceeded' in message.lower():
                print 'maximum tolerance exceeded: {} - {}'.format(integrate_name, distance)
                return None
            else:
                raise

        cluster_features = arcpy.CollectEvents_stats(integrate_name, cluster_name)
        add_attributes(integrate_name, cluster_name, output_name, )

        arcpy.Delete_management(integrate_name)
        arcpy.Delete_management(cluster_name)

        arcpy.AddField_management(output_name, 'zoom_level', "SHORT")
        expression = "int({})".format(levelInfo['level'])
        arcpy.CalculateField_management(output_name, 'zoom_level', expression, "PYTHON")

    return output_name

def add_attributes(integratedFeatureClass, clustersFeatureClass, outputFeatureClass, mergeNumericFields=[], mergeTextFields=False):
    '''
    TODO: Add join option for string fields; 
    TOOD: Add customizable fields for summary

    summarizes attributes for clusters
    First - The first input value is used.
    Last -The last input value is used.
    Join -Merge the values together using the delimiter value to separate the values (only valid if the output field type is text).
    Min -The smallest input value is used (only valid if the output field type is numeric).
    Max -The largest input value is used (only valid if the output field type is numeric).
    Mean -The mean is calculated using the input values (only valid if the output field type is numeric).
    Median -The median is calculated using the input values (only valid if the output field type is numeric).
    Sum -The sum is calculated using the input values (only valid if the output field type is numeric).
    StDev -The standard deviation is calculated using the input values (only valid if the output field type is numeric).
    Count -The number of values included in statistical calculations. This counts each value except null values.
    '''
    fieldmappings = arcpy.FieldMappings()
    fieldmappings.addTable(clustersFeatureClass)

    integrated_fields = arcpy.ListFields(integratedFeatureClass)
    for f in integrated_fields:
        if f.type in ['Integer', 'SmallInteger', 'Double', 'Single']:
            template = f.name + '_{}'
            for m in ['Min', 'Max', 'Mean', 'Sum', 'Count', 'StDev', 'Median']:
                fieldmap = arcpy.FieldMap()
                fieldmap.addInputField(integratedFeatureClass, f.name)
                field = fieldmap.outputField
                field.name = template.format(m.lower())
                field.aliasName = template.format(m.lower())
                fieldmap.mergeRule = m
                fieldmap.outputField = field
                fieldmappings.addFieldMap(fieldmap)

    arcpy.SpatialJoin_analysis(clustersFeatureClass, integratedFeatureClass, outputFeatureClass, "#", "#", fieldmappings)

def create_membership_where(featureClass, fieldName, values):
    where_membership = ','.join(["'{}'".format(v) for v in values])
    return arcpy.AddFieldDelimiters(featureClass, fieldName) + " IN ({})".format(where_membership)

def get_unique_field_values(featureClass, fieldName, where=None, getValueFunction=None):
    unique_values = []
    with arcpy.da.SearchCursor(featureClass, [fieldName], where) as cursor:
        for r in cursor:
            if getValueFunction:
                unique_values.append(getValueFunction(r[0]))
            else:
                unique_values.append(str(r[0]))
    final_list = list(set(unique_values))
    return final_list

all_results = []
delete_files = []
def cluster_complete(results):
    clustered_feature_classes, workerName, delete_file = results
    print 'worker complete({})'.format(workerName)

    all_results += clustered_feature_classes
    delete_files.append(delete_file)

def run_clustering(featureClass, clusterField, outputDirectory, includeLevels=None):
    global all_results
    global delete_files
    arcpy.gp.overwriteOutput = True
    pool = multiprocessing.Pool(multiprocessing.cpu_count() - 1)
    for v in get_unique_field_values(featureClass, clusterField):
        where = create_membership_where(featureClass, clusterField, [v])
        worker_arguments = (featureClass, outputDirectory, None, where, v, includeLevels)
        pool.apply_async(cluster_worker, worker_arguments, callback=cluster_complete)
    pool.close()
    pool.join()

    final_output = featureClass + '_clustered'
    arcpy.CopyFeatures_management(all_results.pop(), final_output)
    arcpy.Append_management(all_results, final_output, "TEST")

    for d in delete_files:
        arcpy.Delete_management(d)

#=============================================================================================================
# TESTING
#=============================================================================================================
class TestClustering(unittest.TestCase):
    def test_run_clustering(self):
        feature_class = r'D:\data\EPA\FRS_INTERESTS.gdb\FACILITY_INTERESTS_WM'
        clusterField = 'STATE_CODE'
        outputDirectory = r'D:\data\EPA'
        levels = range(4, 16)
        run_clustering( feature_class, clusterField, outputDirectory, includeLevels=levels )

    def _add_attributes(self):
        arcpy.gp.overwriteOutput = True
        integrated_feature_class = r'C:\gis_working_directory\tempDE.gdb\integrate_10'
        clusters = r'C:\gis_working_directory\tempDE.gdb\clusters_10'
        outputDirectory = r'C:\gis_working_directory\tempDE.gdb\clusters_10_joined'
        add_attributes(integrated_feature_class, clusters, outputDirectory)

if __name__ == '__main__':
    unittest.main()