[GRASS-user] Problems Running r.flow

GRASS GRASS user
1

   In -6.4svn from last Friday's weekly snapshot.

   After setting "g.region rast=aber5m" I tried running "r.flow elevin=aber5m
aspin=aber5aspect flout=aberFlowline lgout=aberFlowlength" and ran into
resolution problems:

Reading input files: elevation
ERROR: Elevation file's resolution differs from current region resolution

   g.region -p reports:

nsres: 4.99982001
ewres: 4.99979979

   g.region rast=aber5m -p reports:

nsres: 4.99982001
ewres: 4.99979979

   I also get the same error message about the aspect map (derived from the
elevation map) which has the same resolution as above.

   Any ideas what's going on? More importantly, how do I get his module
running here?

Rich

2

Hi Rich,

try g.region nsres=5 ewres=5, may be you can bypass this.

best regards

milton

2010/1/27 Rich Shepard <rshepard@appl-ecosys.com>

In -6.4svn from last Friday’s weekly snapshot.

After setting “g.region rast=aber5m” I tried running “r.flow elevin=aber5m
aspin=aber5aspect flout=aberFlowline lgout=aberFlowlength” and ran into
resolution problems:

Reading input files: elevation
ERROR: Elevation file’s resolution differs from current region resolution

g.region -p reports:

nsres: 4.99982001
ewres: 4.99979979

g.region rast=aber5m -p reports:

nsres: 4.99982001
ewres: 4.99979979

I also get the same error message about the aspect map (derived from the
elevation map) which has the same resolution as above.

Any ideas what’s going on? More importantly, how do I get his module
running here?

Rich

grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user

3

On Wed, 27 Jan 2010, Milton Cezar Ribeiro wrote:

try g.region nsres=5 ewres=5, may be you can bypass this.

milton,

   That's how I defined the resolution, but it was returned as slightly less
for some reason. See:

GRASS 6.4.0svn (Oregon):/usr4/grassbase > g.region res=5 -p
projection: 99 (Lambert Conformal Conic)
zone: 0
datum: nad83
ellipsoid: grs80
north: 1334419.25160578
south: 1279151.24118496
west: 769192.9282895
east: 819255.92362222
nsres: 4.99982001
ewres: 4.99979979
rows: 11054
cols: 10013
cells: 110683702

Rich

4

On Thu, Jan 28, 2010 at 1:59 AM, Rich Shepard <rshepard@appl-ecosys.com> wrote:

On Wed, 27 Jan 2010, Milton Cezar Ribeiro wrote:

try g.region nsres=5 ewres=5, may be you can bypass this.

milton,

That's how I defined the resolution, but it was returned as slightly less
for some reason. See:

GRASS 6.4.0svn (Oregon):/usr4/grassbase > g.region res=5 -p
projection: 99 (Lambert Conformal Conic)
zone: 0
datum: nad83
ellipsoid: grs80
north: 1334419.25160578
south: 1279151.24118496
west: 769192.9282895
east: 819255.92362222
nsres: 4.99982001
ewres: 4.99979979

If you add -a to g.region, you obtain the desired resolution by
(automatically) a slightly enlarged region which make the pixels fit
into it properly.

Markus

5

On Thu, 28 Jan 2010, Markus Neteler wrote:

If you add -a to g.region, you obtain the desired resolution by
(automatically) a slightly enlarged region which make the pixels fit into
it properly.

Markus,

   Thank you. I'm still learning about the -a option and when to use it.

   Perhaps this is why some output maps fail.

Rich

6

On Thu, Jan 28, 2010 at 2:46 AM, Rich Shepard <rshepard@appl-ecosys.com> wrote:

On Thu, 28 Jan 2010, Markus Neteler wrote:

If you add -a to g.region, you obtain the desired resolution by
(automatically) a slightly enlarged region which make the pixels fit into
it properly.

Markus,

Thank you. I'm still learning about the -a option and when to use it.

Perhaps this is why some output maps fail.

You have to set to the raster map first. If that has been generated without -a
then perhaps consider to re-do it. I personally dislike odd resolution values
like 4.99982001 which may introduce artifacts later (so, I am g.region -p
maniac to always keep an eye on these things).

Markus

7

On Thu, 28 Jan 2010, Markus Neteler wrote:

You have to set to the raster map first.

   That's what I do, Markus, when I start GRASS and when I change from the
10m resolution maps to the 5m resolution maps.

   What I _really_ need to figure out is why so many modules fail for me. I
cannot get proper output from r.flow, r.param.scale (using the example on
the bottom of that Web page), or the curve outputs from r.shade.aspect.
There must be something I'm not doing properly and I've yet to see what that
is.

Rich

8

Hum
Interesting -a flag.

But g.region res=5 and g.region nsres=5 ewres=5 produce the same result? I never tryed.

Rich, another way is you use n= s= w= e= options, that may be will result something like -a, I suppose.

cheers

milton

2010/1/27 Rich Shepard <rshepard@appl-ecosys.com>

On Wed, 27 Jan 2010, Milton Cezar Ribeiro wrote:

try g.region nsres=5 ewres=5, may be you can bypass this.

milton,

That’s how I defined the resolution, but it was returned as slightly less
for some reason. See:

GRASS 6.4.0svn (Oregon):/usr4/grassbase > g.region res=5 -p
projection: 99 (Lambert Conformal Conic)
zone: 0
datum: nad83
ellipsoid: grs80
north: 1334419.25160578
south: 1279151.24118496
west: 769192.9282895
east: 819255.92362222

nsres: 4.99982001
ewres: 4.99979979

rows: 11054
cols: 10013
cells: 110683702

Rich

grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user

9

On Wed, 27 Jan 2010, Milton Cezar Ribeiro wrote:

But g.region res=5 and g.region nsres=5 ewres=5 produce the same result? I
never tryed.

Milton,

   Yes. res is equal to setting both nsres and weres to the same values.

Rich, another way is you use n= s= w= e= options, that may be will result
something like -a, I suppose.

   That changes the bounds. I don't want to do that other than the small
increment necessary to get an integer resolution value.

Rich

10

Rich,

If your are using a resolution of ~5m, I suggest you adjust your region ~1m. I think that the final result will not suffer so much changes. If so, may be you can try finer resolution (4m?).

g.region n=1334420 s=1279150 w=769190 e=819255 res=5

good luck

milton

2010/1/27 Rich Shepard <rshepard@appl-ecosys.com>

On Wed, 27 Jan 2010, Milton Cezar Ribeiro wrote:

But g.region res=5 and g.region nsres=5 ewres=5 produce the same result? I
never tryed.

Milton,

Yes. res is equal to setting both nsres and weres to the same values.

Rich, another way is you use n= s= w= e= options, that may be will result
something like -a, I suppose.

That changes the bounds. I don’t want to do that other than the small
increment necessary to get an integer resolution value.

Rich

grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user

11

Rich Shepard wrote:

   In -6.4svn from last Friday's weekly snapshot.

   After setting "g.region rast=aber5m" I tried running "r.flow elevin=aber5m
aspin=aber5aspect flout=aberFlowline lgout=aberFlowlength" and ran into
resolution problems:

Reading input files: elevation
ERROR: Elevation file's resolution differs from current region resolution

   g.region -p reports:

nsres: 4.99982001
ewres: 4.99979979

   g.region rast=aber5m -p reports:

nsres: 4.99982001
ewres: 4.99979979

That's a bug in r.flow:

    if (!((region.ew_res == hd.ew_res)
    && (region.ns_res == hd.ns_res)))
  G_fatal_error(_("Elevation file's resolution differs from current region resolution"));

Testing floating-point values for exact equality is usually wrong.
E.g. the above code will generate an error even if the values only
differ by one part in one quadrillion (the "epsilon" for IEEE
double-precision floating point is ~2.22e-16, and represents the
smallest value which can be added to 1.0 without the result being
rounded to 1.0).

Even if you use g.region rast=... to set the current region to the
map's region, the current region may not exactly match due to rounding
errors.

g.region *always* recalculates the resolution as (north-south)/rows
and (east-west)/cols. If you tell it to preserve the resolution, it
adjusts the bounds so that the resulting resolution should remain
unchanged, but rounding errors may ultimately cause it to end up
differing in the least significant digit.

Given that IEEE double precision is sufficient to represent e.g. the
circumference of the earth with an error of roughly 9 nanometres,
rounding errors in the least significant digit really shouldn't be an
issue.

   I also get the same error message about the aspect map (derived from the
elevation map) which has the same resolution as above.

   Any ideas what's going on? More importantly, how do I get his module
running here?

Try the attached patch. It replaces the exact comparison with a check
that the resolutions agree to within 1ppm.

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

(attachments)

r.flow.patch (1.37 KB)

12

Hello grass-users,

I have a problem running a r.mapcalc procedure:

given two maps
dir_x,
dir_y

with 0,1,-1 entries

I want to run:

r.mapcalc "the_map=if(isnull(the_map),the_map[dir_x,dir_y],null())"

in order to make a r.water.outlet on many many areas at the same time.

The error I get is:
syntax error, unexpected NAME, expecting INTEGER or '-'

why does dir_x/dir_y does not return an integer here? Using int(dir_x) also does not work.

Until now I use a r.mapcalc command with many "if", checking the direction of every pixel every time. But this is quite slow.

Achim

13

dir_x and dir_y are maps, and as far as I know they can not use in map[x,y] since it expects relative coordinates for the moving window...even though your raster values are in [0,1,-1]

Christian.

--------------------------------------------------
From: "Achim Kisseler" <ak7@jupiter.uni-freiburg.de>
Sent: Thursday, January 28, 2010 2:13 PM
To: <grass-users@lists.osgeo.org>
Subject: [GRASS-user] r.mapcalc question

Hello grass-users,

I have a problem running a r.mapcalc procedure:

given two maps
dir_x,
dir_y

with 0,1,-1 entries

I want to run:

r.mapcalc "the_map=if(isnull(the_map),the_map[dir_x,dir_y],null())"

in order to make a r.water.outlet on many many areas at the same time.

The error I get is:
syntax error, unexpected NAME, expecting INTEGER or '-'

why does dir_x/dir_y does not return an integer here? Using int(dir_x) also does not work.

Until now I use a r.mapcalc command with many "if", checking the direction of every pixel every time. But this is quite slow.

Achim
_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user

14

Hi Christian,

thanks for fast reply!

I also tried r.mapcalc with eval(), setting variables for dir_x/dir_y. But map[ , ] does not take them with the same error.

Looks like I have to find another way to fasten the if-procedure

in order to make a r.water.outlet on many many areas at the same time.

Cheers,
Achim

15

Hi Achim,

the neighborhood modifier doesn't seem to accept functions (you also
could try rand(), ...)

Marco

Achim Kisseler schrieb:

Hello grass-users,

I have a problem running a r.mapcalc procedure:

given two maps
dir_x,
dir_y

with 0,1,-1 entries

I want to run:

r.mapcalc "the_map=if(isnull(the_map),the_map[dir_x,dir_y],null())"

in order to make a r.water.outlet on many many areas at the same time.

The error I get is:
syntax error, unexpected NAME, expecting INTEGER or '-'

why does dir_x/dir_y does not return an integer here? Using int(dir_x)
also does not work.

Until now I use a r.mapcalc command with many "if", checking the
direction of every pixel every time. But this is quite slow.

Achim
_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user

16

I see,

so variables are functions in this case.

I've been searching for a while, but I cannot find any tool like grasss r.mapcalc that is a bit more flexible, eg. that can handle recursive tasks.

Thanks a lot,
Achim

Marco Lechner - FOSSGIS e.V. schrieb:

Hi Achim,

the neighborhood modifier doesn't seem to accept functions (you also
could try rand(), ...)

Marco

Achim Kisseler schrieb:

Hello grass-users,

I have a problem running a r.mapcalc procedure:

given two maps
dir_x,
dir_y

with 0,1,-1 entries

I want to run:

r.mapcalc "the_map=if(isnull(the_map),the_map[dir_x,dir_y],null())"

in order to make a r.water.outlet on many many areas at the same time.

The error I get is:
syntax error, unexpected NAME, expecting INTEGER or '-'

why does dir_x/dir_y does not return an integer here? Using int(dir_x)
also does not work.

Until now I use a r.mapcalc command with many "if", checking the
direction of every pixel every time. But this is quite slow.

Achim
_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user

_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user

17

On Thu, 28 Jan 2010, Glynn Clements wrote:

That's a bug in r.flow:

   if (!((region.ew_res == hd.ew_res)
    && (region.ns_res == hd.ns_res)))
  G_fatal_error(_("Elevation file's resolution differs from current region resolution"));

Try the attached patch. It replaces the exact comparison with a check that
the resolutions agree to within 1ppm.

Glynn,

   Perhaps I've specified the region incorrectly, but the patch makes no
difference. Upon invoking GRASS I run 'g.region rast=aber5m res=5 -ap'
followed by 'r.flow elevin=aber5m aspin=aber5aspect flout=aberFlowline
lgout=aberFlowlength' and immediately see:

ERROR: Elevation file's resolution differs from current region resolution

   I'm surprised this issue's not come up before.

Thanks,

Rich

18

Christian Schwartze wrote:

> I have a problem running a r.mapcalc procedure:
>
> given two maps
> dir_x,
> dir_y
>
> with 0,1,-1 entries
>
> I want to run:
>
> r.mapcalc "the_map=if(isnull(the_map),the_map[dir_x,dir_y],null())"
>
> in order to make a r.water.outlet on many many areas at the same time.
>
> The error I get is:
> syntax error, unexpected NAME, expecting INTEGER or '-'
>
> why does dir_x/dir_y does not return an integer here? Using int(dir_x)
> also does not work.
>
> Until now I use a r.mapcalc command with many "if", checking the direction
> of every pixel every time. But this is quite slow.

dir_x and dir_y are maps, and as far as I know they can not use in map[x,y]
since it expects relative coordinates for the moving window...even though
your raster values are in [0,1,-1]

To elaborate: the parser only recognises integer literals for the
neighbourhood modifier.

Apart from that, allowing expressions would require knowing in advance
the range of the data, so that it knows how many rows need to be
cached.

I've been searching for a while, but I cannot find any tool like grasss
r.mapcalc that is a bit more flexible, eg. that can handle recursive tasks.

For tasks such as this, the alternatives are either to use a
general-purpose language such as MatLab/Octave or Python+NumPy/SciPy
(which requires reading the entire map into memory), or write a GRASS
module in C/C++ (you could use the SWIG bindings to write a complete
module in Python, but it would be slow, and the bindings are rather
primitive at the moment).

The main problem with extending r.mapcalc is that it's a slippery
slope; the more features we add, the more features we'll be asked to
add and the harder it will be to add such features (as each new
feature has to consider interactions with all existing feature).

As it stands, r.mapcalc provides just enough "glue" (or possibly duct
tape) to connect other modules together. E.g. by allowing simple pre-
and post-processing of the data used or produced by other modules, it
saves us from needing to add many relatively simple special-purpose
modules as well as adding many features to many existing modules.

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

19

Rich Shepard wrote:

> That's a bug in r.flow:
>
> if (!((region.ew_res == hd.ew_res)
> && (region.ns_res == hd.ns_res)))
> G_fatal_error(_("Elevation file's resolution differs from current region resolution"));

> Try the attached patch. It replaces the exact comparison with a check that
> the resolutions agree to within 1ppm.

   Perhaps I've specified the region incorrectly, but the patch makes no
difference. Upon invoking GRASS I run 'g.region rast=aber5m res=5 -ap'
followed by 'r.flow elevin=aber5m aspin=aber5aspect flout=aberFlowline
lgout=aberFlowlength' and immediately see:

ERROR: Elevation file's resolution differs from current region resolution

Did you try "g.region rast=aber5m" without the other options?

If that doesn't work, you can try changing the tolerance value to
something larger than 1e-6, or simply disabling the test.

The test exists to ensure that the data is read cell-for-cell; if any
duplicates occurred due to resampling, it would introduce plateaux
into the elevation surface. For that purpose, it only matters that the
absolute difference between the resolutions is less than
ns_res/(2*rows) or ew_res/(2*cols), i.e. a cumulative error of less
than 1/2 cell.

   I'm surprised this issue's not come up before.

Same here. It's possible that any rounding error is usually masked by
truncation of the value written to the WIND and cellhd files. If it
doesn't affect the decimal value written to the file, it won't matter,
but if it causes the last decimal digit to change, it will.

The files are written using %.15g (i.e. 15 significant figures), while
double-precision floating-point has just under 17. There are 11
distinct values which are equal to 5 when written using %.15g, so
there's roughly a 9% chance that an error in the least significant bit
would affect the decimal representation.

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

20

On Thu, 28 Jan 2010, Glynn Clements wrote:

Did you try "g.region rast=aber5m" without the other options?

Glynn,

   Not until just now. It's busy calculating maps now.

   Why does not specifying res or '-ap' make such a difference? Perhaps I can
solve future situations like this myself if I better understand what's going
on under the hood.

Thanks very much,

Rich