[GRASS-SVN] r56946 - grass/branches/develbranch_6/scripts/i.spectral

svn_grass at osgeo.org svn_grass at osgeo.org
Fri Jun 28 04:57:05 PDT 2013


Author: hamish
Date: 2013-06-28 04:57:05 -0700 (Fri, 28 Jun 2013)
New Revision: 56946

Modified:
   grass/branches/develbranch_6/scripts/i.spectral/description.html
   grass/branches/develbranch_6/scripts/i.spectral/i.spectral
   grass/branches/develbranch_6/scripts/i.spectral/i_spectral.png
Log:
major overhaul: new non-interactive options, additional output formats, fix group map name parsing, lighter use of TMP files through pipes and awk, TMP files with g.tempfile instead of pwd, catch unclean exit for TMP file cleanup, bail out early if interactive mode and no Xmon is open, avoid overflow of enviro variables if many coord picks were made, various quoting, removed unused variables & rename others, reduce x-range by half a unit at either end, cleaner/improved legend text, and a new screenshot to match.


Modified: grass/branches/develbranch_6/scripts/i.spectral/description.html
===================================================================
--- grass/branches/develbranch_6/scripts/i.spectral/description.html	2013-06-28 11:26:00 UTC (rev 56945)
+++ grass/branches/develbranch_6/scripts/i.spectral/description.html	2013-06-28 11:57:05 UTC (rev 56946)
@@ -2,10 +2,19 @@
 
 <em>i.spectral</em> displays spectral response at user specified 
 locations in images.
+It requires the user to either interactively select positions on the
+Xmonitor (the default) or can function non-interactively with sampling
+coordinates specified with the <b>east_north</b> option.
 
+
 <h2>NOTES</h2>
 
-This script needs gnuplot to be installed.
+This script needs <a href="http://www.gnuplot.info">gnuplot</a> to be installed.
+Output may be to a new graphics window on your monitor (the default), or
+to an output file. Available output file formats are PNG, EPS, and SVG.
+In interactive mode use the <b>-m</b> flag to select multiple sampling
+points. Each coordinate pick will appear as a different colored line
+graph.
 
 
 <h2>EXAMPLES</h2>
@@ -21,7 +30,8 @@
 
 <center>
 <img src="i_spectral.png" border=1><br>
-Spectral plot of 3 different land cover types: (1) water, (2) green vegetation, and (3) highway
+Spectral plot of 3 different land cover types: (1) water, (2) green
+vegetation, and (3) highway
 </center>
 
 
@@ -30,7 +40,8 @@
 
 <div class="code"><pre>
 d.rast map_1
-LIST=`g.mlist type=rast mapset=timeseries pat="map_*" | sort -t '_' -k 2 -n | tr '\n' ','| sed 's+,$++g'`
+LIST=`g.mlist type=rast mapset=timeseries pat="map_*" | \
+   sort -t '_' -k 2 -n | tr '\n' ','| sed 's+,$++g'`
 i.spectral -i rast=$LIST
 </pre></div>
 
@@ -40,15 +51,19 @@
 resulting pixel values of all matching maps are drawn in the gnuplot
 output.
 
+
 <h2>SEE ALSO</h2>
 
 <em><a href="d.what.rast.html">d.what.rast</a></em><br>
 <em><a href="d.where.html">d.where</a></em><br>
 <em><a href="r.what.html">r.what</a></em><br>
 
+
 <h2>AUTHOR</h2>
 
 Markus Neteler<br>
-Francesco Pirotti
+Francesco Pirotti<br>
+Hamish Bowman
 
-<p><i>Last changed: $Date$</i>
+<p>
+<i>Last changed: $Date$</i>

Modified: grass/branches/develbranch_6/scripts/i.spectral/i.spectral
===================================================================
--- grass/branches/develbranch_6/scripts/i.spectral/i.spectral	2013-06-28 11:26:00 UTC (rev 56945)
+++ grass/branches/develbranch_6/scripts/i.spectral/i.spectral	2013-06-28 11:57:05 UTC (rev 56946)
@@ -5,31 +5,34 @@
 # MODULE:       i.spectral
 # AUTHOR(S):    Markus Neteler, 18. August 1998
 # PURPOSE:      displays spectral response at user specified locations in group or images
-# COPYRIGHT:    (C) 1999 by the GRASS Development Team
+# COPYRIGHT:    (C) 1999-2013 by 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.
 #
 #############################################################################
-
+#
+#  This script needs Gnuplot:  http://www.gnuplot.info
+#
 # written by Markus Neteler 18. August 1998
 #            neteler geog.uni-hannover.de
 # 
 # bugfix: 25. Nov.98/20. Jan. 1999
 # 3 March 2006: Added multiple images and group support by Francesco Pirotti - CIRGEO
-# this script needs gnuplot
+# June 2013: Scripting improvements and additional output formats by Hamish Bowman
 
 #%Module
-#%  description: Displays spectral response at user specified locations in group or images.
-#%  keywords: imagery, raster, multispectral
+#% description: Displays spectral response at user specified locations for an imagery group or a list of raster maps.
+#% keywords: imagery, raster, multispectral
 #%End
 #%option
 #% key: group
 #% type: string
-#% gisprompt: old,group
-#% description: Group input
+#% gisprompt: old,group,group
+#% description: Imagery group name
 #% required : no
+#% guisection: Input
 #%end
 #%option
 #% key: raster
@@ -38,37 +41,58 @@
 #% description: Raster input maps
 #% multiple : yes
 #% required : no
+#% guisection: Input
 #%end
 #%option
 #% key: output
 #% type: string
-#% gisprompt: output
-#% description: Write output to PNG image
+#% gisprompt: new,file,file
+#% label: Write output to image
+#% description: The default is to open a new window
 #% multiple : no
 #% required : no
+#% guisection: Output
 #%end
-#%flag 
-#%key: i
-#%description: Use image list and not group
+#%Option
+#% key: format
+#% type: string
+#% description: Graphics format for output file
+#% options: png,eps,svg
+#% answer: png
+#% multiple: no
+#% guisection: Output
+#%End
+#%Option
+#% key: east_north
+#% type: double
+#% required: no
+#% multiple: yes
+#% key_desc: east,north
+#% description: Coordinates for query (non-interactive mode)
+#%End
+#% flag 
+#% key: i
+#% description: Use a list of raster images instead of a named group
+#% guisection: Input
 #%end
-#%flag 
-#%key: m
-#%description: Select multiple points
+#% flag 
+#% key: m
+#% description: Select multiple points
 #%end
-#%flag 
-#%key: c
-#%description: Label with coordinates instead of numbering
+#% flag 
+#% key: c
+#% description: Show pick coordinates instead of numbering in the legend
 #%end
 
 if test "$GISBASE" = ""; then
- echo "You must be in GRASS GIS to run this program." >&2
- exit 1
+    echo "You must be in GRASS GIS to run this program." >&2
+    exit 1
 fi
 
 PARAM_NUM=$#
 
 if [ "$1" != "@ARGS_PARSED@" ] ; then
-  exec g.parser "$0" "$@"
+    exec g.parser "$0" "$@"
 fi
 
 #check if present
@@ -83,209 +107,222 @@
     exit 1
 fi
 
-
-# setting environment, so that awk works properly in all languages
+# set environment so that awk works properly in all languages
 unset LC_ALL
 LC_NUMERIC=C
 export LC_NUMERIC
 
+if [ `d.mon -L | grep -w '^x[0-9]' | grep -vc 'not running'` -eq 0 ] \
+   && [ -z "$GIS_OPT_EAST_NORTH" ] ; then
+    g.message -e "No graphics device open (requires d.mon X-windows display)"
+    exit 1
+fi
+
+if [ -z "$GIS_OPT_GROUP" ] && [ "$GIS_FLAG_I" -eq 0 ] ; then
+    g.message -e "Please use either the 'group' parameter or the '-i' flag"
+    exit 1
+fi
+
+#### set up the temp files
 TMP1="`g.tempfile pid=$$`"
+if [ $? -ne 0 ] || [ -z "$TMP1" ] ; then
+    g.message -e "Unable to create temporary file"
+    exit 1
+fi
+
 TMP2="`g.tempfile pid=$$`"
+if [ $? -ne 0 ] || [ -z "$TMP2" ] ; then
+    g.message -e "Unable to create temporary file"
+    exit 1
+fi
 
-# get y-data for gnuplot-data file
-# get data from group files and set the x-axis labels
+TMP_gnuplot="`g.tempfile pid=$$`"
+if [ $? -ne 0 ] || [ -z "$TMP_gnuplot" ] ; then
+    g.message -e "Unable to create temporary file"
+    exit 1
+fi
 
-xlabels=""
+cleanup()
+{
+   rm -f "$TMP1" "$TMP1.raw"
+   rm -f "$TMP2" "$TMP2.pick"[0-9]*
+   rm -f "$TMP_gnuplot"
+}
+trap "cleanup" 2 3 15
 
-if [ ! "$GIS_OPT_GROUP" ] ; then
-   if [ $GIS_FLAG_I -eq 0 ] ; then
-      g.message -e "Please use -i if you don't use the 'group' parameter"
-      exit 1
-   fi
-fi
 
+#### get map names from raster group and set the x-axis labels
+
 if [ $GIS_FLAG_I -eq 0 ] ; then
+    # Parse the group listing and set the x-axis labels as the map names
+    RASTERMAPS=`i.group -g group="$GIS_OPT_GROUP" | tr '\n' ',' | sed -e 's/,$//'`
+    NUMBANDS=`i.group -g group="$GIS_OPT_GROUP" | wc -l`
 
-	# ## PARSES THE GROUP FILES - gets rid of ugly header info from group list output
-	i.group  group=$GIS_OPT_GROUP -l | sed 's+ in + +g' | grep "group" -v | grep "\-\-\-\-$" -v > "$TMP2"
-	as=0
-	lists=""
-	xlabels=""
-	ass=0
-	for i in $( cat "$TMP2" ); do
-		er=`expr $as % 2`
-		as=`expr $as + 1`
-		if [ "$er" -eq "0" ]; then
-			ass=`expr $ass + 1`
-			if [ "$as" -eq "1" ]; then
-  	                  lists=$i
-			  xlabels="'$i' 1"
-			else
-   	                  lists=$lists,$i
-			  xlabels="$xlabels, '$i' $ass"
-		        fi
- 		fi       
-	 done
+    xlabels=`i.group -g group="$GIS_OPT_GROUP" | cut -f1 -d'@' | \
+        sed -e "s/^\|$/'/g" | awk '{print $1,NR ","}' | tr '\n' ' ' | \
+        sed -e 's/, $//'`
 
-NUMBANDS=$ass
-RASTERMAPS="$lists"
-
-# ## get data from list of files and set the x-axis labels
 else
-  RASTERMAPS="$GIS_OPT_RASTER"
-  news=`echo $GIS_OPT_RASTER | sed 's+,+ +g'`
-  as=0
-  for g in $news; do 
-    as=`expr $as + 1`
-    if [ "$as" -eq "1" ]; then
-       xlabels="'$g' $as"
-    else
-       xlabels="$xlabels, '$g' $as"
-    fi
-   done
+    # get list of maps and set the x-axis labels as the map names
+    RASTERMAPS="$GIS_OPT_RASTER"
+    NUMBANDS=`echo "$RASTERMAPS" | tr ',' '\n' | wc -l`
 
-   NUMBANDS=$as
+    xlabels=`echo "$RASTERMAPS" | tr ',' '\n' | cut -f1 -d'@' | \
+        sed -e "s/^\|$/'/g" | awk '{print $1,NR ","}' | tr '\n' ' ' | \
+        sed -e 's/, $//'`
 fi
 
-rm -f "$TMP2"
 
-#RASTERMAPS=`echo $RASTERMAPS | sed 's+ $++g' | sed 's+ +,+g'`
+#### pick sampling positions and extract DNs for gnuplot y-data
 
-if [ "$GIS_FLAG_M" -eq "0" ]; then
-  if [ `d.mon -L | grep -v 'not running' | wc -l | awk '{print $1}'` -eq 2 ] ; then
-	  g.message -e "No graphic device open (d.mon)"
-	  exit 1
-  fi
-  d.where -1 | r.what input=$RASTERMAPS > "$TMP1"
-else 
-  asss=0
-  if [ `d.mon -L | grep -v 'not running' | wc -l | awk '{print $1}'` -eq 2 ] ; then
-	  g.message -e  "No graphic device open (d.mon)"
-	  exit 1
-  fi
-  d.where | r.what input=$RASTERMAPS > "$TMP1"
-  cat "$TMP1" > tmper
-  rm -f "$TMP1"
-  for asas in $(cat "tmper"); do
-    asss=`expr $asss + 1`
-    echo $asas | sed 's+||+|'$asss'|+g' >> "$TMP1"
-  done
-fi
+if [ -z "$GIS_OPT_EAST_NORTH" ] ; then
+   ## Pick the points in the Xmonitor
+   
+   if [ "$GIS_FLAG_M" -eq 0 ] ; then
+      # single pick
+      d.where -1 | r.what input="$RASTERMAPS" null=0 > "$TMP1"
+   else
+      # multi-pick
+      d.where | r.what input="$RASTERMAPS" null=0 > "$TMP1.raw"
+   
+      n=0
+      for line in $( cat "$TMP1.raw" ) ; do
+         n=`expr $n + 1`
+         echo "$line" | sed "s+||+|$n|+g" >> "$TMP1"
+      done
+   fi
 
-cc=0
+else
+   ## Take input coordinates from the command line
+   r.what input="$RASTERMAPS" east_north="$GIS_OPT_EAST_NORTH" null=0 > "$TMP1.raw"
 
-TT="`cat "$TMP1" | sed 's+\*+0+g' | sed 's+|+ +g'`"
-
-NUMCLICKS="`cat "$TMP1" | wc -l`"
-START=0
-
-# tell data to outpu in different files to accept multiple plots - cache coordinates for output
-for hell in $( cat "$TMP1"); do
-   START=`expr $START + 1`
-   #cc=`expr $cc + 1`
-   COORD[$START]="`echo $hell | cut -d'|' -f1,2 | tr '|' ' '`"
-   #LABELS[$START]="`cat "$TMP1" | cut -d'|' -f3 `"
-   TOT="`echo $hell | sed 's+\*+0+g' | sed 's+|+ +g'`"
-   cc=0
-   for ix in $TOT; do
-    cc=`expr $cc + 1`
-    if [ "$cc" -gt 3 ]; then
-      echo $ix >> "$TMP2"_$START
-    fi
+   n=0
+   for line in $( cat "$TMP1.raw" ) ; do
+      n=`expr $n + 1`
+      echo "$line" | sed "s+||+|$n|+g" >> "$TMP1"
    done
+fi
 
-done
 
-rm -f "$TMP1"
+NUMCLICKS=`wc -l < "$TMP1"`
 
-NUM=$NUMBANDS
-#echo $NUM
-#echo $NUMCLICKS RRRRRRRRRRRR
-NUM2=`expr $NUM / $NUMCLICKS`
-NUM2=`expr $NUM2 + 1`
+if [ "$NUMCLICKS" -eq 0 ] ; then
+   g.message -e "No points selected"
+   exit 1
+fi
 
 
-# build data file
-# the x-axis... depending on number of bands
-rm -f data.dum
+### tell data to output in different files to create multiple lines on the
+#  plot, and cache pick coordinates for (optional) use in the legend.
+PICK_NUM=0
+for PICK in $( cat "$TMP1" ) ; do
+   PICK_NUM=`expr "$PICK_NUM" + 1`
 
-  i=0
-  while [ $i != $NUMBANDS ]
-   do
-    i=`expr $i + 1`
-    echo $i >> data.dum
-   done 
+   COORD[$PICK_NUM]=`echo $PICK | cut -d'|' -f1,2 | tr '|' ' '`
+   #LABELS[$PICK_NUM]=`cat "$TMP1" | cut -d'|' -f3 `
 
-r=0
-while [ $r != $NUMCLICKS ]
-do
-  r=`expr $r + 1`
-  paste -d' ' data.dum "$TMP2"_$r > data_$r
-  rm -f "$TMP2"_$r
+   echo "$PICK" | cut -f4- -d'|' | tr '|' '\n' | \
+      awk '{print NR,$0}' > "$TMP2.pick$PICK_NUM"
 done
 
-# paste to data file
-#paste -d' ' data.dum "$TMP2" > data
-xrange=`expr $NUMBANDS + 1`
 
-# build gnuplot script
+
+#### build up the Gnuplot script
+
 if [ -n "$GIS_OPT_OUTPUT" ] ; then
-  echo "set term png large" >> spectrum.gnuplot
-  echo "set output '$GIS_OPT_OUTPUT'" >> spectrum.gnuplot
+   # interestingly enough, a "grass" terminal exists for gnuplot
+   #  which will render directly to the Xmonitor. It wasn't built
+   #  into my copy though so I couldn't test it.
+   case "$GIS_OPT_FORMAT" in
+      png)
+	term_opts="png truecolor large size 825,550"
+	;;
+      eps)
+	term_opts="postscript eps color solid size 6,4"
+	;;
+      svg)
+	term_opts="svg size 825,550 dynamic solid"
+	;;
+   esac
+
+   echo "set term $term_opts" > "$TMP_gnuplot"
+   echo "set output '$GIS_OPT_OUTPUT'" >> "$TMP_gnuplot"
 fi
-echo "set xtics ($xlabels)" >> spectrum.gnuplot
-echo "set grid" >> spectrum.gnuplot
-echo "set title 'Spectral signatures'" >> spectrum.gnuplot
-##echo "set yrange [0:255]" >> spectrum.gnuplot
-echo "set xrange [0:$xrange]" >> spectrum.gnuplot
-echo "set noclabel" >> spectrum.gnuplot
-echo "set xlabel 'Bands'" >> spectrum.gnuplot
-echo "set ylabel 'DN Value'" >> spectrum.gnuplot
 
+xrange=`expr "$NUMBANDS" + 1`
+
+cat << EOF >> "$TMP_gnuplot"
+set xtics ($xlabels)
+set grid
+set title 'Spectral signatures'
+set xrange [0.5 : $xrange - 0.5]
+set noclabel
+set xlabel 'Bands'
+set xtics rotate by -40
+set ylabel 'DN Value'
+EOF
+
+
+## alternate plot options:
+# expand y-axis range to full:
+#echo "set yrange [0 : 255]" >> "$TMP_gnuplot"
+
+
 ## if more then 2 points we can draw an interpolated spline:
 #if [ $PARAM_NUM -gt 2 ]
 #then
-#   echo "set style data linespoints" >> spectrum.gnuplot
-#   echo "plot 'data' with points pt 779, '' smooth csplines t 'spline interpolated'" >> spectrum.gnuplot
+#   echo "set style data linespoints" >> "$TMP_gnuplot"
+#   echo "plot 'data' with points pt 779, '' smooth csplines t 'spline interpolated'" >> "$TMP_gnuplot"
 #else
-   echo "set style data lines" >> spectrum.gnuplot
-  
 
- rm -f data.dum
+# draw markers in all pick lines:
+#echo "set style data linespoints" >> "$TMP_gnuplot"
+# only draw markers on the final pick line: (aesthetic choice)
+echo "set style data lines" >> "$TMP_gnuplot"
 
- i=0
-str=""
-  while [ $i != $NUMCLICKS ]
-   do
-    i=`expr $i + 1`
-    if [ "$i" -eq "1" ]
-     then
-      if [ $GIS_FLAG_C -eq 0 ] ; then
-        str="plot 'data_$i' title '$i'" 
+
+# check if we are in a lat/long location
+if [ `g.region -p | grep -w '^projection' | cut -f2 -d' '` = "3" ] ; then
+  LL=true
+else
+  LL=false
+fi
+
+
+i=0
+while [ "$i" -lt "$NUMCLICKS" ] ; do
+   i=`expr "$i" + 1`
+
+   # less noisy coordinates for legend
+   if [ "$GIS_FLAG_C" -eq 1 ] ; then
+      if [ "$LL" = "true" ] ; then
+         COORDS="${COORD[$i]}"
       else
-        str="plot 'data_$i' title '${COORD[$i]}'" 
+         COORDS=`echo "${COORD[$i]}" | awk '{printf("%.3f, %.3f", $1, $2)}'`
       fi
-     else
-      if [ $GIS_FLAG_C -eq 0 ] ; then
-        str=$str",'data_$i' title '$i'" 
+   fi
+
+   # use 'printf' instead of 'echo -n' to stay POSIX/SUS compatible with BSD & MacOSX
+   if [ "$i" -eq 1 ] ; then
+      if [ "$GIS_FLAG_C" -eq 0 ] ; then
+        printf "plot '$TMP2.pick$i' title 'Pick $i'" >> "$TMP_gnuplot"
       else
-        str=$str",'data_$i' title '${COORD[$i]}'" 
+        printf "plot '$TMP2.pick$i' title '$COORDS'" >> "$TMP_gnuplot"
       fi
-     fi
+   else
+      if [ "$GIS_FLAG_C" -eq 0 ] ; then
+         printf "$str, '$TMP2.pick$i' title 'Pick $i'" >> "$TMP_gnuplot"
+      else
+         printf "$str, '$TMP2.pick$i' title '$COORDS'" >> "$TMP_gnuplot"
+      fi
+   fi
+done
 
-   done 
-str=$str" with linespoints pt 779"
-echo $str >> spectrum.gnuplot
+echo " with linespoints pt 779" >> "$TMP_gnuplot"
 
 
-$GRASS_GNUPLOT spectrum.gnuplot 
- 
+# we don't quote the executable since it could contain the '-persist' command line arg
+$GRASS_GNUPLOT "$TMP_gnuplot" 
 
-i=0
-while [ $i != $NUMCLICKS ]
-   do
- i=`expr $i + 1`
- rm -f data_$i
-done
 
-rm -f spectrum.gnuplot
+cleanup

Modified: grass/branches/develbranch_6/scripts/i.spectral/i_spectral.png
===================================================================
(Binary files differ)



More information about the grass-commit mailing list