Clip Point Features to Polygons with ModelBuilder/ArcPy

Question

I have point features for about 2000 individual trees spread throughout 50 different tree districts. The tree districts are polygons in a polygon feature class and the trees are point features in a point feature class.

I would like to create a mapbook showcasing each individual tree district and the corresponding trees in each district. I could use Data Driven Pages, but I need closeups of each tree district (polygon) removed from its surrounding districts.

I could accomplish this by running a definition query to display a single tree district, clip the trees to that district, and then create a mapbook page. But I'd have to do this 50 times.

How can I use ModelBuilder or ArcPy to perform the aforementioned process?

I've tried using How to Batch Clip in ArcGIS Desktop 10 using Python/ArcPy? as a guide but it doesn't exactly fit my needs.

Answer 1

This is a great example of how relatively simple arcpy.mapping scripts can offer more functionality than Data Driven Pages.

First, if you haven't already, run a Spatial Join with your trees as your target layer and your districts as your point layer using the INTERSECT match option. It's better practice to perform a spatial join and apply a definition query rather than running the Clip tool in each iteration. This is because since geometry operations are usually more performance-intensive and and prone to problems than table operations like definition queries.

Second, edit the following Python snippet to match your data. This code iterates through districts, applies definition queries to both the trees points and districts polygons, zooms the data frame to the current district polygon, and exports a JPG map to the output folder.

import arcpy, os

#Input data
Trees = "Trees" #Point feature class of trees with district attribute
Districts = "Districts" #Polygon feature class of districts
arcpy.env.workspace = r"C:\Temp" #Workspace, replace as needed.
DistrictField = "District" #Name of district ID field in Districts and Trees classes

#MapDocument, DataFrame, and Layer objects
CurrentMXD = arcpy.mapping.MapDocument("CURRENT") #Current map doc
CurrentDF = CurrentMXD.activeDataFrame #Current data frame
Trees_lyr = arcpy.mapping.ListLayers(CurrentMXD, Trees, CurrentDF)[0] #Trees layer object
Districts_lyr = arcpy.mapping.ListLayers(CurrentMXD, Districts, CurrentDF)[0] #Districts layer object

#Iterate through districts
ListOfDistricts = [Row[0] for Row in arcpy.da.SearchCursor(Districts, [DistrictField])]
for District in ListOfDistricts:

    #Apply definition queries. Assumes district name is a string, remove quotes if not string.
    Trees_lyr.definitionQuery = "%s = '%s'" % (arcpy.AddFieldDelimiters(Trees,DistrictField), District)
    Districts_lyr.definitionQuery = "%s = '%s'" % (arcpy.AddFieldDelimiters(Districts,DistrictField), District)

    #Adjust zoom
    CurrentDF.extent = Districts_lyr.getExtent(False) #Extent of queried features

    #Export output jpg
    Output_jpg = os.path.join(arcpy.env.workspace,"District_%s.jpg" % District)
    arcpy.mapping.ExportToJPEG(CurrentMXD, Output_jpg, resolution=300, jpeg_quality=95)

#Reset definition queries
Trees_lyr.definitionQuery = ""
Districts_lyr.definitionQuery = ""

Here is my example data, with district A on the right and District B on the left:

Here is the District_A.jpg output file:

Here is the District_B.jpg output file: