r.mapcalc(1grass)

๐Ÿ’– NAME ๐Ÿ’–

r.mapcalc โ€” the Raster map calculator. Yes bestie, it's r.mapcalc(1grass), straight outta the GRASS User's Manual, and it's about to do math on your raster maps like it's nothing. โœจ

๐Ÿท๏ธ KEYWORDS ๐Ÿท๏ธ

raster, algebra โ€” that's the vibe. ๐ŸŽ€


โœจ SYNOPSIS โœจ

r.mapcalc
r.mapcalc --help
r.mapcalc [-sl]  [expression=string]   [region=string]   [file=name]
[seed=integer]   [nprocs=integer]   [--overwrite]  [--help]  [--verbose]
[--quiet]  [--ui]

๐Ÿšฉ Flags ๐Ÿšฉ

  • -s โ€” Generate random seed (result is non-deterministic) [deprecated]. This flag is deprecated and will be removed in a future release. Seeding is automatic or use parameter seed.
โš ๏ธ Deprecated alert!! -s is on its way out the door bestie. Use seed instead, or just let it seed itself automatically. ๐Ÿ’…
  • -l โ€” List input and output maps ๐Ÿ‘€
  • --overwrite โ€” Allow output files to overwrite existing files
  • --help โ€” Print usage summary
  • --verbose โ€” Verbose module output
  • --quiet โ€” Quiet module output
  • --ui โ€” Force launching GUI dialog

๐ŸŽ€ Parameters ๐ŸŽ€

  • expression=string โ€” Expression to evaluate
  • region=string โ€” The computational region that should be used.
    • Options: current, intersect, union
    • Default: current
    • current โ€” current uses the current region of the mapset
    • intersect โ€” intersect computes the intersection region between all input maps and uses the smallest resolution
    • union โ€” union computes the union extent of all map regions and uses the smallest resolution
  • file=name โ€” File containing expression(s) to evaluate
  • seed=integer โ€” Seed value for the random number generator. Using the same seed ensures identical results, while a randomly generated seed produces different outcomes in each run.
  • nprocs=integer โ€” Number of threads for parallel computing
    • 0: use OpenMP default; >0: use nprocs; <0: use MAX-nprocs
    • Default: 0

๐ŸŒธ DESCRIPTION ๐ŸŒธ

r.mapcalc does arithmetic on raster map layers. It can create brand-new raster map layers that are arithmetic expressions built out of existing raster map layers, integer or floating point constants, and functions. It's giving calculator, but make it geography. ๐Ÿ—บ๏ธโœจ

๐Ÿงฎ Program use ๐Ÿงฎ

An r.mapcalc expression has this form:

result = expression

Where result is the name of a raster map layer that will hold the answer, and expression is any legal arithmetic expression involving existing raster map layers (except result itself, no recursive shenanigans ๐Ÿšซ), integer or floating point constants, and functions the calculator knows about. Parentheses are allowed in the expression and may be nested to any depth. And result gets created in the user's current mapset. ๐Ÿ’–

Since expression= is the first option, it is the default. That means you can pass an expression right on the command line as long as the expression is quoted and there's a space before the first = sign. Example ('foo' is the resulting map):

r.mapcalc "foo = 1"

or:

r.mapcalc โ€™foo = 1โ€™

An unquoted expression (i.e. split over multiple arguments) won't work, and neither will omitting the space before the = sign:

r.mapcalc โ€™foo=1โ€™
Sorry, <foo> is not a valid parameter

To read the command from a file, use file= explicitly, like:

r.mapcalc file=file

or:

r.mapcalc file=- < file

or:

r.mapcalc file=- <<EOF
foo = 1
EOF

The formula you feed to r.mapcalc gets recorded both in the result map title (which shows up in the category file for result) and in the history file for result. It remembers everything. ๐Ÿ™ˆ

Some characters have special meaning to the command shell. If you're entering input to r.mapcalc on the command line, expressions should be enclosed within single quotes. See NOTES, below. ๐Ÿ‘‘

๐ŸŒ Computational regions in r.mapcalc ๐ŸŒ

By default r.mapcalc uses the current region as the computational region โ€” the one you set with g.region โ€” for processing. But sometimes you need a region that's derived from the raster maps in the expression to set the computational region. This is highkey important for modules that use r.mapcalc internally to process time series of satellite images that all have different spatial extents. A module that requires this feature is t.rast.algebra. The region option of r.mapcalc was implemented to address this requirement. It lets you compute and use a region based on all raster maps in an expression. Three modes are supported:

  • Setting the region parameter to current will result in the use of the current region as computational region. This is the default. The current region can be set with g.region.
  • The parameter union will force r.mapcalc to compute the disjoint union of all regions from raster maps specified in the expression. This computed region will then be used as computational region at runtime. The region of the mapset will not be modified. The smallest spatial resolution of all raster maps will be used for processing.
  • The parameter intersect will force r.mapcalc to compute the intersection of all regions from raster maps specified in the expression. This computed region will then be used as computational region at runtime. The region of the mapset will not be modified. The smallest spatial resolution of all raster maps will be used for processing.

๐Ÿ”ฅ Operators and order of precedence ๐Ÿ”ฅ

The following operators are supported:

Operator   Meaning                    Type        Precedence
--------------------------------------------------------------
-          negation                   Arithmetic  12
~          oneโ€™s complement           Bitwise     12
!          not                        Logical     12
^          exponentiation             Arithmetic  11
%          modulus                    Arithmetic  10
/          division                   Arithmetic  10
*          multiplication             Arithmetic  10
+          addition                   Arithmetic   9
-          subtraction                Arithmetic   9
<<         left shift                 Bitwise      8
>>         right shift                Bitwise      8
>>>        right shift (unsigned)     Bitwise      8
>          greater than               Logical      7
>=         greater than or equal      Logical      7
<          less than                  Logical      7
<=         less than or equal         Logical      7
==         equal                      Logical      6
!=         not equal                  Logical      6
&          bitwise and                Bitwise      5
|          bitwise or                 Bitwise      4
&&         logical and                Logical      3
&&&        logical and[1]             Logical      3
||         logical or                 Logical      2
|||        logical or[1]              Logical      2
?:         conditional                Logical      1

(modulus is the remainder upon division) ๐Ÿ’…

[1] The &&& and ||| operators handle null values differently to other operators. See the section entitled NULL support below for more details.

The operators are applied from left to right, with those of higher precedence applied before those with lower precedence. Division by 0 and modulus by 0 are acceptable and give a NULL result. The logical operators give a 1 result if the comparison is true, 0 otherwise. Simple as that, periodt. โญ

๐Ÿ—‚๏ธ Raster map layer names ๐Ÿ—‚๏ธ

Anything in the expression which is not a number, operator, or function name is taken to be a raster map layer name. Examples:

elevation
x3
3d.his

Most GRASS raster map layers meet this naming convention. But if a raster map layer has a name that conflicts with the above rule, it should be quoted. For example, the expression

x = a-b

would be interpreted as: x equals a minus b, whereas

x = "a-b"

would be interpreted as: x equals the raster map layer named a-b. ๐Ÿฆ‹

Also

x = 3107

would create x filled with the number 3107, while

x = "3107"

would copy the raster map layer 3107 to the raster map layer x.

Quotes are not required unless the raster map layer names look like numbers or contain operators, OR unless the program is run non-interactively. Examples given here assume the program is run interactively. See NOTES, below. ๐Ÿ‘‘

r.mapcalc will look for the raster map layers according to the user's current mapset search path. It is possible to override the search path and specify the mapset from which to select the raster map layer. This is done by specifying the raster map layer name in the form:

name@mapset

For example, the following is a legal expression:

result = x@PERMANENT / y@SOILS

The mapset specified does not have to be in the mapset search path. (This method of overriding the mapset search path is common to all GRASS commands, not just r.mapcalc.) ๐Ÿ’œ

๐Ÿ” The neighborhood modifier ๐Ÿ”

Maps and images are data base files stored in raster format, i.e., two-dimensional matrices of integer values. In r.mapcalc, maps may be followed by a neighborhood modifier that specifies a relative offset from the current cell being evaluated. The format is map[r,c], where r is the row offset and c is the column offset. For example, map[1,2] refers to the cell one row below and two columns to the right of the current cell, map[-2,-1] refers to the cell two rows above and one column to the left of the current cell, and map[0,1] refers to the cell one column to the right of the current cell. This syntax permits the development of neighborhood-type filters within a single map or across multiple maps. ๐Ÿ”ฅ

โš ๏ธ The neighborhood modifier cannot be used on maps generated within the same r.mapcalc command run (see "KNOWN ISSUES" section). Noted, bestie, don't try it. ๐Ÿ’€

๐Ÿท๏ธ Raster map layer values from the category file ๐Ÿท๏ธ

Sometimes it's desirable to use a value associated with a category's label instead of the category value itself. If a raster map layer name is preceded by the @ operator, then the labels in the category file for the raster map layer are used in the expression instead of the category value. ๐ŸŽ€

For example, suppose that the raster map layer soil.ph (representing soil pH values) has a category file with labels as follows:

cat     label
------------------
0       no data
1       1.4
2       2.4
3       3.5
4       5.8
5       7.2
6       8.8
7       9.4

Then the expression:

result = @soils.ph

would produce a result with category values 0, 1.4, 2.4, 3.5, 5.8, 7.2, 8.8 and 9.4. Ate. ๐Ÿ’…๐ŸŒˆ

Note that this operator may only be applied to raster map layers, and it always spits out a floating point value in the expression. So bestie, the category label has to start with a valid number, no exceptions. If the category label is an integer, it gets represented as a floating point number. And if the category label doesn't start with a number, or it's justโ€ฆ missing? It gets represented by NULL (no data) in the resulting raster map. ๐Ÿซ 

๐ŸŒˆ Grey scale equivalents and color separates ๐ŸŒˆ

Okay so sometimes you really do need to mess with the colors assigned to map categories, and it's actually super helpful. Like, especially when the spectral properties of cells actually mean something (imagery data energy โœจ), or when the map category values are standing in for real quantities (like when category values are literally true elevation values). Map color manipulation also helps with visual recognition and with map printing. Iconic behavior tbh.

The # operator can either convert map category values into their grey scale equivalents, or extract the red, green, or blue components of a raster map layer into separate raster map layers. Obsessed.

result = #map

This converts each category value in map to a value in the range 0-255, which represents the grey scale level implied by the color for that category. If the map has a grey scale color table, then the grey level is literally what #map evaluates to. Otherwise, it gets computed as:

0.10 * red + 0.81 * green + 0.01 * blue

Alternatively, you can use:

result = y#map

to use the NTSC weightings:

0.30 * red + 0.59 * green + 0.11 * blue

Or, you can use:

result = i#map

to use equal weightings:

0.33 * red + 0.33 * green + 0.33 * blue

The # operator has three other forms: r#map, g#map, b#map. These extract the red, green, or blue components in the named raster map, respectively. The GRASS shell script r.blend extracts each of these components from two raster map layers, and combines them by a user-specified percentage. These forms let you make color separates. For example, to extract the red component from map and store it in the new 0-255 map layer red, you could type:

red = r#map

To assign this map grey colors, type:

r.colors map=red color=rules
black
white

To assign this map red colors, type:

r.colors map=red color=rules
black
red

๐Ÿงฎ Functions ๐Ÿงฎ

Here's the whole lineup of functions currently supported, bestie. The type of the result is in the last column. F means the function always results in a floating point value, I means the function gives an integer result, and * means the result is float if any of the arguments to the function are floating point values, and integer if all arguments are integer.

  • abs(x) โ€” return absolute value of x โ€” type *
  • acos(x) โ€” inverse cosine of x (result is in degrees) โ€” type F
  • asin(x) โ€” inverse sine of x (result is in degrees) โ€” type F
  • atan(x) โ€” inverse tangent of x (result is in degrees) โ€” type F
  • atan(x,y) โ€” inverse tangent of y/x (result is in degrees) โ€” type F
  • ceil(x) โ€” the smallest integral value not less than x โ€” type *
  • cos(x) โ€” cosine of x (x is in degrees) โ€” type F
  • double(x) โ€” convert x to double-precision floating point โ€” type F
  • eval([x,y,...,]z) โ€” evaluate values of listed expr, pass results to z
  • exp(x) โ€” exponential function of x โ€” type F
  • exp(x,y) โ€” x to the power y โ€” type F
  • float(x) โ€” convert x to single-precision floating point โ€” type F
  • floor(x) โ€” the largest integral value not greater than x โ€” type *
  • graph(x,x1,y1[x2,y2..]) โ€” convert the x to a y based on points in a graph โ€” type F
  • graph2(x,x1[,x2,..],y1[,y2..]) โ€” alternative form of graph() โ€” type F
  • if โ€” decision options: โ€” type *
    • if(x) โ€” 1 if x not zero, 0 otherwise
    • if(x,a) โ€” a if x not zero, 0 otherwise
    • if(x,a,b) โ€” a if x not zero, b otherwise
    • if(x,a,b,c) โ€” a if x > 0, b if x is zero, c if x < 0
  • int(x) โ€” convert x to integer [ truncates ] โ€” type I
  • isnull(x) โ€” check if x = NULL
  • log(x) โ€” natural log of x โ€” type F
  • log(x,b) โ€” log of x base b โ€” type F
  • max(x,y[,z...]) โ€” largest value of those listed โ€” type *
  • median(x,y[,z...]) โ€” median value of those listed โ€” type *
  • min(x,y[,z...]) โ€” smallest value of those listed โ€” type *
  • mod(x,y) โ€” return the modulus (the remainder) of x/y โ€” type *
  • mode(x,y[,z...]) โ€” mode value of those listed โ€” type *
  • nmax(x,y[,z...]) โ€” largest value of those listed, excluding NULLs โ€” type *
  • nmedian(x,y[,z...]) โ€” median value of those listed, excluding NULLs โ€” type *
  • nmin(x,y[,z...]) โ€” smallest value of those listed, excluding NULLs โ€” type *
  • nmode(x,y[,z...]) โ€” mode value of those listed, excluding NULLs โ€” type *
  • not(x) โ€” 1 if x is zero, 0 otherwise
  • pow(x,y) โ€” x to the power y โ€” type *
  • rand(a,b) โ€” random value x : a <= x < b โ€” type *
  • round(x) โ€” round x to nearest integer โ€” type I
  • round(x,y) โ€” round x to nearest multiple of y
  • round(x,y,z) โ€” round x to nearest y*i+z for some integer i
  • sin(x) โ€” sine of x (x is in degrees) โ€” type F
  • sqrt(x) โ€” square root of x โ€” type F
  • tan(x) โ€” tangent of x (x is in degrees) โ€” type F
  • xor(x,y) โ€” exclusive-or (XOR) of x and y โ€” type I

Internal variables, the VIPs ๐Ÿ’…:

  • row() โ€” current row of moving window โ€” type I
  • col() โ€” current col of moving window โ€” type I
  • nrows() โ€” number of rows in computation region โ€” type I
  • ncols() โ€” number of columns in computation region โ€” type I
  • x() โ€” current x-coordinate of moving window โ€” type F
  • y() โ€” current y-coordinate of moving window โ€” type F
  • ewres() โ€” current east-west resolution โ€” type F
  • nsres() โ€” current north-south resolution โ€” type F
  • area() โ€” area of current cell in square meters โ€” type F
  • null() โ€” NULL value
Heads up bestie: the row() and col() indexing starts with 1. Not zero. One. ๐Ÿซถ

๐Ÿฆ‹ Data types and their precision ๐Ÿฆ‹

There are three data types you can use inside a mapcalc expression. Each one is defined with a certain precision and a certain value range. The table below reports precision and value range info on AMD64 (x86-64) architectures compiled with GCC/CLANG. And yes, even though this setting is the most frequent one, the ranges could totally differ for different architectures. No cap.

  • int โ€” a 32-bit integer with a range from -2,147,483,647 to +2,147,483,647. The value -2,147,483,648 is reserved for NODATA.
  • float โ€” a 32-bit float (Float32) with a range from -3.4E38 to 3.4E38. However, the integer precision can be only ensured between -16,777,216 and 16,777,216. If your raster overpasses this range, it is strongly suggested to use the type double instead.
  • double โ€” a 64-bit float (Float64) with a range from -1.79E308 to 1.79E308. It is 8 bytes, 15-17 digits precision.
  • null โ€” NULL value. Refer to section "NULL support" below.
๐Ÿ‘€ Note that the value counter wraps around when the value overflows its range. E.g., if your expression is a = int(2147483648), you will get NULL value. For expression a = int(2147483649), you will reach the lowest value possible instead, i.e. -2147483647.

๐Ÿ’œ Floating point values in the expression ๐Ÿ’œ

Floating point numbers are allowed in the expression. A floating point number is a number which contains a decimal point:

2.3   12.0   12.   .81

Floating point values in the expression are handled in a special way. With arithmetic and logical operators, if either operand is float, the other is converted to float and the result of the operation is float. This means, in particular, that division of integers results in a (truncated) integer, while division of floats results in an accurate floating point value. With functions of type * (see table above), the result is float if any argument is float, integer otherwise.

Note: If you calculate with integer numbers, the resulting map will be integer. If you want to get a float result, add the decimal point to integer number(s).

If you want floating point division, at least one of the arguments has to be a floating point value. Multiplying one of them by 1.0 will produce a floating-point result, as will using float():

r.mapcalc "ndvi = float(lsat.4 - lsat.3) / (lsat.4 + lsat.3)"

๐ŸŒธ NULL support ๐ŸŒธ

  • Division by zero should result in NULL.
  • Modulus by zero should result in NULL.
  • NULL-values in any arithmetic or logical operation should result in NULL. (however, &&& and ||| are treated specially, as described below).
  • The &&& and ||| operators observe the following axioms even when x is NULL:
x &&& false == false
false &&& x == false
x ||| true == true
true ||| x == true
  • NULL-values in function arguments should result in NULL (however, if(), eval() and isnull() are treated specially, as described below).
  • The eval() function always returns its last argument
  • The situation for if() is:
if(x)
    NULL if x is NULL; 0 if x is zero; 1 otherwise
if(x,a)
    NULL if x is NULL; a if x is non-zero; 0 otherwise
if(x,a,b)
    NULL if x is NULL; a if x is non-zero; b otherwise
if(x,n,z,p)
    NULL if x is NULL; n if x is negative;
z if x is zero; p if x is positive
  • The (new) function isnull(x) returns: 1 if x is NULL; 0 otherwise. The (new) function null() (which has no arguments) returns an integer NULL.
  • Non-NULL, but invalid, arguments to functions should result in NULL. Like, don't get weird results, get NULL. Examples bestie:
log(-2)
sqrt(-2)
pow(a,b) where a is negative and b is not an integer

NULL support: ok so heads up ๐Ÿ’– โ€” literally any math you do with NULL cells always gives you back a NULL value for those cells. Periodt. If you wanna swap out a NULL cell on the fly, you gotta use the isnull() test function inside an if-statement. That's the move.

Example: say the user wants NULL-valued cells treated like zeros ๐Ÿซถ. To add maps A and B (where B is out here full of NULLs) and get a map C, the user can use a construction like:

C = A + if(isnull(B),0,B)

NULL and conditions:

For the one argument form:

if(x) = NULL        if x is NULL
if(x) = 0        if x = 0
if(x) = 1        otherwise (i.e. x is neither NULL nor 0).

For the two argument form:

if(x,a) = NULL        if x is NULL
if(x,a) = 0        if x = 0
if(x,a) = a        otherwise (i.e. x is neither NULL nor 0).

For the three argument form:

if(x,a,b) = NULL    if x is NULL
if(x,a,b) = b        if x = 0
if(x,a,b) = a        otherwise (i.e. x is neither NULL nor 0).

For the four argument form:

if(x,a,b,c) = NULL    if x is NULL
if(x,a,b,c) = a        if x > 0
if(x,a,b,c) = b        if x = 0
if(x,a,b,c) = c        if x < 0

More generally, all operators and most functions return NULL if any of their arguments are NULL. The functions if(), isnull() and eval() are the exceptions, the special ones, the main characters โœจ. The function isnull() returns 1 if its argument is NULL and 0 otherwise. If the user wants the opposite, the ! operator, e.g. "!isnull(x)" must be used.

All forms of if() return NULL if the first argument is NULL. The 2, 3 and 4 argument forms of if() return NULL if the "selected" argument is NULL, e.g.:

if(0,a,b) = b    regardless of whether a is NULL
if(1,a,b) = a    regardless of whether b is NULL

eval() always returns its last argument, so it only returns NULL if the last argument is NULL.

Note: the user cannot test for NULL using the == operator, because that returns NULL if either or both arguments are NULL ๐Ÿ˜ญ. Like, if x and y are both NULL, then "x == y" and "x != y" are both NULL rather than 1 and 0 respectively.

The behaviour makes sense if the user considers NULL as representing an unknown quantity ๐Ÿ’€. E.g. if x and y are both unknown, then the values of "x == y" and "x != y" are also unknown; if they both have unknown values, the user doesn't know whether or not they both have the same value. iykyk.


๐Ÿ’… NOTES ๐Ÿ’…

Usage from command line

Extra care must be taken if the expression is given on the command line, bestie ๐Ÿซถ. Some characters have special meaning to the UNIX shell. These include, among others:

  • *
  • (
  • )
  • >
  • &
  • |

It is advisable to put single quotes around the expression; e.g.:

'result = elevation * 2'

Without the quotes, the *, which has special meaning to the UNIX shell, would be altered and r.mapcalc would see something other than the *. ๐Ÿ™ˆ

Multiple computations

In general, it's preferable to do as much as possible in each r.mapcalc command. E.g. rather than:

r.mapcalc "$GIS_OPT_OUTPUT.r = r#$GIS_OPT_FIRST * .$GIS_OPT_PERCENT + (1.0 - .$GIS_OPT_PERCENT) * r#$GIS_OPT_SECOND"
r.mapcalc "$GIS_OPT_OUTPUT.g = g#$GIS_OPT_FIRST * .$GIS_OPT_PERCENT + (1.0 - .$GIS_OPT_PERCENT) * g#$GIS_OPT_SECOND"
r.mapcalc "$GIS_OPT_OUTPUT.b = b#$GIS_OPT_FIRST * .$GIS_OPT_PERCENT + (1.0 - .$GIS_OPT_PERCENT) * b#$GIS_OPT_SECOND"

use:

r.mapcalc <<EOF
    $GIS_OPT_OUTPUT.r = r#$GIS_OPT_FIRST * .$GIS_OPT_PERCENT + (1.0 - .$GIS_OPT_PERCENT) * r#$GIS_OPT_SECOND
    $GIS_OPT_OUTPUT.g = g#$GIS_OPT_FIRST * .$GIS_OPT_PERCENT + (1.0 - .$GIS_OPT_PERCENT) * g#$GIS_OPT_SECOND
    $GIS_OPT_OUTPUT.b = b#$GIS_OPT_FIRST * .$GIS_OPT_PERCENT + (1.0 - .$GIS_OPT_PERCENT) * b#$GIS_OPT_SECOND
    EOF

as the latter will read each input map only once ๐Ÿ”ฅ. Efficiency slay.

Backwards compatibility

For the backwards compatibility with GRASS 6, if no options are given, it manufactures file=- (which reads from stdin), so you can continue to use e.g.:

r.mapcalc < file

or:

r.mapcalc <<EOF
foo = 1
EOF

But unless you need compatibility with previous GRASS versions, use file= explicitly, as stated above. ๐Ÿ‘‘

When the map name contains uppercase letter(s) or a dot which are not allowed to be in module option names, the r.mapcalc command will be valid also without quotes:

r.mapcalc elevation_A=1
r.mapcalc elevation.1=1

However, this syntax is not recommended as quotes as stated above more safe. Using quotes is both backwards compatible and valid in future. โœจ

Interactive input in command line

For formulas that the user enters from standard input (rather than from the command line), a line continuation feature now exists ๐ŸŽ€. If the user adds a backslash to the end of an input line, r.mapcalc assumes that the formula being entered by the user continues on to the next input line. There is no limit to the possible number of input lines or to the length of a formula. So go off, write your novel-length expression.

If the r.mapcalc formula entered by the user is very long, the map title will contain only some of it, but most (if not all) of the formula will be placed into the history file for the result map. ๐Ÿ’–

Raster mask handling

r.mapcalc follows the common GRASS behavior of raster mask handling, so the mask is only applied when reading an existing GRASS raster map. This implies that, for example, the command:

r.mapcalc "elevation_exaggerated = elevation * 3"

create a map with NULL cells for the masked-out cells if raster mask is active.

However, when creating a map which is not based on any map, e.g. a map from a constant:

r.mapcalc "base_height = 200.0"

the created raster map is limited only by a computation region but it is not affected by an active raster mask. This is expected because, as mentioned above, the mask is only applied when reading, not when writing a raster map. If the raster mask should be applied in this case, an if() function including the mask raster should be used, e.g. (assuming the mask is called MASK):

r.mapcalc "base_height = if(MASK, 200.0, null())"

When testing expressions related to mask handling keep in mind that when the raster mask is active, you don't see data in masked areas even if they are not NULL ๐Ÿ‘€. See r.mask for details.

eval function

If the output of the computation should be only one map but the expression is so complex that it is better to split it to several expressions, the eval function can be used:

r.mapcalc << EOF
eval(elev_200 = elevation - 200, \
     elev_5 = 5 * elevation, \
     elev_p = pow(elev_5, 2))
elevation_result = (0.5 * elev_200) + 0.8 * elev_p
EOF

This example uses unix-like << EOF syntax to provide input to r.mapcalc.

Note that the temporary variables (maps) are not created and thus it does not matter whether they exists or not ๐Ÿฆ‹. In the example above, if map elev_200 exists it will not be overwritten and no error will be generated. The reason is that the name elev_200 now denotes the temporary variable (map) and not the existing map. The following parts of the expression will use the temporary elev_200 and the existing elev_200 will be left intact and will not be used. If a user want to use the existing map, the name of the temporary variable (map) must be changed.

Using the same map for input and output results

A map cannot be used both as an input and as an output as in this invalid expression oldmap = oldmap + 1, instead a subsequent rename using g.rename is needed when the same name is desired:

r.mapcalc "newmap = oldmap + 1"
g.rename raster=newmap,oldmap

Random number generator initialization

The pseudo-random number generator used by the rand() function can be initialised to a specific value using the seed option. This can be used to replicate a previous calculation. ๐Ÿ”

Alternatively, it can be initialised from the system time and the PID using the -s flag. This should result in a different seed being used each time.

In either case, the seed will be written to the map's history, and can be seen using r.info.

If you want other people to be able to verify your results, it's preferable to use the seed option to supply a seed which is either specified in the script or generated from a deterministic process such as a pseudo-random number generator given an explicit seed. ๐Ÿซถ

Note that the rand() function will generate a fatal error if neither the seed option nor the -s flag are given. ๐Ÿ’€

๐ŸŒˆ EXAMPLES ๐ŸŒˆ

To compute the average of two raster map layers a and b:

ave = (a + b)/2

To form a weighted average:

ave = (5*a + 3*b)/8.0

To produce a binary representation of the raster map layer a so that category 0 remains 0 and all other categories become 1:

mapmask = a != 0

This could also be accomplished by:

mapmask = if(a)

To mask raster map layer b by raster map layer a:

result = if(a,b)

To change all values below 5 to NULL, keep 5 otherwise:

newmap = if(map<5, null(), 5)

To create a map with random values in a defined range (needs either the usage of -s flag or the seed parameter). The precision of the input values determines the output precision (the resulting raster map type):

# write result as integer map (CELL)
random_int   = rand(-100,100)
# write result as double precision floating point map (DCELL)
random_dcell = rand(-100.0,100.0)
# write result as single precision floating point map (FCELL)
random_fcell = float(rand(-100.0,100.0))

The graph() function allows users to specify a x-y conversion using pairs of x,y coordinates. In some situations a transformation from one value to another is not easily established mathematically, but can be represented by a 2-D graph and then linearly interpolated. The graph() function provides the opportunity to accomplish this โœจ. An x-axis value is provided to the graph function along with the associated graph represented by a series of x,y pairs. The x values must be monotonically increasing (each larger than or equal to the previous). The graph function linearly interpolates between pairs. Any x value lower the lowest x value (i.e. first) will have the associated y value returned. Any x value higher than the last will similarly have the associated y value returned. Consider the request:

newmap = graph(map, 1,10, 2,25, 3,50)

X (map) values supplied and y (newmap) values returned:

0, 10
1, 10
1.5, 17.5
2.9, 47.5
4, 50
100, 50

๐Ÿ’€ KNOWN ISSUES ๐Ÿ’€

Okay bestie, gather round, because this is the part where we spill the tea on everything r.mapcalc is not gonna do for you. No cap. ๐Ÿซ–

  • The result variable on the left hand side of the equation should not appear in the expression on the right hand side. Like, the glow-up needs a fresh face. Don't do this:
mymap = if( mymap > 0, mymap, 0)
  • Any maps generated by a r.mapcalc command only exist after the entire command has completed. โœจ All maps are generated concurrently, row-by-row โ€” meaning there's basically an implicit for row in rows {...} wrapped around the whole expression. So the #, @, and [ ] operators cannot be used on a map that was generated within the same r.mapcalc command run. Which means this (the strikethrough code, i.e. this straight-up does not work ๐Ÿ˜ญ) is a no:
newmap = oldmap * 3.14
othermap = newmap[-1, 0] / newmap[1, 0]
  • Continuation lines must end with a \ and have no trailing white space (blanks or tabs). If you do leave white space at the end of a continuation line, the error messages r.mapcalc gives you will be meaningless and the equation will not do what you meant it to do. This is especially important for the eval() function. ๐Ÿšจ
  • Currently, there is no comment mechanism in r.mapcalc. The vibe suggests maybe someday adding a thing where the whole line gets ignored if the user puts a # at the start of a line, as if it weren't there โ€” that would do the trick. ๐Ÿ™ˆ
  • The function should really require the user to type "end" or "exit" instead of just a blank line. That would make separating multiple scripts doable by white space. (As it stands: bestie, it just doesn't.)
  • r.mapcalc does not print a warning in case of operations on NULL cells. It is left to you, the user, to use the isnull() function. ๐Ÿ’…

๐Ÿ“š REFERENCES ๐Ÿ“š

The academic lore, if you will ๐Ÿ‘‘

  • r.mapcalc: An Algebra for GIS and Image Processing, by Michael Shapiro and Jim Westervelt, U.S. Army Construction Engineering Research Laboratory (March/1991).
  • Performing Map Calculations on GRASS Data: r.mapcalc Program Tutorial, by Marji Larson, Michael Shapiro and Scott Tweddale, U.S. Army Construction Engineering Research Laboratory (December 1991)
  • Grey scale conversion is based on the C.I.E. x,y,z system where y represents luminance. See "Fundamentals of Digital Image Processing," by Anil K. Jain (Prentice Hall, NJ, 1989; p 67).

๐Ÿ‘€ SEE ALSO ๐Ÿ‘€

Go be nosy, bestie, these are the siblings:

  • g.region, r.bitpattern, r.blend, r.colors, r.fillnulls, r.mapcalc.simple

๐Ÿ‘‘ AUTHORS ๐Ÿ‘‘

  • Michael Shapiro, U.S.Army Construction Engineering Research Laboratory
  • Glynn Clements

๐Ÿ”— SOURCE CODE ๐Ÿ”—

  • Available at: r.mapcalc source code (history)
  • Latest change: Friday Apr 17 01:04:26 2026 in commit: 32a143358df26b2e518b03c7e0d79541e8d1e5b1
  • Main index | Raster index | Topics index | Keywords index | Graphical index | Full index
  • ยฉ 2003-2026 GRASS Development Team, GRASS 8.5.0 Reference Manual

and that's the whole page, bestie ๐ŸŽ€ this fine specimen of documentation is from GRASS 8.5.0, is literally titled r.mapcalc(1grass) in the footer, and dates back to that version and its commit history as noted above. Read it, love it, go map some cells ๐Ÿ’–โœจ๐Ÿฆ‹๐ŸŒธ๐Ÿซถ

๐Ÿ“ฆ where this page comes from