[GRASS-SVN] r34229 - grass-addons/raster/r.isoregions
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Nov 10 09:41:12 EST 2008
Author: mathieug
Date: 2008-11-10 09:41:12 -0500 (Mon, 10 Nov 2008)
New Revision: 34229
Added:
grass-addons/raster/r.isoregions/r.isoregions
Log:
isoregions creation from grass raster
Added: grass-addons/raster/r.isoregions/r.isoregions
===================================================================
--- grass-addons/raster/r.isoregions/r.isoregions (rev 0)
+++ grass-addons/raster/r.isoregions/r.isoregions 2008-11-10 14:41:12 UTC (rev 34229)
@@ -0,0 +1,268 @@
+#!/bin/ksh
+#
+############################################################################
+#
+# MODULE: r.isoregions
+# AUTHOR(S): Mathieu Grelier (greliermathieu at gmail.com)
+# PURPOSE: isoregions creation from grass raster
+# REQUIREMENTS:
+# - unix utilities : bc and funcs.bc file (additionnal functions for bc : http://docs.google.com/Doc?id=dgqnzvwg_6hh3xn5cg&hl=en)
+# COPYRIGHT: (C) 2007 Mathieu Grelier
+#
+# This program is free software under the GNU General Public
+# License (>=v2). Read the file COPYING that comes with GRASS
+# for details.
+#
+#############################################################################
+
+
+#%Module
+#% description: Create isoregions from a raster
+#% keywords: raster, isoregions, recode
+#%End
+#%option
+#% key: rastername
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: input raster
+#% required : yes
+#%end
+#%option
+#% key: output
+#% type: string
+#% description: Name of the output recoded raster. Required if o flag is not set
+#% required : no
+#%end
+#%option
+#% key: width
+#% type: double
+#% description: width of the interval between classes.
+#% required : yes
+#%end
+#%option
+#% key: scale
+#% type: integer
+#% answer: 2
+#% description: number of digits to be retained after the decimal point for classes
+#% required : no
+#%end
+#%option
+#% key: rules
+#% type: string
+#% answer: rainbow
+#% description: color rules file
+#% required : no
+#%end
+#%flag
+#% key: r
+#% description: replace original raster keeping the same name.
+#%end
+#%flag
+#% key: k
+#% description: keep recode rules file.
+#%end
+
+if [ -z "$GISBASE" ] ; then
+ echo "You must be in GRASS GIS to run this program."
+ exit 1
+fi
+
+if [ "$1" != "@ARGS_PARSED@" ] ; then
+ exec g.parser "$0" "$@"
+fi
+
+## GRASS team recommandations
+#####################################
+
+PROG=`basename $0`
+
+# check if we have awk
+if [ ! -x "`which awk`" ] ; then
+ echo "$PROG: 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
+
+# what to do in case of user break:
+exitprocedure()
+{
+ echo "User break!"
+ cleanup
+ exit 1
+}
+# shell check for user break (signal list: trap -l)
+trap "exitprocedure" 2 3 15
+
+##Config and general procedures
+####################################
+
+#linux
+homedir="$HOME"
+
+#cleanup procedure
+cleanup()
+{
+ if [ "$GIS_FLAG_R" -eq 1 ] ; then
+ g.remove rast="$rastertorecode" >> $LOGFILE 2>&1
+ fi
+ \rm -f $RECODE
+ \rm -f $CATS
+ \rm -f $RASTERRANGE
+}
+
+#fix this path
+if [ -z "$PROCESSDIR" ] ; then
+ PROCESSDIR="$HOME"
+fi
+
+#fix this path
+if [ -z "$LOGDIR" ] ; then
+ LOGDIR="$HOME"
+fi
+LOGFILE="$LOGDIR/r.isoregion.log"
+
+echo "r.colors.isoregions :" >> "$LOGFILE"
+
+
+## necessary checks
+#####################################
+
+
+#####FUNCTIONS
+##########################################################################
+
+#create recode file for r.recode with intervals being given width multiples, between given min and max
+function create_recode_file {
+
+ RECODE="`g.tempfile pid=$$`"
+ if [ $? -ne 0 ] || [ -z "$RECODE" ] ; then
+ echo "ERROR: unable to create temporary file for recode rules" 1>&2
+ exit 1
+ fi
+ CATS="`g.tempfile pid=$$`"
+ if [ $? -ne 0 ] || [ -z "$CATS" ] ; then
+ echo "ERROR: unable to create temporary file for categories" 1>&2
+ exit 1
+ fi
+
+ #min and max scaling
+ bcscale=$4
+ mini=$(echo "scale="$bcscale"; ("$1"/1)" | bc )
+ maxi=$(echo "scale="$bcscale"; ("$2"/1)" | bc )
+ width=$(echo "scale="$bcscale"; ("$3"/1)" | bc )
+ range=$(echo "$maxi - ($mini)" | bc )
+ #debug :
+ echo "mini=$mini"; echo "maxi=$maxi"; echo "width=$width"; echo "scale=$bcscale"; echo "range=$range"
+ if [ "$(echo "if (${width} > ${range}) 1" | bc)" -eq 1 ] ; then
+ echo "ERROR: given width is greater than range (min=$mini, max=$maxi, range=$range" 1>&2
+ exit 1
+ fi
+
+ #first interval starts from min and goes to intial level
+ #intervals width is $width so we use rup function to round $mini up to next multiple of $width
+ levelini=$(echo "rup($mini,$width)" | funcs.bc )
+ #how many regular intervals (with width=$width)?
+ rangefromlevelini=$(echo "$maxi - ($levelini)" | bc )
+ nbregints=$(echo "(floor($rangefromlevelini/$width))" | funcs.bc )
+ #now we know how many categories there will be, we can start to write the cats file for raster legend
+ let "nbcats = $nbregints + 2"
+ echo "$nbcats categories" > $CATS ; echo "recoded by r.isoregions" >> $CATS ; echo "no data" >> $CATS
+ #we can write the first recode rule
+ let "levelnext = $levelini"
+ let "classno = 1"
+ echo *:"$levelini":"$classno" > $RECODE
+ echo "$mini-$levelini" >> $CATS
+
+ #other intervals excepted the last one have the same length ; we start from the second class
+ for ((i = 0; i < "$nbregints" ; i++))
+ do
+ let "levelnext += $width"
+ let "classno++"
+ echo "$levelini":"$levelnext":"$classno" >> $RECODE
+ echo "$levelini-$levelnext" >> $CATS
+ let "levelini = $levelnext"
+ done
+
+ #last interval :
+ let "classno++"
+ echo "$levelini":*:"$classno" >> $RECODE
+ echo "$levelini-$maxi" >> $CATS
+}
+
+##### SCRIPT
+###########################################################################
+
+#isoregions creation
+#script is based on r.recode
+
+#preliminary operations
+echo "prepare raster" >> $LOGFILE
+if [ "$GIS_FLAG_R" -eq 1 ] ; then
+ rastertorecode=tmpraster
+ outputraster="$GIS_OPT_RASTERNAME"
+ g.copy "$GIS_OPT_RASTERNAME","$rastertorecode" >> $LOGFILE 2>&1
+ g.remove rast="$GIS_OPT_RASTERNAME" >> $LOGFILE 2>&1
+else
+ if [ -z "$GIS_OPT_OUTPUT" ] ; then
+ echo "ERROR: output raster name is needed when r flag is not set." 1>&2
+ cleanup
+ exit 1
+ fi
+ rastertorecode="$GIS_OPT_RASTERNAME"
+ outputraster="$GIS_OPT_OUTPUT"
+fi
+
+#retrieve min and max for recode rules file creation
+#region parameters are used to deal with grid cells size and range
+echo "query raster" >> $LOGFILE
+RASTERRANGE="`g.tempfile pid=$$`"
+if [ $? -ne 0 ] || [ -z "$RASTERRANGE" ] ; then
+ echo "ERROR: unable to create temporary file for range parameters" 1>&2
+ cleanup
+ exit 1
+fi
+r.info -r map=$rastertorecode > $RASTERRANGE
+min=$(awk -F "=" '$1=="min" {print $2}' $RASTERRANGE)
+max=$(awk -F "=" '$1=="max" {print $2}' $RASTERRANGE)
+echo "create recode file" >> $LOGFILE
+create_recode_file $min $max "$GIS_OPT_WIDTH" "$GIS_OPT_SCALE"
+if [ $? -ne 0 ] ; then
+ echo "ERROR: unable to create recode file" 1>&2
+ cleanup
+ exit 1
+fi
+
+#recode raster
+echo "recode raster" >> $LOGFILE
+cat "$RECODE" | r.recode input="$rastertorecode" output="$outputraster" --overwrite
+if [ $? -ne 0 ] ; then
+ echo "ERROR: while executing r.recode" 1>&2
+ cleanup
+ exit 1
+fi
+
+echo "finalize" >> $LOGFILE
+#create categories file for legend
+eval `g.gisenv`
+cat "$CATS" > "$GISDBASE"/"$LOCATION_NAME"/"$MAPSET"/cats/"$outputraster"
+if [ $? -ne 0 ] ; then
+ echo "WARNING: unable to write categories file" 1>&2
+fi
+
+#coloring is necessary because mapcalc initialize colors
+r.colors map="$outputraster" rules="$GIS_OPT_RULES" >> $LOGFILE 2>&1
+if [ $? -ne 0 ] ; then
+ echo "ERROR: while executing r.colors" 1>&2
+ cleanup
+ exit 1
+fi
+
+#script end
+if [ "$GIS_FLAG_K" -eq 1 ] ; then
+ cat "$RECODE" > "$PROCESSDIR"/r.isoregions.recode.rules
+fi
+cleanup
Property changes on: grass-addons/raster/r.isoregions/r.isoregions
___________________________________________________________________
Name: svn:executable
+ *
More information about the grass-commit
mailing list