#!/bin/bash
#
# Dylan Beaudette 2008-03-26
#
# 
# example:
# r.angle  elev=elev10m point=p out=ang
# 
# it would be better to compute the percent of the LOS region >= alpha at each cell
# ... but this would also be based on the larger adjacent terrain features
# 
# 
# 
#%Module
#% description: Compute the vertical angle from some (x,y,z) coordinate to all cells in the current region.
#%End
#%option
#% key: elev
#% type: string
#% gisprompt: old,cell,raster
#% description: Name of input elevation raster
#% required : yes
#%END


#%option
#% key: point
#% type: string
#% gisprompt: vector
#% description: Name of a new vector point file to store the reference point.
#% required : yes
#%END

#%option
#% key: output
#% type: string
#% gisprompt: old,cell,raster
#% description: Name of output raster
#% required : yes
#%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

# check if we have awk
AWK=`which awk`
if [ "$AWK" = "" ] ; then
    echo "$PROG: awk required, please install awk/gawk first" 1>&2
    exit 1
fi

# setting environment, so that awk works properly in all languages
unset LC_ALL
LC_NUMERIC=C


# #borrowed from Hamish's r.out.gmt
# # what to do in case of user break:
# exitprocedure()
# {
#     echo "User break!" 1>&2
#     cleanup
#     exit 1
# }
# # shell check for user break (signal list: trap -l)
# trap "exitprocedure" 2 3 15
# 
# cleanup()
# {
#     if [ "$MAP_TYPE" = "DCELL" ] ; then
# 	g.remove $GIS_OPT_output
#     fi
# }


d.where --q -1 | awk '{print $1"|"$2}' | v.in.ascii --q out=$GIS_OPT_point --o

# get x and y coordinates
x=`v.out.ascii --q $GIS_OPT_point | awk -F "|" '{print $1}'`
y=`v.out.ascii --q $GIS_OPT_point | awk -F "|" '{print $2}'`


# get the elevation at that point:
z=`r.what --q in=$GIS_OPT_elev east_north=$x,$y  | awk -F "|" '{print $4}'`

# compute the vertical angle of all surrounding cells
# r.mapcalc "$GIS_OPT_output = if( atan( ($GIS_OPT_elev - $z) / sqrt( (x() - $x)^2 + (y() - $y)^2) ) >= 10, 1, null() )"
r.mapcalc "$GIS_OPT_output = atan( ($GIS_OPT_elev - $z) / sqrt( (x() - $x)^2 + (y() - $y)^2) )"

# give a useful color scale
# pos angles (obstructing the sky) are red
r.colors --q map=ang color=differences

d.erase
d.his --q i=shade h=$GIS_OPT_output
d.legend --q $GIS_OPT_output col=white

d.vect $GIS_OPT_point fcol=yellow col=black icon=basic/circle size=12
d.vect w_smooth type=boundary