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")

• cbunn: I do have access to Pro so I will try it in there. And is it possible to write the output to an temporary file since I don't want to create a new feature class, but rather replace the old COUNT values with the new ones on an existing feature class?
• Lovette: That sounds like a case where you may want to preform the Spatial Join, then read the table as an array in python so that you can count the instances of each unique value and then use a cursor on your existing feature class to update the count values. You could write the spatial join output to a temporary location and delete it after your script is done if you don't need it.
• cbunn: Could you include in your code how to write it out to a temporary file?
• cbunn: Ok I got the temporary files to work but the summarized table is only giving me one row of data and the count is wrong. How do I get the table to give me a summary of the count of points within each polygon?
• Lovette: You may want to add a line of code after the summary statistics and/or join fields that sets all null values to 0. does the "Frequency" field join to your original table?

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.