Selaa lähdekoodia

bring in v.colors, r.colors.stddev from addons (by way of devbr6)
-still shell scrips-

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@34706 15284696-431f-4ddb-bdfa-cd5b030d7da7

Hamish Bowman 16 vuotta sitten
vanhempi
commit
df68e7540f

+ 2 - 0
scripts/Makefile

@@ -22,6 +22,7 @@ SUBDIRS = \
 	m.proj \
 	r.blend \
 	r.buffer \
+	r.colors.stddev \
 	r.fillnulls \
 	r.grow \
 	r.in.aster \
@@ -36,6 +37,7 @@ SUBDIRS = \
 	r.tileset \
 	v.build.all \
 	v.centroids \
+	v.colors \
 	v.convert.all \
 	v.db.addcol \
 	v.db.addtable \

+ 7 - 0
scripts/r.colors.stddev/Makefile

@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../..
+
+PGM = r.colors.stddev
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script

+ 38 - 0
scripts/r.colors.stddev/description.html

@@ -0,0 +1,38 @@
+<H2>DESCRIPTION</H2>
+
+<EM>r.colors.stddev</EM> set raster map color rules based on standard
+deviations from a map's mean value, either as a continuous color gradient
+or in color bands per standard deviation from the mean.
+
+<p>
+With the color band option values less that 1 S.D. from the mean are
+colored green, within 1-2 S.D. are colored yellow, within 2-3 S.D. are
+colored red, and beyond 3 S.D. are colored black.
+<p>
+For a differences map there is an option to lock the center of the color
+table at zero. Values more than two S.D. below the mean will be colored blue;
+values below the mean but less than 2 S.D. away will transition to white,
+and above the mean the colors will simularly transition to full red at +2 S.D.
+
+
+<!--
+<h2>EXAMPLES</h2>
+-->
+
+
+<H2>SEE ALSO</H2>
+
+<EM>
+<A HREF="r.colors.html">r.colors</A><BR>
+<A HREF="r.univar.html">r.univar</A><BR>
+<A HREF="v.colors.html">v.colors</A><BR>
+</EM>
+
+
+<H2>AUTHOR</H2>
+
+Hamish Bowman<BR>
+<i>Dunedin, New Zealand</i>
+
+<p>
+<i>Last changed: $Date$</i>

+ 145 - 0
scripts/r.colors.stddev/r.colors.stddev

@@ -0,0 +1,145 @@
+#!/bin/sh
+############################################################################
+#
+# MODULE:       r.colors.stddev
+# AUTHOR:       M. Hamish Bowman, Dept. Marine Science, Otago Univeristy,
+#                 New Zealand
+# PURPOSE:      Set color rules based on stddev from a map's mean value.
+#
+# COPYRIGHT:    (c) 2007 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: Set color rules based on stddev from a map's mean value.
+#% keywords: raster
+#%End
+#% option
+#% key: input
+#% type: string
+#% gisprompt: old,cell,raster
+#% key_desc: name
+#% description: Name of input raster map 
+#% required: yes
+#%end
+#%flag
+#% key: b
+#% description: Color using standard deviation bands
+#%end
+#%flag
+#% key: z
+#% description: Force center at zero
+#%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
+
+
+eval `r.univar -g "$GIS_OPT_INPUT"`
+# $? is result of the eval not r.univar (???)
+#if [ $? -ne 0 ] ; then
+#   g.region -e "Problem running r.univar"
+#   exit 1
+#fi
+
+
+if [ $GIS_FLAG_Z -eq 0 ] ; then
+
+   MEAN_MINUS_2STDEV=`echo "$mean $stddev" | awk '{print $1 - 2*$2}'`
+   MEAN_PLUS_2STDEV=`echo "$mean $stddev" | awk '{print $1 + 2*$2}'`
+
+   if [ $GIS_FLAG_B -eq 0 ] ; then
+     # smooth free floating blue/white/red
+     r.colors "$GIS_OPT_INPUT" color=rules << EOF
+       0% blue
+       $MEAN_MINUS_2STDEV blue
+       $mean white
+       $MEAN_PLUS_2STDEV red
+       100% red
+EOF
+   else
+     # banded free floating  black/red/yellow/green/yellow/red/black
+     MEAN_MINUS_1STDEV=`echo "$mean $stddev" | awk '{print $1 - $2}'`
+     MEAN_MINUS_3STDEV=`echo "$mean $stddev" | awk '{print $1 - 3*$2}'`
+     MEAN_PLUS_1STDEV=`echo "$mean $stddev" | awk '{print $1 + $2}'`
+     MEAN_PLUS_3STDEV=`echo "$mean $stddev" | awk '{print $1 + 3*$2}'`
+
+     # reclass with labels only works for category (integer) based maps
+     #r.reclass input="$GIS_OPT_INPUT" output="${GIS_OPT_INPUT}.stdevs" << EOF
+
+     # >3 S.D. outliers colored black so they show up in d.histogram w/ white background
+     r.colors "$GIS_OPT_INPUT" color=rules << EOF
+       0% black
+       $MEAN_MINUS_3STDEV black
+       $MEAN_MINUS_3STDEV red
+       $MEAN_MINUS_2STDEV red
+       $MEAN_MINUS_2STDEV yellow
+       $MEAN_MINUS_1STDEV yellow
+       $MEAN_MINUS_1STDEV green
+       $MEAN_PLUS_1STDEV green
+       $MEAN_PLUS_1STDEV yellow
+       $MEAN_PLUS_2STDEV yellow
+       $MEAN_PLUS_2STDEV red
+       $MEAN_PLUS_3STDEV red
+       $MEAN_PLUS_3STDEV black
+       100% black
+EOF
+   fi
+
+
+else 
+   # data centered on 0  (e.g. map of deviations)
+   r.mapcalc "r_col_stdev_abs_$$ = abs($GIS_OPT_INPUT)"
+   eval `r.info -r "r_col_stdev_abs_$$"`
+
+   # current r.univar truncates percentage to the base integer
+   STDDEV2=`r.univar -eg "r_col_stdev_abs_$$" perc=95.4500 | grep ^percentile | cut -f2 -d=`
+
+   if [ $GIS_FLAG_B -eq 0 ] ; then
+     # zero centered smooth blue/white/red
+     r.colors "$GIS_OPT_INPUT" color=rules << EOF
+       -$max blue
+       -$STDDEV2 blue
+       0 white
+       $STDDEV2 red
+       $max red
+EOF
+   else
+     # zero centered banded  black/red/yellow/green/yellow/red/black
+
+     # current r.univar truncates percentage to the base integer
+     STDDEV1=`r.univar -eg "r_col_stdev_abs_$$" perc=68.2689 | grep ^percentile | cut -f2 -d=`
+     STDDEV3=`r.univar -eg "r_col_stdev_abs_$$" perc=99.7300 | grep ^percentile | cut -f2 -d=`
+
+     # >3 S.D. outliers colored black so they show up in d.histogram w/ white background
+     r.colors "$GIS_OPT_INPUT" color=rules << EOF
+       -$max black
+       -$STDDEV3 black
+       -$STDDEV3 red
+       -$STDDEV2 red
+       -$STDDEV2 yellow
+       -$STDDEV1 yellow
+       -$STDDEV1 green
+       $STDDEV1 green
+       $STDDEV1 yellow
+       $STDDEV2 yellow
+       $STDDEV2 red
+       $STDDEV3 red
+       $STDDEV3 black
+       $max black
+EOF
+   fi
+
+   g.remove rast="r_col_stdev_abs_$$" --quiet
+fi
+

+ 7 - 0
scripts/v.colors/Makefile

@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../..
+
+PGM = v.colors
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script

+ 80 - 0
scripts/v.colors/description.html

@@ -0,0 +1,80 @@
+<H2>DESCRIPTION</H2>
+
+<EM>v.colors</EM> is much like <em>r.colors</em>, but may be used for vector maps.
+You give it a vector map and numeric data column, together with color rules
+like you would do for a raster. It creates a new column in the database with
+R:G:B values suitable for use with '<tt>d.vect -a</tt>'.
+
+<P>
+How it works: it creates a dummy raster map with the same data range as
+the vector's column then runs <em>r.colors</em> for that temporary map.
+It then uses <em>r.what.colors</em> for each value found by <em>v.db.select</em>
+and uploads it a new column in the vector map's attribute database.
+
+<P>
+It is planned that this script will be replaced with a C display module which
+renders thematic vector maps directly instead requiring the overhead of
+saving the colors into the DB.
+
+<P>
+If the target column name given by the <b>rgb_column</b> option does
+not exist, it will be created. The default name is "<tt>GRASSRGB</tt>".
+
+<h2>EXAMPLES</h2>
+
+Create a random sample point map, query raster map values for those points,
+and colorize output.
+
+<div class="code"><pre>
+# Spearfish dataset
+g.region -d
+v.random out=rand5k_elev n=5000
+v.db.addtable map=rand5k_elev column='elevation double precision'
+v.what.rast vector=rand5k_elev raster=elevation.10m column=elevation
+v.colors map=rand5k_elev column=elevation color=bcyr
+d.mon x0
+d.vect -a rand5k_elev
+</pre></div>
+
+
+<P>
+
+Colorizing a TIN generated by <em>v.delaunay</em>:
+
+<div class="code"><pre>
+# new columns for x,y,z of centroids
+v.db.addtable map=tin \
+   columns="east double precision, north double precision, height double precision, GRASSRGB varchar(11)"
+
+# transfer geometry for colorizing (we need the centroid height)
+v.to.db tin opt=coor col="east,north,height"
+v.db.select tin
+
+v.colors tin column=height rgb_column=GRASSRGB color=rainbow
+
+# display colorized triangles
+d.mon x0
+d.vect -a tin
+</pre></div>
+
+
+<H2>SEE ALSO</H2>
+
+<EM>
+<A HREF="r.colors.html">r.colors</A><BR>
+<A HREF="r.colors.stddev.html">r.colors.stddev</A><BR>
+<A HREF="r.what.color.html">r.what.color</A><BR>
+<A HREF="v.db.addcol">v.db.addcol</A><BR>
+<A HREF="v.db.select.html">v.db.select</A><BR>
+<A HREF="db.execute.html">db.execute</A>
+</EM>
+
+
+
+<H2>AUTHOR</H2>
+
+Hamish Bowman<BR>
+<i>Dunedin, New Zealand</i>
+
+<p>
+<i>Last changed: $Date$</i>

+ 335 - 0
scripts/v.colors/v.colors

@@ -0,0 +1,335 @@
+#!/bin/sh
+############################################################################
+#
+# MODULE:       v.colors
+# AUTHOR:       M. Hamish Bowman, Dept. Marine Science, Otago Univeristy,
+#                 New Zealand
+# PURPOSE:      Populate a GRASSRGB column with a color map and data column
+#		Helper script for thematic mapping tasks
+#
+# COPYRIGHT:    (c) 2008 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: Set color rules for features in a vector using a numeric attribute column.
+#% keywords: vector
+#%End
+#% option
+#% key: map
+#% type: string
+#% gisprompt: old,vector,vector
+#% key_desc: name
+#% description: Name of input vector map 
+#% required: yes
+#%end
+#%option
+#% key: column
+#% type: string
+#% description: Name of column containg numeric data
+#% required : yes
+#%end
+#%option
+#% key: layer
+#% type: integer
+#% description: Layer number of data column
+#% answer: 1
+#% required: no
+#%end
+#%Option
+#% key: rgb_column
+#% type: string
+#% required: no
+#% description: Name of color column to populate
+#% answer: GRASSRGB
+#% guisection: Colors
+#%End
+#% option
+#% key: range
+#% type: double
+#% required: no
+#% multiple: no
+#% key_desc: min,max
+#% description: Manually set range (min,max)
+#%End
+#% option
+#% key: color
+#% type: string
+#% key_desc: style
+#% description: Type of color table
+#% required: no
+#% guisection: Colors
+#%end
+#%Option
+#% key: raster
+#% type: string
+#% required: no
+#% description: Raster map name from which to copy color table
+#% gisprompt: old,cell,raster
+#% guisection: Colors
+#%End
+#%Option
+#% key: rules
+#% type: string
+#% required: no
+#% description: Path to rules file
+#% gisprompt: old_file,file,input
+#% guisection: Colors
+#%End
+#%Flag
+#% key: s
+#% description: Save placeholder raster map for use with d.legend
+#%End
+#%Flag
+#% key: n
+#% description: Invert colors
+#% guisection: Colors
+#%End
+
+
+## TODO: implement -e (equalized) and -g (logarithmic) methods in r.colors
+##   'v.db.select column= | wc -l' to set region size (1xLength)
+##   then create r.in.ascii 1xLength matrix with data (WITHOUT uniq)
+##   and run r.colors on that raster map.
+##      or
+##   v.to.rast, r.colors -g, then parse colr/ file. But that's resolution dependent
+
+
+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
+
+
+GRASS_VERBOSE=0
+export GRASS_VERBOSE
+
+cleanup()
+{
+   if [ -n "$TMP" ] ; then
+      \rm "$TMP" "${TMP}_vcol.sql"
+   fi
+   eval `g.findfile elem=windows file="tmp_vcolors.$$" | grep '^name='`
+   if [ -n "$name" ] ; then
+      unset WIND_OVERRIDE
+      g.remove region="tmp_vcolors.$$"
+   fi
+   eval `g.findfile elem=cell file="tmp_colr_$$" | grep '^name='`
+   if [ -n "$name" ] ; then
+      g.remove rast="tmp_colr_$$"
+   fi
+}
+trap "cleanup" 2 3 15
+
+
+### setup enviro vars ###
+eval `g.gisenv`
+: ${GISBASE?} ${GISDBASE?} ${LOCATION_NAME?} ${MAPSET?}
+
+#### does map exist in CURRENT mapset?
+eval `g.findfile element=vector file="$GIS_OPT_MAP" mapset="$MAPSET"`
+VECT_MAPSET=`echo "$GIS_OPT_MAP" | grep '@' | cut -f2 -d'@'`
+if [ -z "$VECT_MAPSET" ] ; then
+   VECT_MAPSET="$MAPSET"
+fi
+if [ ! "$file" ] || [ "$VECT_MAPSET" != "$MAPSET" ] ; then
+   g.message -e "Vector map <$GIS_OPT_MAP> not found in current mapset"
+   exit 1
+else
+   VECTOR=`echo "$GIS_OPT_MAP" | cut -f1 -d'@'`
+fi
+
+
+
+# sanity check mutually exclusive color options
+CTEST=0
+for OPT in COLOR RASTER RULES ; do
+   eval "OPTNAME=\$GIS_OPT_${OPT}"
+   if [ -n "$OPTNAME" ] ; then
+      CTEST=`expr $CTEST + 1`
+   fi
+done
+if [ "$CTEST" -ne 1 ] ; then
+    g.message -e "Pick one of color, rules, or raster options"
+    exit 1
+fi
+
+if [ -n "$GIS_OPT_COLOR" ] ; then
+   #### check the color rule is valid
+   COLOR_OPTS=`ls $GISBASE/etc/colors/`
+   COLOR_OPTS="`echo $COLOR_OPTS | tr '\n' ' '` random grey.eq grey.log rules"
+
+   HAVE_COLR=0
+   for COLR in $COLOR_OPTS ; do
+       if [ "$COLR" = "$GIS_OPT_COLOR" ] ; then
+          HAVE_COLR=1
+          break
+       fi
+   done
+   if [ "$HAVE_COLR" -eq 0 ] ; then
+       g.message -e "Invalid color rule <$GIS_OPT_COLOR>"
+       echo "       Valid options are: $COLOR_OPTS"
+       exit 1
+   fi
+elif [ -n "$GIS_OPT_RASTER" ] ; then
+   r.info -r "$GIS_OPT_RASTER" > /dev/null
+   if [ $? -ne 0 ] ; then
+       g.message -e "Unable to open raster map <$GIS_OPT_RASTER>"
+       exit 1
+   fi
+elif [ -n "$GIS_OPT_RULES" ] ; then
+   if [ ! -r "$GIS_OPT_RULES" ] ; then
+       g.message -e "Unable to read color rules file <$GIS_OPT_RULES>"
+       exit 1
+   fi
+fi
+
+
+#### column checks
+# check input data column 
+NCOLUMN_TYPE=`v.info -c map="$GIS_OPT_MAP" layer="$GIS_OPT_LAYER" | grep -i "|$GIS_OPT_COLUMN$" | cut -f1 -d'|'`
+if [ -z "$NCOLUMN_TYPE" ] ; then
+    g.message -e "Column <$GIS_OPT_COLUMN> not found"
+    exit 1
+elif [ "$NCOLUMN_TYPE" != "INTEGER" ] && [ "$NCOLUMN_TYPE" != "DOUBLE PRECISION" ] ; then
+    g.message -e "Column <$GIS_OPT_COLUMN> is not numeric"
+    exit 1
+fi
+#echo "column <$GIS_OPT_COLUMN> is type [$NCOLUMN_TYPE]"
+
+# check if GRASSRGB column exists, make it if it doesn't
+TABLE=`v.db.connect -g map="$GIS_OPT_MAP" layer="$GIS_OPT_LAYER" fs=";" | cut -f2 -d';'`
+COLUMN_TYPE=`v.info -c map="$GIS_OPT_MAP" layer="$GIS_OPT_LAYER" | grep -i "|$GIS_OPT_RGB_COLUMN$" | cut -f1 -d'|'`
+if [ -z "$COLUMN_TYPE" ] ; then
+    # RGB Column not found, create it
+    echo "Creating column <$GIS_OPT_RGB_COLUMN> ..."
+    v.db.addcol map="$GIS_OPT_MAP" layer="$GIS_OPT_LAYER" column="$GIS_OPT_RGB_COLUMN varchar(11)"
+    if [ $? -ne 0 ] ; then
+       g.message -e "Creating color column"
+       exit 1
+    fi
+elif [ "$COLUMN_TYPE" != "CHARACTER" -a "$COLUMN_TYPE" != "TEXT" ] ; then
+    g.message -e "Column <$GIS_OPT_RGB_COLUMN> is not of compatible type (found $COLUMN_TYPE)"
+    exit 1
+else
+    NUM_CHARS=`db.describe -c "$TABLE" | grep -i " $GIS_OPT_RGB_COLUMN:" | cut -f4 -d':'`
+    if [ "$NUM_CHARS" -lt 11 ] ; then
+       g.message -e "Color column <$GIS_OPT_RGB_COLUMN> is not wide enough (needs 11 characters)"
+       exit 1
+    fi
+fi
+
+
+# find data range
+if [ -n "$GIS_OPT_RANGE" ] ; then
+   #order doesn't matter
+   MINVAL=`echo "$GIS_OPT_RANGE" | grep '[[:digit:]]' | grep ',' | cut -f1 -d','`
+   MAXVAL=`echo "$GIS_OPT_RANGE" | grep '[[:digit:]]' | grep ',' | cut -f2 -d','`
+else
+   echo "Scanning values ..."
+   MINVAL=`v.db.select map="$GIS_OPT_MAP" column="$GIS_OPT_COLUMN" layer="$GIS_OPT_LAYER" | sort -n | grep '^[-0-9]' | head -n 1`
+   MAXVAL=`v.db.select map="$GIS_OPT_MAP" column="$GIS_OPT_COLUMN" layer="$GIS_OPT_LAYER" | sort -n | grep '^[-0-9]' | tail -n 1`
+fi
+echo " min=[$MINVAL]  max=[$MAXVAL]"
+if [ -z "$MINVAL" ] || [ -z "$MAXVAL" ] ; then
+   g.message -e "Scanning data range"
+   exit 1
+fi
+
+
+# setup internal region
+g.region save="tmp_vcolors.$$"
+WIND_OVERRIDE="tmp_vcolors.$$"
+export WIND_OVERRIDE
+g.region rows=2 cols=2
+
+# create dummy raster map
+if [ "$NCOLUMN_TYPE" = "INTEGER" ] ; then
+   r.mapcalc "tmp_colr_$$ = if( row() == 1, $MINVAL, $MAXVAL)"
+else
+   r.mapcalc "tmp_colr_$$ = double( if( row() == 1, $MINVAL, $MAXVAL) )"
+fi
+
+
+if [ -n "$GIS_OPT_COLOR" ] ; then
+   COLOR_CMD="color=$GIS_OPT_COLOR"
+elif [ -n "$GIS_OPT_RASTER" ] ; then
+   COLOR_CMD="raster=$GIS_OPT_RASTER"
+elif [ -n "$GIS_OPT_RULES" ] ; then
+   COLOR_CMD="rules=$GIS_OPT_RULES"
+fi
+if [ $GIS_FLAG_N -eq 1 ] ; then
+   FLIP_FLAG="-n"
+else
+   FLIP_FLAG=""
+fi
+r.colors map="tmp_colr_$$" "$COLOR_CMD" $FLIP_FLAG
+
+
+
+TMP="`g.tempfile pid=$$`"
+if [ $? -ne 0 ] || [ -z "${TMP}" ] ; then
+    g.message -e "Unable to create temporary file"
+    cleanup
+    exit 1
+fi
+
+
+# calculate colors and write SQL command file
+echo "Looking up colors ..."
+v.db.select map="$GIS_OPT_MAP" layer="$GIS_OPT_LAYER" column="$GIS_OPT_COLUMN" | \
+  sort -n | grep '^[-0-9]' | uniq | \
+  r.what.color -i in="tmp_colr_$$" | sed -e 's/: /|/' | grep -v '|\*$' | \
+  ( while read LINE ; do
+      #echo "LINE=[$LINE]"
+      VALUE=`echo $LINE | cut -f1 -d'|'`
+      COLR=`echo $LINE | cut -f2 -d'|'`
+      echo "UPDATE $TABLE SET $GIS_OPT_RGB_COLUMN = '$COLR' WHERE $GIS_OPT_COLUMN = $VALUE;" >> "${TMP}_vcol.sql"
+  done )
+
+
+if [ ! -e "${TMP}_vcol.sql" ] ; then
+    g.message -e "No values found in color range"
+    cleanup
+    exit 1
+fi
+
+# apply SQL commands to update the table with values
+NUM_RULES=`wc -l < "${TMP}_vcol.sql"`
+echo "Writing $NUM_RULES colors ..."
+#less "$TMP"
+db.execute input="${TMP}_vcol.sql"
+if [ $? -ne 0 ] ; then
+   g.message -e "Processing SQL transaction"
+   cleanup
+   exit 1
+fi
+
+
+if [ $GIS_FLAG_S -eq 1 ] ; then
+   g.rename rast="tmp_colr_$$","vcolors_$$"
+   echo "Raster map containing color rules saved to <vcolors_$$>"
+   # TODO save full v.colors command line history
+   r.support map="vcolors_$$" history="" \
+	source1="vector map=$GIS_OPT_MAP" \
+	source2="column=$GIS_OPT_COLUMN" \
+	title="Dummy raster to use as thematic vector legend" \
+	description="generated by v.colors using r.mapcalc"
+   r.support map="vcolors_$$" \
+	history="RGB saved into <$GIS_OPT_RGB_COLUMN> using <$GIS_OPT_COLOR$GIS_OPT_RASTER$GIS_OPT_RULES>"
+else
+   g.remove rast="tmp_colr_$$"
+fi
+cleanup
+
+exit 0
+#v.db.dropcol map=vcol_test col=GRASSRGB
+#d.vect -a vcol_test icon=basic/circle color=none size=8
+