1r.mapcalc(1) Grass User's Manual r.mapcalc(1)
2
6 r.mapcalc
7
9 r.mapcalc performs arithmetic on raster map layers. New raster map
10 layers can be created which are arithmetic expressions involving exist‐
11 ing raster map layers, integer or floating point constants, and func‐
12 tions.
13
15 If used without command line arguments, r.mapcalc will read its input,
16 one line at a time, from standard input (which is the keyboard, unless
17 redirected from a file or across a pipe). Otherwise, the expression on
18 the command line is evaluated. r.mapcalc expects its input to have the
19 form:
20 result=expression
22 where result is the name of a raster map layer to contain the result of
24 the calculation and expression is any legal arithmetic expression
25 involving existing raster map layers, integer or floating point con‐
26 stants, and functions known to the calculator. Parentheses are allowed
27 in the expression and may be nested to any depth. result will be cre‐
28 ated in the user's current mapset.
29 The formula entered to r.mapcalc by the user is recorded both in the
31 result map title (which appears in the category file for result) and in
32 the history file for result.
33 Some characters have special meaning to the command shell. If the user
35 is entering input to r.mapcalc on the command line, expressions should
36 be enclosed within single quotes. See NOTES, below.
37
39 The following operators are supported:
40 Operator Meaning Type Precedence
41 --------------------------------------------------------------
42 ! not Logical 6
43 ^ exponentiation Arithmetic 5
44 % modulus Arithmetic 4
45 / division Arithmetic 4
46 * multiplication Arithmetic 4
47 + addition Arithmetic 3
48 - subtraction Arithmetic 3
49 == equal Logical 2
50 != not equal Logical 2
51 > greater than Logical 2
52 >= greater than or equal Logical 2
53 < less than Logical 2
54 <= less than or equal Logical 2
55 && and Logical 1
56 || or Logical 1
57 &&& and[1] Logical 1
58 ||| or[1] Logical 1
59 # color separator operator Arithmetic
60 (modulus is the remainder upon division)
61 [1] The &&& and ||| operators handle null values differently to
63 other operators. See the section entitled NULL support below for more
64 details.
65 The operators are applied from left to right, with those of higher
67 precedence applied before those with lower precedence. Division by 0
68 and modulus by 0 are acceptable and give a NULL result. The logical
69 operators give a 1 result if the comparison is true, 0 otherwise.
70
72 Anything in the expression which is not a number, operator, or function
73 name is taken to be a raster map layer name. Examples:
74 elevation
77 x3
78 3d.his
79 Most GRASS raster map layers meet this naming convention. However, if
82 a raster map layer has a name which conflicts with the above rule, it
83 should be quoted. For example, the expression
84 x = a-b
87 would be interpreted as: x equals a minus b, whereas
90 x = "a-b"
93 would be interpreted as: x equals the raster map layer named a-b
96 Also
98 x = 3107
101 would create x filled with the number 3107, while
104 x = "3107"
107 would copy the raster map layer 3107 to the raster map layer x.
110 Quotes are not required unless the raster map layer names look like
112 numbers or contain operators, OR unless the program is run non-interac‐
113 tively. Examples given here assume the program is run interactively.
114 See NOTES, below.
115 r.mapcalc will look for the raster map layers according to the user's
117 current mapset search path. It is possible to override the search path
118 and specify the mapset from which to select the raster map layer. This
119 is done by specifying the raster map layer name in the form:
120 name@mapset
123 For example, the following is a legal expression:
126 result = x@PERMANENT / y@SOILS
129 The mapset specified does not have to be in the mapset search path.
132 (This method of overriding the mapset search path is common to all
133 GRASS commands, not just r.mapcalc.)
134
136 Maps and images are data base files stored in raster format, i.e., two-
137 dimensional matrices of integer values. In r.mapcalc, maps may be fol‐
138 lowed by a neighborhood modifier that specifies a relative offset from
139 the current cell being evaluated. The format is map[r,c], where r is
140 the row offset and c is the column offset. For example, map[1,2]
141 refers to the cell one row below and two columns to the right of the
142 current cell, map[-2,-1] refers to the cell two rows above and one col‐
143 umn to the left of the current cell, and map[0,1] refers to the cell
144 one column to the right of the current cell. This syntax permits the
145 development of neighborhood-type filters within a single map or across
146 multiple maps.
147
149 Sometimes it is desirable to use a value associated with a category's
150 label instead of the category value itself. If a raster map layer name
151 is preceded by the @ operator, then the labels in the category file for
152 the raster map layer are used in the expression instead of the category
153 value.
154 For example, suppose that the raster map layer soil.ph (representing
156 soil pH values) has a category file with labels as follows:
157 cat label
160 ------------------
161 0 no data
162 1 1.4
163 2 2.4
164 3 3.5
165 4 5.8
166 5 7.2
167 6 8.8
168 7 9.4
169 Then the expression:
172 result = @soils.ph
175 would produce a result with category values 0, 1.4, 2.4, 3.5, 5.8, 7.2,
178 8.8 and 9.4.
179 Note that this operator may only be applied to raster map layers and
181 produces a floating point value in the expression. Therefore, the cat‐
182 egory label must start with a valid number. If the category label is
183 integer, it will be represented by a floating point number. I the cate‐
184 gory label does not start with a number or is missing, it will be rep‐
185 resented by NULL (no data) in the resulting raster map.
186
188 It is often helpful to manipulate the colors assigned to map cate‐
189 gories. This is particularly useful when the spectral properties of
190 cells have meaning (as with imagery data), or when the map category
191 values represent real quantities (as when category values reflect true
192 elevation values). Map color manipulation can also aid visual recogni‐
193 tion, and map printing.
194 The # operator can be used to either convert map category values to
196 their grey scale equivalents or to extract the red, green, or blue com‐
197 ponents of a raster map layer into separate raster map layers.
198 result = #map
201 converts each category value in map to a value in the range 0-255 which
204 represents the grey scale level implied by the color for the category.
205 If the map has a grey scale color table, then the grey level is what
206 #map evaluates to. Otherwise, it is computed as:
207 0.10 * red + 0.81 * green + 0.01 * blue
210 Alternatively, you can use:
213 result = y#map
216 to use the NTSC weightings:
219 0.30 * red + 0.59 * green + 0.11 * blue
222 Or, you can use:
225 result = i#map
228 to use equal weightings:
231 0.33 * red + 0.33 * green + 0.33 * blue
234 The # operator has three other forms: r#map, g#map, b#map. These
237 extract the red, green, or blue components in the named raster map,
238 respectively. The GRASS shell script r.blend extracts each of these
239 components from two raster map layers, and combines them by a user-
240 specified percentage. These forms allow color separates to be made.
241 For example, to extract the red component from map and store it in the
242 new 0-255 map layer red, the user could type:
243 red = r#map
246 To assign this map grey colors type:
249 r.colors map=red color=rules
252 black
253 white
254 To assign this map red colors type:
257 r.colors map=red color=rules
260 black
261 red
262
265 The functions currently supported are listed in the table below. The
266 type of the result is indicated in the last column. F means that the
267 functions always results in a floating point value, I means that the
268 function gives an integer result, and * indicates that the result is
269 float if any of the arguments to the function are floating point values
270 and integer if all arguments are integer.
271 function description type
274 ---------------------------------------------------------------------------
275 abs(x) return absolute value of x *
276 acos(x) inverse cosine of x (result is in degrees) F
277 asin(x) inverse sine of x (result is in degrees) F
278 atan(x) inverse tangent of x (result is in degrees) F
279 atan(x,y) inverse tangent of y/x (result is in degrees) F
280 cos(x) cosine of x (x is in degrees) F
281 double(x) convert x to double-precision floating point F
282 eval([x,y,...,]z) evaluate values of listed expr, pass results to z
283 exp(x) exponential function of x F
284 exp(x,y) x to the power y F
285 float(x) convert x to single-precision floating point F
286 graph(x,x1,y1[x2,y2..]) convert the x to a y based on points in a
287 graph F
288 if decision options: *
289 if(x) 1 if x not zero, 0 otherwise
290 if(x,a) a if x not zero, 0 otherwise
291 if(x,a,b) a if x not zero, b otherwise
292 if(x,a,b,c) a if x > 0, b if x is zero, c if x < 0
293 int(x) convert x to integer [ truncates ] I
294 isnull(x) check if x = NULL
295 log(x) natural log of x F
296 log(x,b) log of x base b F
297 max(x,y[,z...]) largest value of those listed *
298 median(x,y[,z...]) median value of those listed *
299 min(x,y[,z...]) smallest value of those listed *
300 mode(x,y[,z...]) mode value of those listed
301 *
302 not(x) 1 if x is zero, 0 otherwise
303 pow(x,y) x to the power y *
304 rand(a,b) random value x : a < x < b
305 round(x) round x to nearest integer I
306 sin(x) sine of x (x is in degrees) F
307 sqrt(x) square root of x F
308 tan(x) tangent of x (x is in degrees) F
309 Internal variables:
311 row() current row of moving window
312 col() current col of moving window
313 x() current x-coordinate of moving window
314 y() current y-coordinate of moving window
315 ewres() current east-west resolution
316 nsres() current north-south resolution
317 null() NULL value
318 Note, that the row() and col() indexing starts with 1.
319
321 Floating point numbers are allowed in the expression. A floating point
322 number is a number which contains a decimal point:
323 2.3 12.0 12. .81
324 Floating point values in the expression are handled in a special way.
325 With arithmetic and logical operators, if either operand is float, the
326 other is converted to float and the result of the operation is float.
327 This means, in particular that division of integers results in a (trun‐
328 cated) integer, while division of floats results in an accurate float‐
329 ing point value. With functions of type * (see table above), the
330 result is float if any argument is float, integer otherwise.
331 Note: If you calculate with integer numbers, the resulting map will be
333 integer. If you want to get a float result, add the decimal point to
334 integer number(s).
335 If you want floating point division, at least one of the arguments has
337 to be a floating point value. Multiplying one of them by 1.0 will pro‐
338 duce a floating-point result, as will using float():
339 r.mapcalc "ndvi=float(lsat.4 - lsat.3) / (lsat.4 + lsat.3)"
340
343 Division by zero should result in NULL. Modulus by zero should result
344 in NULL. NULL-values in any arithmetic or logical operation should
345 result in NULL. (however, &&& and ||| are treated specially, as
346 described below). The &&& and ||| operators observe the following
347 axioms even when x is NULL:
348 x &&& false == false
349 false &&& x == false
350 x ||| true == true
351 true ||| x == true
352 NULL-values in function arguments should result in NULL (however,
353 if(), eval() and isnull() are treated specially, as described below).
354 The eval() function always returns its last argument The situation for
355 if() is:
356 if(x)
357 NULL if x is NULL; 0 if x is zero; 1 otherwise
358 if(x,a)
359 NULL if x is NULL; a if x is non-zero; 0 otherwise
360 if(x,a,b)
361 NULL if x is NULL; a if x is non-zero; b otherwise
362 if(x,n,z,p)
363 NULL if x is NULL; n if x is negative;
364 z if x is zero; p if x is positive
365 The (new) function isnull(x) returns: 1 if x is NULL; 0 otherwise. The
366 (new) function null() (which has no arguments) returns an integer NULL.
367 Non-NULL, but invalid, arguments to functions should result in NULL.
368 Examples:
369 log(-2)
370 sqrt(-2)
371 pow(a,b) where a is negative and b is not an integer
372 NULL support: Please note that any math performed with NULL cells
375 always results in a NULL value for these cells. If you want to replace
376 a NULL cell on-the-fly, use the isnull() test function in a if-state‐
377 ment.
378 Example: The users wants the NULL-valued cells to be treated like
380 zeros. To add maps A and B (where B contains NULLs) to get a map C the
381 user can use a construction like:
382 C=A + if(isnull(B),0,B)
385 NULL and conditions:
388 For the one argument form:
390 if(x) = NULL if x is NULL
391 if(x) = 0 if x = 0
392 if(x) = 1 otherwise (i.e. x is neither NULL nor 0).
393 For the two argument form:
396 if(x,a) = NULL if x is NULL
397 if(x,a) = 0 if x = 0
398 if(x,a) = a otherwise (i.e. x is neither NULL nor 0).
399 For the three argument form:
402 if(x,a,b) = NULL if x is NULL
403 if(x,a,b) = b if x = 0
404 if(x,a,b) = a otherwise (i.e. x is neither NULL nor 0).
405 For the four argument form:
408 if(x,a,b,c) = NULL if x is NULL
409 if(x,a,b,c) = a if x > 0
410 if(x,a,b,c) = b if x = 0
411 if(x,a,b,c) = c if x < 0
412 More generally, all operators and most functions return NULL if *any*
413 of their arguments are NULL.
414 The functions if(), isnull() and eval() are exceptions.
415 The function isnull() returns 1 if its argument is NULL and 0 other‐
416 wise. If the user wants the opposite, the ! operator, e.g.
417 "!isnull(x)" must be used.
418 All forms of if() return NULL if the first argument is NULL. The 2, 3
420 and 4 argument forms of if() return NULL if the "selected" argument is
421 NULL, e.g.:
422 if(0,a,b) = b regardless of whether a is NULL
423 if(1,a,b) = a regardless of whether b is NULL
424 eval() always returns its last argument, so it only returns NULL if
425 the last argument is NULL.
426 Note: The user cannot test for NULL using the == operator, as that
428 returns NULL if either or both arguments are NULL, i.e. if x and y are
429 both NULL, then "x == y" and "x != y" are both NULL rather than 1 and 0
430 respectively.
431 The behaviour makes sense if the user considers NULL as representing an
432 unknown quantity. E.g. if x and y are both unknown, then the values of
433 "x == y" and "x != y" are also unknown; if they both have unknown val‐
434 ues, the user doesn't know whether or not they both have the same
435 value.
436
438 To compute the average of two raster map layers a and b:
439 ave = (a + b)/2
440 To form a weighted average:
441 ave = (5*a + 3*b)/8.0
442 To produce a binary representation of the raster map layer a so that
443 category 0 remains 0 and all other categories become 1:
444 mask = a != 0
445 This could also be accomplished by:
446 mask = if(a)
447 To mask raster map layer b by raster map layer a:
448 result = if(a,b)
449 To change all values below 5 to NULL:
450 newmap = if(map<5, null(), 5)
451 The graph function allows users to specify a x-y conversion using
452 pairs of x,y coordinates. In some situations a transformation from one
453 value to another is not easily established mathematically, but can be
454 represented by a 2-D graph. The graph() function provides the opportu‐
455 nity to accomplish this. An x-axis value is provided to the graph
456 function along with the associated graph represented by a series of x,y
457 pairs. The x values must be monotonically increasing (each larger than
458 or equal to the previous). The graph function linearly interpolates
459 between pairs. Any x value lower the lowest x value (i.e. first) will
460 have the associated y value returned. Any x value higher than the last
461 will similarly have the associated y value returned. Consider the
462 request:
463 newmap = graph(map, 1,10, 2,25, 3,50)
464 X (map) values supplied and y (newmap) values returned:
465 0, 10
466 1, 10,
467 1.5, 16.5
468 2.9, 47.5
469 4, 50
470 100, 50
471
474 Extra care must be taken if the expression is given on the command
475 line. Some characters have special meaning to the UNIX shell. These
476 include, among others:
477 * ( ) > & |
479 It is advisable to put single quotes around the expression; e.g.:
481 result = 'elevation * 2'
482 Without the quotes, the *, which has special meaning to the UNIX
483 shell, would be altered and r.mapcalc would see something other than
484 the *.
485 If the input comes directly from the keyboard and the result raster map
487 layer exists, the user will be asked if it can be overwritten. Other‐
488 wise, the result raster map layer will automatically be overwritten if
489 it exists.
490 Quoting result is not allowed. However, it is never necessary to quote
492 result since it is always taken to be a raster map layer name.
493 For formulas that the user enters from standard input (rather than from
495 the command line), a line continuation feature now exists. If the user
496 adds \e to the end of an input line, r.mapcalc assumes that the formula
497 being entered by the user continues on to the next input line. There
498 is no limit to the possible number of input lines or to the length of a
499 formula.
500 If the r.mapcalc formula entered by the user is very long, the map
502 title will contain only some of it, but most (if not all) of the for‐
503 mula will be placed into the history file for the result map.
504 When the user enters input to r.mapcalc non-interactively on the com‐
506 mand line, the program will not warn the user not to overwrite existing
507 map layers. Users should therefore take care to assign program outputs
508 raster map names that do not yet exist in their current mapsets.
509
511 Continuation lines must end with a \ and have NO trailing white space
512 (blanks or tabs). If the user does leave white space at the end of
513 continuation lines, the error messages produced by r.mapcalc will be
514 meaningless and the equation will not work as the user intended. This
515 is important for the eval() function.
516 Error messages produced by r.mapcalc are almost useless. In future,
518 r.mapcalc should make some attempt to point the user to the offending
519 section of the equation, e.g.:
520 x = a * b ++ c
521 ERROR: somewhere in line 1: ... b ++ c ...
522 Currently, there is no comment mechanism in r.mapcalc. Perhaps adding
525 a capability that would cause the entire line to be ignored when the
526 user inserted a # at the start of a line as if it were not present,
527 would do the trick.
528 The function should require the user to type "end" or "exit" instead of
530 simply a blank line. This would make separation of multiple scripts
531 separable by white space.
532 r.mapcalc does not print a warning in case of operations on NULL cells.
534 It is left to the user to utilize the isnull() function.
535
537 r.mapcalc: An Algebra for GIS and Image Processing, by Michael Shapiro
538 and Jim Westervelt, U.S. Army Construction Engineering Research Labora‐
539 tory (March/1991).
540 Performing Map Calculations on GRASS Data: r.mapcalc Program Tutorial,
542 by Marji Larson, Michael Shapiro and Scott Tweddale, U.S. Army Con‐
543 struction Engineering Research Laboratory (December 1991)
544 Grey scale conversion is based on the C.I.E. x,y,z system where y rep‐
546 resents luminance. See "Fundamentals of Digital Image Processing," by
547 Anil K. Jain (Prentice Hall, NJ, 1989; p 67).
548 g.region, r.bitpattern, r.blend, r.colors, r.fillnulls
550
552 Michael Shapiro, U.S.Army Construction Engineering Research Laboratory
553 Glynn Clements
555 Last changed: $Date: 2007/01/22 13:05:21 $
557GRASS 6.2.2 r.mapcalc(1)