Problems if "Extract by Mask script loop"

deleted-user-22pwUMwu80-d · ‎01-10-2011

I'm a beginner using Python and I have some doubts.
My first script is write for extract some areas of interest (vector) of many topographic charts (raster) using "Extract By Mask"....and put a output name similar to the input name of the charts...

First I set the input folder "E:/GIS_Mestrado/Corumbatai/3_Cartas_um_metro"...there is about 80 charts...
The areas of interest were set in "E:/GIS_Mestrado/Corumbatai/Geodados/Grid_Cartas_sem_Grade.shp"...

And a have to loop the extract by mask in all the rasters in input folder...
Anyone can help me ??

Thanks...
PS: I'm using ArcGIS 9.2 yet...

Below is the script...
---------------------------
# Extracts the cells of a raster that correspond with the areas
# defined by a mask.
# Author: Michel
# Date: 10 de Janeiro de 2011
# -----------------------------

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

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

# Check out any necessary licenses
gp.CheckOutExtension("spatial")

# Load required toolboxes...
gp.AddToolbox("C:/Arquivos de programas/ArcGIS/ArcToolbox/Toolboxes/Spatial Analyst Tools.tbx")

GP.workspace = "E:/GIS_Mestrado/Corumbatai/3_Cartas_um_metro"
out_workspace = "E:/GIS_Mestrado/Corumbatai/4_Cartas_Con/"
mascara = "E:/GIS_Mestrado/Corumbatai/Geodados/Grid_Cartas_sem_Grade.shp"

# Get a list of the rasters in the workspace
rasters = GP.ListRasters()

# Loop through the list of rasters
rasters.reset()
raster = rasters.next()

while raster:
    # Set the outputname for each output to be the same as the input
    output = out_workspace + raster

    # Process: Extract by Mask...
    gp.ExtractByMask_sa(raster, mascara, output)

    # Loop da função...
    raster = rasters.next()

except:
    # If an error occurred while running a tool, then print the messages.
    print gp.GetMessages()

BradPosthumus · ‎01-11-2011

There are two errors in your script that I notice:

1. You have an "except" clause but there is no corresponding "try" statement. I've added it in the code below to encompass all geoprocessing statements.
2. Python is case-sensitive, therefore the variable "gp" is not the same as "GP".

It's also a good idea to include a generic error statement in the "except" clause. If there's a Python error in your script that doesn't show up in "gp.GetMessages()" your script will finish with no indication that it actually crashed.

One more thing: in this forum you can wrap your code in code tags (by selecting the code and hitting the "#" symbol above the text box) to preserve your formatting.

# Extracts the cells of a raster that correspond with the areas
# defined by a mask.
# Author: Michel
# Date: 10 de Janeiro de 2011
# -----------------------------

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

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

try:
    # Check out any necessary licenses
    gp.CheckOutExtension("spatial")

    # Load required toolboxes...
    gp.AddToolbox("C:/Arquivos de programas/ArcGIS/ArcToolbox/Toolboxes/Spatial Analyst Tools.tbx")

    gp.workspace = "E:/GIS_Mestrado/Corumbatai/3_Cartas_um_metro"
    out_workspace = "E:/GIS_Mestrado/Corumbatai/4_Cartas_Con/"
    mascara = "E:/GIS_Mestrado/Corumbatai/Geodados/Grid_Cartas_sem_Grade.shp"

    # Get a list of the rasters in the workspace
    rasters = gp.ListRasters()

    # Loop through the list of rasters
    rasters.reset()
    raster = rasters.next()

    while raster:
        # Set the outputname for each output to be the same as the input
        output = out_workspace + raster

        # Process: Extract by Mask...
        gp.ExtractByMask_sa(raster, mascara, output)

        # Loop da função...
        raster = rasters.next()

except:
    # If an error occurred while running a tool, then print the messages.
    print "Error in script"
    print gp.GetMessages()

ChrisSnyder · ‎01-11-2011

You might want to use the Clip_management tool (clips rasters w/ option to use a polygons). Note this is different from the Clip_analysis tool which is for vectors.

If your charts are in a color mapped multiband image format (.tif, .jpg, etc), you will find that the the Clip_management tool is much easier to deal with.

Also, be sure to use the "gp.snapraster" environment setting to ensure that the output rasters have the same cell alignment as the input rasters.

So it might look something like this:

import sys, string, os, arcgisscripting
gp = arcgisscripting.create()
try:
    gp.CheckOutExtension("spatial")
    gp.AddToolbox("C:/Arquivos de programas/ArcGIS/ArcToolbox/Toolboxes/Spatial Analyst Tools.tbx")
    gp.workspace = "E:/GIS_Mestrado/Corumbatai/3_Cartas_um_metro"
    out_workspace = "E:/GIS_Mestrado/Corumbatai/4_Cartas_Con/"
    mascara = "E:/GIS_Mestrado/Corumbatai/Geodados/Grid_Cartas_sem_Grade.shp"
    rasters = gp.ListRasters()
    raster = rasters.next()
    while raster:
        gp.extent = gp.describe(raster).extent
        gp.snapraster = raster
        output = out_workspace + raster + ".tif" #reformat as a .tif or whatever - it's that easy!
        gp.Clip_management(raster, gp.extent, output, mascara, "", "CLIPPINGGEOMETRY")
        gp.extent = "" #clear the extent
        raster = rasters.next()
except:
    print "Error in script"
    print gp.GetMessages().

ChrisSnyder · ‎01-11-2011

You might want to use the Clip_management tool (clips rasters w/ option to use a polygons). Note this is different from the Clip_analysis tool which is for vectors.

If your charts are in a color mapped multiband image format (.tif, .jpg, etc), you will find that the the Clip_management tool is much easier to deal with.

Also, be sure to use the "gp.snapraster" environment setting to ensure that the output rasters have the same cell alignment as the input rasters.

So it might look something like this:

import sys, string, os, arcgisscripting
gp = arcgisscripting.create()
try:
    gp.CheckOutExtension("spatial")
    gp.AddToolbox("C:/Arquivos de programas/ArcGIS/ArcToolbox/Toolboxes/Spatial Analyst Tools.tbx")
    gp.workspace = "E:/GIS_Mestrado/Corumbatai/3_Cartas_um_metro"
    out_workspace = "E:/GIS_Mestrado/Corumbatai/4_Cartas_Con/"
    mascara = "E:/GIS_Mestrado/Corumbatai/Geodados/Grid_Cartas_sem_Grade.shp"
    rasters = gp.ListRasters()
    raster = rasters.next()
    while raster:
        gp.extent = gp.describe(raster).extent
        gp.snapraster = raster
        output = out_workspace + raster + ".tif" #reformat as a .tif or whatever - it's that easy!
        gp.Clip_management(raster, gp.extent, output, mascara, "", "CLIPPINGGEOMETRY")
        gp.extent = "" #clear the extent
        raster = rasters.next()
except:
    print "Error in script"
    print gp.GetMessages()

BradPosthumus · ‎01-11-2011

In 9.2 you can only clip rasters with rectangular coordinates using the Clip_management tool.

I had no idea they changed it so you can use polygons in 9.3 and beyond. I gave up on using that tool back in the 9.2 days, but this makes it much more useful and you don't need to the Spatial Analyst extension. Good to know!

deleted-user-22pwUMwu80-d · ‎01-11-2011

Thanks for all,
The script posted by bposthumus works fine...

The charts are single band..anyway, I try the second script and it's doesn't work...I search for error, but I don't fine nothing...I think that is because is single band...

In next time I use "#" for the codes..sorry...

Thanks!

SergeRafanoharana · ‎01-11-2012

Great script from you guys... I tried the Python one and it worked perfectly... Thumbs up...

Problems if &quot;Extract by Mask script loop&quot;

Problems if "Extract by Mask script loop"