Need to Extract by Mask on Features in a Feature Class

Question

I have a polygon feature class of a country that has water basin polygons as individual features. I am trying iterate each feature and Extract by Mask using a set of rasters - for instance on Basin_ID I need to extract by mask from crop_a, crop_b, crop_c. I was able to create this sequence in Model Builder using the "Iterate Feature Selection" tool but I want to script this so that I have more control and because I need to collect a sum from zonal statistics down the line.

I have "tried" opening a cursor on the country, creating a for loop (row in cursor:) and then running a select query on the feature to make a new feature layer to use as a mask for the extract by mask tool. I assume that is what the model builder is doing but I can't seem to recreate this myself.

The beginning of my model essentially looks like this:

Answer 1

Here's a script that should work with your data, if I got you correctly. I tested it with my own data and it creates the desired output rasters. Note that depending on the number of features in your feature class and the number of input rasters you pass to the tool, it could be running for a while. Make sure to change the lines that name your output rasters if needed.

import arcpy, os
from arcpy.sa import *


# Set input parameters
workspace = r"C:\TESTWORKSPACE.gdb"
tomatoProductionRaster = r"C:\TESTWORKSPACE.gdb\tomatos"
appleProductionRaster = r"C:\TESTWORKSPACE.gdb\apples"
maskFC = r"C:\TESTWORKSPACE.gdb\canada"
maskField = u"Basin_ID"


# Set environments
arcpy.env.workspace = workspace
arcpy.env.scratchWorkspace = workspace
arcpy.env.overwriteOutput = True


# Get Spatial Analyst Extension
if arcpy.CheckExtension("Spatial") == "Available":
    arcpy.CheckOutExtension("Spatial")

    # ----------------------
    # START PROCESSING DATA
    # ----------------------

    # Create lists
    rasterList = [tomatoProductionRaster, appleProductionRaster]
    # you could replace the line above with the following line in order to automatically chose the rasters in your workspace:
    ## rasterList = arcpy.ListRasters()

    # Start Iteration
    # Create cursor on your maskField (i.e. Basin_ID) and the shape of the feature:
    with arcpy.da.SearchCursor(maskFC, [maskField, "SHAPE@"]) as cursor:
        for row in cursor:
            # Create Feature Layer of single Basin_ID
            arcpy.MakeFeatureLayer_management(row[1], "tempFeature")
            # Create name part with Basin ID
            newName1 = u"basin{}".format(row[0])
            # Iterate through input rasters
            for raster in rasterList:
                # Create name part with raster being processed
                newName2 = u"{}".format(os.path.split(raster)[1])
                outName = u"{}_{}".format(newName1, newName2)
                # Run extract by mask for all rasters in your list
                extractOut = ExtractByMask(raster, "tempFeature")
                extractOut.save("{}\\{}".format(workspace, outName))
    # Return extension:
    arcpy.CheckInExtension("Spatial")

else:
    print(u"Necessary spatial analyst extention unavailable. Script stopped")

There would of course be nicer ways to check for the extension (using try and except statments) but i thought this would be sufficient for the job.

In order to use ZonalStatistics() after the extraction, you could for instance use the command arcpy.ListRasters("basin*"), which would list all your output rasters and then continue with the ZonalStatistics() in a similar manner as shown for the ExtractByMask() tool.