3

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: enter image description here

1 Answer 1

3

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.

4
  • @Trojan1887 I'm curious to know if this worked well for you.
    – dru87
    Aug 11, 2017 at 9:22
  • I am still having issues with my code, especially when adding the new steps for zonal statistics and so forth
    – Trojan1887
    Aug 12, 2017 at 13:56
  • I am testing the make feature layer output by saving it to the gdb but I get an empty feature class
    – Trojan1887
    Aug 12, 2017 at 18:33
  • maybe add a screenshot of the attribute table of your feature class containing the basin_id. if all is right with your feature class it shouldn't produce empty results. at least it didn't when i tested my code. in case you're stuck with the ZonalStatistics() you should possibly ask a new question.
    – dru87
    Aug 14, 2017 at 5:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.