How to clip raster image by the extent of each row in a shapefile with ArcGIS 9.3

Question

I'm using Python to clip a raster image by each row in a shapefile with ArcGIS 9.3. I am hoping to end up with one raster image for each polygon row, clipped to the extent of the polygon row. However, my current script clips the raster to the entire shapefile, instead of setting the extent by each polygon row. My script is as follows:

# Import system modules
import sys
import string
import os
import arcgisscripting

# Create the Geoprocessor object
gp = arcgisscripting.create(9.3)

# Load required toolboxes...
gp.AddToolbox("C:/Program Files (x86)/ArcGIS/ArcToolbox/Toolboxes/Data Management Tools.tbx")

# Designate a workspace
gp.workspace = "L:\Clip"

# Assign variables...
clipto = "test.shp"

# grids to clip to
clipit = "vfcm.tif" # features to be clipped

# Cursor time
rows = gp.SearchCursor(clipto)
row = rows.Next()

# clip management
while row:
    n = str(row.VAHU6) # assign a variable for the processing message based on a field
print "clipping to: " +n # tells what grid was clipped

# Clip_management <in_raster> <rectangle> <out_raster> {in_template_dataset}   {nodata_value} {NONE | ClippingGeometry}
outputname = str(row.VAHU6) + "_c"
desc = gp.Describe(clipto)
extent = str(desc.Extent)
gp.Clip_management(clipit, extent, outputname, clipto, "", ClippingGeometry) #update str if necessary
print str(row.VAHU6)
row=rows.next()

# reset the array...
del rows
del gp

Any ideas how to fix it?

Thanks!

Answer 1

You need the shape field of the row to clarify what geometry you're clipping by. So somewhere near your while loop (either inside or just before), specify the geometry like so:

clipgeo = row.Shape

And then use the object (in this case, clipgeo) as your feature to clip by instead of the entire shapefile.

A similar problem can be seen here: How to batch clip to single-part shapefile using ArcGIS 9.3 and Python 2.5?

Answer 2

For those who want to do this in 10.2, this modification of above worked for me:

# Import system modules
import sys
import string
import os
import arcgisscripting

# Create the Geoprocessor object
gp = arcgisscripting.create(10.2)

# Load required toolboxes...
gp.AddToolbox("C:/Program Files (x86)/arcgis/Desktop10.2/ArcToolbox/Toolboxes/Data Management Tools.tbx")

# Designate a workspace
gp.workspace = "D:/Clip"

# Assign variables...
clipto = "test.shp"

# grids to clip to
clipit = "250K.ecw" # features to be clipped

# Cursor time
rows = gp.SearchCursor(clipto)



# clip management
for row in rows:
    n = str(row.Name[:9]) # assign a variable for the processing message based on a field
    print "clipping to: " +n # tells what grid was clipped

    # Clip_management <in_raster> <rectangle> <out_raster> {in_template_dataset}   {nodata_value} {NONE | ClippingGeometry}
    clipgeo = row.Shape.Extent
    outputname = "%s.tif" % (n)
    desc = gp.Describe(clipto)
    extent = str(desc.Extent)
    gp.Clip_management(clipit, extent, outputname, clipto, "", None) #update str if necessary
    print n

    # reset the array...
    del rows
    del gp

Answer 3

cross-posting this solution to this thread AND user3397's other thread

@lpinner was on the right track, by stating that the extent being used is the entire extent of the 'clipto' shapefile. This is an underlying attribute of that shapefile, and doesn't change according to the current feature being processed by a cursor. So, every time you're feeding the extent object into the Clip tool, you're giving it this same overall extent, over and over.

You can manipulate the output by changing the Extent environment setting - this is done by setting the gp.Extent property. In my experience, there are two different ways of doing this:

1. the legacy method involved getting and setting the extent as a space-delimited string,
2. the newer (9.3+ I think) method involves instantiating an 'Extent' object, and getting and setting its Xmax, YMax, etc properties

Option 1: legacy method

# Set extent property with a space-delimited string
# extents should be in this order: XMin YMin XMax YMax 
gp.Extent = "-180 -90 180 90"

Option 2: modern method

# Create Extent object, and set its attributes
ext = gp.CreateObject('Extent')
ext.XMin, ext.YMin, ext.XMax, ext.YMax = -180, -90, 180, 90
gp.Extent = ext

I'm not sure which ArcGIS versions restrict you to which method (off the top of my head, I think you can do both in 9.3+, but you can't do Option 2 in 9.2-).

So, to solve your problem, you need to extract the XMin, YMin, XMax and YMax lat/longs from each polygon, and change gp.Extent every time. I.e.:

while row:
    # do your other processing and extract extent co-ords from current polygon        
    ext = gp.CreateObject('Extent')
    ext.XMin, ext.YMin, ext.XMax, ext.YMax = currentXMin, currentYMin, currentXMax, currentYMax
    gp.Extent = ext
    gp.Clip_management #...continue with Clip operation

A third option, I guess, would be to write each polygon to a new temporary feature class. This temporary FC's Extent property would then be the same as the individual polygon's extent. This would be a slow operation if you wrote these temp files to disk, but you can write them to ArcGIS's 'in_memory' workspace to keep it nippy. I.e.:

# Get the primary key field name for the polygon shapefile
oidField = gp.Describe(clipit).OIDFieldName

while row: # Iterate through each polygon
    oid = row.GetValue(oidField) # get the unique id for this polygon

    # Use Select tool to extract individual polygon to temporary in_memory file
    # use the unique ID as the SQL query to extract this record
    gp.Select_management(clipit, 'in_memory/temp', '"%s" = %s' % (oidField, oid))  

    # get the polygon's extent
    ext = gp.Describe('in_memory/temp').Extent

From here, you can either use this in_memory/temp polygon as ClippingGeometry for the raster Clip tool, or you can extract the Extent object's values and feed them into the Clip tool's 'extent' parameters as text.

I'll cut my answer off here, and wait and see if it helps. I can clarify any points some more my referring to my scripts at work - I'm writing this off the top off my head at the moment.

Answer 4

One can easily access extent the same way as geometry is accesed, through row obj.

clipextent=str(row.Shape.Extent)

In arcgis 9.3 you can access each coordinate as Extent.XMin like timmy says. To make it more simple you can just make it as a string and you exactly get 2 coordinate pairs in string as required for clipping.

And when running python script all goes well. Just when added to arcToolbox, it might have errors with dataset in memory.. (as I have it)