#!/bin/bash

############################################################################
#
# 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-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
# June 2013: Scripting improvements and additional output formats by Hamish Bowman

#%Module
#% 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,group
#% description: Imagery group name
#% required : no
#% guisection: Input
#%end
#%option
#% key: raster
#% type: string
#% gisprompt: old,cell,raster
#% description: Raster input maps
#% multiple : yes
#% required : no
#% guisection: Input
#%end
#%option
#% key: output
#% type: string
#% 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
#%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
#%end
#% 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
fi

PARAM_NUM=$#

if [ "$1" != "@ARGS_PARSED@" ] ; then
    exec g.parser "$0" "$@"
fi

#check if present
if [ ! -x "`which gnuplot`" ] ; then
    g.message -e "gnuplot required, please install first"
    exit 1
fi

#### check if we have awk
if [ ! -x "`which awk`" ] ; then
    g.message -e "awk required, please install awk or gawk first"
    exit 1
fi

# 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

TMP_gnuplot="`g.tempfile pid=$$`"
if [ $? -ne 0 ] || [ -z "$TMP_gnuplot" ] ; then
    g.message -e "Unable to create temporary file"
    exit 1
fi

cleanup()
{
   rm -f "$TMP1" "$TMP1.raw"
   rm -f "$TMP2" "$TMP2.pick"[0-9]*
   rm -f "$TMP_gnuplot"
}
trap "cleanup" 2 3 15


#### 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`

    xlabels=`i.group -g group="$GIS_OPT_GROUP" | cut -f1 -d'@' | \
        sed -e "s/^\|$/'/g" | awk '{print $1,NR ","}' | tr '\n' ' ' | \
        sed -e 's/, $//'`

else
    # 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`

    xlabels=`echo "$RASTERMAPS" | tr ',' '\n' | cut -f1 -d'@' | \
        sed -e "s/^\|$/'/g" | awk '{print $1,NR ","}' | tr '\n' ' ' | \
        sed -e 's/, $//'`
fi


#### pick sampling positions and extract DNs for gnuplot y-data

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

else
   ## Take input coordinates from the command line
   r.what input="$RASTERMAPS" east_north="$GIS_OPT_EAST_NORTH" 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


NUMCLICKS=`wc -l < "$TMP1"`

if [ "$NUMCLICKS" -eq 0 ] ; then
   g.message -e "No points selected"
   exit 1
fi


### 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`

   COORD[$PICK_NUM]=`echo $PICK | cut -d'|' -f1,2 | tr '|' ' '`
   #LABELS[$PICK_NUM]=`cat "$TMP1" | cut -d'|' -f3 `

   echo "$PICK" | cut -f4- -d'|' | tr '|' '\n' | \
      awk '{print NR,$0}' > "$TMP2.pick$PICK_NUM"
done



#### build up the Gnuplot script

if [ -n "$GIS_OPT_OUTPUT" ] ; then
   # 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

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" >> "$TMP_gnuplot"
#   echo "plot 'data' with points pt 779, '' smooth csplines t 'spline interpolated'" >> "$TMP_gnuplot"
#else

# 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"


# check if we are in a lat/long location
if [ `g.region -pu | 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
         COORDS=`echo "${COORD[$i]}" | awk '{printf("%.3f, %.3f", $1, $2)}'`
      fi
   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
        printf "plot '$TMP2.pick$i' title '$COORDS'" >> "$TMP_gnuplot"
      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

echo " with linespoints pt 779" >> "$TMP_gnuplot"


# we don't quote the executable since it could contain the '-persist' command line arg
$GRASS_GNUPLOT "$TMP_gnuplot" 


cleanup