#!/bin/sh # # This program is Free Software under the GNU GPL (>=v2). # Calculate centroid of raster area (center of gravity) # Cited from the book: # Markus Neteler and Helena Mitasova: # Open Source GIS: A GRASS GIS Approach. # Kluwer Academic Publishers, Boston, Dordrecht, 464 pp, # ISBN: 1-4020-7088-8, http://www.grassbook.org if test "$GISBASE" = ""; then echo "You must be in GRASS to run this program." exit 1 fi usage() { echo " r.centroid:" echo " Calculates center of gravity (centroid) of raster area" echo " USAGE: r.centroid raster_map area_number" exit 1 } # shell check for user break (signal list: trap -l) trap "rm -f $TMPFILE > /dev/null;\ echo \"User break!\"; exit 1" 2 3 9 15 TMP=$$ TMPFILE=`g.tempfile pid=$TMP` if [ $# = 0 ] then g.ask type=old element=cell unixfile=$TMPFILE . $TMPFILE MAP="${fullname}" echo $MAP elif [ $1 = "help" -o $1 = "-help" -o $1 = "-h" ] ; then usage else MAP=$1 fi rm -f $TMPFILE #error check: if [ ! "$MAP" ] then exit 1 fi #read area number from command line: AREANO=$2 if [ ! $AREANO ] ; then echo -n "Specify area for calculations: " read AREANO fi if [ $AREANO -lt 0 ] ; then echo "ERROR: area parameter must be a number." usage exit 1 fi # 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.report -qhen $MAP u=me | sed -e 's/ //' |\ grep "^|$AREANO|." |\ cut -d'|' -f4| awk '{printf "%.2f", $1}'` #is area not present in map? if [ ! $AREA ] ; then echo "ERROR: Selected area $AREANO not found in map $MAP." exit 1 fi # determine current resolution and LOCATION units: EWRES=`awk ' /e-w/ { print $3}' $LOCATION/WIND` NSRES=`awk ' /n-s/ { print $3}' $LOCATION/WIND` if [ -f $LOCATION/../PERMANENT/PROJ_UNITS ] ; then UNITS=`cat $LOCATION/../PERMANENT/PROJ_UNITS |\ awk '/units:/ {print $2}'` else UNITS="cellunits" fi echo "Area of basin $AREANO: $AREA meters^2" echo "Current cell resol. [$UNITS]: EW: $EWRES, NS: $NSRES" #set MASK to get only selected area: g.rename rast=MASK,$TMP.MASK 2> /dev/null r.mapcalc MASK="if($MAP == $AREANO)" echo "Showing selected area..." d.rast $MAP echo "Calculating x_min and y_min of area..." #calculate x_min XMIN=`r.stats -1gnq $MAP |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 -1gnq $MAP |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 -1gnq $MAP |cut -d' ' -f1 | awk 'BEGIN{ sum = 0.0 ; calc = 0.0 ; xmin2 = 0.0 ewres = '$EWRES' ; nsres = '$NSRES' xmin = '$XMIN' ; area = '$AREA'} NR == 1{xmin2 = xmin * 1.0} {calc = ($1 - xmin2) * ewres * nsres} {sum = sum + calc} END{printf "Center of gravity x_c: %.2f\n", sum/area + xmin2}' # calculate y_c: r.stats -1gnq $MAP |cut -d' ' -f2 | awk '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: %.2f\n", sum/area+ymin2}' #restore old MASK, if present g.remove MASK > /dev/null g.rename rast=$TMP.MASK,MASK 2> /dev/null