topic Re: Map Algebra statement to interpolate between two rasters in Python Questions

Map Algebra statement to interpolate between two rasters

BenRufenacht — Tue, 15 Apr 2014 22:54:54 GMT

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?

Re: Map Algebra statement to interpolate between two rasters

NeilAyres — Wed, 16 Apr 2014 03:58:41 GMT

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

Re: Map Algebra statement to interpolate between two rasters

JamesCrandall — Sat, 11 Dec 2021 10:43:38 GMT

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)

Re: Map Algebra statement to interpolate between two rasters

BenRufenacht — Wed, 16 Apr 2014 16:16:54 GMT

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.

Re: Map Algebra statement to interpolate between two rasters

BenRufenacht — Sat, 11 Dec 2021 10:43:41 GMT

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

Re: Map Algebra statement to interpolate between two rasters

BenRufenacht — Sat, 11 Dec 2021 10:43:44 GMT

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?

Re: Map Algebra statement to interpolate between two rasters

BenRufenacht — Sat, 11 Dec 2021 10:43:46 GMT

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