Raster Calculator/Map Algebra and NumPy?

[MatthewBrown1]

Hi,

Can anyone suggest an efficient method for replacing the above tools? I have three map algebra equations that expand to over a dozen Math tools (FYI I'm implementing the Universal Soil Loss Equation). I could put these together with a Python script, but I'm not looking forward to that task.

Would Numpy be able to accomplish something like this:

0.0012* Power(float(%Mean Annual Rainfall%),2) * "%soil erodibility factor raster%" * "%LS%" *"%landcover raster%"

Without creating spaghetti code?

[KevinHibma · Accepted solution]

Hey Matthew,

You should be able to put together a script which does the same functionality using some of the math tools or SA math/raster functions to achieve this result. After talking with a colleague on the Spatial Analyst team, she put together a few lines of Python that might get you started on this.

import arcpy from arcpy.sa import *  #Checkout SA extension arcpy.CheckOutExtension("Spatial")  # May want to set environments Cellsize, extent, snap, and workspaces  #Cast datasets as Rasters MAR = Raster("C:/data/meanannurain") #this doesn't actually need to be cast b/c it's used in a GP tool not with operators SEFR = Raster("C:/data/soilerodfact") LS = Raster("C:/data/ls") LCR = Raster("C:/data/landcover")  #Universal Soil Loss Equation Result = 0.0012* (Float("MeanAnnualRain") ** 2) * SEFR * LS * LCR  #Save the temperary result Result.save("C:/output/result")

For reference, heres the operators you can use in SA on once you've converted to a raster object (as shown in the script):
http://resources.arcgis.com/en/help/main/10.1/index.html#/An_overview_of_the_Map_Algebra_Operators/0...
Or if you'd perfer, you can continue using the tools which are supported in the Runtime. The ones you're after live in this toolset here: http://resources.arcgis.com/en/help/main/10.1/index.html#/An_overview_of_the_Math_toolset/009z000000...

Hopefully this moves you forward.

[MatthewBrown1]

Hi Kevin,

That looks a lot easier than I thought it would be!

ArcGIS + Python = 😄

Thanks!

I had to use a couple of intermediate steps for operations like Sin. This didn't work:

Z = (math.sin(slope * (3.14 / 180)) / 0.0896) ** 1.3

(I guess that's where arcpy.RasterToNumPyArray would be useful?)

But this did:

sinSlope = Sin((slope * (3.14 / 180))) 
Z = (sinSlope / 0.0896) ** 1.3