Select nearest polygons from point that sum to a value in the point

Question

I have a polygon tessellation representing 5 acres, and is populated with an estimated value of children between the ages of 0-10. I have a point feature class that represents child care locations and it contains a capacity number. I want to somehow represent how many of the nearest polygons would be needed to fill the capacity of the point value. So if I have a point with a capacity of 50, I would want to know how many of the closest polygons that point would cover.

VinceE · Answer

Here's an attempt. Hexagons have "populations" between 1 and 5 (light to dark). Black points are the child care centers and their capacities are shown in the adjacent black labels. Regions made up of hexagons that have aggregated enough "population" to reach each capacity for each center are shown in the colors, with the "actual" population gathered labeled with the callouts.

The caveats:

* Regions can overlap if centers are too close (see red/blue). Lots of fixes for this, depending on what you're after.
* Childcare centers can be "over" capacity. See green region, where the max capacity is 80, and 81 have been "gathered" from adjacent hexagons. Also lots of fixes, depending on needs.

Code below, in simple terms: for each center, add the closest hexagon to a region, followed by the next closest, etc., until enough hexagons have been added so the summed population reaches the capacity of the child care center. This runs extremely quickly with ~1,300 hexagons, but is definitely not the most optimized.

You'll have to edit this code if you want to use it with your own feature classes. Happy to provide more details if you need them.

import arcpy
import math
import os
import sys

arcpy.env.workspace = r"YOUR_GEODATABASE_PATH"
arcpy.env.overwriteOutput = True

# Names of hexagon bins and the associated points.
hex_bin = "HexBin"
hex_pnt = "HexAggPoint"

# Create an info dictionary about the hexagons
# {oid: {HEX_GEOM: geometry, POP: integer}}
with arcpy.da.SearchCursor(hex_bin, ["OID@", "SHAPE@", "POP"]) as scurs:
    hex_pop_dict = {oid: {"HEX_GEOM": geom, "POP": pop} for oid, geom, pop in scurs}

# Empty dictionary to store info about our output regions.
clusters = {}

# Loop over points, logging assembled region information in 'clusters' dict.
with arcpy.da.SearchCursor(hex_pnt, ["OID@", "SHAPE@XY", "CAPACITY"]) as scurs:
    for oid, centroid, capacity in scurs:

# Get list of (hex_id, distance_to_point) tuples.
        distances = []
        for hex_id, hex_record in hex_pop_dict.items():
            hex_centroid_coor = hex_record["HEX_GEOM"].centroid
            dist = math.dist(centroid, (hex_centroid_coor.X, hex_centroid_coor.Y))
            distances.append((hex_id, dist), )

# Sort tuples ascending distance from the point in question.
        distances.sort(key=lambda x: x[1])

# Loop distance tuples, using the hex ID to accumulate population per hex.
        # Log the hex IDs that are used. Stop adding more hexes onces the capacity
        #   of the point is reached.
        gathered_hexes = []
        accumulated_pop = 0
        for hex_id, _ in distances:
            gathered_hexes.append(hex_id)
            accumulated_pop += hex_pop_dict[hex_id]["POP"]

if accumulated_pop >= capacity:
                break

# For this Point ID, associated a list of hexes that will create region.
        clusters[oid] = gathered_hexes
        print(f"HEX ID: {oid}
CAPACITY: {capacity}"
              f"
ACTUAL ACCUM: {accumulated_pop}"
              f"
HEX IDS: {gathered_hexes}
")

# CREATE OUTPUT FC, ADD FIELDS
region_fc = arcpy.management.CreateFeatureclass(out_path=arcpy.env.workspace,
    out_name="HEX_REGION", geometry_type="POLYGON", spatial_reference=2232)
for fld_name in ["POINT_ID", "HEX_ID", "POP"]:
    arcpy.management.AddField(region_fc, fld_name, "SHORT")

# WRITE OUTPUT RECORDS
with arcpy.da.InsertCursor(region_fc, ["SHAPE@", "POINT_ID", "HEX_ID", "POP"]) as icurs:
    # For this point and list of hexes that should create the region...
    for point_id, hex_id_list in clusters.items():

# ...and for each hex in the list...
        for hex_id in hex_id_list:

# Using hex's ID, insert Geometry and Population,
            #   plus the ID of both the parent point and the hex.
            icurs.insertRow([hex_pop_dict[hex_id]["HEX_GEOM"],
                             point_id,
                             hex_id,
                             hex_pop_dict[hex_id]["POP"]])

# Dissolve the hex polygons based on the associated point. Sum the pop of the region.
dissolve_fc_path = os.path.join(arcpy.env.workspace, "HEX_REGION_DISSOLVE")
arcpy.management.Dissolve(in_features="HEX_REGION", out_feature_class=dissolve_fc_path,
            dissolve_field="POINT_ID", statistics_fields=[["POP", "SUM"], ])

LRoberts_RDM · Answer

Hi @SeaRM,
I feel like that would be basically the same as symbolizing my points by proportional symbol, using the capacity as my value. I need it to somehow be affected by the polygons below it. So if I have a point with a capacity of 50, and it's surrounded by polygons with a value of 1, it would cover more. But if I had a point with a capacity of 50 and it sits on top of a polygon of 50, it would just cover that polygon. Maybe this is more a spatial analyst function?

SeaRM · Answer

I think in this case the Weighted Thiessen Polygons built for your child care locations (weighted by the capacity field) can be useful.