Counting points within polygon with ArcPy?

Question

I am trying to script a process where I need to make a selection on a point feature class, count the number of points that fall within each polygon on another feature class, and then update a COUNT field with the new value. I can preform the selection and creating a new layer from the selection via ArcPy but I cannot figure out how to count the points within the polygon via ArcPy. I tried

import arcpy
try:  
   #select all trip hazards not marked 'Closed'  
   arcpy.SelectLayerByAttribute_management ("Sidewalk_Trip_Hazards",   "NEW_SELECTION", "\"Status\" <> \'Closed\'") 

   #make layer of selected hazards
   arcpy.MakeFeatureLayer_management("Sidewalk_Trip_Hazards","notClosedHazards")

    #spatially join points to hexbins
    arcpy.SpatialJoin_analysis("notClosedHazards", "LaytonHexBins", "in_memory/SpatialJoin")

    #summarize points based off of HexID  
    arcpy.Statistics_analysis("in_memory/SpatialJoin", "in_memory/SummarizedTable", [["HexID", "Count"]])

except:
    print(arcpy.GetMessages())  

Is there a way via ArcPy to do a points in polygon count?

Answer 1

SummarizeWithin looks to only be a ArcGIS Online and ArcGIS Pro tool. However, the documentation states --

Summarize Within performs the functions of the Spatial Join and Summary Statistics tools.

Therefore, you would first compute a spatial join between your points and polygons, attributing each point with the unique id of their containing polygon, and then you can summarize the count of points for each unique id. Something like...

arcpy.SpatialJoin_analysis(target_features (Points), join_features (Polys), out_feature_class)
arcpy.Statistics_analysis (out_feature_class, out_table, [["UID", "Count"]])

UPDATED with slightly modified methods...

import arcpy

points =
polygons =
polygonID = ## unique polygon field name
countField = ## field to be updated
expression = "recalc(!FREQUENCY!)"
codeblock = """def recalc(freq):
    if freq > -1:
        return freq
    else:
        return 0"""

arcpy.SpatialJoin_analysis(points, polygons, "in_memory/PointsInPolys")
## case field returns count per unique UID
arcpy.Statistics_analysis ("in_memory/PointsInPolys", "in_memory/SS_PointsInPolys", [[polygonID, "Count"]], polygonID)
arcpy.JoinField_management(polygons, polygonID, "in_memory/SS_PointsInPolys", polygonID, "FREQUENCY")
arcpy.CalculateField_management(polygons, countField, expression, "PYTHON", codeblock)
arcpy.DeleteField_management(polygons, "FREQUENCY")

Answer 2

You don't mention ArcGIS Pro in your post/tags, so I would assume you don't have Pro. That would explain why you don't have the SummarizeWithin tool. If you have the Business Analyst extension, you could use Summarize Points; see here for checking out an extension within ArcPy.

If, however, you have none of those options, you can use some relatively simple Python. Follow these steps...each linked page has the syntax for the tools:

1. Spatial Join, passing in your points as the target features and the polygons as the join features. Note that if your polygons overlap each other, then you should select the option for a one-to-many join. Use either Intersect or Within as the match_option, whichever fits your needs.
2. Run a Dissolve, passing the spatially joined layer as the input, passing the OID or some other unique ID field that originated in the polygons layer as the dissolve field, and finally passing joined layer's unique ID field (i.e., the one that refers to unique points) as a statistics field with the statistic COUNT.
3. Finally, run a Join, with your polygons as the input and the result of the dissolve as the join_table.