[GRASS-user] How to find out an angle between to points on the map (for r.plane, r.lake) - shellscript

[Philipp Steigenberger]

Dylan Beaudette schrieb:
> On Thursday 28 February 2008, Philipp Steigenberger wrote:
>> Hi
>> I want to work with r.plane and r.lake to get simulate floodings in DTMs.
>>
>> In the section of interest there is a river which runs from about SSW to
>> NNE. For r.plane I need the
>> azimuth. Is there a tool in GRASS which tells me the angle between two
>> points I mark on the map like I can measure in between two points with
>> d.measure?
>>
>>
> 
> Simple trigonometry is a start. Or on a computer the atan2() function is 
> useful:
> 
> atan2( (end_y - start_y) , (end_x - start_x) ) * 180/pi


Dylan Beaudette schrieb:

> 
> attached is a patch (against todays SVN) to optionally print the bearing 
> between clicks.
> 
> Dylan
> 
> 
> 

Dylan,
sorry, I have no idea how to incllude this patch to my system.
I'm neither good in R or C or Fortran or Python.
> attached is a patch (against todays SVN
(Also I have problems to install the actual SVN)

But I'm used to write shell scrpits, so I wrote this one for bc.

#!/bin/sh
LIST_1=`d.what.rast -t1`
E1=`echo $LIST_1 | awk -F : '{print $1}'`
N1=`echo $LIST_1 | awk -F : '{print $2}'`

LIST_2=`d.what.rast -t1 --q`
E2=`echo $LIST_2 | awk -F : '{print $1}'`
N2=`echo $LIST_2 | awk -F : '{print $2}'`

SCALE="scale=8"
pi=$(echo "scale=10; 4 * a ( 1 )" | bc -l)

echo "epsilon = `echo "$SCALE;  360-( a ( ( $E1 - $E2 ) / ( $N1 - $N2 )
) ) * 180 / $pi" | bc -l`°"
echo "if epsilon > 360"
echo "epsilon' = `echo "$SCALE;  ( a ( ( $E1 - $E2 ) / ( $N1 - $N2 ) ) )
* 180 / $pi * -1" | bc -l`°"

In GRASS I work in Transversal Mercator (Gauß-Krüger)
I know the script isn't perfect, cause it doesn't work if N1 = N2, but
for me this doesn't matter.

Philipp