#!/bin/bash
#
############################################################################
#
# 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.6, developed over Grass 6.1 CVS
#
# COPYRIGHT:    (c) 2006 Maciej Sieczka
#
#               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. Follow the instruction in html manual
###	page to obtain it.
###
### 2.	When creating the input elevation raster, make sure it's extents cover
###	your whole region of interest, the no-data cells are NULL, and the
###	resolution is correct. Same as most Grass raster modules this one is
###	region sensitive too.

#%Module
#%  description: Interpolate raster using the nnbathy natural neighbor interpolation program.
#%END
#%option
#% key: input
#% type: string
#% gisprompt: old,cell,raster
#% description: The raster to interpolate, with NULL assigned to no-data cells
#% required : yes
#%END
#%option
#% key: output
#% gisprompt: new,cell,raster
#% type: string
#% description: Name of the output raster
#% required : yes
#%END

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

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

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

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

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

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

PROG=`basename $0`

# cleanup temp files
cleanup()
{
 \rm -f $TMP $TMP.${PROG}.input_xyz $TMP.${PROG}.output_xyz $TMP.${PROG}.output_grd $TMP.${PROG}.region
}

# 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

INPUT="$GIS_OPT_INPUT"
OUTPUT="$GIS_OPT_OUTPUT"

### DO IT ###
echo
# grab the current region settings
g.region -g > $TMP.${PROG}.region

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

# extract variables from the grabbed region
north=`awk 'BEGIN {FS="="} (NR==1) {print $2}' $TMP.${PROG}.region`
south=`awk 'BEGIN {FS="="} (NR==2) {print $2}' $TMP.${PROG}.region`
west=`awk 'BEGIN {FS="="} (NR==3) {print $2}' $TMP.${PROG}.region`
east=`awk 'BEGIN {FS="="} (NR==4) {print $2}' $TMP.${PROG}.region`
nsres=`awk 'BEGIN {FS="="} (NR==5) {print $2}' $TMP.${PROG}.region`
ewres=`awk 'BEGIN {FS="="} (NR==6) {print $2}' $TMP.${PROG}.region`
rows=`awk 'BEGIN {FS="="} (NR==7) {print $2}' $TMP.${PROG}.region`
cols=`awk 'BEGIN {FS="="} (NR==8) {print $2}' $TMP.${PROG}.region`

nn_north=`awk -v res="$nsres" 'BEGIN {FS="="} (NR==1) {printf "%.8f",$2-res/2}' $TMP.${PROG}.region`
nn_south=`awk -v res="$nsres" 'BEGIN {FS="="} (NR==2) {printf "%.8f",$2+res/2}' $TMP.${PROG}.region`
nn_west=`awk -v res="$ewres" 'BEGIN {FS="="} (NR==3) {printf "%.8f",$2+res/2}' $TMP.${PROG}.region`
nn_east=`awk -v res="$ewres" 'BEGIN {FS="="} (NR==4) {printf "%.8f",$2-res/2}' $TMP.${PROG}.region`

null=NaN
type=double

# interpolate
echo
echo "nnbathy IN ACTION - MAY TAKE SOME TIME" 1>&2
echo
echo "PLEASE STAND BY UNTIL 'ALL DONE' IS PRINTED" 1>&2
echo

nnbathy -W 0 -P alg=nn -n "$cols"x"$rows" -x $nn_west $nn_east -y $nn_north $nn_south -i $TMP.${PROG}.input_xyz > $TMP.${PROG}.output_xyz
# Y in r.stats -1gn output is in descending order, thus -y must be MAX MIN 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: $north
south: $south
east: $east
west: $west
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

### ALL DONE ###

cleanup

echo
echo "ALL DONE" 1>&2
echo