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

GRASS GRASS user
1

I’m working on the preparation of the statement of Agricultural Capacity in general terms the final map has classes ranging 1-8.

To get this map I use three themes with their respective notes:

  • Map of Slope: with classes ranging 1 to 8;
  • Soil Map: with classes ranging 1 to 8;
  • Map of flood risk: with classes ranging 1 to 8;

I - After converting all to a vector format and “merge” the three layers, I have multiple polygons with various combinations, eg:
a) 1,3,1
b) 5,2,1
c) 4,2,1
d) 1,1,8

II - To determine the final value of the class is chosen the highest value:
a) 3
b) 5
c) 4
d) 8

I need to automate this task because the last step (II) is very slow (actually manually), which would be the best way to do this? I think the whole process can be done with the data in raster format, but have not found a formula that I can apply in r.mapcalc to this problem.

Thanks in advance,

Marcello Benigno B. de Barros Filho
Prof. do Curso Superior de Tecnologia em Geoprocessamento - IFPB
Mestre em Ciências Geodésicas e Tecnologias da Geoinformação - UFPE
Doutorando em Tecnologia Ambiental e Recursos Hídricos - UFPE
http://profmarcello.blogspot.com
http://about.me/marcello.benigno

2

Not sure I understand the logic of taking the highest value, but you can do this simply with the original rasters using: r.series input=slope,soils,floods out=agri_risk method=max Then convert to vectors only at the end.

···

On 27/05/2013 16:41, Marcello Benigno wrote:

I’m working on the preparation of the statement of Agricultural Capacity in general terms the final map has classes ranging 1-8.

To get this map I use three themes with their respective notes:

  • Map of Slope: with classes ranging 1 to 8;
  • Soil Map: with classes ranging 1 to 8;
  • Map of flood risk: with classes ranging 1 to 8;

I - After converting all to a vector format and “merge” the three layers, I have multiple polygons with various combinations, eg:
a) 1,3,1
b) 5,2,1
c) 4,2,1
d) 1,1,8

II - To determine the final value of the class is chosen the highest value:
a) 3
b) 5
c) 4
d) 8

I need to automate this task because the last step (II) is very slow (actually manually), which would be the best way to do this? I think the whole process can be done with the data in raster format, but have not found a formula that I can apply in r.mapcalc to this problem.

Thanks in advance,

Marcello Benigno B. de Barros Filho
Prof. do Curso Superior de Tecnologia em Geoprocessamento - IFPB
Mestre em Ciências Geodésicas e Tecnologias da Geoinformação - UFPE
Doutorando em Tecnologia Ambiental e Recursos Hídricos - UFPE
http://profmarcello.blogspot.com
http://about.me/marcello.benigno

This mail was received via Mail-SeCure System.

_______________________________________________
grass-user mailing list
[grass-user@lists.osgeo.org](mailto:grass-user@lists.osgeo.org)
[http://lists.osgeo.org/mailman/listinfo/grass-user](http://lists.osgeo.org/mailman/listinfo/grass-user)

This mail was received via Mail-SeCure System.


-- 
Micha Silver
GIS Consulting
052-3665918
[http://www.surfaces.co.il](http://www.surfaces.co.il)
3

Hi Micha,

About:

Not sure I understand the logic of taking the highest value

This was defined by an agronomist, I also do not understand why =/

However, your answer is exactly what I need, did not know r.series, thank you!

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())’

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

How do I create a single expression for this problem?

Regards

Marcello Benigno B. de Barros Filho
Prof. do Curso Superior de Tecnologia em Geoprocessamento - IFPB
Mestre em Ciências Geodésicas e Tecnologias da Geoinformação - UFPE
Doutorando em Tecnologia Ambiental e Recursos Hídricos - UFPE
http://profmarcello.blogspot.com
http://about.me/marcello.benigno

4

How about using r.reclass ? Create a text file “slope_reclass.txt” which contains rows: 0 thru 2 = 1 2 thru 5 = 2 5 thru 10 = 3 … 70 thru 100 = 7 then run r.reclass in=slope out=slope_rc rule=slope_reclass.txt Be aware that the reclass map is not a “real” raster. If you need to use it for other things, then a second r.mapcalc expression will be necessary to create on on-disk raster map. Such as: r.mapcalc slope_reclass=slope_rc Regards, Micha

···

On 27/05/2013 22:36, Marcello Benigno wrote:

Hi Micha,

About:

Not sure I understand the logic of taking the highest value

This was defined by an agronomist, I also do not understand why =/

However, your answer is exactly what I need, did not know r.series, thank you!

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())’

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

How do I create a single expression for this problem?

Regards

Marcello Benigno B. de Barros Filho
Prof. do Curso Superior de Tecnologia em Geoprocessamento - IFPB
Mestre em Ciências Geodésicas e Tecnologias da Geoinformação - UFPE
Doutorando em Tecnologia Ambiental e Recursos Hídricos - UFPE
http://profmarcello.blogspot.com
http://about.me/marcello.benigno

This mail was received via Mail-SeCure System.

-- 
Micha Silver
GIS Consulting
052-3665918
[http://www.surfaces.co.il](http://www.surfaces.co.il)
5

Micha,

Thank you very for helping me.

Regards,

Marcello Benigno B. de Barros Filho
Prof. do Curso Superior de Tecnologia em Geoprocessamento - IFPB
Mestre em Ciências Geodésicas e Tecnologias da Geoinformação - UFPE
Doutorando em Tecnologia Ambiental e Recursos Hídricos - UFPE
http://profmarcello.blogspot.com
http://about.me/marcello.benigno

6

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@gclements.plus.com>

7

On 28/05/13 21:33, Glynn Clements wrote:

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

Or keep the FCELL or DCELL and just use r.recode:

r.recode input=slope output=slope.recode rules=-<<EOF
*:2:1
2:5:2
5:10:3
10:15:4
15:45:5
45:70:6
70:*:7
EOF

?

Moritz

8

Glynn wrote:

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=\- &lt;&lt;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

note that r.reclass can read floating point maps too, but
it's a bit buggy/weird,
  https://trac.osgeo.org/grass/ticket/1525#comment:10

Hamish

9

Hamish wrote:

> 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

note that r.reclass can read floating point maps too, but
it's a bit buggy/weird,
  https://trac.osgeo.org/grass/ticket/1525#comment:10

Right; r.reclass just generates the table, the actuall reclass is
performed by lib/raster. If the base map is floating-point, it will be
coerced to integer using its associated quantisation rules.

Note that r.reclass doesn't pay any attention to the map's
quantisation rules (even if it did, they could be changed between the
reclass being created and being used, and the rules at the time of use
would have effect).

The buggy rounding code should probably be replaced with e.g.

  *v = (CELL) floor(fv + 0.5);
or:
  *v = (CELL) ceil(fv - 0.5);
or:
  *v = (CELL) (fv > 0 ? floor(fv + 0.5) : ceil(fv - 0.5));

These differ on whether values lying exactly on a 0.5 boundary are
rounded toward positive infinity, toward negative infinity, or away
from zero respectively. The last one mimics the behaviour of C99's
round() function.

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