r.mask with NULL map doesn't work

[ferrouswheel]

When presenting r.mask with map containing all NULL values, it accepts the map and sets it as a mask.

However, despite setting it as a mask, it does nothing. Instead of preventing all data from being available to read, it allows all data to be read.

This is counter-intuitive. r.mask should either complain when setting a mask of all NULL values, or it should make it impossible to access any data (because there are no non-NULL cells in the source map).

While one might question the use of a mask of all NULL cells, this can easily occur when scripting GRASS commands.

[hamish]

Hi,

I can reproduce it,

setup:

# fully null input map for r.mask:
G65> g.region -d
G65> r.mapcalc "allnull = null()"

G65> r.info -r allnull 
min=NULL
max=NULL

btw r.univar is broken for that, with --quiet you get no output at all!

what r.mask does:

G65> echo "* = 1" | r.reclass input=allnull output=MASK.test

G65> r.info -r MASK.test
min=1
max=1

again, this time with partially null input map for r.mask:

G65> r.mapcalc "partnull = if(row() > 100, row(), null())"

r.info -r partnull
min=101
max=477

echo "* = 1" | r.reclass input=partnull output=partial.MASK.test

G65> r.univar partial.MASK.test
total null and non-null cells: 302418
total null cells: 63400
n: 239018
minimum: 1
maximum: 1

so it is a quirk of how r.reclass works. I'm not sure, but it seems like it could be an intended feature. One which should be documented in the man page of course...

Keeping r.reclass as-is, here's a possible solution:

eval `r.info -r allnull`
if [ "$min" = "NULL" -a "$max" = "NULL" ] ; then
   RECLASS_TO="NULL"
else
   RECLASS_TO="1"
fi
echo "* = $RECLASS_TO" | r.reclass input="$GIS_OPT_INPUT" output=MASK

Hamish

[hamish]

patch vs devbr6 for r.mask on an all-null file

[glynn]

Replying to hamish:

what r.mask does:

G65> echo "* = 1" | r.reclass input=allnull output=MASK.test

G65> r.info -r MASK.test
min=1
max=1

The range file generated by r.reclass is based upon the range of values in the reclass table, regardless of the values present in the base map. E.g.:

$ r.mapcalc --o 'foo = 2'
$ r.reclass --o in=foo out=bar rules=-
Enter rule(s), "end" when done, "help" if you need it
CELL: Data range is 2 to 2
> 1 = 1
> 2 = 2
> 3 = 3
> end
$ r.info -r bar
min=1
max=3

However, that isn't the real problem here.

1. A map which is a reclass of an all-null map is read as all-null:

   $ echo "* = 1" | r.reclass input=allnull output=testmask rules=-
   $ r.mapcalc --o 'foo = testmask'
   $ r.info -r foo
   min=NULL
   max=NULL
   

2. A MASK map which is a reclass of an all-null map doesn't mask anything out:

   $ echo "* = 1" | r.reclass input=allnull output=MASK rules=-
   $ r.mapcalc --o 'foo = elevation.dem'
   $ g.remove rast=MASK
   $ r.info -r foo
   min=1066
   max=1840
   

3. A (non-reclass) all-null MASK map masks everything out:

   $ r.mapcalc 'MASK = allnull'
   $ r.mapcalc --o 'foo = elevation.dem'
   $ g.remove rast=MASK
   $ r.info -r foo
   min=NULL
   max=NULL
   

Clearly, the code in lib/raster/get_row.c which reads a MASK map doesn't handle reclass maps correctly.

Specifically, embed_mask() calls get_map_row_nomask() then, if it's a reclass map, do_reclass_int(). Finally, any cell whose value is zero is marked as null.

However, get_map_row_nomask() doesn't embed nulls (so any nulls will be read as zero). When the mask isn't a reclass, this works okay, as nulls are stored with zero as the cell value. But do_reclass_int() expects nulls to have been embedded; a zero will be treated as such and reclassed.

[glynn]

Replying to glynn:

Clearly, the code in lib/raster/get_row.c which reads a MASK map doesn't handle reclass maps correctly.

Possible fix:

--- lib/raster/get_row.c	(revision 56625)
+++ lib/raster/get_row.c	(working copy)
@@ -919,11 +919,13 @@
 	return;
     }
 
-    if (R__.fileinfo[R__.mask_fd].reclass_flag)
+    if (R__.fileinfo[R__.mask_fd].reclass_flag) {
+	embed_nulls(R__.mask_fd, mask_buf, row, CELL_TYPE, 0, 0);
 	do_reclass_int(R__.mask_fd, mask_buf, 1);
+    }
 
     for (i = 0; i < R__.rd_window.cols; i++)
-	if (mask_buf[i] == 0)
+	if (mask_buf[i] == 0 || Rast_is_c_null_value(&mask_buf[i]))
 	    flags[i] = 1;
 
     G__freea(mask_buf);

[neteler]

Replying to glynn:

Replying to glynn:

Clearly, the code in lib/raster/get_row.c which reads a MASK map doesn't handle reclass maps correctly.

Possible fix:

--- lib/raster/get_row.c	(revision 56625)
+++ lib/raster/get_row.c	(working copy)
@@ -919,11 +919,13 @@

...

Glynn, do you plan to apply that to G7? And is a backport possible/recommended?

[glynn]

Replying to neteler:

Glynn, do you plan to apply that to G7?

Try r58723.

And is a backport possible/recommended?

I have no idea.