[GRASS-user] Trying to reclassify raster using non-linear functions in GRASS

GRASS GRASS user
1

Dear all,

I am using ArcInfo Workstation and Grid in order to to produce Gaussian functions of elevation from a DEM with minimum value at 800m, as follow:

elev_x800 = elev90 - 800
elev_x800sq = pow(elev_x800, 2)
elev_x800sqn = -1 * elev_x800sq
elev_x800sqnd = elev_x800sqn / 500000
elev_x800exp = exp(elev_x800sqnd)
elev_x800gaus = elev_x800exp * 0.000798
elev_x800gi = (elev_x800gaus * -11330) + 10.044

I was wondering how to apply that in GRASS (because I want to move towards Open Software). Can you help me?

Cheers,

Dr. Gabriel Zorello Laporta
Instituto Oscar Freire, Informática Médica, Faculdade de Medicina
Universidade de São Paulo
Avenida Doutor Arnaldo, 455 - 3° andar (corredor central)
Bairro Pacaembu, São Paulo/SP
CEP 01246-903
Tel (55) 11 3061-7386

2

Hi Gabriel,

On Fri, Nov 8, 2013 at 1:54 AM, Gabriel Zorello Laporta
<gabriellaporta@gmail.com> wrote:

I am using ArcInfo Workstation and Grid in order to to produce Gaussian
functions of elevation from a DEM with minimum value at 800m, as follow:

elev_x800 = elev90 - 800
elev_x800sq = pow(elev_x800, 2)
elev_x800sqn = -1 * elev_x800sq
elev_x800sqnd = elev_x800sqn / 500000
elev_x800exp = exp(elev_x800sqnd)
elev_x800gaus = elev_x800exp * 0.000798
elev_x800gi = (elev_x800gaus * -11330) + 10.044

I was wondering how to apply that in GRASS (because I want to move towards
Open Software). Can you help me?

You can use r.mapcalc:

r.mapcalc expression="elev_x800 = elev90 - 800"

You can see more info about raster algebra, here:

http://grass.osgeo.org/grass64/manuals/r.mapcalc.html

Have a nice day.

Pietro

3

Gabriel Zorello Laporta wrote:

I am using ArcInfo Workstation and Grid in order to to produce Gaussian
functions of elevation from a DEM with minimum value at 800m, as follow:

elev_x800 = elev90 - 800
elev_x800sq = pow(elev_x800, 2)
elev_x800sqn = -1 * elev_x800sq
elev_x800sqnd = elev_x800sqn / 500000
elev_x800exp = exp(elev_x800sqnd)
elev_x800gaus = elev_x800exp * 0.000798
elev_x800gi = (elev_x800gaus * -11330) + 10.044

I was wondering how to apply that in GRASS (because I want to move towards
Open Software). Can you help me?

You should be able to use those exact same expressions in r.mapcalc,
e.g.

  r.mapcalc <<EOF
  elev_x800 = elev90 - 800
  elev_x800sq = pow(elev_x800, 2)
  elev_x800sqn = -1 * elev_x800sq
  elev_x800sqnd = elev_x800sqn / 500000
  elev_x800exp = exp(elev_x800sqnd)
  elev_x800gaus = elev_x800exp * 0.000798
  elev_x800gi = (elev_x800gaus * -11330) + 10.044
  EOF

The above will write out each step as a map; if you only need the
final result (elev_x800gi), you can omit all of the intermediate steps
and use e.g.:

  r.mapcalc "elev_x800gi = (exp(-((elev90 - 800)^2) / 500000) * 0.000798 * -11330) + 10.044"

Or, if you want to keep the expressions in the above form but only
want the last one written out as a map, use eval(), e.g.:

  r.mapcalc <<EOF
  eval(elev_x800 = elev90 - 800,\
       elev_x800sq = pow(elev_x800, 2),\
       elev_x800sqn = -1 * elev_x800sq,\
       elev_x800sqnd = elev_x800sqn / 500000,\
       elev_x800exp = exp(elev_x800sqnd),\
       elev_x800gaus = elev_x800exp * 0.000798)
  elev_x800gi = (elev_x800gaus * -11330) + 10.044
  EOF

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

4

On Fri, Nov 8, 2013 at 8:14 AM, Glynn Clements <glynn@gclements.plus.com>wrote:

Gabriel Zorello Laporta wrote:

> I am using ArcInfo Workstation and Grid in order to to produce Gaussian
> functions of elevation from a DEM with minimum value at 800m, as follow:
>
> elev_x800 = elev90 - 800
> elev_x800sq = pow(elev_x800, 2)
> elev_x800sqn = -1 * elev_x800sq
> elev_x800sqnd = elev_x800sqn / 500000
> elev_x800exp = exp(elev_x800sqnd)
> elev_x800gaus = elev_x800exp * 0.000798
> elev_x800gi = (elev_x800gaus * -11330) + 10.044
>
> I was wondering how to apply that in GRASS (because I want to move
towards
> Open Software). Can you help me?

You should be able to use those exact same expressions in r.mapcalc,
e.g.

        r.mapcalc <<EOF
        elev_x800 = elev90 - 800
        elev_x800sq = pow(elev_x800, 2)
        elev_x800sqn = -1 * elev_x800sq
        elev_x800sqnd = elev_x800sqn / 500000
        elev_x800exp = exp(elev_x800sqnd)
        elev_x800gaus = elev_x800exp * 0.000798
        elev_x800gi = (elev_x800gaus * -11330) + 10.044
        EOF

The above will write out each step as a map; if you only need the
final result (elev_x800gi), you can omit all of the intermediate steps
and use e.g.:

        r.mapcalc "elev_x800gi = (exp(-((elev90 - 800)^2) / 500000) *
0.000798 * -11330) + 10.044"

Or, if you want to keep the expressions in the above form but only
want the last one written out as a map, use eval(), e.g.:

        r.mapcalc <<EOF
        eval(elev_x800 = elev90 - 800,\
             elev_x800sq = pow(elev_x800, 2),\
             elev_x800sqn = -1 * elev_x800sq,\
             elev_x800sqnd = elev_x800sqn / 500000,\
             elev_x800exp = exp(elev_x800sqnd),\
             elev_x800gaus = elev_x800exp * 0.000798)
        elev_x800gi = (elev_x800gaus * -11330) + 10.044
        EOF

This is an interesting feature. I added this example to r.mapcalc manual

in G7 [r58173].

It is in Notes section. I was not sure where to put it and it is only a
dummy example. Maybe, it would be helpful to have this "Gaussian functions"
example in examples but I don't know what it is, so I don't want to write
about it.

It is bit unfortunate that this example does not work in mapcalc in wxGUI.
I'm not sure if it is a serious issue and where to mention it
(manual/gui/...).

Vaclav

https://trac.osgeo.org/grass/changeset/58173

PS: MarkusN, I've tried to do it I passive voice this time.

--

Glynn Clements <glynn@gclements.plus.com>
_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user