nav[aria-label="Primary Navigation"] { padding: 0; & ul { list-style: none; width: 100%; display: flex; flex-direction: row; justify-content: start; align-items: start; gap: 30px; padding: 0; & li { margin: 0; } & ul li { list-style: none; } } }

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

Trouver plus de publications taguées avec

Accepted answers

All comments

JamieWilliams

I had the same problem. I solved it by adding the full file path to the feature name.
eg. pathFeature = "R:\\folderContainingShpFile\\" + str(feature)

gp.clip_management(raster,"446716,298473 4425996,437405 499153,348855 4492682,68626",outRaster,pathFeature,"255","ClippingGeometry")

ChrisSnyder

Not sure what's wrong... Take look at the sample code: http://forums.arcgis.com/threads/9512-Rectangle-for-Clip_management-tool-in-arcpy (my post Aug 2nd @10:25).

Are you sure that your raster and shapefile are using the same projection?

bogdanpalade1

Tried changing the shapefile name but still the same error....

ChrisSnyder

Try naming your shapefile something else than a number (1.shp). Generally it is a bad idea to name a layer or field with a number as the 1st character.

bogdanpalade1

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

DanPatterson_Retired

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

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

error.JPG

bogdanpalade1

Thank you Chris,

I will check it out now...

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

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

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

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.

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.

Rejoindre la conversation

Les membres connectés peuvent publier, suivre les mises à jour, et plus encore. Nouveau ici ? Inscrivez-vous gratuitement.

Liens rapides

Community Help
Find useful guides, FAQs, and documents to help you navigate and make the most of Esri Community.