r.mapcalc and the neighborhood modifier

Lars Schylberg larss at fmi.kth.se
Thu Sep 24 12:41:34 EDT 1992


Maybe this is a little complex example but this a shell script that I
wrote a year ago when was testing different neighbourhood operators.
This script does a distance weighting on the surrounding pixels.  The
same methodolgy can be used with 8 or more neighbors.  I wrote another
version that used 8 neighbors as well, but that script is still in
Grass 3.1 version.  Well I will include that script as well in the
end.  But be aware of that the second script is a 3.1 version.

The setting of the discrimination value is a little tricky.  Please
read the original article if you want to use this script.
I have not used this script very much so I would be glad for 
comments if you find bugs.  As input you should use maps in nominal
or ordinal scale (i.e. classed data ).  Do not run this on a map
with too many categories since this will create a very big rulefile
for r.mapcalc.

Have fun !!!

Lars

Lars Schylberg                 Email: larss at fmi.kth.se    (Internet)
Department of Photogrammetry          larss at sekth.bitnet  (Bitnet)
Royal Institute of Technology  Tel.   +46 8 790 86 33     (office)
S-100 44  STOCKHOLM, SWEDEN    Fax.   +46 8 790 66 10


#########################################################################
r.prox4.sh
#########################################################################
#!/bin/sh
#
# r.prox4.sh ( version 1.0 )
# 
# Neighbourhood filter based on selecting the level of influence
# exerted on the central pixel by a predetermined set of nearest 
# neighbors, uses a proximity function derived by analogy from 
# a gravitational attractive force model
#  
# This version is based on the article: 
# I.L. Thomas (1980), Spatial Postprocessing of Spectrally
# Classified Landsat Data, Photogrammetric Engineering and
# Remote Sensing, Volume 46, No 9, pp 1201-1206
#
# This version uses a 4 neighbourhood around the center pixel as
# in the original paper.
#
# Author: Lars Schylberg  910522
#         Department of Photogrammetry, KTH, Stockholm
#
#  email: larss at fmi.kth.se 
#
#
#--------------------------------------------------------------------------
# Check if GRASS is running
#
test "$GISRC" || echo "GRASS is not running" || exit 2
#----------------------------------------------------------------------------
#
#  Evaluate arguments
#
if [ $# != 3 ]
then
    echo
    echo Usage: `basename $0` 
    echo '       input=mapname '
    echo '       output=mapname '
    echo '       dis=level of discrimination ( read the article for help ) '
    echo
    echo '  This program does a neighborhood clean up operation'
    echo '  based on a proximity function.  The discrimination value'
    echo '  determines how much that is removed'
    echo
    exit 1
fi
#
#  parse input arguments
#
for i do
	case $i in
		input=*)
			IN=`echo $i | sed s/input=//` ;;
		output=*)
			OUT=`echo $i | sed s/output=//` ;;
                dis=*)
                        DIS=`echo $i | sed s/dis=//` ;;
		*)
			echo ""
			echo "Unrecognized option: $i"
			echo 'Options: input=mapname '
                        echo '         output=mapname '
                        echo '         dis=value'
			exit 1
	esac
done
#-------------------------------------------------------------------------
#
#  Check the input arguments
#
eval `g.findfile element=cell file=$IN`
if [ ! "$file" ]
then
    echo "$IN - cell file not found"
    exit 2
fi
eval `g.findfile element=cell file=$OUT`
if [ "$file" ]
then
    echo "$OUT - exists already"
    exit 2
fi
#--------------------------------------------------------------------------
echo
echo 'Working ...'
#
nsres=`g.region -p raster=$IN | awk -f: '{ if (NR==7) print $2 }'`
ewres=`g.region -p raster=$IN | awk -f: '{ if (NR==8 ) print $2 }'`
d25=$nsres
d85=$nsres
d45=$ewres
d65=$ewres
#
r.stats -cqz $IN | awk '{print $1}' > categories
#
# Create rule for mapcalc
#
echo $OUT '= \
eval( \' > rulefile
for i in  `cat categories`
do
echo '\
q2 = eval( qt = eval( if( '$IN'[-1,0] == '$i' )) , qt * 2.0 ) ,\
q25 = eval( if( '$IN'[0,0] == '$i' && '$IN'[-1,0] == '$i', 2.0, 1.0)) ,\
q4 = eval( qt = eval( if( '$IN'[0,-1] == '$i' )) , qt * 2.0 ) ,\
q45 = eval( if( '$IN'[0,0] == '$i' && '$IN'[0,-1] == '$i', 2.0, 1.0)) ,\
q6 = eval( qt = eval( if( '$IN'[0,1] == '$i' )) , qt * 2.0 ) ,\
q65 = eval( if( '$IN'[0,0] == '$i' && '$IN'[0,1] == '$i', 2.0, 1.0)) ,\
q8 = eval( qt = eval( if( '$IN'[1,0] == '$i' )) , qt * 2.0 ) ,\
q85 = eval( if( '$IN'[0,0] == '$i' && '$IN'[1,0] == '$i', 2.0, 1.0)) ,\
f'$i' = eval(((q2 * q25)/('$d25'.0 *'$d25'.0)) + ((q4 * q45)/('$d45'.0 *'$d45'.0)) + \
((q6 * q65)/('$d65'.0 *'$d65'.0)) + ((q8 * q85)/('$d85'.0 *'$d85'.0)) ) , \' \
 >> rulefile
done
echo 'fmax = max( \'>> rulefile
for i in  `cat categories`
do
echo 'f'$i' ,\'>> rulefile
done
echo '0 ) ,\'>> rulefile
echo 'fc = 0 ,\'>> rulefile
echo 'fclass = max( \'>> rulefile
for i in  `cat categories`
do
echo 'eval( if( fmax == f'$i', '$i' )) ,\'>> rulefile
done 
echo '0), \'>> rulefile
echo 'if( fmax < '$DIS', fclass, '$IN'[0,0] )\
)
'>> rulefile
#
# Run mapcalc
#
cat rulefile | r.mapcalc
#
# set a nice colortable to the outputfile, the same as the input file
#
cp $LOCATION/colr/$IN $LOCATION/colr/$OUT
cp $LOCATION/cats/$IN $LOCATION/cats/$OUT
#
# Clean up
#
/bin/rm rulefile
/bin/rm categories


#########################################################################
Gprox8.sh
#########################################################################
#!/bin/sh
#
# Gprox8.sh ( version 1.0 )
# 
# Neighbourhood filter based on selecting the level of influence
# exerted on the central pixel by a predetermined set of nearest 
# neighbors, uses a proximity function derived by analogy from 
# a gravitational attractive force model
#  
# This version is based on the article: 
# I.L. Thomas (1980), Spatial Postprocessing of Spectrally
# Classified Landsat Data, Photogrammetric Engineering and
# Remote Sensing, Volume 46, No 9, pp 1201-1206
#
# This version uses a 8 neighbourhood around the center pixel
#
# Author: Lars Schylberg  910522
#         Department of Photogrammetry, KTH, Sweden
#
#  email: larss at fmi.kth.se
#
# Input arguments:
# $1 : filename of input cell file
# $2 : filename of output cell file
# $3 : level of discrimination ( try 0.0012 if you don't know, for 5 m cells)
#
case $# in
0) echo 'Usage: '$0' input result discrim '; exit 2;;
1) echo 'Usage: '$0' input result discrim'; exit 2;;
2) echo 'Usage: '$0' input result discrim'; exit 2;;
3) ;;
*) echo 'Usage: '$0' input result discrim '; exit 2;;
esac
#
# Check if GRASS is running
#
test "$GISRC" || echo "GRASS is not running" || exit 2
#
# Test if input file is exists
#
if test ! -f $LOCATION/cell/$1
then
    echo 'File ' $LOCATION/cell/$1 ' not found'
    exit 2
fi
#
# Test if result file is exists
#
if test -f $LOCATION/cell/$2
then
    echo 'File ' $LOCATION/cell/$2 ' exists already'
    exit 2
fi
#
echo
echo 'Working ...'
#
nsres=`Gwindow - layer=$1 print | awk -f: '{ if (NR==7) print $2 }'`
ewres=`Gwindow - layer=$1 print | awk -f: '{ if (NR==8 ) print $2 }'`
d25=$nsres
d85=$nsres
d45=$ewres
d65=$ewres
d15=`echo $ewres $nsres | awk '{printf "%.2f\n",exp(log(($1*$1)+($2*$2))*0.5)}'`
d35=$d15
d75=$d15
d95=$d15
dis=$3
#
# Create rule for mapcalc
#
echo $2 '= \
eval( \' > rulefile
for i in  `Gstats -cz $1 | awk -f: '{print $1}'`
do
echo '\
q5 = eval( qt = eval( if( '$1'[0,0] == '$i' )) , qt * 2 ) ,\
q2 = eval( qt = eval( if( '$1'[-1,0] == '$i' )) , qt * 2 ) ,\
q4 = eval( qt = eval( if( '$1'[0,-1] == '$i' )) , qt * 2 ) ,\
q6 = eval( qt = eval( if( '$1'[0,1] == '$i' )) , qt * 2 ) ,\
q8 = eval( qt = eval( if( '$1'[1,0] == '$i' )) , qt * 2 ) ,\
q1 = eval( qt = eval( if( '$1'[-1,-1] == '$i' )) , qt * 2 ) ,\
q3 = eval( qt = eval( if( '$1'[-1,1] == '$i' )) , qt * 2 ) ,\
q7 = eval( qt = eval( if( '$1'[1,-1] == '$i' )) , qt * 2 ) ,\
q9 = eval( qt = eval( if( '$1'[1,1] == '$i' )) , qt * 2 ) ,\
f'$i' = eval( ((q2 * q5)/('$d25' *'$d25')) + ((q4 * q5)/('$d45' *'$d45')) + \
((q6 * q5)/('$d65' *'$d65')) + ((q8 * q5)/('$d85' *'$d85')) + \
((q1 * q5)/('$d15' * '$d15')) + ((q3 * q5)/('$d35' * '$d35')) + \
((q7 * q5)/('$d75' * '$d75')) + ((q9 * q5)/('$d95' * '$d95')) ) , \' \
 >> rulefile
done
echo 'fmax = max( \'>> rulefile
for i in  `Gstats -cz $1 | awk -f: '{print $1}'`
do
echo 'f'$i' ,\'>> rulefile
done
echo '0 ) ,\'>> rulefile
echo 'fc = 0 ,\'>> rulefile
echo 'fclass = max( \'>> rulefile
for i in  `Gstats -cz $1 | awk -f: '{print $1}'`
do
echo 'eval( if( fmax == f'$i', '$i' )) ,\'>> rulefile
done 
echo '0), \'>> rulefile
echo 'if( fmax > '$dis', fclass ) \
)
'>> rulefile
#
# Run mapcalc
#
Gmapcalc2 < rulefile
#
# set a nice colortable to the outputfile, the same as the input file
#
eval `gisenv`
LOCATION=$GISDBASE/$LOCATION_NAME/$MAPSET
cp $LOCATION/colr/$1 $LOCATION/colr/$2
cp $LOCATION/cats/$1 $LOCATION/cats/$2
#
# Clean up
#
#/bin/rm rulefile



More information about the grass-user mailing list