Raster Calculator for Python Script

SAC_PhySAC_Phy · ‎06-27-2013

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.

ChrisSnyder · ‎06-27-2013

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)

SAC_PhySAC_Phy · ‎06-27-2013

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 · ‎06-27-2013

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 · ‎06-27-2013

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.

RhettZufelt · ‎06-27-2013

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 · ‎06-27-2013

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.

MattSayler · ‎06-27-2013

(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 😉

RhettZufelt · ‎06-27-2013

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()

SAC_PhySAC_Phy · ‎07-01-2013

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()