Mapping Center

In a blog entry posted last week called Filling and clipping a raster, Aileen described how to fill in holes in a "bad DEM" using data from an existing "good DEM", then clip the filled DEM to the outline of a feature.

The blog post suggested using some ArcGIS geoprocessing tools that are available with the Spatial Analyst extension. As with most GIS operations, there is more than one way to get to the final answer! In this blog post, I describe how Map Algebra can be used to achieve the same results.

We can break this process down into three steps:

• create a "clip raster" from an existing polygon feature,
• fill in the missing cell values in the bad DEM with values in the good DEM to create a filled DEM, and 
• clip the filled DEM to the extent of the clip raster.

To demonstate, I used the following datasets:

• WashingtonState – a polygon feature of the state of Washington,
• badDEM – the DEM that needs to be filled and clipped, and
• goodDEM – the DEM from the ESRI Data and Maps DVD (world_elevation\etopo2\etopo2 GRID) that I used to fill the badDEM.

You can accomplish the entire process using three Map Algebra statements:

So what exactly are these statements doing?

-- Shapegrid is the map algebra equivalent of Feature to Raster, and defines:

• the feature you want to convert (in this case, WashingtonState),
• the name of the column in the <shapefile>.dbf file used to assign grid cell values (I’ve used # to set this to no field), and
• the cell size of the output grid (I’ve set this to 50) – note that the units are the same as for your input features (in this case, feet).

The output is a raster that closely resembles the Washington State boundary (how closely it resembles the polygon is a function of the cell size of the raster). All cells have a value of either 0 or the value from the attribute you selected to assign grid cell values.

-- here I’ve combined two steps into one using the Con and Is Null statements:

1. The Con statement sets up a conditional test: if a cell in the raster passes the test, it is given a certain value (in this case, the value of the cell in the good DEM). If it fails the test, it is given a different value (in this case, the value of the cell in the bad DEM.)
2. The conditional statement uses Is Null to determine whether the cell being processed in the bad DEM is NoData or not. This statement returns a value of 1 if the input value is NoData and 0 for cells that are not.

The conditional statement is then: set up a condition that tests badDEM for NoData cells; if the cell has a NoData value, give it the value of the cell in the goodDEM; otherwise give it the original value in the badDEM. The output is the badDEM with new values for NoData cells from the goodDEM.

-- Selectmask is the Map Algebra equivalent of the Extract by Mask Spatial Analyst geoprocessing tool.  It asks: what is the raster to extract from (DEMfill), and what is the mask to use as the extent (land_grid)?  The output is a raster that contains elevation values only for the cells that fall completely within the clip raster. All other cells will have NoData values.

Note that another way to clip the filled DEM to the WashingtonState clip raster would be to use another Con statement:

-- this yields the same result, but you have to know the grid cell values that the Shapegrid conversion will produce, in this example 0. The conditional statement is: set up a condition that tests if the cells in land_grid equal 0 (note that you need to use "==" or EQ symbol for the "equal to" condition); if they do, give them the value from DEMfill; otherwise, keep the original value (which in the case of land_grid is NoData because anything outside of land_grid has no elevation value).

So your next question might be, "How you can create a ModelBuilder model using this Math Algebra approach?" You can use the Single Output Map Algebra (SOMA) and Multi Output Map Algebra (MOMA) geoprocessing tools. Below is what the ModelBuilder model looked like when I used the three SOMA expressions to fill and clip a raster. Each line of text in a yellow box in the model is a SOMA expression.

This is what the SOMA tool interface looks like with a Map Algebra statement filled in:

Note that the Input raster data shown in ModelBuilder includes DEMfill and land_grid. These are shown as inputs to the SOMA expression in the ModelBuilder interface. The output raster is also shown in ModelBuilder.

You can download this .zip file which contains the the model described above so you can try out the methods for yourself.

Perhaps the best resource to learn more about Map Algebra is the book by Dr. C. Dana Tomlin called Geographic Information Systems and Cartographic Modeling. This book, published in 1990, was based on the PhD dissertation that Dr. Tomlin wrote in 1983. In his book, Dr. Tomlin introduced the the concepts of a simple and an elegant set-based algebraic approach to manipulating geographic data that has since become the basis of much of the work in raster analysis. Read more about Dr. Tomlin and his book here: http://www.design.upenn.edu/new/larp/facultybio.php?fid=121.