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 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.

Find more posts tagged with

Accepted answers

All comments

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!

rzufelt

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

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

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_

rzufelt

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

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

SAC_PhySAC_Phy

It doesn't report anything.
I'm not sure how I see command line messages (I'm very new to this.) The only messages I can see are the ones I code in by arcpy.AddMessage("Message")

So the first option didn't give me any message (though I tried it with the arcpy.AddMessage() as well) and I tried elevation.save()

I have since changed it to elevation.save("elevation") because I figured it needed the name to save the raster to since the workspace was already set and knew where to save.

I assume I have writing privileges because before this I have a terrain to raster method, and the terrains are being saved to the same location as I'm trying to save this raster (FD_addr), but this is before I set the workspace. I only set the workspace because I wanted a list created with the rasters that are saved in the FD_addr location. So unless setting the workspace alters the writing privilege, I assume I can write files to that location.

I'm still at a loss. I feel like I've tried so many options and still don't know what's wrong or how to fix this. It seems as if it should work!

Interesting. Are you sure you hard-coded it to an "existing" FGDB? If so, what do you mean it didn't work? Didn't create the raster with your hard-coded name, created a blank raster? Error messages?

What does it report if make the modifications in red?
import arcpy, os
arcpy.CheckOutExtension('Spatial')
arcpy.AddMessage("Creating Elevation Raster")
arcpy.env.workspace = FD_addr
list = arcpy.ListRasters("*", "ALL")
raster0 = arcpy.Raster(list[0])
raster1 = arcpy.Raster(list[1])
elevation = arcpy.sa.Con(raster1 - raster0 >= 0, raster1 - rasters0, 0)
print "temporary: ",elevation.isTemporary
elevation.save(FD_addr)
print "temporary: ",elevation.isTemporary
R_

Also, it looks like elevation.save(FD_addr) is telling it to save with the filename of the workspace. Workspace is normally a folder or FGDB, not a file.
It's my understanding that elevation.save() would save it in the current workspace as raster dataset "elevation".

Now that I think about it, before trying the above code, what happens if you try this (of course, make sure the FD_addr is set to a FGDB that the current user has write access to):
import arcpy, os
arcpy.CheckOutExtension('Spatial')
arcpy.AddMessage("Creating Elevation Raster")
arcpy.env.workspace = FD_addr
list = arcpy.ListRasters("*", "ALL")
raster0 = arcpy.Raster(list[0])
raster1 = arcpy.Raster(list[1])
elevation = arcpy.sa.Con(raster1 - raster0 >= 0, raster1 - rasters0, 0)
elevation.save()

rzufelt

Still didn't work 😕

I'm setting the workspace here:

env.workspace = FD_addr
list = arcpy.ListRasters("*", "ALL")

(not sure how you make it so that code appears in gray?)

I feel that has to be working because the rasters being inserted into list are the rasters in that location (I double checked.)
FD_addr is set way in the beginning of the code.

Interesting. Are you sure you hard-coded it to an "existing" FGDB? If so, what do you mean it didn't work? Didn't create the raster with your hard-coded name, created a blank raster? Error messages?

What does it report if make the modifications in red?

import arcpy, os
arcpy.CheckOutExtension('Spatial')
arcpy.AddMessage("Creating Elevation Raster")
arcpy.env.workspace = FD_addr
list = arcpy.ListRasters("*", "ALL")
raster0 = arcpy.Raster(list[0])
raster1 = arcpy.Raster(list[1])
elevation = arcpy.sa.Con(raster1 - raster0 >= 0, raster1 - rasters0, 0)
print "temporary: ",elevation.isTemporary
elevation.save(FD_addr)
print "temporary: ",elevation.isTemporary

R_

Also, it looks like elevation.save(FD_addr) is telling it to save with the filename of the workspace. Workspace is normally a folder or FGDB, not a file.
It's my understanding that elevation.save() would save it in the current workspace as raster dataset "elevation".

Now that I think about it, before trying the above code, what happens if you try this (of course, make sure the FD_addr is set to a FGDB that the current user has write access to):

import arcpy, os
arcpy.CheckOutExtension('Spatial')
arcpy.AddMessage("Creating Elevation Raster")
arcpy.env.workspace = FD_addr
list = arcpy.ListRasters("*", "ALL")
raster0 = arcpy.Raster(list[0])
raster1 = arcpy.Raster(list[1])
elevation = arcpy.sa.Con(raster1 - raster0 >= 0, raster1 - rasters0, 0)
elevation.save()

MattSayler

(not sure how you make it so that code appears in gray?)

http://forums.arcgis.com/threads/48475-Please-read-How-to-post-Python-code 😉

SAC_PhySAC_Phy

I don't see in the code anywhere that you are setting the variable 'FD_addr' (maybe it is above the imports somewhere?) to actually set the env.workspace.

Could try somehting like:
elevation.save("C:/output/file_gdb.gdb/elevation")      # or whatever the path/filename you want for your out raster
   
If this works, at least then you know it's an issue with the workspace/path being set properly.

R_

Still didn't work 😕

I'm setting the workspace here:

env.workspace = FD_addr
list = arcpy.ListRasters("*", "ALL")

(not sure how you make it so that code appears in gray?)

I feel that has to be working because the rasters being inserted into list are the rasters in that location (I double checked.)
FD_addr is set way in the beginning of the code.

rzufelt

I don't see in the code anywhere that you are setting the variable 'FD_addr' (maybe it is above the imports somewhere?) to actually set the env.workspace.

Could try somehting like:

elevation.save("C:/output/file_gdb.gdb/elevation")      # or whatever the path/filename you want for your out raster

If this works, at least then you know it's an issue with the workspace/path being set properly.

R_

SAC_PhySAC_Phy

Yep - on the last line try:
elevation.save("my_output")
The raster "my_output" should then be saved to the 'FD_addr' workspace.

Question though: Does 'FD_addr' point to a feature dataset within a geodatabase? Last I checked, I thought that you can't write rasters to a feature dataset (has to be directly to the GDB).

Tried both
elevation.save("my_output")
my_output.save(FD_addr)
"my_output".save(FD_addr)

But again, neither worked/saved.

FD_addr is really just bad naming on my behalf. This is part of something much larger, and when I first began this, I had intended for FD_addr to be the location inside the (file) geodatabase where the feature dataset would be created..

So FD_addr points to within the geodatabase and not within the feature dataset.
Bad naming, I know...
But I have other rasters being created with FD_addr as their output address and they do indeed go to the correct spot.

ChrisSnyder

Yep - on the last line try:

elevation.save("my_output")

The raster "my_output" should then be saved to the 'FD_addr' workspace.

Question though: Does 'FD_addr' point to a feature dataset within a geodatabase? Last I checked, I thought that you can't write rasters to a feature dataset (has to be directly to the GDB).

SAC_PhySAC_Phy

The inputs need to be "raster objects" to work in teh new map algebra. Try this:

import arcpy, os
arcpy.CheckOutExtension('Spatial')
arcpy.AddMessage("Creating Elevation Raster")
arcpy.env.workspace = FD_addr
list = arcpy.ListRasters("*", "ALL")
raster0 = arcpy.Raster(list[0])
raster1 = arcpy.Raster(list[1])
elevation = arcpy.sa.Con(raster1 - raster0 >= 0, raster1 - rasters0, 0)
elevation.save(FD_addr)

While this went through without giving any error message, no raster was created (or saved?).

ChrisSnyder

The inputs need to be "raster objects" to work in teh new map algebra. Try this:

import arcpy, os
arcpy.CheckOutExtension('Spatial')
arcpy.AddMessage("Creating Elevation Raster")
arcpy.env.workspace = FD_addr
list = arcpy.ListRasters("*", "ALL")
raster0 = arcpy.Raster(list[0])
raster1 = arcpy.Raster(list[1])
elevation = arcpy.sa.Con(raster1 - raster0 >= 0, raster1 - rasters0, 0)
elevation.save(FD_addr)

Join the Conversation

Signed in members can post, follow updates, and more. New here? Register a free account.

Quick Links

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