[GRASS-user] Problems Running r.flow

[Glynn Clements]

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

-------------- next part --------------
Index: raster/r.flow/io.c
===================================================================
--- raster/r.flow/io.c	(revision 40694)
+++ raster/r.flow/io.c	(working copy)
@@ -27,6 +27,7 @@
 #include <stdio.h>
 #include <stdlib.h>
 #include <unistd.h>
+#include <math.h>
 #include <grass/glocale.h>
 #include "r.flow.h"
 #include "mem.h"
@@ -68,6 +69,12 @@
     return G_open_cell_old(fname, mapset);
 }
 
+static int compare_regions(const struct Cell_head *a, const struct Cell_head *b)
+{
+    return (fabs(a->ew_res - b->ew_res) < 1e-6 * b->ew_res &&
+	    fabs(a->ns_res - b->ns_res) < 1e-6 * b->ns_res);
+}
+
 void read_input_files(void)
 {
     DCELL *barc;
@@ -77,8 +84,7 @@
     G_message(_("Reading input files: elevation"));
 
     fd = open_existing_cell_file(parm.elevin, &hd);
-    if (!((region.ew_res == hd.ew_res)
-	  && (region.ns_res == hd.ns_res)))
+    if (!compare_regions(&region, &hd))
 	G_fatal_error(_("Elevation file's resolution differs from current region resolution"));
 
     for (row = 0; row < region.rows; row++) {
@@ -93,8 +99,7 @@
     if (parm.aspin) {
 	G_message(_("Reading input files: aspect"));
 	fd = open_existing_cell_file(parm.aspin, &hd);
-	if (!((region.ew_res == hd.ew_res)
-	      && (region.ns_res == hd.ns_res)))
+	if (!compare_regions(&region, &hd))
 	    G_fatal_error(_("Resolution of aspect file differs from "
 			    "current region resolution"));