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.