Storing Rasters in an Array

[Kevin_Smith]

As I am new to this, please forgive me for any glaring errors.  I am trying to loop through a series of folders find the "wsgrid004" raster, clip it using a polygon, store the results in an array. Once all the rasters have been processed, mosaic them together.

All I get is a general syntax error, any help would be appreciated!

# Import modules import arcpy  # Set enviroment settings env.workspace = Input_Model  # Check out the ArcGIS Spatial Analyst extension license arcpy.CheckOutExtension("Spatial")  # Load Required Toolboxes arcpy.AddToolBox (r"C:\Program Files (x86)\ArcGIS\Desktop10.0\ArcToolbox\Toolboxes\Data Management Tools.tbx")  #Enable the Overwriting of Existing GeoProcessing Results gp.overwrite = "True"  #Local Variables Count = 0 array = arcpy.Array() wse_grid = "wsgridp004 + Count" Folders = arcpy.ListWorspaces("*", "Folder")  # Argument 1 is the GeoRAS Directory Input_Model = arcpy.GetParameterAsText(0) # Argument 2 is the 100-Year Floodplain FP_Poly = arcpy.GetParameterAsText(1)  # Lists all folders in the Input_Model ListWorkspaces("*","Folder")  print "Looking for Sub-Directories" for Folder in Folders:     #Find the wsgridp004 Raster     print "Finding 100 YR WSEL Raster"     rasters = arcpy.ListRasters("wsgridp004", "Grid")                                               for "wsgridp004" in rasters:         Count = Count + 1         #Clip Raster with 100-Yr Floodplain Polygon         print "Clipping Raster"         arcpy.Clipmanagement("rasters", "#", "wse_grid", "FP_Poly", "0", "ClippingGeometry")         print arcpy.GetMessages()         array.add(wse_grid) #Mosaics Clipped 100-Yr Floodplain Grids print "Mosaicing Clipped Water Surface Grids" arcpy.MosaictoNewRaster_Management("array.next", "Input_Model", "wse_100_Yr","NAD_1983_StatePlane_Texas_Central_FIPS_4203_Feet", "32_BIT_FLOAT", "3", "1", " ", " " ) print arcpy.GetMessages() print "Complete"

[JimCousins]

2 syntax errors:
for "sswgridp004" in rasters
should be
for "sswgridp004" in rasters:
Print "Complete"
should be
print "Complete

Regards,
Jim

[XanderBakker · Accepted solution]

Hi Kevin,

I hate to break the news to you, but there is a lot more wrong with the script then what Jim is mentioning. I'll mention a few points:

• You access the environment setting "env" before they have been initialized ("from arcpy import env" is missing)


• You set the workspace property of the environment setting using a variable "Input_Model", which is not available at that point, since it is obtained from the first parameter later in the script


• You use "gp.overwrite", but gp is not set (you should use the "env" for this purpose and set it properly). The overwrite property does not exist this should be overwriteOutput


• You check out a Spatial Analyst extension, but the tools you're using do not require the sa extension.


• You load the Data Management Tools toolbox, but this is really not necessary (was done in previous versions)


• You set a variable wse_grid equal to a hard coded string including " + Count", when you probably want to add the count to the string


• As mentioned before your first argument is read after its first use "Input_Model = arcpy.GetParameterAsText(0)" should move to the top


• You preform an "arcpy.ListRasters". This works on the current workspace. To do this for each folder, you should set the env.workspace for each folder


• Next you  loop through the rasters with the statement for "wsgridp004" in rasters:  This will not work. Normally you would use something like for raster in rasters: and then access the raster object in the following steps


• In the arcpy.Clipmanagement statement, you specify strings, when you should specify the variables.


• arcpy.Clipmanagement does not exits, it should be arcpy.Clip_management


• "FP_Poly" is specified as hard coded string, but should be specified as variable (without double quotes)


• Next you add a variable to the array with "array.add", but this variable is never set correctly (the initial wrong value is assigned multiple times)


• In the arcpy.MosaicToNewRaster_management you again specify hard coded strings, when you should specify variables (for instance Input_Model should be without quotes.


• The string "array.next" should be a list of rasters separated by semicolons. Since you obtain it from different workspaces, you should include the path for each raster

To get you started I've included some code. The code however has not been tested...

import arcpy import os from arcpy import env env.overwriteOutput = True  # Argument 1 is the GeoRAS Directory Input_Model = arcpy.GetParameterAsText(0) # Argument 2 is the 100-Year Floodplain FP_Poly = arcpy.GetParameterAsText(1)   Count = 0 lstRas = '' rasName = 'wsgridp004'  # set workspace for list folders env.workspace = Input_Model Folders = arcpy.ListWorspaces("*", "Folder")  for Folder in Folders:     env.workspace = Folder     #Find the wsgridp004 Raster     rasters = arcpy.ListRasters(rasName , "GRID")     for raster in rasters:         Count += 1         wse_grid = Input_Model + os.sep + "{0}{1}".format(rasName ,Count)          #Clip Raster with 100-Yr Floodplain Polygon         arcpy.Clip_management(raster, "#", wse_grid, FP_Poly, "0", "ClippingGeometry", "NO_MAINTAIN_EXTENT")         print arcpy.GetMessages()         if lstRas == '':             lstRas = wse_grid         else:             lstRas = lstRas + ";" + wse_grid  #Mosaics Clipped 100-Yr Floodplain Grids print "Mosaicing Clipped Water Surface Grids" arcpy.MosaicToNewRaster_management(lstRas, Input_Model, "wse_100_Yr","NAD_1983_StatePlane_Texas_Central_FIPS_4203_Feet", "32_BIT_FLOAT", "3", "1") print arcpy.GetMessages() print "Complete"

Assuming from the line where you load the toolbox that you are using ArcGIS 10.0 you should look into this topic:
A quick tour of ArcPy

The training.esri.com site offers some Python/arcoy courses that might be interesting:
http://training.esri.com/gateway/index.cfm?fa=search.results&searchterm=arcpy&search=Search&CourseTy...


Kind regards,

Xander

[Kevin_Smith]

Xander,

Thanks for the suggestions, the first round of code was a rough. After Jim's help, I went back and spent a fair amount of time going back through the code and making changes that addressed a fair amount your points. However, I was having a rough time with the list rasters and array but you gave excellent tips.

Thanks!

[Kevin_Smith]

Xander,

Thanks for you help, the advice about environments and looping was extremely helpful. I have run into another problem. I keep on getting the below error:
<class 'arcgisscripting.ExecuteError'>: Failed to execute. Parameters are not valid. ERROR 000157: Input and target dataset should have the same number of bands Failed to execute (MosaicToNewRaster).

When I run the tool in ArcMap with the exact same settings, the tool completes no problems. Any thoughts?

# Import modules
import arcpy
import os

# Set enviroment settings
from arcpy import env
env.overwriteOutput = True


# Argument 1 is the GeoRAS Directory
Input_Model = arcpy.GetParameterAsText(0)
# Argument 2 is the 100-Year Floodplain
FP_Poly = arcpy.GetParameterAsText(1)

#Local Variables
Count = 0
P004 = 'wsgridp004'
lstRstr = ' ' 

#Set Workspace for list folders
env.workspace = Input_Model
Folders = arcpy.ListWorkspaces("*", "Folder")

print "Looking for Sub-Directories"
for Folder in Folders:
    env.workspace = Folder
    #Find the wsgridp004 Raster
    print "Finding 100 YR WSEL Raster"
    rasters = arcpy.ListRasters(P004, "Grid")                                          
    for raster in rasters:
        Count = Count + 1
        wse_grid = Input_Model + os.sep + "{0}{1}".format(P004 ,Count)
        #Clip Raster with 100-Yr Floodplain Polygon
        print "Clipping Raster"
        arcpy.Clip_management(raster, "#", "wse_grid", FP_Poly, "0", "ClippingGeometry")
        print arcpy.GetMessages()
        if lstRstr == '':
            lstRstr = wse_grid
        else:
            lstRstr = lstRstr + ";" + wse_grid
#Mosaics Clipped 100-Yr Floodplain Grids
print "Mosaicing Clipped Water Surface Grids"
arcpy.MosaicToNewRaster_management(lstRstr, Input_Model,"wse_100_Yr"," ", "32_BIT_FLOAT", "3", "1", "MAXIUM", "FIRST")
print arcpy.GetMessages()
print "Complete"

[TimBarnes]

It is telling you the problem here: Input and target dataset should have the same number of bands 🙂

When you run a mosaic through ArcCatalog, it will default to the number of bands in your input datasets. In your script you have defined it explicitly to 3 bands, such as in an RGB image, even though it appears you are using grid rasters which have only one band.

[XanderBakker]

Hi Kevin,

Tim is right about the number of bands which should be the same for all raster you want to mosaic together. The syntax of Mosaic To New Raster (Data Management) is this:

MosaicToNewRaster_management (input_rasters, output_location, raster_dataset_name_with_extension, {coordinate_system_for_the_raster}, {pixel_type}, {cellsize}, number_of_bands, {mosaic_method}, {mosaic_colormap_mode})

In my script this is what is assigned:

input_rasters = lstRas
output_location = Input_Model
raster_dataset_name_with_extension = "wse_100_Yr"
{coordinate_system_for_the_raster} = "NAD_1983_StatePlane_Texas_Central_FIPS_4203_Feet"
{pixel_type} = "32_BIT_FLOAT"
{cellsize} = "3"
number_of_bands = "1"
{mosaic_method} = omitted
{mosaic_colormap_mode} = omitted

This was taken from your original settings. So you are executing a mosaic using a cell size of 3 and 1 band. Verify your input rasters. If the settings are not correct, change them in the script. If the number of bands in the rasters is not the same, you will not be able to mosaic them together.

Kind regards,

Xander

[Kevin_Smith]

Sorry for the delay in my response, you know how it goes always a thousand different things to do. 🙂

I originally set the number of bands by looking at raster data set properties. I thought maybe I had a rouge raster with a different amount of bands, so I went back and checked every input raster and they are all the same. (See below)
[ATTACH=CONFIG]28511[/ATTACH]

Any thoughts??

[XanderBakker]

Any thoughts??

Hi Kevin,

I had a closer look at the script, and noticed you changed a few things.

The line:

lstRstr = ' '

should be:

lstRstr = ''

You initially set the variable to single space, when it should be an empty string. Later the code checks for an empty string and will not find it. This causes the list of rasters to contain a "space" for the first raster name.


The line:

arcpy.Clip_management(raster, "#", "wse_grid", FP_Poly, "0", "ClippingGeometry")

should be:

arcpy.Clip_management(raster, "#", wse_grid, FP_Poly, "0", "ClippingGeometry")

In this case the output raster will have a unique name for each raster and not always "wse_grid".


The line:

arcpy.MosaicToNewRaster_management(lstRstr, Input_Model,"wse_100_Yr"," ", "32_BIT_FLOAT", "3", "1", "MAXIUM", "FIRST") 

should be:

arcpy.MosaicToNewRaster_management(lstRstr, Input_Model,"wse_100_Yr"," ", "32_BIT_FLOAT", "3", "1", "MAXIMUM", "FIRST") 

The provided mosaic_method was not valid due to typo. I also notice that you blanked the output coordinate system. This may not be a problem if you environment settings are correct.

Try again and see what happens.

Kind regards,

Xander

[Kevin_Smith]

D'oh, I cannot believe that I missed spelling and variable errors. Although, I would have missed the spacing in the quotes. I appreciate everyone's help, the script is up and running as we speak! Once again, many thanks!