Batch Raster Clipping

[Leo_KrisPalao]

Hi ArcGIS Users,

Good Day!

I have been searching for answers in the web but haven't find any plausible answers yet. What I want to do is automate the clipping of my raster files using a set of shapefiles. Basically, what I want to do is:

raster 1 clip with shapefile 1
raster 2 clip with shapefile 2
raster 3 clip with shapefile 3
.
.
.
.
raster 37 clip with shapefile 37

My raster files and shapefiles are in different directory.

My initial code is like the one below which I got from exporting the model builder to a python script.

How can I automate my clipping process in such a way that it will loop to all my raster data and match each raster data with a specified shapefile. I am currently using ArcGIS 10.2. Can you help me build the code using the initial code that I have below?

Thanks for any help.

-Leo

# Import arcpy module
import arcpy


# Local variables:
Shapefile = r'E:\folder\folder\folder\folder\shapefile1.shp'
Raster_Dataset = r'E:\folder\folder\folder\raster1.tif'
Output_Raster_Dataset = r'E:\folder\test.tif'

# Process: Clip
arcpy.Clip_management(Raster_Dataset, "", Output_Raster_Dataset, Shapefile, "", "ClippingGeometry", "NO_MAINTAIN_EXTENT")

[DuncanHornby]

This can all be done in model builder no need to create a python script. You need to read the help file about iterators and inline substitution. Read those pages to understand the model below which will do what you want.

[ATTACH=CONFIG]33864[/ATTACH]

[Leo_KrisPalao]

If I use the model that you suggested, can I input multiple raster layers (*.tif)? and when I run the model, it will automatically match the first shapefile in my folders with my first raster?

I have tried something similar before, but what it does is the first shapefile in my folder reiterates through all of my rasters in my Raster List ( i have used raster layer and enable "List of values" option) like:

SHP 1 clip with raster 1
SHP 1 clip with raster 2
SHP 1 clip with raster 3
.
.
.
SHP 1 clip with raster 37

what I would love to do is to match each SHP with my 1st raster in my list, then 2nd shapefile with my 2nd raster, and so on.

But I will try what I can do in the model builder as you suggested.

Thanks,
-Leo

[DuncanHornby]

This is why I parse the output of the iterator to get the shapefile name so c:\temp\fred.shp will parse to "fred". I then use inline substitution to grab the equivalent raster from c:\raster\fred.tif.

[Anonymous User]

If your file names really are named like "raster1" and "shapefile1" and there is a 1:1 relationship, this code should work:

import arcpy, os

raster_ws = r'E:\path_to_raster_fold'
shp_ws = r'E:\path_to_shp_fold'

# output workspace
out_ws = r'E:\path_to_output_fold'

# get lists
arcpy.env.workspace = raster_ws
rasters = [os.path.join(raster_ws, r) for r in arcpy.ListRasters()]

arcpy.env.workspace = shp_ws
shapes = [os.path.join(shp_ws, s) for s in arcpy.ListFeatureClasses()]

# get dictionary {raster1 : shapefile1.shp, ...}
clip_dict = dict(zip(rasters, shapes))

# clip all reasters
for raster, shapefile in clip_dict.iteritems():
    out_raster = os.path.join(out_ws, os.path.basename(raster))
    arcpy.Clip_management(raster, "", out_raster, shapefile, "", "ClippingGeometry", "NO_MAINTAIN_EXTENT")
    print 'Clipped {0}'.format(os.path.basename(raster))

[Leo_KrisPalao]

Hi Caleb,

Thanks for your reply and sharing the code. I would like to modify the code to incorporate calculation of cell statistics after the clipping process. Would you suggest incorporating it in the code that you gave me, or it will be much simpler to write another code just for calculating cell statistics.

Thanks,
-Leo

[Leo_KrisPalao]

Hi Caleb and other ArcGIS Python users,

I modified the code a bit to incorporate computation of cell statistics in the code. Basically what I have done so far is to create a container >>> outputlist = [] that will store the output of the clipping. Then used this list to compute for cell statistics. I think I got the code correctly until Line 28 >>> outputList.append(out_raster) but I have a problem in computing the cell statistics. The error says:

Traceback (most recent call last):
  File "C:\Python27\ArcGIS10.2\Lib\site-packages\pythonwin\pywin\framework\scriptutils.py", line 326, in RunScript
    exec codeObject in __main__.__dict__
  File "E:\Python_scripts\Clip\rasterClip_with_cellstats.py", line 29, in <module>
    for raster in outputList():
TypeError: 'list' object is not callable

The code with a little modification is below:

import arcpy, os, sys
from arcpy import env
from arcpy.sa import *
arcpy.env.overwriteOutput = True

raster_ws = r'E:\Test_awd\Criteria_SP_2mm_v2'
shp_ws = r'E:\Test_awd\Season1'

# output workspace
out_ws = r'E:\Test_awd'

# get lists
arcpy.env.workspace = raster_ws
rasters = [os.path.join(raster_ws, r) for r in arcpy.ListRasters()]

arcpy.env.workspace = shp_ws
shapes = [os.path.join(shp_ws, s) for s in arcpy.ListFeatureClasses()]
         
# get dictionary {raster1 : shapefile1.shp, ...}
clip_dict = dict(zip(rasters, shapes))

# clip all rasters
outputList = []
for raster, shapefile in clip_dict.iteritems():
    out_raster = os.path.join(out_ws, os.path.basename(raster))
    arcpy.Clip_management(raster, "", out_raster, shapefile, "", "ClippingGeometry", "NO_MAINTAIN_EXTENT")
    print 'Clipped {0}'.format(os.path.basename(raster))
    outputList.append(out_raster)
for raster in outputList():
    sumRaster = CellStatistics(["outputList"], "SUM", "DATA")
    sumRaster.save(r'E:\Test_awd')

Can you help me solve the error that I am getting.

Thanks,
-Leo

[Anonymous User]

Leo,

You were close...You do not need the () after outputList. Try modifying the last portion of the code:

import arcpy, os, sys
from arcpy import env
from arcpy.sa import *
arcpy.env.overwriteOutput = True
arcpy.CheckOutExtension('Spatial')

raster_ws = r'E:\Test_awd\Criteria_SP_2mm_v2'
shp_ws = r'E:\Test_awd\Season1'

# output workspace
out_ws = r'E:\Test_awd'

# get lists
arcpy.env.workspace = raster_ws
rasters = [os.path.join(raster_ws, r) for r in arcpy.ListRasters()]

arcpy.env.workspace = shp_ws
shapes = [os.path.join(shp_ws, s) for s in arcpy.ListFeatureClasses()]
         
# get dictionary {raster1 : shapefile1.shp, ...}
clip_dict = dict(zip(rasters, shapes))

# clip all rasters
outputList = []
for raster, shapefile in clip_dict.iteritems():
    out_raster = os.path.join(out_ws, os.path.basename(raster))
    arcpy.Clip_management(raster, "", out_raster, shapefile, "", "ClippingGeometry", "NO_MAINTAIN_EXTENT")
    print 'Clipped {0}'.format(os.path.basename(raster))
    outputList.append(out_raster)

# output list is already a list, so no need to wrap it in [ ]
sumRaster = CellStatistics(outputList, "SUM", "DATA")
sumRaster.save(r'E:\Test_awd')
print "Created Sum Raster"
arcpy.CheckInExtension('Spatial')

[Leo_KrisPalao]

Hi Caleb,

I tried revising the code this morning and run it again, but the process is stuck on sumRaster = CellStatistics(ouputList, "SUM", "DATA"). Is this an indication that we cannot use list of raster values in Cellstatistics function?

Below is the revised code:

import arcpy, os, sys
from arcpy import env
from arcpy.sa import *
arcpy.env.overwriteOutput = True
arcpy.CheckInExtension("Spatial")

raster_ws = r'E:\Test_awd\Criteria_SP_2mm_v2'
shp_ws = r'E:\Test_awd\Season1'

# output workspace
out_ws = r'E:\Test_awd'

# get lists
arcpy.env.workspace = raster_ws
rasters = [os.path.join(raster_ws, r) for r in arcpy.ListRasters()]

arcpy.env.workspace = shp_ws
shapes = [os.path.join(shp_ws, s) for s in arcpy.ListFeatureClasses()]
         
# get dictionary {raster1 : shapefile1.shp, ...}
clip_dict = dict(zip(rasters, shapes))

# clip all rasters
outputList = []
for raster, shapefile in clip_dict.iteritems():
    out_raster = os.path.join(out_ws, os.path.basename(raster))
    arcpy.Clip_management(raster, "", out_raster, shapefile, "", "ClippingGeometry", "NO_MAINTAIN_EXTENT")
    print 'Clipped {0}'.format(os.path.basename(raster))
    outputList.append(out_raster)

for raster in outputList:
    sumRaster = CellStatistics(outputList, "SUM", "DATA")
    sumRaster.save(r'E:\Test_awd\test.tif')
print "Created Sum Raster"

This is the code that I am getting:

Traceback (most recent call last):
  File "C:\Python27\ArcGIS10.2\Lib\site-packages\pythonwin\pywin\framework\scriptutils.py", line 326, in RunScript
    exec codeObject in __main__.__dict__
  File "E:\Python_scripts\Clip\rasterClip_with_cellstats.py", line 32, in <module>
    sumRaster = CellStatistics(outputList, "SUM", "DATA")
  File "C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy\arcpy\sa\Functions.py", line 2931, in CellStatistics
    ignore_nodata)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy\arcpy\sa\Utils.py", line 47, in swapper
    result = wrapper(*args, **kwargs)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy\arcpy\sa\Functions.py", line 2927, in Wrapper
    [function] + Utils.flattenLists(in_rasters_or_constants))
RuntimeError: ERROR 000824: The tool is not licensed.

I tried to use cell statistics in a separate code, and it works. Hence, my question on using list of rasters as input in Cellstatistics. The code looks like this:

import arcpy
from arcpy import env
from arcpy.sa import *

env.workspace = r'E:\Test_awd\output'

Raster01 = "RE_suitable_dek001.tif"
Raster02 = "RE_suitable_dek002.tif"
Raster03 = "RE_suitable_dek003.tif"
Raster04 = "RE_suitable_dek004.tif"
Raster05 = "RE_suitable_dek005.tif"

arcpy.CheckOutExtension("Spatial")

outCellStatistics = CellStatistics([Raster01, Raster02, Raster03, Raster04, Raster05], "SUM", "DATA")
outCellStatistics.save(r'E:\Test_awd\output.tif')

I think the code is working okay, but get stuck with the line >>> sumRaster = CellStatistics(outputList, "SUM", "DATA")

Thanks,
-Leo

[Anonymous User]

Leo,

When you say "stuck on the the sum raster process" what do you mean? If it just looks unresponsive that is probably because you are running this analysis on 37 rasters, so this could take a long time depending on the resolution of your rasters.

I tried revising the code this morning and run it again, but the process is stuck on sumRaster = CellStatistics(ouputList, "SUM", "DATA"). Is this an indication that we cannot use list of raster values in Cellstatistics function?

I do not think this is the case...From the tool" rel="nofollow" target="_blank">http://resources.arcgis.com/en/help/main/10.1/index.html#//009z... help documents:

[in_raster_or_constant,...] A list of input rasters for which a statistic will be calculated for each cell within the Analysis window.

This code:

# output list is already a list, so no need to wrap it in [ ]
sumRaster = CellStatistics(outputList, "SUM", "DATA")
sumRaster.save(r'E:\Test_awd')

Should be correct because it is passing in a list of rasters for the sum analysis.

This part of your code is not correct:

# you do not want to have this in a loop, remove the for loop
for raster in outputList: # take this out
    sumRaster = CellStatistics(outputList, "SUM", "DATA")
    sumRaster.save(r'E:\Test_awd\test.tif')

because you are running the the cell statistics on the entire list of rasters for every iteration of your list (which would be 37 times if there are still 37 rasters in your list). You are also overwriting the test.tif every time as well.

Also, judging from your error:

Traceback (most recent call last): 
File "C:\Python27\ArcGIS10.2\Lib\site-packages\pythonwin\pywin\framework\scriptutils.py", line 326, in RunScript 
exec codeObject in __main__.__dict__ 
File "E:\Python_scripts\Clip\rasterClip_with_cellstats.py", line 32, in <module> 
sumRaster = CellStatistics(outputList, "SUM", "DATA") 
File "C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy\arcpy\sa\Functions.py", line 2931, in CellStatistics 
ignore_nodata) 
File "C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy\arcpy\sa\Utils.py", line 47, in swapper 
result = wrapper(*args, **kwargs) 
File "C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy\arcpy\sa\Functions.py", line 2927, in Wrapper 
[function] + Utils.flattenLists(in_rasters_or_constants)) 
RuntimeError: ERROR 000824: The tool is not licensed.

import arcpy, os, sys
from arcpy import env
from arcpy.sa import *
arcpy.env.overwriteOutput = True
arcpy.CheckInExtension("Spatial")  

The line highlighted in red needs to be changed to:

arcpy.CheckOutExtension("Spatial")

Because you were checking the license in before running the tool, therefore you were now allowing access to spatial analyst tools. The full code I supplied in my last post should do the trick. I would suggest trying that and if it is hanging up on the sum statistics part, try running that as a separate script.