Raster Calculator for Python Script

[SAC_PhySAC_Phy]

I cannot get a conditional evaluation calculation to be done on two rasters.

I want an ouput raster to be created from the positive difference of two rasters (so no negative numbers)

The part of the code that isn't working

import arcpy
import os
from arcpy import env
from arcpy.sa import *
import arcgisscripting
arcpy.CheckOutExtension('Spatial')

I've tried two different ways, both failing

Method one

        arcpy.AddMessage("Creating Elevation Raster")
env.workspace = FD_addr
list = arcpy.ListRasters("*", "ALL")
elevation = Con(list[1]-list[0]>=0, list[1]-list[0],0)
elevation.save(FD_addr)

Method two
arcpy.AddMessage("Creating Elevation Raster")
expression = "Con({}-{}>=0, {}-{},0)".format(list[1], list[0], list[1], list[0])
OR = "elevation3"
OutRaster = os.join.path(FD_addr,OR)
gp = arcgisscripting.create()
gp.SingleOutputMapAlgebra_sa(expression, OutRaster)



my list has two rasters (I've checked to make sure they are there): a DSM and a DEM, and I want to find the positive elevation difference between the two.

Any help would be appreciated! I've been stuck on this for ages.

[SAC_PhySAC_Phy]

So this is what the part of the code that isn't creating the raster looks like (libraries and extensions included)

import arcpy
import os
from arcpy import env
from arcpy.sa import *
arcpy.CheckOutExtension("Spatial") # for arcpy.sa

 env.workspace = FD_addr  # DEM/DSM raster locations
 raster_list = arcpy.ListRasters("*", "ALL") # adds DEM/DSM raster
 size = len(raster_list) #### Make some exception here if != 2
 raster0 = arcpy.Raster(raster_list[0])
 raster1 = arcpy.Raster(raster_list[1])
 elevation = arcpy.sa.Con(raster1-raster0 >= 0, raster1-raster0,0)
 elevation.save("elevation") # saved in FD_addr
 

[RhettZufelt]

Not really sure how to see Arcmessages as I do everything in the command line, and not really familiar with the Arc command window/python tools.

How are you running this?  In the command window, or making a script tool or what?

If running in console window, after it runs, if you type:

print elevation
print raster0
print raster1

do you get any results?

also, can you create the elevation raster that you are after within arcmap using the tools?  I am wondering if it is creating an empty output, so is not saving it.

I would make sure that it is possible to do what you are trying using arcmap tools/model and get the desired output.  then work on the python code (export from the model builder and modify).

Might also look into the DEM format.  I have tried to use this in the past for certain tasks, and found that it does not support state plane coordinates.  Perhaps one of the input rasters is not supported or valid for this operation?  If you can do the raster calculations in ArcMap, then that would rule this out as well.

R_

[SAC_PhySAC_Phy]

Yes, I can create the elevation raster from the two rasters produced by the code by using the Raster Calculator.
I tried copying the code as a python snippet, but I read on the documentation for the Raster Calculator that it is only supported inside arcmap and will not work in a python script. The equivalent is supposedly how I'm trying to implement it.

In Raster Calculator, I type in the following:

elevation = Con("Surface_RAS"-"Ground_RAS">=0,"Surface_RAS"-"Ground_RAS",0)

The python snippet is:

arcpy.gp.RasterCalculator_sa("""elevation = Con("Surface_RAS"-"Ground_RAS">=0,"Surface_RAS"-"Ground_RAS",0)""","C:/Summer_13/Geo.gdb/elevation")

(Surface_RAS and Ground_Ras are the two raster inputs that the code creates.)
Works perfectly.

I'm creating a new Python toolbox, and running the script as a tool.

DEM format? I'm not sure I know what you're referencing.

Not really sure how to see Arcmessages as I do everything in the command line, and not really familiar with the Arc command window/python tools.

How are you running this?  In the command window, or making a script tool or what?

If running in console window, after it runs, if you type:

print elevation
print raster0
print raster1

do you get any results?

also, can you create the elevation raster that you are after within arcmap using the tools?  I am wondering if it is creating an empty output, so is not saving it.

I would make sure that it is possible to do what you are trying using arcmap tools/model and get the desired output.  then work on the python code (export from the model builder and modify).

Might also look into the DEM format.  I have tried to use this in the past for certain tasks, and found that it does not support state plane coordinates.  Perhaps one of the input rasters is not supported or valid for this operation?  If you can do the raster calculations in ArcMap, then that would rule this out as well.

R_

[SAC_PhySAC_Phy]

This works:

arcpy.gp.RasterCalculator_sa("""elevation = Con("Surface_RAS"-"Ground_RAS">=0,"Surface_RAS"-"Ground_RAS",0)""","C:/Summer_13/Geo.gdb/elevation")

But trying to edit it so that it takes the variables doesn't work.

DEM = arcpy.Raster(raster_list[0])
 DSM = arcpy.Raster(raster_list[1])
 expression = """elevation = Con({}-{}>=0,{}-{},0)""".format(DSM,DEM,DSM,DEM)
 raster_out = os.path.join(FD_addr,"elevation")
 arcpy.gp.RasterCalculator_sa(expression,raster_out)

Probably something is wrong with the expression variable.
I have to look into figuring out how to convert the """....""" to make my expression

[RhettZufelt]

env.workspace = FD_addr  # DEM/DSM raster locations
raster_list = arcpy.ListRasters("*", "ALL") # adds DEM/DSM raster

is where I got the information that you were using a DEM (Digital Elevation Model by USGS).

What I meant by using ArcMap tools, was not actually using the raster calculator, but the spatial analys "Con" tool.

What I would do is open ArcMap, open the catalog window, right-click on Toolbox (or My Toolboxes) folder and choose "New" Model.

Then, drag the sa.Con tool onto the model (if you can't find it, open the search window in ArcMap, select Tools, and type in Con, hit enter)

double click the Con icon and fill in the parameters, click save and run the model.

If you get the results you want, open the model, choose file, export python script.  This will give you a starting point, might have to clean up the syntax, etc. (i.e. Arc accepts "-" in path, python does not) and should give you an idea of how to code the "tool".

R_

Noticed this uses the gp.Con_sa tool not the sa.Con tool.  however, it is 10.1 that dumps this code out, and don't see ANY documentation for it.

[SAC_PhySAC_Phy]

What I meant by using ArcMap tools, was not actually using the raster calculator, but the spatial analys "Con" tool.

What I would do is open ArcMap, open the catalog window, right-click on Toolbox (or My Toolboxes) folder and choose "New" Model.

Perfect, this worked!!!
I separated the process into two different models. I performed a Minus first, and used that as the input raster which had the conditional evaluated on.

Thanks so much for your help and patience!