topic Batch Raster Clipping in Python Questions

Batch Raster Clipping

Leo_KrisPalao — Fri, 16 May 2014 08:36:30 GMT

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")

Re: Batch Raster Clipping

DuncanHornby — Fri, 16 May 2014 09:38:30 GMT

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]

Re: Batch Raster Clipping

Leo_KrisPalao — Fri, 16 May 2014 11:05:23 GMT

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

Re: Batch Raster Clipping

DuncanHornby — Fri, 16 May 2014 11:58:57 GMT

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.

Re: Batch Raster Clipping

Anonymous User — Sat, 11 Dec 2021 09:44:13 GMT

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))

Re: Batch Raster Clipping

Leo_KrisPalao — Tue, 20 May 2014 08:24:35 GMT

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

Re: Batch Raster Clipping

Leo_KrisPalao — Sat, 11 Dec 2021 09:44:15 GMT

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

Re: Batch Raster Clipping

Anonymous User — Sun, 12 Dec 2021 16:10:03 GMT

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')

Re: Batch Raster Clipping

Leo_KrisPalao — Sat, 11 Dec 2021 09:44:18 GMT

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

Re: Batch Raster Clipping

Anonymous User — Sat, 11 Dec 2021 09:44:20 GMT

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#//009z0000007q000000]tool 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.

Re: Batch Raster Clipping

Leo_KrisPalao — Sat, 11 Dec 2021 09:44:23 GMT

Hi Caleb,

Thanks for your great help. The problem with the code is the loop in the cell statistics portion and the arcpy.CheckInExtension("Spatial). I changed the code as you suggested and it worked.

I am sharing the final code so other users can use it:

Again, thanks for your help.
-Leo

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

raster_ws = r'C:\raster\folderpath'
shp_ws = r'E:\shapefile\folderpath'

# output workspace
out_ws = r'C:\outputpath'

# 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 = [] # container
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) # aopend all the output to the container

# compute the sum of all rasters with ignore no data option
sumRaster = CellStatistics(outputList, "SUM", "DATA")
sumRaster.save(r'C\outputpath\ofsumrasters.tif') # remove .tif if GRID format is desired
print "Created Sum Raster"

Re: Batch Raster Clipping

Leo_KrisPalao — Sat, 11 Dec 2021 09:44:26 GMT

Hi Caleb,

Sorry again for this request. I am again needed for some assistance in completing the python code that I need.

In the code that we completed before, I am trying to reiterate the resulting raster (one raster only) to multiple shapefiles (provinces).

I added some code below, but I am having an error, i.e.,

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\clip_raster_with_multifeatures_v2.py", line 46, in <module>
    arcpy.Clip_management(sumRaster, "", out_prov, fc, "", "clippingGeometry", "NO_MAINTAIN_EXTENT")
File "C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy\arcpy\management.py", line 12803, in Clip
    raise e
ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000445: Extension is invalid for the output raster format.
Failed to execute (Clip).

I think the error is something to do with the filename.

Below is the code the modified code that I am trying to run:

import arcpy, os, sys
from arcpy import env
from arcpy.sa import *
arcpy.env.overwriteOutput = True
arcpy.CheckOutExtension("Spatial")
arcpy.env.workspace = r'E:\output_workspace'
arcpy.env.extent = "MAXOF"

raster_ws = r'E:\raster_workspace'
shp_ws = r'E:\shape_workspace'

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

# 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 = [] # container
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) # append all the output to the container

# compute the sum of all rasters with ignore no data option
sumRaster = CellStatistics(outputList, "SUM", "DATA")
sumRaster.save(r'E:\output_workspace\output.tif')
print "Created Sum Raster"

# clip sumRaster with shapefiles
province = r'E:\workspace_of_shapefiles_per_province'

# Process: Iterate Feature Classes
arcpy.env.workspace = province
fcList = [os.path.join(province, p) for p in arcpy.ListFeatureClasses()] # Get list
for fc in fcList:
    out_prov = os.path.join(out_ws, os.path.basename(fc))
    arcpy.Clip_management(sumRaster, "", out_prov, fc, "", "clippingGeometry", "NO_MAINTAIN_EXTENT")

Thanks very much,
-Leo

Re: Batch Raster Clipping

Leo_KrisPalao — Sat, 11 Dec 2021 09:44:29 GMT

Hi Caleb,

I am almost corrected the code. It seems to work properly. But my problem is the naming of the file. I am having trouble removing the .shp extension in my features list. For instance my filename input is out_prov + "clip.tif". When it saved the file it carries the .shp extension name of my feature list.

This is the output file name that I get: province.shpclip.tif

What I like to have is provinceclip.tif

Here is the last 10 lines of the code that I am using:

# clip sumRaster with shapefiles
province = r'E:\path\province_shapefiles'

# Process: Iterate Feature Classes
arcpy.env.workspace = province
fcList = [os.path.join(province, p) for p in arcpy.ListFeatureClasses()] # Get list
for fc in fcList: # clip the sumRaster with all the shapefiles in the folder
    out_prov = os.path.join(out_ws, os.path.basename(fc))
    arcpy.Clip_management(sumRaster, "", out_prov + "clip.tif", fc, "", "clippingGeometry", "NO_MAINTAIN_EXTENT")
    print "scoring maps clipped by province"

Thanks,
-Leo

Re: Batch Raster Clipping

Anonymous User — Sat, 11 Dec 2021 09:44:31 GMT

Leo,

Try changing the line in red:

# clip sumRaster with shapefiles
province = r'E:\path\province_shapefiles'

# Process: Iterate Feature Classes
arcpy.env.workspace = province
fcList = [os.path.join(province, p) for p in arcpy.ListFeatureClasses()] # Get list
for fc in fcList: # clip the sumRaster with all the shapefiles in the folder
    out_prov = os.path.join(out_ws, os.path.splitext(os.path.basename(fc))[0])
    arcpy.Clip_management(sumRaster, "", out_prov + "clip.tif", fc, "", "clippingGeometry", "NO_MAINTAIN_EXTENT")
    print "scoring maps clipped by province"

Also, you should mark this thread as answered to let others know your issue has been resolved. If you do not know how, see this [url=http://resources.arcgis.com/en/help/forums-mvp/]page[/url].

Re: Batch Raster Clipping

Leo_KrisPalao — Tue, 10 Jun 2014 00:58:35 GMT

Hi Caleb,

Sorry for this late reply. Why is it in the thread there is no option to mark the reply as best (check mark). I can only vote your reply but not mark it as best to signify that my query has been resolved.

Thanks,
-Leo

Re: Batch Raster Clipping

CharlottePeters95540 — Tue, 21 May 2024 19:38:22 GMT

Even though your suggestion was posted 10 years ago, this worked great for me! I was pulling my hair out trying to get a sub-model with an iterator to work within a model that also had an iterator, but your suggestion was so much easier to implement. Brilliant workaround.