#!/bin/sh

############################################################################
#
# MODULE:       r.surf.nnbathy
#
# AUTHOR(S):    Maciej Sieczka, sieczka@biol.uni.wroc.pl
#		Intitute of Plant Biology, Wroclaw University, Poland
#
# PURPOSE:	Interpolate raster surface using the "nnbathy" natural
#		neighbor interpolation program.
#
# VERSION:	1.95, developed over GRASS 6.3 CVS
#
# COPYRIGHT:    (c) 2006-2007 Maciej Sieczka
#
# NOTES:	Developed within the course and for the purpose of the CAVES
#		project, funded under EU 6FP NEST programme.
#
#               This program is free software under the GNU General Public
#               License (>=v2). Read the file COPYING that comes with GRASS
#               for details.
#
#############################################################################

# NOTES:
#
# 1. Requires nnbathy executable v 1.75 or later. Follow the instruction in
#    html manual page to obtain it.
#
# 2. When creating the input raster map, make sure it's extents cover
#    your whole region of interest, the cells you wish to interplate on are
#    NULL, and the resolution is correct. Same as most GRASS raster modules
#    this one is region sensitive too.

# CHANGELOG:

# 1.6, 2006.06.15:
# - first public release

# 1.9, 2007.01.02:
# - parse g.region -g with eval, not awk
# - support all interpolation methods nnbathy provides
# - create history for output raster with r.support
# - require nnbathy 1.69 (major bug was fixed)
# - try to detect if nnbathy failed and exit cleanly
# - documentation extended
# - minor cleanups
# - todo: there is a progress indicator switch (-%) in nnbathy now, but using it
#   slows down the interpolation 3-4 times on my machine, while works OK on nnbathy
#   author's; if sorted out, I'll use -%

# 1.95, 2007.11.12:
# - require nnbathy 1.75 - bugfixes and speed improvemenets for large grids;
#   refer to file CHANGELOG in nn sorce code for details

# 1.96, 2008.25.03:
# - handle spaces in pathnames
# - get rid of Bashims
# - better run-time error trapping
# - require nnbathy 1.76 - contains a bugfix; refer to file CHANGELOG in nn
#   sorce code for details
# - cosmetics
# - manual cleaned up

#%Module
#%  description: Interpolate raster using the nnbathy natural neighbor interpolation program.
#%END

#%option
#% key: input
#% type: string
#% gisprompt: old,cell,raster
#% description: The raster map to interpolate on
#% required : yes
#%END

#%option
#% key: output
#% gisprompt: new,cell,raster
#% type: string
#% description: Name of the output raster map
#% required : yes
#%END

#%option
#% key: alg
#% type: string
#% options: l,nn,ns
#% answer: nn
#% description: Interpolation algorithm for nnbathy to use
#% required : yes
#%END


# called from GRASS?
if test "$GISBASE" = ""; then
    echo
    echo "ERROR: You must be in GRASS GIS to run this program." 1>&2
    echo
    exit 1
fi

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

# check if we have awk
if [ ! -x "`which awk`" ] ; then
    echo
    echo "ERROR: awk required. Please install awk or gawk first." 1>&2
    echo
    exit 1
fi

# check if we have nnbathy
if [ ! -x "`which nnbathy`" ] ; then
    echo
    echo "ERROR: nnbathy required. To obtain it please follow the instruction in r.surf.nnbathy manual". 1>&2
    echo
    exit 1
fi

# check nnbathy version
nnv=`nnbathy -v | sed 's/ /\n/g' | sort -nr | head -n1`
nnv_ok=`echo $nnv | awk '{ if ($0<1.76) {print 0} else {print 1} }'`

    if [ $nnv_ok -eq 0 ] ; then
	echo "ERROR: nnbathy version >= 1.76 is required."
        exit 1
    fi

# set up temporary files
TMP="`g.tempfile pid=$$`"
if [ $? -ne 0 ] || [ -z "$TMP" ] ; then
    echo
    echo "ERROR: Unable to create temporary files." 1>&2
    echo
    exit 1
fi

# set environment so that awk works properly in all languages
unset LC_ALL
export LC_NUMERIC=C

PROG=`basename $0`

# define cleanup procedure
proc_cleanup()
{
 # Reset traps before normal script termination to avoid bogus ERROR message, as
 # we put a trap on signal 0.
 trap - 0 2 3 15 
 rm -f "$TMP" "$TMP.${PROG}.input_xyz" "$TMP.${PROG}.output_xyz" "$TMP.${PROG}.output_grd"
}

# define run-time error handling procedure
proc_runtime_error()
{
 echo
 echo "ERROR: There was an error at the script's run time. Please try to debug the problem if you can and let me know by email." 1>&2
 echo
 exit 1
}

# define user-break procedure
proc_user_break()
{
 echo
 echo "User break!"
 echo
 proc_cleanup
 exit 1
}

# set trap for when script terminates
trap "proc_runtime_error" 0

# trap user break
trap "proc_user_break" 2 3 15

# assign main variables from user's input
INPUT="$GIS_OPT_INPUT"
OUTPUT="$GIS_OPT_OUTPUT"
ALG="$GIS_OPT_ALG"



### DO IT ###

set -e	# Make script terminate (ie. emit signal 0) if any statement returns a
	# non-0 value. Then we trap signal 0, which lets handle such errors.
	# However, the trap on signal 0 must be reset before normal script
	# termination to avoid a bogus ERROR message then - this is done in the
	# cleanup procedure.

# grab the current region settings
eval `g.region -gp`

# spit out non-null (-n) raster coords + values to be interpolated
r.stats -1gn input="${INPUT}" > "$TMP.${PROG}.input_xyz"

# set the working region for nnbathy (it's cell-center oriented)
nn_n=`echo $n | awk -v res="$nsres" '{printf "%.8f",$1-res/2.0}'`
nn_s=`echo $s | awk -v res="$nsres" '{printf "%.8f",$1+res/2.0}'`
nn_w=`echo $w | awk -v res="$ewres" '{printf "%.8f",$1+res/2.0}'`
nn_e=`echo $e | awk -v res="$ewres" '{printf "%.8f",$1-res/2.0}'`

null=NaN
type=double

# interpolate
echo
echo "WARNING: nnbathy is performing the interpolation now. This may take some time." 1>&2
echo
echo "WARNING: Once it completes 'All done.' message is printed." 1>&2
echo

nnbathy -W 0 -P alg=$ALG -n ${cols}x${rows} -x $nn_w $nn_e -y $nn_n $nn_s -i "$TMP.${PROG}.input_xyz" > "$TMP.${PROG}.output_xyz"
# Y in "r.stats -1gn" output is in descending order, thus -y must be in MAX MIN order, not MIN MAX, for nnbathy not to produce a grid upside-down

# convert the X,Y,Z nnbathy output into a GRASS ASCII grid, then import with r.in.ascii:

# 1 create header
cat << EOF > "$TMP.${PROG}.output_grd"
north: $n
south: $s
east: $e
west: $w
rows: $rows
cols: $cols
null: $null
type: $type
EOF

# 2 do the conversion
echo "Converting nnbathy output to GRASS raster." 1>&2
echo

awk -v cols="$cols" '
BEGIN {col_cur=1; ORS=" "}
{
 if (col_cur==cols) {ORS="\n"; col_cur=0; print $3; ORS=" "}
		    else {print $3}
 col_cur++
}' "$TMP.${PROG}.output_xyz" >> "$TMP.${PROG}.output_grd"

# 3 import
r.in.ascii input="$TMP.${PROG}.output_grd" output="${OUTPUT}" > /dev/null

# store comand history in raster's metadata

r.support map=${OUTPUT} history=""
r.support map=${OUTPUT} history="script syntax:"
r.support map=${OUTPUT} history=""
r.support map=${OUTPUT} history="r.surf.nnbathy alg=$ALG input=${INPUT} output=${OUTPUT}"
r.support map=${OUTPUT} history=""
r.support map=${OUTPUT} history="nnbathy syntax:"
r.support map=${OUTPUT} history=""
r.support map=${OUTPUT} history="nnbathy -W 0 -P alg=$ALG -n ${cols}x${rows} "
r.support map=${OUTPUT} history="-x $nn_w $nn_e "
r.support map=${OUTPUT} history="-y $nn_n $nn_s "
r.support map=${OUTPUT} history="-i tmp_in > tmp_out"
r.support map=${OUTPUT} history=""

### ALL DONE ###

proc_cleanup

echo
echo "All done." 1>&2
echo