#!/bin/sh

############################################################################
#
# MODULE:       r.reclass.area
# AUTHOR(S):    NRCS
# PURPOSE:      Reclasses a raster map greater or less than user specified area size (in hectares)
# COPYRIGHT:    (C) 1999-2008 by the GRASS Development Team
#
#               This program is free software under the GNU General Public
#               License (>=v2). Read the file COPYING that comes with GRASS
#               for details.
#
#############################################################################
# 3/2007: added label support MN
# 3/2004: added parser support MN
# 11/2001 added mapset support markus
# 2/2001 fixes markus
# 2000: updated to GRASS 5
# 1998 from NRCS, slightly modified for GRASS 4.2.1

#%Module
#%  description: Reclasses a raster map greater or less than user specified area size (in hectares).
#%  keywords: raster, statistics, aggregation
#%End

#%option
#% key: input
#% type: string
#% gisprompt: old,cell,raster
#% description: Name of input raster map
#% required : yes
#%END

#%option
#% key: output
#% type: string
#% gisprompt: new,cell,raster
#% description: Name for output raster map
#% required : yes
#%END

#%option
#% key: lesser
#% type: double
#% description: Lesser value option that sets the <= area size limit [hectares]
#% guisection: Area
#%END

#%option
#% key: greater
#% type: double
#% description: Greater value option that sets the >= area size limit [hectares]
#% guisection: Area
#%END

if  [ -z "$GISBASE" ] ; then
    echo "You must be in GRASS GIS to run this program." >&2
    exit 1
fi

if [ "$1" != "@ARGS_PARSED@" ] ; then
  exec g.parser "$0" "$@"
fi

PROG=`basename "$0"`

#### check if we have awk
if [ ! -x "`which awk`" ] ; then
    g.message -e "awk required, please install awk or gawk first"
    exit 1
fi

# setting environment, so that awk works properly in all languages
unset LC_ALL
LC_NUMERIC=C
export LC_NUMERIC

infile=`echo "$GIS_OPT_INPUT" | cut -f1 -d'@'`
outfile="$GIS_OPT_OUTPUT"

g.region -pu | head -n 1 |grep 0 > /dev/null
if [ $? -eq 0 ] ; then
 g.message -e "xy-locations are not supported"
 g.message -e "Need projected data with grids in meter"
 exit 1
fi

if [ -z "$GIS_OPT_GREATER" -a -z "$GIS_OPT_LESSER" ] ; then
  g.message -e message="You have to specify either lesser= or greater="
  exit 1
fi

if [ -n "$GIS_OPT_GREATER" -a -n "$GIS_OPT_LESSER" ] ; then
  g.message -e message="lesser= and greater= are mutually exclusive"
  exit 1
fi

if [ -n "$GIS_OPT_LESSER" ] ; then
  op=0
  limit="$GIS_OPT_LESSER"
fi
if [ -n "$GIS_OPT_GREATER" ] ; then
  op=1
  limit="$GIS_OPT_GREATER"
fi

file2="$infile.clump.$outfile"

eval `g.findfile element=cell file="$infile"`
filename="$fullname"
BASE="$name"
if [ -z "$filename" ] ; then
  g.message -e "Raster map <$infile> not found"
  exit 1
else
  infile="$filename"
fi

eval `g.findfile element=cell file="$file2"`
filename2="${fullname}"
BASE="$name"
if test "$filename2" ; then
  g.message -e "Temporal raster map <$filename2> exists"
  exit 1
else
 g.message "Generating a clumped raster map..."
 r.clump input="$infile" output="$file2"
fi

## calculation in acres
#if test "$limit" = ""; then
#echo 
#echo "Generating a reclass rules file by acres"
#r.stats -az in=$file2,$file |awk '{acre=$3 * 0.0002471; printf("%d = %.0f\n",$1,acre)}' >$infile.rules
#else
#  if test $op = 0; then
#echo 
#echo "Generating a reclass rules file by acres less than or equal to $limit"
#  r.stats -az in=$file2,$infile | awk '{limit='$limit'; acre=$3 * 0.0002471; 
#{if (acre <= limit) printf("%d = %d\n",$1,$2)}}' >$infile.rules
#  else
#echo 
#echo "Generating a reclass rules file by acres greater than or equal to $limit"
#  r.stats -az in=$file2,$infile | awk '{limit='$limit'; acre=$3 * 0.0002471; 
#{if (acre >= limit) printf("%d = %d\n",$1,$2)}}' >$infile.rules
#  fi
#fi

## calculation in hectares
#if test "$limit" = ""; then
# echo 
# echo "Generating a reclass rules file categorized by hectares"
# r.stats -an in=$file2,$infile |awk '{hectares=$3 * 0.0001; 
# printf("%d = %d %.4f\n",$1,hectares,hectares)}' > "$infile.rules"
# else
if test $op = 0; then
   g.message "Generating a reclass rules file with area size less than or equal to $limit hectares..."
   r.stats -aln in="$file2","$infile" fs='|' | \
      awk -F'|' -v LIM="$limit" \
        '{hectares=$5 * 0.0001; 
          {if (hectares <= LIM)
	     printf("%d = %d %s\n",$1,$3,$4) } }' > "$infile.rules"
else
   g.message  "Generating a reclass rules file with area size greater than or equal to $limit hectares..."
   r.stats -aln in="$file2","$infile" fs='|' | \
      awk -F'|' -v LIM="$limit" \
        '{hectares=$5 * 0.0001;
          {if (hectares >= LIM)
	     printf("%d = %d %s\n",$1,$3,$4) } }' > "$infile.rules"
fi
#fi

if test "$outfile" = ""; then
   outfile="${infile}_${limit}"
fi
g.message "Generating output raster map <$outfile>..."
cat "$infile.rules" | r.reclass i="$file2" o="$outfile.recl"
r.mapcalc "$outfile = $outfile.recl"
r.colors map=$outfile raster=$infile --quiet

g.remove rast="$outfile.recl","$file2" --quiet

#####cleanup
rm -f "$infile.rules"

exit 0