#!/bin/sh
#
# Markus Neteler
# V 1.0 - 15. Jan. 2000, 2012
# This program is free software under the GNU GPL (>=v2).
# Calculate centroid of raster area (center of gravity)

#%Module
#% description: calculates centroid of a GRASS raster area map
#% keywords: raster, statistics
#%End
#%option
#% key: map
#% type: string
#% gisprompt: old,cell,raster
#% description: Name of raster map
#% required : yes
#%End

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

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

TMP=$$

eval `g.gisenv`
: ${GISBASE?} ${GISDBASE?} ${LOCATION_NAME?} ${MAPSET?}
LOCATION=$GISDBASE/$LOCATION_NAME/$MAPSET

grep 'proj: ll' $LOCATION/../PERMANENT/PROJ_INFO > /dev/null
if [ $? -eq 0 ] ; then
 g.message -e "This modules does not work in LatLong locations"
 exit 1
fi

#### 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

name=$GIS_OPT_MAP

# example: calculate watershed (minimum size: 1000 cell units)
#  r.watershed elevation=dgm25 basin=basin threshold=1000
# now select your watershed by masking:
#  r.mapcalc MASK="if(basin == 6)"
# check it
#   d.rast basin   

# here we go for centroid calculation:
# centroid is defined as
#               N
#  x_c = 1/A * SUM (x_i * a_i)
#              i=1 
#
#               M
#  y_c = 1/A * SUM (y_i * a_i)
#              i=1
# with 
#  N: total number of cells in x direction
#  M: total number of cells in y direction
#  x_i: distance of cell center from left boundary
#  y_i: distance of cell center from upper boundary
#  a_i: area of ith cell

# calculate area in square meters:
AREA=`r.stats --q -an $name | cut -d' ' -f2`
export AREA
MORETHANONE=`echo $AREA| cut -d' ' -f2 | wc -w`
if [ $MORETHANONE -gt 1 ]
then
 echo "ERROR: more than one area in this map!"
 echo "Use r.mask to masking area of interest"
 exit
fi

# determine current resolution
eval `g.region -g`
EWRES=$ewres
NSRES=$nsres

if [ -f $LOCATION/../PERMANENT/PROJ_UNITS ] ; then
  UNITS=`cat $LOCATION/../PERMANENT/PROJ_UNITS |grep units |cut -d' ' -f2`
else
  UNITS="cellunits"
fi

# loop over areas
ALIST=`r.category $name | cut -f1`

for i in $ALIST ; do
 AREANO=$i

 echo "Area of basin $AREANO: $AREA meters^2"
 echo "Current cell resolution [$UNITS]: EW: $EWRES, NS: $NSRES"

 #set MASK to get only selected area:
 g.rename --q rast=MASK,$TMP.MASK 2> /dev/null
 r.mapcalc MASK="if($name == $AREANO)"
 # d.rast $name

 echo "Calculating x_min and x_min of area..."
 #calculate x_min
 XMIN=`r.stats --q -1gn $name |cut -d ' ' -f1 | awk 'BEGIN{min = 0.0}
 NR == 1{min = $1}
        {if ($1 < min) {min = $1}}
 END{print min}'`

 #calculate y_min
 YMIN=`r.stats --q -1gn $name |cut -d ' ' -f2 | awk 'BEGIN{min = 0.0}
 NR == 1{min = $1}
        {if ($1 < min) {min = $1}}
 END{print min}'`

 echo "Calculating centroid..."

 # calculate x_c:
 r.stats --q -1gn $name |cut -d ' ' -f1 | gawk 'BEGIN{
     sum = 0.0 ; calc = 0.0 ; xmin2 = 0.0
     ewres = '$EWRES' ; nsres = '$NSRES'
     xmin  = '$XMIN'  ; area  = '$AREA'}
  NR == 1{xmin2 = xmin * 1.0 ; ewres2 = ewres * 1.0 ; nsres2 = nsres * 1.0}
         {calc = ($1 - xmin2) * ewres2 * nsres2}
         {sum = sum + calc}
  END{printf "Center of gravity x_c: %.8f\n", sum/area + xmin2}'

 # calculate y_c:
 r.stats --q -1gn $name |cut -d ' ' -f2 | gawk 'BEGIN{
     sum = 0.0 ; calc = 0.0 ; ymin2 = 0.0
     ewres = '$EWRES' ; nsres = '$NSRES'
     ymin  = '$YMIN'  ; area  = '$AREA'}
  NR == 1{ymin2 = ymin * 1.0 }
         {calc = ($1 - ymin2) * ewres * nsres}
         {sum = sum + calc}
  END{printf "Center of gravity y_c: %.8f\n", sum/area+ymin2}'

done

#restore eventually old MASK
r.mask --q -r MASK
g.rename --q rast=$TMP.MASK,MASK 2> /dev/null

echo ""

rm -f $TMP

exit 0