Labeling each polygon based on the number of points inside it using Python

Question

I have a parcel layer and a trees layer, I want to label each parcel with the number of trees it contains and the sum of DBH. I wrote a script that calculates this and it works stand-alone, but I can't figure out how to get it to work with a Label Expression. Here is the code I've tried (ArcMap 10.2):

import arcpy
import numpy
def FindLabel():
  mxd = arcpy.mapping.MapDocument("CURRENT")
  parcels = arcpy.mapping.ListLayers(mxd, "Parcel*")[0]
  trees = arcpy.mapping.ListLayers(mxd, "Surveyed Trees*")[0]
  arcpy.SelectLayerByLocation_management(parcels, "CONTAINS", trees, "#", "NEW_SELECTION")
  arcpy.SelectLayerByLocation_management(trees, "WITHIN", parcels, "#", "NEW_SELECTION") 
  with arcpy.da.SearchCursor(parcels, ["OBJECTID"]) as cur:
    for row in cur:
      arcpy.SelectLayerByLocation_management(trees, "WITHIN", parcels, "#", "NEW_SELECTION")
      treeCount = int(arcpy.GetCount_management(trees).getOutput(0))
      field = arcpy.da.TableToNumPyArray(trees, ("DBH"), skip_nulls = True)
      sum = field["DBH"].sum()
      message = "Trees = " + str(treeCount) + "\n" + "DBH = " + str(sum)
      return message

This sends ArcMap into a long loop and then seems to return the last calculated value for each parcel, so everything is labeled "1 Tree 7 DBH".

Maybe I shouldn't be using a SearchCursor?

Solution

Here's the code I ended up using. It's still processor intensive, but this fits our workflow fine.

import arcpy
def FindLabel ( [Loc_Name] ):
  mxd = arcpy.mapping.MapDocument("CURRENT")
  trees = arcpy.mapping.ListLayers(mxd, "ProgramSurveyGround*")[0]
  arcpy.MakeTableView_management(trees, "treesTable")
  parcels = arcpy.mapping.ListLayers(mxd, "ParcelProcessing*")[0]
  Query = ' "Loc_Name"' + " = '" + [Loc_Name] + "'"
  treeSelect = arcpy.SelectLayerByAttribute_management("treesTable", "NEW_SELECTION", Query)
  treeCount = int(arcpy.GetCount_management("treesTable").getOutput(0))
  dbhArray = arcpy.da.TableToNumPyArray("treesTable", ('DBH'), skip_nulls = True)
  dbhCount = dbhArray['DBH'].sum()
  myText = "Trees : " + str(treeCount) + "\n" + "DBH : " + str(dbhCount)
  return myText

Answer 1

Three calls to arcpy.SelectLayerByLocation_management and one call to arcpy.da.TableToNumPyArray!

Suppose you have solved the bug in your code but calling 4 process costly geoprocessing functions in a label expression is not a rational choice and doesn't follow the best practice.

1. You can either pre-process the data and use spatial join (Merge Rule = Count)
2. Create relationship classes between the two layer and count the number of related features and show as label. This is both dynamic and process friendly. However you should take care of point feature deletion and addition and auto relate the point feature to the polygon upon different update events.