GEOG 485:
GIS Programming and Software Development

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 import *

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

# Check out the Spatial Analyst extension

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

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:

  • Notice the lines of code that check out the Spatial Analyst extension before doing any map algebra and check it back in after finishing. 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 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, open ArcMap, click Customize > Extensions and ensure Spatial Analyst is checked;
  • you do not have any of the datasets open in ArcMap;
  • 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).


Read the sections of Chapter 5 that talk about environment variables and licenses (5.9 & 5.11) which we covered in this part of the lesson.