topic Re: Slow Count Points in Polygons Script in Python Questions

Slow Count Points in Polygons Script

HollyCopeland — Thu, 06 Oct 2011 21:44:25 GMT

I wrote this script to iterate through a polygon fc and count the number of points in each polygon and put that information in a new field. The script works, but it is REALLY slow. If I use Hawth's Tools to do the same thing, it takes 20-30 minutes. With this script, it was set to take 4 days. Any ideas what I've done wrong--I have a feeling it's related to the search cursor. Should I not be using the select by location/select by attribute ESRI functions? Thanks, Holly

#-------------------------------------------------------------------------------
# Name:       Count Points in Polygon
# Purpose:
#
# Author:      hcopeland
#
# Created:     22/08/2011
# Copyright:   (c) hcopeland 2011
# Licence:     <your licence>
#-------------------------------------------------------------------------------
#!/usr/bin/env python
import arcpy
from arcpy import env

#getcontext().prec = 8

ws = r'D:\Work\Projects\SageGrouseEasementScenarios2011\Data\FinalModelData.gdb'
arcpy.env.workspace = ws

# set overwrite outputs to true
arcpy.env.overwriteOutput = 1

POINT = 'D:\Work\Projects\SageGrouseEasementScenarios2011\Data\Final_Development_Buildouts\Cumulative_buildouts.gdb\Cumulative_Unconstrained'
POLYGON = 'LekParcelFC'

print 'Making Multipart into Single Part....'
#Convert multipart to singlepart
#arcpy.MultipartToSinglepart_management(POLYGON,"LekParcelFCSingle")

# Converts unselectable feature class into a selectable feature layer
arcpy.MakeFeatureLayer_management(POINT,"point")
arcpy.MakeFeatureLayer_management(POLYGON,"polygon")

rows = arcpy.SearchCursor("polygon")

print 'Calculating Points in Polygons...'

# Loop through each row and count the points in each polygon record
try:
    for row in rows:
        # Select each record inside of the polygon feature class
        #expression = "\"OBJECTID\" =" + str(row.OBJECTID) + ' AND "PARCELNB" <> ' + "''"
        expression = "\"OBJECTID\" = " + str(row.OBJECTID) + ""
        print expression
        SelPoly = arcpy.SelectLayerByAttribute_management("polygon", "NEW_SELECTION",expression)

        # Select all the point that are inside of the polygon record
        SelPts = arcpy.SelectLayerByLocation_management("point", "WITHIN", SelPoly, 0, "NEW_SELECTION")

        # Count the points that are in each polygon
        GetCount = arcpy.GetCount_management(SelPts)

        # Calculate the ASSOC_PTS field with the counted points
        arcpy.CalculateField_management("polygon", "CEFeaturesPreventedHigh", GetCount, "VB", "")

except:
    arcpy.GetMessages(2)
    arcpy.AddError("Script Bombed")

print 'Count Points in Polygon Script Finished!'

Re: Slow Count Points in Polygons Script

ChrisSnyder — Thu, 06 Oct 2011 22:20:47 GMT

Yes, that would indeed be a slow boat to China... A faster method: Use the SpatialJoin or Identity tool to assign the polygons OBJECTID onto the points. Then use the Frequency tool to count the occurances (points per polygon OBJECTID). Then use a simple tabular join and field calc to put the count (aka "FREQUENCY") values into a new field in the polygon FC.

Re: Slow Count Points in Polygons Script

HollyCopeland — Sat, 11 Dec 2021 16:57:11 GMT

Spatial join might work, but does it work for overlapping polygons? It is critical that the points are accurately counted in my case. I adjusted a few parts of my script and it has sped it up considerably. It's still not fast enough to process 300,000 records in under a few hours, but I thought I'd post this improved version in case it was helpful to someone else out there. I still don't understand why this simple procedure has to take so long! :confused:

#-------------------------------------------------------------------------------
# Name:       Count Points in Polygon
# Purpose:
#
# Author:      hcopeland
#
# Created:     22/08/2011
# Copyright:   (c) hcopeland 2011
# Licence:     <your licence>
#-------------------------------------------------------------------------------
#!/usr/bin/env python
import arcpy
from arcpy import env

#getcontext().prec = 8

ws = r'D:\Work\Projects\SageGrouseEasementScenarios2011\Data\FinalModelData.gdb'
arcpy.env.workspace = ws

# set overwrite outputs to true
arcpy.env.overwriteOutput = 1

POINT = 'D:\Work\Projects\SageGrouseEasementScenarios2011\Data\Final_Development_Buildouts\Cumulative_buildouts.gdb\Cumulative_Baseline'
POLYGON = 'LekParcelFC'

# Converts unselectable feature class into a selectable feature layer
arcpy.MakeFeatureLayer_management(POINT,"point")

rows = arcpy.SearchCursor(POLYGON)

print 'Calculating Points in Polygons...'

### Loop through each row and count the points in each polygon record
try:
    for row in rows:
        expression2 = "\"OBJECTID\" = " + str(row.OBJECTID) + ""
        #expression = "\"OBJECTID\" = " + str(row.OBJECTID) + ""' AND "PARCELNB" <> ' + "'""'"
        print expression2
        # Select each record inside of the polygon feature class
        arcpy.MakeFeatureLayer_management(POLYGON,"polygon2",expression2)

        # Select all the points that are inside of the polygon record
        arcpy.SelectLayerByLocation_management("point", "intersect", "polygon2", 0, "NEW_SELECTION")

        # Count the points that are in each polygon
        GetCount = arcpy.GetCount_management("point")

        # Calculate the ASSOC_PTS field with the counted points
        arcpy.CalculateField_management("polygon2", "CEFeaturesPreventedBaseline", GetCount, "VB", "")

except:
    arcpy.GetMessages(2)
    arcpy.AddError("Script Bombed")

print 'Count Points in Polygon Script Finished!'

Re: Slow Count Points in Polygons Script

ChrisSnyder — Tue, 11 Oct 2011 15:57:46 GMT

but does it work for overlapping polygons?

There's only one way to find out...

Using the Spatial Join tool, you will need to specify a Join operation of "ONE_TO_MANY" to deal with the overlapping polys. In addition to the afore mentioned Identity tool, the Intersect tool will also work. Both the Identity and Intersect correctly handle the one-to-many point/poly relationship by default.

Don't reinvent the wheel... There are (at least three) out of the box tools that are designed to do exactly what you want!