[GRASS-user] nulls left by r.surf.nnbathy

[Maciej Sieczka]

Maciej Sieczka wrote:

> Try using the attached earlier version of r.surf.nnbathy (1.2).

Now attaching it. Really.

Maciek
-------------- next part --------------
#!/bin/bash
#
############################################################################
#
# MODULE:       r.surf.nnbathy
#
# AUTHOR(S):    Maciej Sieczka, sieczka at biol.uni.wroc.pl
#		Intitute of Plant Biology, Wroclaw University, Poland
#
# PURPOSE:	Interpolate raster surface using the "nnbathy" natural
#		neighbor interpolation program.
#
# VERSION:	1.2, developed over Grass 6.1 CVS 2006.05.18
#
# 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, download from
###	http://www.marine.csiro.au/~sakov/. At the time of writing the script
###	1.57 was the most recent version.
###
### 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. 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." >&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, please download from http://www.marine.csiro.au/~sakov/ and put the nnbathy executable in your $PATH." 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}.region $TMP.${PROG}.canvas_xyz $TMP $TMP.${PROG}.input_xyz $TMP.${PROG}.output_xyz $TMP.${PROG}.output_grd
}

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

# extract the current region settings
g.region -g > $TMP.${PROG}.region

# set variables
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`
rows=`awk 'BEGIN {FS="="} (NR==7) {print $2}' $TMP.${PROG}.region`
cols=`awk 'BEGIN {FS="="} (NR==8) {print $2}' $TMP.${PROG}.region`
null=nan
type=double

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

# spit out the "canvas" coords for nnbathy to interpolate in
r.stats -1g input="${INPUT}" | awk '{print $1" "$2}' > $TMP.${PROG}.canvas_xyz

# interpolate
echo
echo "nnbathy in action. May take some time. Please stand by." 1>&2
echo
nnbathy -i $TMP.${PROG}.input_xyz -o $TMP.${PROG}.canvas_xyz > $TMP.${PROG}.output_xyz

# 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. Please stand by." 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"

#### TESTING ###
#cp $TMP.${PROG}.output_grd dem.asc


# 3 import
r.in.ascii input=$TMP.${PROG}.output_grd output=${OUTPUT}

### ALL DONE ###

echo
echo "Done." 1>&2
echo

cleanup