Raster ExtractByMask with python

[bogdanpalade1]

Hello,

I am trying to clip a raster with a polygon shapefile.The shapefile contains several polygons representing administrative regions. The method is to temporarily export each polygon to a new shapefile and to extract the coresponding piece of raster.

The problem is that the output of

gp.ExtractByMask_sa(raster, shapefile, newRaster) 

is always just one band of the original RGB raster. Do you know how to extract all the 3 bands?

The original raster is .img

Thank you

[ChrisSnyder]

Check out the Clip_managment tool (not the Clip_analysis one). Be sure to set the snap raster in the Envr. Settings, otherwise the pixel allignment may differ from the original.

[bogdanpalade1]

Hello chris,

Sorry for the late reply. Clip_management asks for a rectangle to do the clip...instead I want to use a shapefile to do the clip with.

[DanPatterson_Retired]

http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//009z0000002n000000.htm
indicates that a grid stack is the default output, I am not using 10 so I don't know what other output options there, so check those options

[bogdanpalade1]

Dan,

Thank you for the reply but this works only in ArcGIS 10 environment (by the way, congrats to esri team for that although I still don't understand why only 3 bands can be extracted and not all of them since I guess that not all the people are working with RGB or pan images). I'm using 9.1 and 9.3 at the time so...I have to deal with it for the time being.

In the end I think I will have to do it separatly band by band and run some kind of layer stack afterwards....It will be soooo....slow!

Thanks,
Bogdan

[ChrisSnyder]

In v9.3.1 ArcInfo level, the Clip_managment tool does indeed have an option for using a vector featureclass as a psudo-mask (see the "in_template_dataset" parameter). This is tool to use when you have a multi-band image, since the output can be any ESRI-supported format (not just GRID like the ExtractByMask tool). Check it out...

[bogdanpalade1]

Thank you Chris,

I will check it out now...

[bogdanpalade1]

Ok, so the clip_management allows me to use a feature class extent to clip the raster but it's still necessary to introduce the rectangle parameters (I thought that I could find something to substitute the rectangle by a feature class).
Well, the initial purpose was to cut some rasters with a polygon feature class (administrative boundaries). The approach was first to extract single polygons to new feature classes and to use these to cut the rasters located in a specific folder. When I thought all was great I ran into an error: after extracting the single polygons to new shapefiles the clipping process is interupted by the following message:
ERROR 000732: Output Extent: Dataset 1.shp does not exist or is not supported
Failed to execute (Clip).
I attached a print screen of the error message.

The script is the following:

import arcgisscripting, os, string
gp = arcgisscripting.create()
gp.AddToolbox("C:/Program Files/ArcGIS/ArcToolbox/Toolboxes/Data Management Tools.tbx")
gp.SetProduct("ArcInfo")


clip = "D:/hyperspectralImages/torrejon/torrejon_ALI/clip/clip.shp"   #the polygon feature class
rows = gp.SearchCursor(clip)
row = rows.Next()
while row:  # here I am extracting each polygon to an individual feature class
    FID = row.GetValue("FID")
    hoja = row.GetValue("hoja")
    Expression = '\"FID\" = ' + str(FID)
    outFeatures = "D://hyperspectralImages//torrejon//torrejon_ALI//outClips//"+str(hoja)+".shp" #the single fc's will be named after the field "hoja"
    print "Extracting sheet number "+str(hoja)                                                   #like 1.shp, 2.shp...etc
    gp.Select_analysis(clip,outFeatures,Expression)
    row=rows.Next()
print "Finished extracting individual sheets.Begining to clip..."

inputFolder= "D:/hyperspectralImages/torrejon/torrejon_ALI/rasters" #the folder where all the rasters to be clipped are located
gp.Workspace = inputFolder
rasters = gp.ListRasters()
rasters.Reset()
raster = rasters.Next()

clipFeatures = "D://hyperspectralImages//torrejon//torrejon_ALI//outClips//" #the folder where the polygon feature classes used for clipping are located

while raster: # the loop to clip each raster by each polygon feature class
    for feature in os.listdir(clipFeatures):
        if feature.endswith(".shp"):
            outRaster = "D://hyperspectralImages//torrejon//torrejon_ALI//output//"+feature[0:-3]+".tif"
            print "Clipping "+raster+" with sheet number "+feature[0:-3]
            gp.clip_management(raster,"446716,298473 4425996,437405 499153,348855 4492682,68626",outRaster,feature,"255","ClippingGeometry")
                
print "Done"

The 1.shp, 2.shp ..and so on do exist in the right folder...

Do you have any idea of what's happening?

Thanks,
Bogdan

[DanPatterson_Retired]

Python paths use a single forward slash /
or double back slash \\
not double forward slash //
as you have in several places

[bogdanpalade1]

Dan,
Thank you for your advice. I changed the paths accordingly but the problem remains the same.