[GRASS-user] Suggestion for creating a script to create a map of agricultural capacity

[Glynn Clements]

Marcello Benigno wrote:

> If you can help me also in this tip, I would be very grateful, I'm
> reclassifying the slope as follows:
> 
> r.mapcalc class_1 = 'if(slope <= 2.0, 1, null())'
> r.mapcalc class_2 = 'if(2.0 < slope <= 5, 2, null())'
> r.mapcalc class_3 = 'if(5 < slope <= 10, 3, null())'
> r.mapcalc class_4 = 'if(10 < slope <= 15, 4, null())'
> r.mapcalc class_5 = 'if(15 < slope <= 45, 5, null())'
> r.mapcalc class_6 = 'if(45 < slope <= 70, 6, null())'
> r.mapcalc class_7 = 'if(slope > 70, 7, null())'

1. Note that "a < b <= c" is wrong; you need to use "a < b && b <= c".

2. You could optimise the above by using a single r.mapcalc command,
e.g. (assuming Unix shell syntax):

	r.mapcalc <<EOF
	class_1 = if(slope <= 2.0, 1, null())
	class_2 = if(2.0 < slope && slope <= 5, 2, null())
	class_3 = if(5 < slope && slope <= 10, 3, null())
	class_4 = if(10 < slope && slope <= 15, 4, null())
	class_5 = if(15 < slope && slope <= 45, 5, null())
	class_6 = if(45 < slope && slope <= 70, 6, null())
	class_7 = if(slope && slope > 70, 7, null())
	EOF

This will generate all 7 output maps in one go, reading the input map
only once (most r.mapcalc calculations spend more time reading and
writing raster maps than performing calculations).

3. This:

> r.mapcalc slope.reclass = class_1 + class_2 + class_3 + class_4 + class_5 +
> class_6 + class_7

will result in an all-null map, as arithmetic operators return null if
either input is null. Changing the "null()" with "0" in the
calculation of the class_* maps would fix this, or you could just use
r.patch instead of r.mapcalc.

Sticking with r.mapcalc, it would be simpler (and more efficient) to
generate the output map directly from the input, without generating
any intermediate maps, e.g.:

	r.mapcalc <<EOF
	slope.reclass = \
	if(slope <= 2,  1,\
	if(slope <= 5,  2,\
	if(slope <= 10, 3,\
	if(slope <= 15, 4,\
	if(slope <= 45, 5,\
	if(slope <= 70, 6,\
	7))))))
	EOF

[The above could be written on one line, but splitting the expression
with backslash-newline improves legibility.]

But for this specific case, it's probably simpler to just generate an
integer version of the slope map (if it isn't already) and use
r.reclass, e.g.

	r.reclass input=slope output=slope.reclass rules=- <<EOF
	0 thru 2 = 1
	3 thru 5 = 2
	6 thru 10 = 3
	11 thru 15 = 4
	16 thru 45 = 5
	46 thru 70 = 6
	* = 7
	EOF

-- 
Glynn Clements <glynn at gclements.plus.com>