[GRASS-de] Re: [GRASS-de] Rasterfunktion zum Finden des Schwerpunktes einer Fläche?

[Markus Neteler]

On Fri, Sep 06, 2002 at 12:21:29PM +0200, Till Francke wrote:
> Hallo,
> ich suche eine Rasterfunktion, die mir von attributgleichen Flächen
> jeweils den Schwerpunkt (oder einen ähnlichen representativen Punkt,
> Medianpunkt o.ä.) liefert.

Hallo Till,

ich habe mal so ein Skript geschrieben (anbei). Es sollte
hoffentlich das tun, was du brauchst.

Benutzung:
 r.centroid ammprv 22

Die Nummer gibt die Category-number der abzufragenden Rasterflaeche
an (evtl. mit d.what.rast herausfinden).

Ergebnis z.B.
 Area of basin 22: 6208000000.00 meters^2
 Current cell resolution [meters]: EW: 500, NS: 500
 MASK = (if(double(eq(ammprv[0,0],22))))
  100%
  100%
 Calculating x_min and x_min of area...
 Calculating centroid...
 Center of gravity x_c: 1663789.04
 Center of gravity y_c: 5111224.04
 
Ciao

 Markus
-------------- nächster Teil --------------
#!/bin/sh
#
# Markus Neteler
#   neteler at itc.it
# This program is free software under the GNU GPL (>=v2).
# Calculate centroid of raster area (center of gravity)

if test "$GISBASE" = ""; then
  echo "You must be in GRASS to run this program."
  exit 1
fi

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

function 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 resolution [$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)"
d.rast $MAP

echo "Calculating x_min and x_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 eventual old MASK
g.remove MASK > /dev/null
g.rename rast=$TMP.MASK,MASK 2> /dev/null