[GRASS-SVN] r53159 - in grass-addons/grass6/raster: . r.cog
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Sep 10 22:32:16 PDT 2012
Author: hamish
Date: 2012-09-10 22:32:16 -0700 (Mon, 10 Sep 2012)
New Revision: 53159
Added:
grass-addons/grass6/raster/r.cog/
grass-addons/grass6/raster/r.cog/r.cog
Removed:
grass-addons/grass6/raster/r.cog/v.points.cog
Modified:
grass-addons/grass6/raster/r.cog/Makefile
grass-addons/grass6/raster/r.cog/description.html
Log:
+ new addon module: r.cog for center of gravity and trending
Modified: grass-addons/grass6/raster/r.cog/Makefile
===================================================================
--- grass-addons/grass6/vector/v.points.cog/Makefile 2012-09-10 23:47:11 UTC (rev 53158)
+++ grass-addons/grass6/raster/r.cog/Makefile 2012-09-11 05:32:16 UTC (rev 53159)
@@ -1,6 +1,6 @@
MODULE_TOPDIR = ../..
-PGM = v.points.cog
+PGM = r.cog
include $(MODULE_TOPDIR)/include/Make/Script.make
Modified: grass-addons/grass6/raster/r.cog/description.html
===================================================================
--- grass-addons/grass6/vector/v.points.cog/description.html 2012-09-10 23:47:11 UTC (rev 53158)
+++ grass-addons/grass6/raster/r.cog/description.html 2012-09-11 05:32:16 UTC (rev 53159)
@@ -1,50 +1,67 @@
<h2>DESCRIPTION</h2>
-<em>v.points.cog</em> condenses points or centroids sharing
-a common attribute into a single point in a new vector map at their
-average position (center of gravity).
+<em>r.cog</em> finds the center of gravity of a raster map, with respect
+to occupied cells in the x,y plane, and the mean z-value of the raster
+data. By default it prints out the <tt>easting,northing</tt> of the center
+cell and the mean of the raster map to give the "pivot point". If the <b>-a</b>
+flag is given it will also print out the generalized azimuth and dip angle
+of the surface. Angles are measured in degrees; azimuth is measured CCW from
+the +x axis, and dip angle is measured downwards along the direction of the
+azimuth.
+
<p>
-For this to work well your clusters of points must be gregarious
-(Gaussian distribution) - if two groups of points habitate in two
-corners of the map the output point will fall in the center and
-match niether well.
+These outputs can be printed in shell script format suitable for parsing by
+using the <b>-g</b> flag. In that case the order of terms is
+<tt>easting|northing|elevation|azimuth|dip</tt>.
+
<p>
-If needed you can use <em>v.digit</em> to adjust point positions
-created by this module, <!-- repel 2 points if v.distance < minimum? -->
-or use this module as a preprocessing step before running
-<em>v.label.sa</em> to deal with overlapping labels automatically.
+If the <b>plane</b> option is used, a new raster map is created in the plane
+of the overall azimuth and dip trends.
+<p>
+For a categorical maps you can optionally select the categories to consider.
+This mode is not currently compatible with the angle output or plane output
+options.
+
+<p>
+Optionally it will draw a marker on the Xmonitor using the <em>d.mark</em>
+module from GRASS addons (d.shortcuts).
+
+
<!--<p> appropriate for lat/lon? (weighted avg for y by cos(y)?) -->
+<!--
<h2>EXAMPLE</h2>
-Create single points at the average of all points in the
-<tt>bugsites</tt> map, and place a single label there.
+TODO
<div class="code"><pre>
-v.points.cog in=bugsites out=bug_cog column=str1
-d.vect -c bugsites color=none icon=basic/circle
-
-d.vect bug_cog disp=attr attrcol=str1 lcolor=black \
- lsize=14 xref=center yref=center bgcolor=white
</pre></div>
+-->
+<h2>BUGS</h2>
+<i>The dip angle is currently wrong.</i>
+<!-- see comments in the code -->
+
+
<h2>SEE ALSO</h2>
<em>
-<a href="v.label.sa.html">v.label.sa</a><br>
-<a href="v.label.html">v.label</a><br>
-<a href="v.digit.html">v.digit</a>
+<a href="v.points.cog.html">v.points.cog</a> (addon module)<br>
+<a href="d.mark.html">d.mark</a> (addon module)<br>
+<a href="r.param.scale.html">r.param.scale</a><br>
+<a href="r.plane.html">r.plane</a><br>
+<a href="r.univar.html">r.univar</a>
</em>
<h2>AUTHOR</h2>
Hamish Bowman<br>
-<i>Dept. Marine Science<br>
+<i>Dept. of Geology<br>
University of Otago<br>
Dunedin, New Zealand</i>
<br>
Copied: grass-addons/grass6/raster/r.cog/r.cog (from rev 53158, grass-addons/grass6/vector/v.points.cog/v.points.cog)
===================================================================
--- grass-addons/grass6/raster/r.cog/r.cog (rev 0)
+++ grass-addons/grass6/raster/r.cog/r.cog 2012-09-11 05:32:16 UTC (rev 53159)
@@ -0,0 +1,226 @@
+#!/bin/sh
+#
+############################################################################
+#
+# MODULE: r.cog
+#
+# AUTHOR(S): Hamish Bowman, Dunedin, New Zealand
+#
+# PURPOSE: Find the center of gravity of a raster map, by area cover
+#
+# COPYRIGHT: (c) 2012 Hamish Bowman, and 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.
+#
+#############################################################################
+#%Module
+#% description: Find the center of gravity of a raster map, by area cover.
+#% keywords: raster, cluster
+#%End
+#%option
+#% key: map
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Name of input raster map
+#% required : yes
+#%end
+#%option
+#% key: cat
+#% label: Specify category value(s) to work on
+#% description: Example: 1,3,7-9,13
+#%end
+#%option
+#% key: plane
+#% description: Name for optional planar trend map at the center of gravity
+#% type: string
+#% gisprompt: new,cell,raster
+#%end
+#%flag
+#% key: d
+#% description: Draw a marker on the point using the d.mark addon module
+#%end
+#%flag
+#% key: a
+#% description: Display map trend (average azimuth and dip angle)
+#%end
+#%flag
+#% key: g
+#% description: Print in shell script style
+#%end
+
+##%flag
+##% key: m
+##% description: Use median instead of mean for z-center
+##%end
+## or trimmed mean?
+
+
+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
+
+MAP="$GIS_OPT_MAP"
+CAT=`echo "$GIS_OPT_CAT" | sed -e 's/,/ /g' -e 's/-/ thru /'`
+#g.message -d message="cats = [$CAT]"
+
+# check for input map
+eval `g.findfile element=cell file="$MAP"`
+if [ ! "$file" ] ; then
+ g.message -e "Raster map <$MAP> does not exist."
+ exit 1
+fi
+
+
+if [ -n "$CAT" ] ; then
+ # r.reclass flattens it to all be = 1, so no slope or aspect
+ # using r.mapcalc to extract the cats could fix that, but parsing would be a pain
+ if [ -n "$GIS_OPT_PLANE" ] ; then
+ g.message -e "The cat and plane options are not (currently) compatible"
+ exit 1
+ fi
+ if [ "$GIS_FLAG_A" -eq 1 ] ; then
+ g.message -e "The cat and angle options are not (currently) compatible"
+ exit 1
+ fi
+fi
+
+
+cleanup()
+{
+ if [ -n "$CAT" ] ; then
+ g.remove --quiet rast="$TMP_MAP"
+ fi
+ if [ "$GIS_FLAG_A" -eq 1 ] ; then
+ true; # g.remove --quiet rast="$TMP_MAP.slope,$TMP_MAP.dx,$TMP_MAP.dy"
+ fi
+}
+trap "cleanup" 2 3 15
+
+
+# -p needs -a
+if [ -n "$GIS_OPT_PLANE" ] ; then
+ if [ "$GIS_FLAG_A" -eq 0 ] ; then
+ GIS_FLAG_A=1
+ fi
+fi
+
+
+TMP_MAP="tmp.rcog.$$"
+if [ -n "$CAT" ] ; then
+ r.reclass in="$MAP" out="$TMP_MAP" rules=- << EOF
+$CAT = 1
+* = NULL
+EOF
+ MAP="$TMP_MAP"
+fi
+
+COG_XY=`r.out.xyz in="$MAP" out=- --quiet | cut -f1,2 -d'|' | \
+ awk -F'|' \
+ 'BEGIN { sum_x=0.0; sum_y=0.0; i=0 }
+ { sum_x += $1; sum_y += $2; i++ }
+ END { if(i>0) { printf("%.15g|%.15g\n", sum_x/i, sum_y/i) } }'`
+
+if [ "$GIS_FLAG_D" -eq 1 ] ; then
+ if [ ! -x `which d.mark` ] ; then
+ g.message -e "d.mark addon module not found. Please install."
+ exit 1
+ fi
+ d.mark -m at=`echo "$COG_XY" | tr '|' ','`
+fi
+
+if [ -n "$CAT" ] ; then
+ #d.rast $TMP_MAP
+ eval `r.univar -g map="$GIS_OPT_MAP" zones="$TMP_MAP" | grep '^me'`
+else
+ eval `r.univar -g map="$MAP" | grep '^me'`
+fi
+
+#maybe use median or trimmed mean? += 'r.univar -e'
+
+COG_Z="$mean"
+#g.message -d message="mean=$mean median=$median"
+#g.message -d message=" "
+
+EAST=`echo "$COG_XY" | cut -f1 -d'|'`
+NORTH=`echo "$COG_XY" | cut -f2 -d'|'`
+if [ "$GIS_FLAG_G" -eq 0 ] ; then
+ echo "Center of gravity: $EAST,$NORTH"
+ echo "Average elevation: $COG_Z"
+else
+ echo -n "$COG_XY|$COG_Z"
+fi
+
+
+if [ "$GIS_FLAG_A" -eq 0 ] ; then
+ if [ "$GIS_FLAG_G" -eq 1 ] ; then
+ echo
+ fi
+ cleanup
+ exit 0
+fi
+
+
+# calc mean slope and dip angle
+r.slope.aspect elev="$MAP" \
+ slope="$TMP_MAP.slope" \
+ dx="$TMP_MAP.dx" \
+ dy="$TMP_MAP.dy" \
+ prec=double --quiet
+
+## FIXME:
+## {mean || median} slope is the wrong thing to use here. e.g. if you had some
+## terrain riddled with canyons and pinacles, but no overall gradient, the
+## mean slope would be quite high, but the overall region gradient would be
+## near zero.
+## If all else fails we could use a 2D least-squares on the DEM, or perhaps run
+## v.surf.rst with very high smoothing and low tension ...?
+#
+# maybe the r.univar mean or median of the magnitude of the axes' 1st derivatives:
+# r.mapcalc "tmp.magnitude = sqrt(tmp.rcog.$$.dx*tmp.rcog.$$.dx \
+# + tmp.rcog.$$.dy*tmp.rcog.$$.dy)"
+# r.univar -g tmp.magnitude
+# ??
+
+
+DIP=`r.univar -eg map="$TMP_MAP.slope" | grep '^median=' | cut -f2 -d'='`
+DX=`r.univar -g map="$TMP_MAP.dx" | grep '^mean=' | cut -f2 -d'='`
+DY=`r.univar -g map="$TMP_MAP.dy" | grep '^mean=' | cut -f2 -d'='`
+
+#r.univar map="$TMP_MAP.slope" -e
+#d.rast "$TMP_MAP"
+PI=3.14159265358979323846
+AZ=`echo "$DX $DY" | awk -v PI="$PI" '{print atan2($2, $1) * (180/PI)}'`
+
+if [ "$GIS_FLAG_G" -eq 0 ] ; then
+ echo "Azimuth: $AZ"
+ echo "Dip angle: $DIP"
+else
+ echo "|$AZ|$DIP"
+fi
+
+#g.message -d message=". slope=$DIP dx=$DX dy=$DY"
+
+
+if [ -z "$GIS_OPT_PLANE" ] ; then
+ cleanup
+ exit 0
+fi
+
+
+#### create a plane
+
+# note we need to add 270 since r.plane is wierd
+theta=`echo "$AZ" | awk '{print ($1 + 270.) % 360}'`
+
+r.plane name="$GIS_OPT_PLANE" dip="$DIP" azimuth="$theta" \
+ east="$EAST" north="$NORTH" elev="$COG_Z" type=double --quiet
+#r.colors "$GIS_OPT_PLANE" rast="$MAP"
+
+cleanup
Deleted: grass-addons/grass6/raster/r.cog/v.points.cog
===================================================================
--- grass-addons/grass6/vector/v.points.cog/v.points.cog 2012-09-10 23:47:11 UTC (rev 53158)
+++ grass-addons/grass6/raster/r.cog/v.points.cog 2012-09-11 05:32:16 UTC (rev 53159)
@@ -1,133 +0,0 @@
-#!/bin/sh
-#
-############################################################################
-#
-# MODULE: v.points.cog
-#
-# AUTHOR(S): Hamish Bowman
-#
-# PURPOSE: Condense points or centroids sharing a common attribute into a single point
-#
-# COPYRIGHT: (c) 2010 Hamish Bowman, and 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.
-#
-#############################################################################
-#%Module
-#% description: Condense points or centroids sharing a common attribute into a single point.
-#% keywords: vector, cluster
-#%End
-#%option
-#% key: input
-#% type: string
-#% gisprompt: old,vector,vector
-#% description: Name of input vector map
-#% required : yes
-#%end
-#%option
-#% key: output
-#% type: string
-#% gisprompt: new,vector,vector
-#% description: Name for output vector map
-#% required : yes
-#%end
-#%option
-#% key: column
-#% type: string
-#% description: Column containing common attribute
-#% required : yes
-#%end
-#%option
-#% key: layer
-#% type: integer
-#% answer: 1
-#% description: Layer number
-#% required: no
-#%end
-
-##%option
-##% key: type
-##% type: string
-##% description: Feature type(s)
-##% options: point,centroid
-##% answer: point
-##% required: no
-##% multiple: 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
-
-MAP="$GIS_OPT_INPUT"
-COLUMN="$GIS_OPT_COLUMN"
-LAYER="$GIS_OPT_LAYER"
-
-# check for input map
-eval `g.findfile element=vector file="$MAP"`
-if [ ! "$file" ] ; then
- g.message -e "Vector map <$MAP> does not exist."
- exit 1
-fi
-
-# check for column
-if [ `v.info -c "$MAP" layer="$LAYER" --quiet | cut -f2 -d'|' | grep -c "^$COLUMN$"` -ne 1 ] ; then
- g.message -e "Column <$COLUMN> not found."
- exit 1
-fi
-
-# get column details so we can recreate it
-DB=`v.db.connect "$MAP" -g layer="$LAYER" fs='|' | cut -f4 -d'|'`
-COLUMN_DESC=`db.describe -c table="$MAP" database="$DB" | grep " $COLUMN:" | cut -f3- -d:`
-if [ `echo "$COLUMN_DESC" | grep -c CHARACTER` -eq 1 ] ; then
- COLUMN_TYPE="string"
- COLUMN_LEN=`echo "$COLUMN_DESC" | cut -f2 -d:`
- COLUMN_DEFN="varchar($COLUMN_LEN)"
-else
- COLUMN_TYPE="number"
- COLUMN_DEFN=`echo "$COLUMN_DESC" | cut -f1 -d:`
-fi
-
-# cheap hack to avoid conflict
-if [ $COLUMN = "cat" ] ; then
- OUT_COLUMN="cat_"
-else
- OUT_COLUMN="$COLUMN"
-fi
-
-
-(
-IFS='|'
-for ITEM in `v.db.select "$MAP" -c column="$COLUMN" layer="$LAYER" | sort | uniq | tr '\n' '|'` ; do
- #echo "[$ITEM]"
- if [ "$COLUMN_TYPE" = "string" ] ; then
- WHERE_STR="$COLUMN = '$ITEM'"
- else
- WHERE_STR="$COLUMN = $ITEM"
- fi
-
- v.out.ascii "$MAP" column="$COLUMN" layer="$LAYER" where="$WHERE_STR" | \
- awk -F'|' \
- 'BEGIN { sum_x=0; sum_y=0; i=0 }
- { sum_x += $1; sum_y += $2; i++ }
- END { if(i>0) { printf("%.15g|%.15g|%s\n", sum_x/i, sum_y/i, $4) } }'
-done
-) | v.in.ascii out="$GIS_OPT_OUTPUT" columns="x double, y double, $OUT_COLUMN $COLUMN_DEFN"
-#echo columns="x double, y double, $OUT_COLUMN $COLUMN_DEFN"
-
-retval=$?
-
-# cleanup cheap hack
-if [ $COLUMN = "cat" ] ; then
- v.db.dropcol map="$GIS_OPT_OUTPUT" column="cat_" --quiet
-fi
-
-exit $retval
More information about the grass-commit
mailing list