1

I need to iterate through each of the features of a feature class successively to use each one as the mask to extract from a raster with the ExtractByMask tool. Basically I'm looking for the "Iterate Feature Selection" function from ModelBuilder for Python scripting.

If I understand correctly, using a search cursor will return only the row of values from the attribute table, not the feature shape itself.

Would iteratively creating a new shape with "Select_Analysis" be a feasible way to go about this instead?

I suppose creating a new shapefile for each feature would take quite long.

Is there another, more efficient way?

2

The SearchCursor will return geometries with the SHAPE@ token which can be used as extracting features etc.:

SHAPE@ —A geometry object for the feature.

import arcpy

feature_class = r'C:\test.gdb\polygon'

with arcpy.da.SearchCursor(feature_class,'SHAPE@') as cursor:
    for row in cursor:
        #do something with row[0]

You can of course also return field values and use for example to name the outputs. In that case use a list of tokens and fields:

with arcpy.da.SearchCursor(feature_class,['SHAPE@','field1']) as cursor:
    for row in cursor:
        #do something

row[0] will be the geometry and row[1] the value of field1

  • ok, I think this is very close to what I'm looking for, but when I call the 'SHAPE@' token to define the cursor and then use row[1] as input for the ExtractByPolygon tool I get: "ERROR 000859: The required parameter Polygon is empty, or is not the type of Multiple Value. Failed to execute (ExtractByPolygon)." Using arcpy.Describe(row[1]) and print(desc.dataType) gives me: "RuntimeError: Object: Describe input value is not valid type" I suspect that I have to create an actual feature class with row[1] and use that as input for processing. Any ideas? Many thanks! – Yan Hildebrand Oct 22 '17 at 16:48
  • The indexing starts with zero so if you only have the shape token in your cursor you need to use row[0] – BERA Oct 22 '17 at 16:51
  • ah sorry, I should have mentioned that, no I have the ID as well, so shape is at 1. – Yan Hildebrand Oct 22 '17 at 18:22
  • In your question you are using Extract by Mask, now you say Extract by Polygon. They require different inputs. – BERA Oct 22 '17 at 18:25
  • yes, I realized I should use ExtractByPolygon if I want to use a Polygon Feature Class as a mask. But this doesn't seem to work either. – Yan Hildebrand Oct 23 '17 at 19:18
0

Here's how I handled selecting multiple polygons scattered throughout an area that clips from a single raster image(DEM). My area covered four counties, so I had to change the DEM (hardcoded) each time I ran it. May be a little klunky, but it works. All of my polygons begin with the same naming convention. I use a lot of print statements in my code to verify it's running and that I'm getting the expected outputs. My sanity check....

    import arcpy, os
    from arcpy import env
    from arcpy.sa import *
    arcpy.CheckOutExtension("Spatial")
    env.workspace = "//w63-atlas/Stormwater/Waterways/Structure Flood Severity"
    bslash = "\\"
    Raster = bslash+"\\w63-atlas\Stormwater\AMEC\Platte\TOPO\DEM\platte_dem"
    print " Raster lowkayshunn ",Raster
    tile = raw_input("Enter FEMA TILE : C4, C5, C6, D4, D5, D6, E5, E6, F5, F6, G5, or G6: ")
    outFC = "//w63-atlas/Stormwater/Waterways/Structure Flood Severity"
    prefix = "Export_"
    suffix = "_Culv*_1.shp"
    depiwant = prefix+tile+suffix
    slash = "//"
    aewout = outFC+"/"
    print "aewout = ", aewout
    print "depiwant = ", depiwant
    outFC = outFC+"/"+depiwant
    print "Raster = ",Raster
    print "outFC = ",outFC
    arcpy.env.scratchWorkspace = "//w63-
    atlas/Stormwater/Waterways/Rasterdatasets.gdb"
    gdb = "//w63-atlas/Stormwater/Waterways/Rasterdatasets.gdb"
    print "gdb = ", gdb
    fcList  = arcpy.ListFeatureClasses("Export_"+tile+"_Culv*_1.shp")
    print "Shape files to clip from: " , fcList
    arcpy.AddMessage(Raster)
    for i in fcList:
        print "aewout+i ",aewout+i
        print "Shape file: ",i
        Culvnum = i[14:-4] #strip out the name
        rasnew = gdb+"/"+tile+"rz"+Culvnum[:-2]
        print "Culvnum ", Culvnum
        outExtractByMask = ExtractByMask(Raster,i) #create temp raster file
        arcpy.AddMessage(outExtractByMask)
        print "Raster image = ",Raster
        print "current Culv = ", i
        print "outExtractByMask = ", outExtractByMask
        outExtractByMask.save(rasnew) # save temp raster to final name
        print "Raster CLip NAMe: ", rasnew
        print "All done creating raster files.  NEXT !!!"
        exit
0

Here is a solution using ExtractByMaskto iterate over a polygon feature class and batch clip a raster by each feature. In the polygon feature class grid_id is the unique name of each feature and is simply an integer that is used for naming the raster outputs. This will save each clipped grid in an IMAGINE format though this can be changed by editing the suffix. You will need to update the paths in 3 places. Don't for get the trailing \\ in the first argument of the outmask.save function so the last level in the path is recognized as a folder.

grid = "path\\to\\grid"    
clip_grid = "path\\to\\polygonFeatureClass"

with arcpy.da.SearchCursor(clip_grid,['SHAPE@','grid_id']) as rows:
    for row in rows:
        print("Starting Grid Cell " + str(row[1]))
        outmask = arcpy.sa.ExtractByMask(grid, row[0])
        outmask.save("path\\to\\grid\\folder\\" + str(row[1]) + ".img")

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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