[GRASS-user] eliminate isolated pixel groups

[Glynn Clements]

Florian Kindl wrote:

> I am working on binary (1;null) maps depicting streams.
> What I try to accomplish is a filter which removes isolated groups of
> pixels.
> It is easy to remove a single pixel with r.mapcalc:
> 
> r.mapcalc "$RO= \
> if( \
> ( \
> isnull($RI[-1,-1])&&isnull($RI[-1,0])&&isnull($RI[-1,1])&& \
> isnull($RI[0,-1])&&isnull($RI[0,1])&& \
> isnull($RI[1,-1])&&isnull($RI[1,0])&&isnull($RI[1,1])  \
> ) \
> ,null() \
> ,$RI \
> )"

A couple of quick tips about r.mapcalc expressions:

1. Read complex expressions from stdin using a here-document, e.g.:

	r.mapcalc <<EOF
	$RO=
	if(
	(
	isnull($RI[-1,-1])&&isnull($RI[-1,0])&&isnull($RI[-1,1])&&
	isnull($RI[0,-1])&&isnull($RI[0,1])&&
	isnull($RI[1,-1])&&isnull($RI[1,0])&&isnull($RI[1,1]) 
	)
	,null()
	,$RI
	)
	EOF

2. Binary maps are easier to process if they are zero/one rather than
null/one. Use "r.mapcalc outmap = isnull(inmap)" to create a 0/1 map
first. Then you can omit the isnull() functions from the complex
expressions, e.g.:

	r.mapcalc <<EOF
	$RO=
	if(
	$NULL[-1,-1] && $NULL[-1,0] && $NULL[-1,1] && 
	$NULL[0,-1] && $NULL[0,1] && 
	$NULL[1,-1] && $NULL[1,0] && $NULL[1,1] 
	,null()
	,$RI
	)
	EOF

> - test all 8 neighbours and null() the center pixel or leave it
>   unchanged.
> 
> Now, what if I try to expand this concept:
> 
>   * * * * *
>   * + + + *
>   * + + + *
>   * + + + *
>   * * * * *
> 
> IF the outer ring (cells symbolized by a "*") fullfills a certain
> condition (say, all cells are null)
> THEN _all_ cells of the inner area (the "+" in the above scheme
> are to be nulled.
> 
> In other words, I have to modify more than one cell of the output raster
> at the same time, modify more than one cell at one step of the moving window.
> 
> Has anybody an idea how to accomplish this, either with r.mapcalc or
> some other tool?

My first attempt would be:

1. r.mapcalc "${RI}.null = isnull(${RI})"

2. r.mfilter input=${RI}.null ouput=${RI}.count filter=ring5x5

where ring5x5 is:

	TITLE	5x5 sum
	MATRIX	5
	1 1 1 1 1
	1 0 0 0 1
	1 0 0 0 1
	1 0 0 0 1
	1 1 1 1 1
	DIVISOR	1
	TYPE	P

3. r.mapcalc "${RI}.ring = ${RI}.count == 20"

4. r.mfilter input=${RI}.ring ouput=${RI}.change filter=box3x3

where box3x3 is:

	TITLE	3x3 sum
	MATRIX	3
	1 1 1
	1 1 1
	1 1 1
	DIVISOR	1
	TYPE	P

5. r.mapcalc "${RO} = if(${RI}.change,null(),${RI})"

6. g.remove rast=${RI}.null,${RI}.count,${RI}.ring,${RI}.change

While you could use r.mapcalc instead of r.mfilter for steps 2 and 4,
the r.mfilter approach scales more readily to larger windows.

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