[GRASS-user] r.mapcalc median filters entire map instead of only null values

GRASS GRASS user
1

I'm trying to use r.neighbors and r.mapcalc in a script to fill holes in a raster where there are null values. A temporary raster is created with r.neighbors that contains a 3x3 median-interpolated surface of the primary DEM. r.mapcalc should only use the value in this temporary raster where a null value occurs in the primary DEM.

Here's the relevant section from my script:

# Procedure for filling nulls with multiple median filter passes

if [ "$PASSES" -gt 1 ] ; then

COUNT=1
while [ "$COUNT" -le "$PASSES" ] ; do
          
  r.neighbors --o input=${MAP} output=fill.tmp method=median size=3
  r.mapcalc "${MAP}_fill = if(isnull(${MAP}), fill.tmp, ${MAP})"
        
  MAP=${MAP}_fill
  COUNT=$(($COUNT + 1))
  g.remove rast=fill.tmp

done

else
echo ""
r.neighbors --o input=${MAP} output=fill.tmp method=median size=3
r.mapcalc "${MAP}_fill = if(isnull("$MAP"), fill.tmp, "$MAP")"
MAP=${MAP}_fill
g.remove rast=fill.tmp

fi

My problem is that when my script is run, the output $MAP_fill raster has it's entire surface smoothed with the median-filtered one, not just the null values like I intended. But if I run the r.neighbors and r.mapcalc lines in isolation outside of the script, I get the intended null-filled DEM with no unwanted smoothing. Am I not quoting the r.mapcalc line properly?

FWIW, I've also tried different quoting methods for the r.mapcalc line, with the same problem occurring each time:

r.mapcalc ${MAP}_fill = "if(isnull($"MAP"), fill.tmp, $"MAP")"
r.mapcalc "${MAP}_fill = if(isnull(${MAP}), fill.tmp, ${MAP})"
r.mapcalc ${MAP}_fill = "if(isnull($MAP), fill.tmp, $MAP)"

Any ideas?

~ Eric.

2

Patton, Eric wrote:

I'm trying to use r.neighbors and r.mapcalc in a script to fill holes
in a raster where there are null values. A temporary raster is created
with r.neighbors that contains a 3x3 median-interpolated surface of
the primary DEM. r.mapcalc should only use the value in this temporary
raster where a null value occurs in the primary DEM.

Here's the relevant section from my script:

# Procedure for filling nulls with multiple median filter passes

if [ "$PASSES" -gt 1 ] ; then

COUNT=1
while [ "$COUNT" -le "$PASSES" ] ; do
          
  r.neighbors --o input=${MAP} output=fill.tmp method=median size=3
  r.mapcalc "${MAP}_fill = if(isnull(${MAP}), fill.tmp, ${MAP})"
        
  MAP=${MAP}_fill
  COUNT=$(($COUNT + 1))
  g.remove rast=fill.tmp

done

else
echo ""
r.neighbors --o input=${MAP} output=fill.tmp method=median size=3
r.mapcalc "${MAP}_fill = if(isnull("$MAP"), fill.tmp, "$MAP")"
MAP=${MAP}_fill
g.remove rast=fill.tmp

fi

My problem is that when my script is run, the output $MAP_fill raster
has it's entire surface smoothed with the median-filtered one, not
just the null values like I intended. But if I run the r.neighbors and
r.mapcalc lines in isolation outside of the script, I get the intended
null-filled DEM with no unwanted smoothing.

FWIW, the above script (with PASSES=5 and MAP set to a test map) works
fine for me; the non-null cells remain untouched.

Am I not quoting the r.mapcalc line properly?

No, both r.mapcalc commands are fine.

FWIW, I've also tried different quoting methods for the r.mapcalc
line, with the same problem occurring each time:

r.mapcalc ${MAP}_fill = "if(isnull($"MAP"), fill.tmp, $"MAP")"

This is bogus.

r.mapcalc "${MAP}_fill = if(isnull(${MAP}), fill.tmp, ${MAP})"
r.mapcalc ${MAP}_fill = "if(isnull($MAP), fill.tmp, $MAP)"

Either of these should work.

r.mapcalc simply concatenates all of its arguments, with spaces in
between. Usually, it's best to quote the entire expression with double
quotes.

You only need to use braces around a variable name if the substitution
is immediately followed by a character which is valid as part of a
variable name, e.g. ${MAP}_fill. Adding braces where they aren't
required is harmless.

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

3

Thanks for checking the script code, Glynn. I'm going to have to investigate what's going on locally on my machine.

~ Eric.

-----Original Message-----
From: Glynn Clements
To: Patton, Eric
Cc: grassuser@grass.itc.it
Sent: 1/17/2007 3:23 PM
Subject: Re: [GRASS-user] r.mapcalc median filters entire map instead of onlynull values

Patton, Eric wrote:

I'm trying to use r.neighbors and r.mapcalc in a script to fill holes
in a raster where there are null values. A temporary raster is created
with r.neighbors that contains a 3x3 median-interpolated surface of
the primary DEM. r.mapcalc should only use the value in this temporary
raster where a null value occurs in the primary DEM.

Here's the relevant section from my script:

# Procedure for filling nulls with multiple median filter passes

if [ "$PASSES" -gt 1 ] ; then

COUNT=1
while [ "$COUNT" -le "$PASSES" ] ; do
          
  r.neighbors --o input=${MAP} output=fill.tmp method=median size=3
  r.mapcalc "${MAP}_fill = if(isnull(${MAP}), fill.tmp, ${MAP})"
        
  MAP=${MAP}_fill
  COUNT=$(($COUNT + 1))
  g.remove rast=fill.tmp

done

else
echo ""
r.neighbors --o input=${MAP} output=fill.tmp method=median size=3
r.mapcalc "${MAP}_fill = if(isnull("$MAP"), fill.tmp, "$MAP")"
MAP=${MAP}_fill
g.remove rast=fill.tmp

fi

My problem is that when my script is run, the output $MAP_fill raster
has it's entire surface smoothed with the median-filtered one, not
just the null values like I intended. But if I run the r.neighbors and
r.mapcalc lines in isolation outside of the script, I get the intended
null-filled DEM with no unwanted smoothing.

FWIW, the above script (with PASSES=5 and MAP set to a test map) works
fine for me; the non-null cells remain untouched.

Am I not quoting the r.mapcalc line properly?

No, both r.mapcalc commands are fine.

FWIW, I've also tried different quoting methods for the r.mapcalc
line, with the same problem occurring each time:

r.mapcalc ${MAP}_fill = "if(isnull($"MAP"), fill.tmp, $"MAP")"

This is bogus.

r.mapcalc "${MAP}_fill = if(isnull(${MAP}), fill.tmp, ${MAP})"
r.mapcalc ${MAP}_fill = "if(isnull($MAP), fill.tmp, $MAP)"

Either of these should work.

r.mapcalc simply concatenates all of its arguments, with spaces in
between. Usually, it's best to quote the entire expression with double
quotes.

You only need to use braces around a variable name if the substitution
is immediately followed by a character which is valid as part of a
variable name, e.g. ${MAP}_fill. Adding braces where they aren't
required is harmless.

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

4

Patton, Eric wrote:

I'm trying to use r.neighbors and r.mapcalc in a script to fill holes
in a raster where there are null values. A temporary raster is created
with r.neighbors that contains a 3x3 median-interpolated surface of
the primary DEM. r.mapcalc should only use the value in this temporary
raster where a null value occurs in the primary DEM.

perhaps modify/extend the r.fillnulls script? A flag to use an alternate
interpolation technique might be nice.

Hamish

5

It's ok, Glynn tested the script out and it worked for him, so I can't
explain what was happening. My script even works for me now. User error,
I guess :wink:

I didn't know r.fillnulls was a script, as I've never used it before.
Maybe I could take a look at using other interpolation methods as you
suggest.

~ Eric.

--
Eric Patton
email: epatton@nrcan.gc.ca

-----Original Message-----
From: Hamish [mailto:hamish_nospam@yahoo.com]
Sent: Tuesday, February 06, 2007 4:58 AM
To: Patton, Eric
Cc: grassuser@grass.itc.it
Subject: Re: [GRASS-user] r.mapcalc median filters entire map
instead of only null values

Patton, Eric wrote:
> I'm trying to use r.neighbors and r.mapcalc in a script to
fill holes
> in a raster where there are null values. A temporary raster
is created
> with r.neighbors that contains a 3x3 median-interpolated surface of
> the primary DEM. r.mapcalc should only use the value in
this temporary
> raster where a null value occurs in the primary DEM.

perhaps modify/extend the r.fillnulls script? A flag to use
an alternate interpolation technique might be nice.

Hamish

6

Eric wrote:

I didn't know r.fillnulls was a script, as I've never used it before.
Maybe I could take a look at using other interpolation methods as you
suggest.

see-
M. Neteler. SRTM and VMAP0 data in OGR and GRASS. GRASS Newsletter,
3:2-6, June 2005. ISSN 1614-8746.

Hamish