r.mapcalc

r.mapcalc is GRASS's powerful map algebra utility. The general format of a mapcalc statement is newmap = f(map1,map2,...). It can be run interactively by typing r.mapcalc and then entering statements at the mapcalc> prompt, or noninteractively with statements entered with the mapcalc command or redirected from a command file. mapcalc works within the current region window and resolution, and current mask: masked cells are assigned zeroes, and the cellhd file for any layer created in mapcalc is defined from the current region settings.

Review the available operands and their precedences, and the available functions below. Note that these follow standard C programming conventions, e.g. a double equals sign for logical equality. The expression

        conifer = if(vegcover==3,1,0)
creates the new binary raster map conifer from the existing raster map vegcover with coniferous forest cells equal to one and all other cells equal to zero.  The expression could be written more compactly as

        conifer = if(vegcover==3)

mapcalc statements can operate on multiple raster layers at once, and perform all kinds of complex mathematical or logical operations, e.g.:

        siteclass = if((slope<0.05 && elevation.dem<1400 && geology==5) \
        || (slope<0.1 && vegcover==3),100,if(streams,15,1))
Cells with either slope less than 0.05, elevation greater than 1400 and geology category 5, or with slope less than 0.1 and vegetation category 3, are assigned a value of 100; otherwise, cells in streams are assigned a value of 15; other cells are assigned a value of 1.

mapcalc handles floating-point calculations, but its final output is truncated by GRASS 4 to integer cell values.  GRASS 5 maintains floating-point values.

Note that while GRASS 4 maintained no distinction between zero and missing cell values, GRASS 5 supports this distinction.  (By default, missing data values are coded as *.)  When you do map algebra on missing-data cells, the output cells also have mising data.  To convert missing-data cells to zeroes, use the isnull() function:

        fixmap = 'if(isnull(origmap),0,origmap)'

Long mapcalc statements can be broken up into multiple lines terminated by \ characters (with no trailing blanks).

mapcalc's neighbor operators make it useful for creating all kinds of filters. The expression

        smooth.elev = (elev[-1,-1]+elev[-1,0]+elev[-1,1]+elev[0,-1] \
        +elev[0,0]+elev[0,1]+elev[1,-1]+elev[1,0]+elev[1,1])/9
assigns to each cell the average elevation in its 9-cell neighborhood (equivalent to an r.neighbors operation). mapcalc can also be used on satellite imagery for edge-detection, image constrast stretching with saturation, etc. Shapiro and Westervelt's tutorial (in the grass.docs directory) illustrates use of mapcalc for calculating slope and aspect maps from digital elevation models, performing watershed analyses, image processing, etc.

Be careful about using characters which have special meanings in the UNIX shell:

 * ( ) > & |
It is good practice to use single quotes around expressions containing such characters.

If raster cell categories reference numeric labels in a map's cats file, the @ symbol makes mapcalc operate on these label values directly rather than on the cell values:

        newmap = map1*@map2
Each cell in newmap is the product of the corresponding cell category in map1 times the corresponding cell label in map2.

Alternately, the @ symbol may be used to identify the mapset containing a map:

        newmap2 = conifer@john/(log(elevation.dem@PERMANENT))
The # operator converts map cell values to grayscale equivalents

        graymap = #origmap

The r#, g# and b# operators extract red, green and blue color separates from the map, by reference to its color table:

        ndvi = '(r#map-g#map)/r#map+g#map)'

Terminate interactive use of mapcalc with an exit command or ctrl-D.

Type long mapcalc expressions or sequences of expressions in a text file for redirection to mapcalc:

r.mapcalc < command.file
If an expression will overwrite an existing file, mapcalc in interactive mode (expressions entered at the mapcalc> prompt) will ask you to confirm the overwrite before proceeding. In non-interactive mode (expressions entered at the GRASS command line or redirected from a text file) mapcalc will overwrite existing files without warning.

mapcalc's error messages ("you have confused me") are pretty useless. You have to figure out syntax errors for yourself, and should verify that your mapcalc expressions actually execute the way you want.

By default, maps are titled with the mapcalc expression which created them. You can re-title and label any maps created with mapcalc using r.support.

mapcalc can quickly generate huge numbers of files. Try to give your maps descriptive names so you don't lose track of what they are. At the end of your GRASS session, please remove any unwanted maps you have created.

Operators (order of precedence in parentheses)

Relative offset:

x[m,n] returns the value of the cell in map x offset by m rows and n columns from the reference cell.

Arithmetic:

*       multiplication                  (4)
/       division                        (4)
%       modulus                         (4)
+       addition                        (3)
-       subtraction                     (3)
Comparison:
==      equal                           (2)
!=      not equal                       (2)
<       less than                       (2)
<=      less than or equal              (2)
>       greater than                    (2)
>=      greater than or equal           (2)
Logical:
&&      logical and (intersection)      (1)
||      logical or (union)
(1)

Mapcalc functions:

abs(x)          absolute value of x
exp(x)          e to the power x
exp(x,y)        x to the power y
float(x)        converts integer map x to floating-point
int(x)          converts floating-point map x in integer
log(x)          natural logarithm of x
log(x,y)        logarithm base y of x
sqrt(x)         square root of x
min(x,y,...)    minimum of x, y, ...
max(x,y,...)    maximum of x, y, ...
round(x)        rounds x to integer
if(x)           returns 1 if x!=0; 0 otherwise
if(x,y)         returns y if x!=0; 0 otherwise
if(x,y,z)       returns y if x!=0; z otherwise
if(x,y,z,a)     returns y if x>0; z if x==0; a otherwise
eval(a,b,...,x) returns x (evaluating all maps a, b, ...)

Examples:

upland.forest = (elevation.dem>=1200 && elevation.dem<=1450) \ 
&& (vegcover==3 || vegcover==4 || vegcover==5)
--forest cells within elevation range coded as 1, other cells as zero.
lowpass3x3 = (image[-1,-1]+image[-1,0]+image[-1,1]+image[0,-1] \
+image[0,0]+image[0,1]+image[1,-1]+image[1,0]+image[1,1])/9
--low-pass filter, equivalent to r.neighbors i=image o=lowpass3x3 size=3 method=average
smooth.elev = (elev[-1,-1]+2*elev[-1,0]+elev[-1,1]+2*elev[0,-1] \
+8*elev[0,0]+2*elev[0,1]+elev[1,-1]+2*elev[1,0]+elev[1,1])/20
--creates smoothed elevation map from weighted average of 9-cell neighborhood
grow.area = if(area,area,if(area[-1,0],area[-1,0],if(area[0,-1], \
area[0,-1],if(area[1,0],area[1,0],if(area[0,1],area[0,1])))))
--equivalent to r.grow i=area o=grow.area
tvi = 100 * sqrt(max(0,(tm.4 - tm.3)/(tm.4 + tm.3) + 0.5))
--creates transformed vegetation index from Landsat TM bands 3 (red) and 4 (near infra-red).



course syllabus