Week 4: Raster Analysis & Map Algebra

Back To Week 3 and Vector Overlay | On to Week 5 and Raster <-> Vector Conversion


Many of the functions used for vector analysis are either inefficient or ineffective in dealing with certain spatial questions, processes, or data types. Cell-based processing with Grid introduces a different approach to performing spatial analysis. The functionalities included in Grid are many, so again, this will simply be a cursory treatment of some of the basic things Grid can do.

Before you start, you may want to import the Pack Forest DEM grid. Once you have downloaded the file, unzip it:

trees% gzip -d dem.dem.gz

Then create a lattice (grid) called "dem" from the USGS-format DEM file:

Arc: demlattice dem.dem dem # # int

Map Algebra, Grid Functions, and Grid Commands

Most of the functionality in grid is used in the form of map algebra. Map algebra is a method of treating individual raster, or grid layers, as members of algebraic equations. The basic form of a map algebra equation is

output  assignment-operator  expression

For example, the output is a grid representing slope, the assignment-operator is "=," and the expression is the slope function applied to an elevation grid:

slp = slope (dem, percentrise)

In this example, the new surface grid (dem_new) is the sum of the existing surface grid (dem) and a grid describing where and how deep the proposed landfill is (landfill):

dem_new = dem + landfill

Some functions do not use mathematical expressions at all; in the next example, we create a new grid with zones representing the mean elevation per stand:

mean_elev = zonalmean (stands, dem)


Most of what you do in Grid takes the form of these algebraic functions. The usage for a function always looks like

Usage:  (I) FLOWDIRECTION (<surface_grid>,  {o_drop_grid}, {NORMAL | FORCE})
Usage:  (F) SLOPE (<grid>, {DEGREE | PERCENTRISE})
        (F) SLOPE (<grid>, <z_factor> {DEGREE | PERCENTRISE})
Usage:  (*) ABS (<grid | scalar | number>)

Where the (I) means that the output will be an integer grid or scalar, (F) means that the output will be floating-point (decimal), and (*) means that the output will be integer if all inputs are integer; however, if one or more inputs if floating-point, the output will be floating-point. This is an important distinction, since VAT (value attribute tables) are created only for integer grids.

Some of what you will do in Grid takes the form of commands, such as

Usage:  FILL <in_grid> <out_grid> {SINK | PEAK} {z_limit} {out_dir_grid}
Usage:  SETCELL <MAXOF | MINOF | cellsize | grid>
Usage:  SETMASK <OFF | grid>

Note that these commands do not take the form of algebraic equations, since they do not have an equal sign.

Functions follow standard algebraic rules for precedence of calculation, and functions may be nested within each other, and are worked from the inside out.

Basic Function Types

Grid functions can be classified according to how they deal with spatial extents.

Local functions create outputs in which the output cell values are determined on a cell-by-cell basis without regard for the value of neighboring cells. For example, the calculation of the landfill grid in the above example is a local function. In the output grid, each cell takes a value determined solely by the values of the corresponding cells in the 2 input grids. Although in reality, the values of the surface elevations are affected by their neighboring locations, in the Grid cell processing, the value of neighboring cells has no effect.

Focal functions create outputs in which the value of the output grid is affected by the value of neighboring cells. Filters are commonly used to smooth out data. A common filter, the "low pass" filter, passes a 3 x 3 cell window over the input grid. The values of the 9 cells in this "focus" are averaged, and the resultant value is placed in the central cell of the 3 x 3 focus in the output grid.

Zonal functions create outputs in which the value of output cells are determined in part by the spatial association between cells in the input grids. In the "zonalmean" command example above, the zone grid is timber stands. Every cell with the same value in the timber stands grid is part of a single zone. The mean value per stand is generated from these zones. Once the mean value per zone is determined, this value is placed on the cell in the output grid.

Conditional Cell Processing

Many analytical tasks in Grid rely on Boolean and conditional statements. At the simplest, use the CON function:

CON (<condition>, <true_expression>, 
    {<condition>, <true_expression>},....
    {<condition>, <true_expression>}, 

For example, to create the slp_4 coverage from last week's exercises, I used the following statement to create 4 zones of different slope interval. The CON function is nested within the GRIDPOLY function. First, the slope classes are created and written to a temporary grid. the temporary grid is then converted to a polygon coverage. The temporary grid is deleted and the polygon coverage persists.

Grid: slp_4 = gridpoly (con (slp lt 10, 1, ~
                        slp ge 10 and slp lt 25, 2, ~
                        slp ge 25 and slp lt 40, 3, ~

The conditional statement sets up 4 output classes (1, 2, 3, and 4), based on the slope values (<10, 10 - 25, 25 - 40, and >= 40).

In addition to the CON function, there is the IF statement and the DOCELL block. If you are interested in conditional cell processing, please see the on-line documentation for these.

Displaying Grids

There are several simple commands for displaying grids. As with coverages, you must set your map extent before you can display a grid:

Grid: &stat 9999

Grid: mapextent slp

Once your map extent is set, display the grid with the GRIDPAINT or GRIDSHADES commands. GRIDSHADES uses maps cell values to the shade numbers in the current shadeset. GRIDPAINT uses either a pre-defined set of shades or a colormap file. Linear or equal-area stretches can be applied in either of these commands.

Usage: GRIDPAINT <grid> {item}
                 {IDENTITY | LINEAR | EQUALAREA | remap_table}
                 {WRAP | NOWRAP} {NOMINAL | GRAY | colormap_file}
Usage: GRIDSHADES <grid> {item}
                  {IDENTITY | LINEAR | EQUALAREA | remap_table}
                  {WRAP | NOWRAP}

Grid: shadeset colorrange
Grid: gridshades dem # linear nowrap
Grid: gnds white
Grid: gridpaint dem # linear nowrap gray

Querying Grids

There are 2 basic ways to query grids: visually, and logically.

To answer the question "What value does a grid have at a given location?" use the CELLVALUE command:

Usage: CELLVALUE <grid> <xy | *> {item...item | NONE}

If you use the * option, use the pointer to click on the cell of interest. If you know the x, y coordinate of the location if interest, you can enter these instead. You also have the option of listing 0 or more items from the grid's VAT.

To answer the question "Where are cells of a given value?" use the GRIDQUERY command:

Usage: GRIDQUERY <grid> {item}
                 {IDENTITY | LINEAR | EQUALAREA | remap_table}
                 {WRAP | NOWRAP} <logical_expression>

GRIDQUERY will act on the logical expression you enter, and then display the selected cells using the currently active shadeset. You can perform stretches exactly the same way as you would with the GRIDPAINT and GRIDSHADES commands.

Grid: gridq dem # # # value gt 1000


It will be impossible to cover too many examples, but here are a few representative samples of some Grid processes.

Grid: setcell dem 
Grid: age_grid = polygrid (tty4, age_1998)
0       100   :   1
100     200   :   2
200     300   :   3
300     400   :   4
400     500   :   5
500     600   :   6
600     700   :   7
700     800   :   8
800     900   :   9
900    1000   :   10
1000   1100   :   11
1100   1200   :   12
1200   1300   :   13
1300   1400   :   14
1400   1500   :   15
1500   1600   :   16
1600   1700   :   17
1700   1800   :   18
1800   1900   :   19
1900   2000   :   20
2000   2100   :   21

Grid: dem_100 = reclass (dem, elev_100.txt)
Grid: age_elev = combine (age_grid, dem_100)
Grid: strgrd = linegrid (str2, str2dnrty, #, #, 30)
Grid: str_100 = con (eucdistance (strgrd) lt 100, 1)
Grid: shade = hillshade (dem, #, #, shade)
Grid: shadedelete all
Grid: shadetype color
Grid: shadecolorramp 1 50 green yellow
Grid: shadecolorramp 51 50 yellow red
Grid: $$display = slice (dem, eqinterval, 100)


Back To Week 3 and Vector Overlay | On to Week 5 and Raster <-> Vector Conversion