[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