#!/bin/sh
##################################################################
#
# MODULE:       v.ascii.addrastcol
# AUTHOR(S):    R. Todd Jobe <toddjobe@unc.edu>
# PURPOSE:      Add a raster value field to a csv of points
# COPYRIGHT:    (C) 2009 by R. Todd Jobe
#
#  This program is free software; you can redistribute it and/or modify
#  it under the terms of the GNU General Public License as published by
#  the Free Software Foundation; either version 2 of the License, or
#  (at your option) any later version.
#
#  This program is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#  GNU General Public License for more details.
#
##################################################################
#%Module
#% description: Outputs an ASCII points file that has the intersected values from a raster
#% keywords: vector, ascii, raster
#%End
#%Flag
#% key: h
#% description: File includes a header
#%End
#%Flag
#% key: o
#% description: Overwrite original files
#%End
#%Flag
#% key: r
#% description: Remove vector files created from input on exit
#%End
#%Flag
#% key: m
#% description: For dems, convert feet to meters
#%End
#%Option
#% key: input
#% type: string
#% required: yes
#% multiple: yes
#% label: files to be intersected with the raster
#% description: Need to be in the same coordinate system as the raster
#% gisprompt: old_file,file,input
#%End
#%Option
#% key: rast
#% type: string
#% required: yes
#% multiple: no
#% key_desc: raster
#% label: raster to intersect 
#% description: Need to be in the same coordinate system as the input
#%End
#%Option
#% key: suffix
#% type: string
#% required: no
#% multiple: no
#% key_desc: string
#% label: Suffix for created files
#% answer: _isect.csv
#% description: Overridden by -o
#%End
#%Option
#% key: fs
#% type: string
#% required: no
#% multiple: no
#% key_desc: character
#% description: Field separator
#% answer: ,
#%End
#%Option
#% key: x
#% type: integer
#% required: no
#% multiple: no
#% label: Number of column used as x coordinate (points mode)
#% description: First column is 1
#% answer: 1
#%End
#%Option
#% key: y
#% type: integer
#% required: no
#% multiple: no
#% label: Number of column used as y coordinate (points mode)
#% description: First column is 1
#% answer: 2
#%End

if [ -z "$GISBASE" ] ; then
    echo "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

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

# Error Codes 
E_NOTFOUND=66

#### Main Code ###################################################
# Default csv files
#${GIS_OPT_INPUT="PlantPop05_SP_Feet.csv,PlantPop99_SP_Feet.csv,Pop_Centroids_SP_Feet.csv"}

# Need space-delimited arguments for looping
GIS_OPT_INPUT=${GIS_OPT_INPUT//,/ }

# test for input raster map
eval `g.findfile element=cell file="$GIS_OPT_RAST"`
if [ ! "$file" ] ; then
    g.message -e "Raster map <$GIS_OPT_RAST> not found in mapset search path"
    exit 1
fi

for i in $GIS_OPT_INPUT; do
    # Check that input files exist
    if [ ! "$i" ] ; then 
	echo "file $i does not exist" >&2
	exit $E_NOTFOUND
    fi

    # Create header line if necessary
    skip=0
    if [ ! -z $GIS_FLAG_H ]; then
	skip=1
    fi

    # Intermediate variables
    map=`basename $i .csv`
    tmp=`tempfile -s ".csv"`
    of="${map}${GIS_OPT_SUFFIX}"
    col=`typeVars.awk -F, $i`

    # Get the raster type
    if `r.info -t $GIS_OPT_RAST | grep -q "FCELL"`; then
	rtype="double precision"
    else
	rtype="int"
    fi
    
    # Grass-specific calls
    v.in.ascii --overwrite input="$i" output="$map" fs="$GIS_OPT_FS" \
	skip="$skip" column="$col" x="$GIS_OPT_X" y="$GIS_OPT_Y"
    v.db.addcol map="$map" \
	columns="$GIS_OPT_RAST $rtype"
    v.what.rast vector="$map" column="$GIS_OPT_RAST" \
	raster="$GIS_OPT_RAST"
    v.out.ascii input="$map" fs="," \
	columns="$GIS_OPT_RAST" | \
    awk -v m=$GIS_FLAG_M -F, \
    '{
         # Convert Ft. elevation to meters if necessary
         if($4!="" && m==1) $4=$4 * 0.30480060
         print $4
    }' >"$tmp" 
    # Join raster values to input and write to output
    paste -d$GIS_OPT_FS "$i" "$tmp" >"$of"

    # Cleanup
    rm "$tmp"
    if [ $GIS_FLAG_R -eq 1 ] ; then
	g.remove -f --quiet vect=$map
    fi
done