is there an available example of r.mapcalc that uses the
neighborhood modifier, so that the result of the target
cell depends on the weighted values of the 8 surrounding
cells ?
if there is could i get a copy ?
Tom Nelson
C.O.E.
Fort Worth, Texas
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@fmi.kth.se (Internet)
Department of Photogrammetry larss@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@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@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
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@fmi.kth.se (Internet)
Department of Photogrammetry larss@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@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@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