Select to view content in your preferred language

Map Algebra statement to interpolate between two rasters

792
6
04-15-2014 03:54 PM
BenRufenacht
Occasional Contributor
I am having some trouble converting a script I wrote for 9.3 to 10.0 using Map Algebra.  It is a linear interpolation between two rasters on a log-linear scale.

Here is the 9.3 code:

gp.SingleOutputMapAlgebra_sa("POW(10 , LOG10(4) - %s * (LOG10(4) - LOG10(10)) / (%s - %s))" %(D25, D25, D10), X25, "'';''")


and here is what I have for 10, which does not work:

X25 = Power(((Log10(4)) - Raster("D25") * (Log10(4) - Log10(10)) / (Raster("D25") - Raster("D10")),10)
X25.save(ws + "\X25")


Does anyone have any suggestions on what I may be doing wrong?
Tags (2)
0 Kudos
6 Replies
NeilAyres
MVP Alum
What doesn't work, what's the error message?
You may want to try a double "\" when you save the output.
0 Kudos
JamesCrandall
MVP Frequent Contributor
I am unsure if your Power() method is properly formed. 

http://resources.arcgis.com/en/help/main/10.1/index.html#//009z00000097000000

It takes two parameters but it doesn't look like you are supplying these.


Power (in_raster_or_constant1, in_raster_or_constant2)



Also, I noticed that there is an extra "(" paren char:


X25 = Power(((Log10(4)) - Raster("D25") * (Log10(4) - Log10(10)) / (Raster("D25") - Raster("D10")),10)

0 Kudos
BenRufenacht
Occasional Contributor
What doesn't work, what's the error message?
You may want to try a double "\" when you save the output.


This is the very descriptive error I am getting:

<type 'exceptions.RuntimeError'>: UNKNOWN


I tried the double "\\" and it did not make a difference.

I should note that I simplified my equation just to test the map algebra and I was able to get a simple raster addition to work and save.
0 Kudos
BenRufenacht
Occasional Contributor
I am unsure if your Power() method is properly formed. 

http://resources.arcgis.com/en/help/main/10.1/index.html#//009z00000097000000

It takes two parameters but it doesn't look like you are supplying these.


Power (in_raster_or_constant1, in_raster_or_constant2)



Also, I noticed that there is an extra "(" paren char:



My first parameter is the map algebra statement and the second is the constant 10 to cancel out the log interpolation result.

Thanks for pointing out the extra "(".
0 Kudos
BenRufenacht
Occasional Contributor
Now I am thinking that my map algebra is fine.  I think my issues may be in how I am stetting the environmental variables.  Here is more of the script:

# Import system modules
import arcpy, sys
from arcpy.sa import *

# Check out the ArcGIS Spatial Analyst extension license
arcpy.CheckOutExtension("Spatial")

# Set workspace from input
ws = sys.argv[1]
db = sys.argv[2]
arcpy.env.workspace = db
arcpy.env.scratchWorkspace = db

# Set input rasters
D10 = sys.argv[3]
D25 = sys.argv[4]
D50 = sys.argv[5]
D100 = sys.argv[6]
D500 = sys.argv[7]

# Set environment variables from input rasters
arcpy.env.snapRaster = D10
arcpy.env.extent = D10
arcpy.env.outputCoordinateSystem = D10
cord = arcpy.env.outputCoordinateSystem
arcpy.env.cellSize = D10
cs = arcpy.env.cellSize
arcpy.env.overwriteOutput = True


# 10 and 25 Log-Linear Interpolation
#X25 = ws + "\Calc1.gdb\X25"
X25 = Power((Log10(4)) - Raster(D25) * ((Log10(4) - Log10(10)) / (Raster(D25) - Raster(D10))),10)
X25.save(ws + "\X25")


Anyone see anything wrong here?
0 Kudos
BenRufenacht
Occasional Contributor
I was able to get it working with the following code:

# Import system modules
import arcpy
from arcpy.sa import *

# Check out the ArcGIS Spatial Analyst extension license
arcpy.CheckOutExtension("Spatial")

# Set workspace from input
ws = arcpy.GetParameterAsText(0)
db = arcpy.GetParameterAsText(1)
arcpy.env.workspace = ws
arcpy.env.scratchWorkspace = db

# Set input rasters
D10 = arcpy.GetParameterAsText(2)
D25 = arcpy.GetParameterAsText(3)
D50 = arcpy.GetParameterAsText(4)
D100 = arcpy.GetParameterAsText(5)
D500 = arcpy.GetParameterAsText(6)

# Set environment variables from input rasters
arcpy.env.snapRaster = r"D10"
arcpy.env.extent = r"D10"
arcpy.env.outputCoordinateSystem = r"D10"
#cord = arcpy.env.outputCoordinateSystem
arcpy.env.cellSize = r"D10"
#cs = arcpy.env.cellSize
#arcpy.env.overwriteOutput = True


# 10 and 25 Log-Linear Interpolation Map Algebra
X25 = Power((Log10(4)) - Raster(D25) * ((Log10(4) - Log10(10)) / (Raster(D25) - Raster(D10))),10)
X25.save(ws + "\X25")
0 Kudos