Count number of features inside polygon

[CCWeedcontrol]

I need some help on my code to count the number of features inside a polygon. I have the following but after testing it, it doesn't seem to be correct. I am trying to figure out how many parcel lots are within the subdivision boundary. My test show after running the code that the "Stats" table has 9 Frequencies but in actually there is 10 parcels within the subdivision. I have tried "HAVE_THEIR_CENTER_IN" but all the results are 1.

import arcpy, os

mxd = arcpy.mapping.MapDocument("current")
df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
lyr = arcpy.mapping.ListLayers(mxd, "Taxlots")[0]
subs = arcpy.mapping.ListLayers(mxd, "SubsOutsideCities")[0]
#arcpy.MakeFeatureLayer_management(par, "parLyr")

arcpy.SpatialJoin_analysis(subs, lyr,"in_memory/points_SpatialJoin", "JOIN_ONE_TO_MANY", "KEEP_ALL", "", "CONTAINS")
arcpy.Statistics_analysis("points_SpatialJoin", "in_memory/stats", "Join_Count COUNT","OBJECTID")
‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

[DanPatterson_Retired]

Your shapes are all singlepart?  If not, run the Multipart to singlepart tool

[CCWeedcontrol]

Yes they are all singlepart.

[DanPatterson_Retired]

many to one

one to one

one to many

switch the order of the join

which options did you try?

[CCWeedcontrol]

There is only two options for the spatialJoin_analysis, i am not sure what you mean by the other options.

JOIN_ONE_TO_ONE

JOIN_ONE_TO_MANY

[JoeBorgione]

Looks like you are using ArcMap to perform you analysis.  You might try to do it by hand, and then go to the geoprocessing results button; that'll tell what ArcMap just did for you and you can use that in your python coding.  I do that all the time with tools in Pro.

[CCWeedcontrol]

I have tried that and get the same results.

[RandyBurton · Accepted solution]

Perhaps something like:

countFC = 'polygons'  # feature containing outer polygons
countFld = ['ParcelCount', 'SHAPE@'] # ParcelCount is for count results, SHAPE@ is geometry of polygon

parcels = 'parcels' # parcels to be counted

with arcpy.da.UpdateCursor(countFC, countFld) as cursor: # cycle through outer polygon layer
    for row in cursor:
        arcpy.management.SelectLayerByLocation(parcels, "HAVE_THEIR_CENTER_IN", row[1], "", "NEW_SELECTION") # row[1] is SHAPE@
        row[0] = int(arcpy.GetCount_management(parcels)[0]) # row[0] to be updated with count
        cursor.updateRow(row) # update row‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

[CCWeedcontrol]

Randy I was able to get to work, thanks for the help.

import arcpy

mxd = arcpy.mapping.MapDocument("current")
df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
lyr = arcpy.mapping.ListLayers(mxd, "Taxlots")[0]
subs = arcpy.mapping.ListLayers(mxd, "SubsOutsideCities")[0]
#arcpy.MakeFeatureLayer_management(par, "parLyr")

with arcpy.da.UpdateCursor(subs, ["VxCount1", "Shape@"]) as cursor: # cycle through outer polygon layer
    for row in cursor:
        arcpy.management.SelectLayerByLocation(lyr, "HAVE_THEIR_CENTER_IN", row[1]) # row[1] is SHAPE@
        row[0] = int(arcpy.GetCount_management(lyr)[0]) # row[0] to be updated with count
        cursor.updateRow(row) # update row‍‍‍‍‍‍‍‍‍‍‍‍‍

[DanPatterson_Retired]

save the result to disk, and open the table showing the OBJECTID for the join count.

click on every feature in turn and see if one of the features is actually 2 parts, or one feature is partially on top of the other