1.6.3 Example: Performing map algebra on a raster

Here’s another simple script that finds all cells over 3500 meters in an elevation raster and makes a new raster that codes all those cells as 1. Remaining values in the new raster are coded as 0. This type of “map algebra” operation is common in site selection and other GIS scenarios.

Something you may not recognize below is the expression Raster(inRaster). This function just tells ArcGIS that it needs to treat your inRaster variable as a raster dataset so that you can perform map algebra on it. If you didn't do this, the script would treat inRaster as just a literal string of characters (the path) instead of a raster dataset.

# This script uses map algebra to find values in an
#  elevation raster greater than 3500 (meters).

import arcpy
from arcpy.sa import *

# Specify the input raster
inRaster = "C:/Data/Elevation/foxlake"
cutoffElevation = 3500 

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

# Make a map algebra expression and save the resulting raster
outRaster = Raster(inRaster) > cutoffElevation
outRaster.save("C:/Data/Elevation/foxlake_hi_10")
    
# Check in the Spatial Analyst extension now that you're done
arcpy.CheckInExtension("Spatial")

Begin by examining this script and trying to figure out as much as you can based on what you remember from the previous scripts you’ve seen.

The main points to remember on this script are:

• Unlike other tools/functions we've seen so far, map algebra functions and the Raster class aren't in the arcpy module, but in the arcpy.sa sub-module.  Thus, this script imports both arcpy and arcpy.sa.
• The script uses an alternative syntax for importing modules: from <module> import <class> or *.  When a module is imported with this syntax, any classes imported can be invoked without prefixing the parent module.  In other words, using from arcpy.sa import * allows us to refer to the Raster class like this:
  outRaster = Raster(inRaster) > cutoffElevation
  whereas using import arcpy.sa would require this syntax when referring to the Raster class:
  outRaster = arcpy.sa.Raster(inRaster) > cutoffElevation
• You may notice that PyScripter may flag the from statement and later the statement that invokes the Raster class with a warning icon.  Your script will still run, but the warning appears because it can be poor practice to import all classes/functions from a module.  It is often better to import only those items that are needed for the task at hand.  In this case, changing * to Raster would be the better practice and would eliminate the warning.  We've used the * notation to be consistent with Esri's code samples, but also encourage you to limit what you import from modules when possible.
• Notice the lines of code that check out the Spatial Analyst extension before doing any map algebra and check it back in after finishing.  Checking the extension out before using something from the arcpy.sa module was required in ArcMap.  However, in ArcGIS Pro, the checkout step is required only if you're operating in an environment where extension licenses are shared (i.e., using a Concurrent Use license).  If you were to run this script on an educational copy of the software, the checkout and checkin statements could be omitted.  We've left them in because they don't hurt, and they may be needed in your work environment.  
• Because each line of code takes some time to run, avoid putting unnecessary code between checkout and checkin. This allows others in your organization to use the extension if licenses are limited. The extension automatically gets checked back in when your script ends, thus some of the Esri code examples you will see do not check it in. However, it is a good practice to explicitly check it in, just in case you have some long code that needs to execute afterward, or in case your script crashes and against your intentions "hangs onto" the license.
• inRaster begins as a string, but is then used to create a Raster object once you run Raster(inRaster). A Raster object is a special object used for working with raster datasets in ArcGIS. It's not available in just any Python script: you can use it only if you import the arcpy module at the top of your script.
• cutoffElevation is a number variable that you declare early in your script and then use later on when you build the map algebra expression for your outRaster.
• Your expression outRaster = Raster(inRaster) > cutoffElevation is saying, in plain terms, "Make a new raster and call it outRaster. Do this by taking all the cells of the raster dataset at the path of inRaster that are greater than the number I assigned to the variable cutoffElevation."
• outRaster is also a Raster object, but you have to call the method outRaster.save() in order to make it permanent on disk. The save() method takes one argument, which is the path to which you want to save.

Now try to run the script yourself using the FoxLake digital elevation model (DEM) in your Lesson 1 data folder. If it doesn’t work the first time, verify that:

• you have supplied the correct input and output paths;
• your path name contains forward slashes (/) or double backslashes (\\), not single backslashes (\);
• you have the Spatial Analyst extension installed and enabled. To check this, from ArcPro, click on the Project tab > Licensing and ensure Spatial Analyst is enabled;
• you do not have any of the datasets open in ArcGIS products;
• the output data does not exist yet. If you want to be able to overwrite the output, you need to add the line arcpy.env.overwriteOutput = True. This line can be placed immediately after import arcpy.

You can experiment with this script using different values in the map algebra expression (try 3000 for example).