Map Algebra and Python?

ArvindBhuta · ‎09-02-2011

Hey all,

I'm trying to input two rasters and derive a "mean" raster from both of them. The thing is I need to do this 4500 times. Any help would be appreciated. I've been messing around with a code, but I don't even know if I have it right. Right now, it runs successfully but no output and I'm thinking this is because I forgot to define something?

Also, I don't know even if I can run it like the way I envision it?
What I'm trying to do is
mean (1895_01Tmin, 1895_02Tmax) = 1895_01Tmean
mean (1895_02Timin, 1895_02tmax) = 1895_02Tmean
....
mean (XXX_XXTmin, XXXX_XXTmax) = XXXX_XXTmean

I posted my code below and I really appreciate any help!

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

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

# Parameters
Workspace01 = gp.GetParameterAsText(0)
gp.AddMessage ("input Workspace01: " + Workspace01)
#Workspace01 = "C:\test\ppt_study_counties"

Workspace02 = gp.GetParameterAsText(1)
gp.AddMessage ("input Workspace02: " + Workspace02)
#Workspace02 = "C:\test\ppt_study_counties"

output_Workspace = gp.GetParameterAsText(2)
gp.AddMessage ("output Workspace: " + output_Workspace)
#OutWorkspace01 = "C:\test\zonal_stats_ppt_counties"

try:
    #Get the raster datasets in the input Workspace01 and loop through them from the start
    gp.Workspace01 = Workspace01
    InputRasters01 = gp.ListRasters01()
    gp.AddMessage(InputRasters01)
    InputRasters01.reset()
    InputRaster01 = InputRasters01.next()
    gp.AddMessage(InputRaster01)

    gp.Workspace02 = Workspace02
    InputRasters02 = gp.ListRasters02()
    gp.AddMessage(InputRasters02)
    InputRasters02.reset()
    InputRaster02 = InputRasters02.next()
    gp.AddMessage(InputRaster02)

    while InputRaster:

        #Set the output table name
        outRaster = OutputWorkspace + "\\" + Workspace01 + Workspace02 + "tmen"
        InExpression = "MEAN(Workspace01, Workspace02)"
        gp.AddMessage(outRaster)

        # Check out Spatial Analyst extension license
        #gp.CheckOutExtension("Spatial")

        # Load required toolboxes. You may need to alter the paths if you are using a 64 bit OS
        #gp.AddToolbox("C:\Program Files (x86)\ArcGIS\ArcToolbox\Toolboxes\Spatial Analyst Tools.tbx")

        # Process: MapAlgebra
        gp.SingleOutputMapAlgebra_sa(InExpression, Outfile)
        # print InputRaster
        InputRaster = InputRasters.next()


except:
    gp.AddMessage(gp.GetMessages(2))
    print gp.GetMessages(2)

TonyPagani · ‎09-08-2011

Arvind,

Here is a script which should give you the mean of 2 rasters. The number of rasters in each of the workspaces must be the same or the script will error.

Regards,
Tony

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

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

# Check out Spatial Analyst extension license
gp.CheckOutExtension("Spatial")


# Parameters
Workspace01 = gp.GetParameterAsText(0)
gp.AddMessage ("input Workspace01: " + Workspace01)
#Workspace01 = r"C:\temp\test\Workspace1"

Workspace02 = gp.GetParameterAsText(1)
gp.AddMessage ("input Workspace02: " + Workspace02)
#Workspace02 = r"C:\temp\test\Workspace2"

OutWorkspace = gp.GetParameterAsText(2)
gp.AddMessage ("output Workspace: " + output_Workspace)
#OutWorkspace = r"C:\temp\test\outputwkspc"

try:
#Get the raster datasets in the input Workspace01 and loop through them from the start
    #I removed the messages for easier reading of the code.
    #Set the environment workspace to workspace one. gp.workspace is an environment variable and must be set befor listing data
    gp.workspace = Workspace01
    InputRasters01 = gp.ListRasters()#This will list all the rasters in Workspace01

    gp.workspace = Workspace02#change the workspace to the next workspace you wish to list data for.
    InputRasters02 = gp.ListRasters()
    iterInputRasters02 = iter(InputRasters02)#create an iterater for rasters02 so it may be incremented.
  
    for InputRaster in InputRasters01:#This will loop through your first list.
        #Becasue the workspace was changed to create the second list the path to the first list is no longer valid. So the script does not fail
        #we need to set the full path name to items in InputRasters01.
        fullInputRaster = Workspace01 + os.sep + InputRaster

        InputRaster02 = iterInputRasters02.next()#increment        
        #Set the output table name
        nameRaster1 = gp.describe(fullInputRaster).baseName
        nameRaster2 = gp.describe(InputRaster02).baseName
        
        outRaster = OutWorkspace + os.sep + "mean" + InputRaster + ".tif"#create a tif with first raster name. Use a tif because the name could become too long for esri grid.
        InExpression = "MEAN(" + fullInputRaster + "," + InputRaster02 +")"#use the variable so the raster is updated for each SOMA call

        # Process: MapAlgebra
        gp.SingleOutputMapAlgebra_sa(InExpression, outRaster)
except:
    gp.AddMessage(gp.GetMessages(2))
    print gp.GetMessages(2)

ArvindBhuta · ‎09-27-2011

Hey Tony,

I really appreciate it. I'll check it out and run it through. Looking forward to trying it. Thanks a ton again!

Regards,

Arvind