Calculate compass heading between 2 UTM grid points

[skorchio]

I have been trying unsuccessfully to calculate the "bearing" between 2 Universal Transverse Mercator points. Lets call them head and tail.

head has points $XH = 756979.0 $YH = 7951269.2 (x,y)
tail has point $XT = 765484.9 $YT = 7951229.6 (x,y)

what I have been trying is as follows: (variables already defined...and chomped)

use Math::Trig;
$angle = atan(($YT - $YH) / ($XT - $XH));
$heading = deg2rad($angle);
printf ("heading approx. %.2f degrees\n", $heading);

But I am getting nonsense numbers back from this. I did a fair bit of googleing and installed Math::Trig, which provides the atan function and deg2rad function. atan returns the angle in radians, and the deg2rad I thought returned the heading between 2 points in degrees. Anyone tried this before and got it to work?

I also have tried the following, stolen from another person's script:

$angle = atan2($YL - $YF, $XL - $XF);
$angledeg = (180/3.1415926)*$angle;
if ( $angledeg < -90 ) {
        $headingdeg = abs((180/3.1415926)*$angle) + 90;
}
elsif ( $angledeg < 0 ) {
        $headingdeg = abs((180/3.1415926)*$angle) + 180;
}
elsif ( $angledeg < 90 ) {
        $headingdeg = 90 - (180/3.1415926)*$angle;
}
else {
        $headingdeg = 180 - (180/3.1415926)*$angle;
}
printf ("Line heading approx. %.2f degrees\n", $headingdeg);

but I am still getting the quadrant wrong. Its doing my head in.

[histrungalot]

Shouldn't it be rad2deg and not deg2rad? Also with atan the return value is in the range [-pi/2, pi/2]. Maybe use atan2 which the return value is in the range [-pi, pi], but you'll need to change to return [0,360). Start with an easy test like 0, 45, 180, 270 and 315 to see if you got it right. Test it!

#!/usr/bin/perl
use Math::Trig;
$XH = 1.0;
$YH = 7.2;
$XT = 0.0;
$YT = 6.2;
$angle = atan2(($YT - $YH) , ($XT - $XH));
#$heading = deg2rad($angle);
$heading = rad2deg($angle);
$heading += 360 if ( $heading < 0 );
printf ("heading approx. %.2f degrees\n", $heading);

$XH = 1.0;
$YH = 7.2;
$XT = 2.0;
$YT = 8.2;
$angle = atan2(($YT - $YH) , ($XT - $XH));
#$heading = deg2rad($angle);
$heading = rad2deg($angle);
$heading += 360 if ( $heading < 0 );
printf ("heading approx. %.2f degrees\n", $heading);

$XH = 756979.0;
$YH = 7951269.2;
$XT = 765484.9;
$YT = 7951229.6;
$angle = atan2(($YT - $YH) , ($XT - $XH));
#$heading = deg2rad($angle);
$heading = rad2deg($angle);
$heading += 360 if ( $heading < 0 );
printf ("heading approx. %.2f degrees\n", $heading);

Output:

heading approx. 225.00 degrees
heading approx. 45.00 degrees
heading approx. 359.73 degrees

[histrungalot]

I forgot the translation to north. Nee to add into above!

[histrungalot]

Not that anyone is going to read this (nothing else to do at 11:30pm but sleep, but why do that when I could learn something about UTM which I'll use ... but I digress), but it turns out that I needed to learn a few things about UTM. Below is the latest version that I think would actually work. The issue before is that I was not taking into account the zone and not handling the zero that could be possible with the denominator of the atan2. I think the easiest approach is to convert from UTM to lat/lon and then find the angle using the lat/lon. I pick Clarke 1866 but I don't think it matters, any conversion should work. I think (I'll spend more time thinking about that one). But change it if you think it will make a difference or change it and run some tests to see if it matters (as I typed it I figured, well check it out). I was correct it only changes the answer in the 1.0e-6.

#!/usr/bin/perl

# This script prints out the bearing in degrees between to point given in UTM format.
#
#
use strict;
use warnings;
use Geo::Coordinates::UTM;
use Math::Trig;

my ($x1,$y1,$x2,$y2) = (0,0,0,0);
my @names = ellipsoid_names;
# I don't think the ellipsoid makes a difference
($y1,$x1)=utm_to_latlon('clarke 1866','30V',512543.0,6406592.0);
($y2,$x2)=utm_to_latlon('clarke 1866','30V',514543.0,6406592.0);
my $ans1 = getAngleBetweenPoints($x1,$y1,$x2,$y2);
print "Bearing using \'$names[5]\':  $ans1\n";
($y1,$x1)=utm_to_latlon('Australian National','30V',512543.0,6406592.0);
($y2,$x2)=utm_to_latlon('Australian National','30V',514543.0,6406592.0);
my $ans2 = getAngleBetweenPoints($x1,$y1,$x2,$y2);
print "Bearing using \'$names[1]\':  $ans2\n";
print "Difference between \'$names[5]\' and \'$names[1]\':  ".($ans1-$ans2)."\n";

#
# Ported from:  http://stackoverflow.com/questions/642555/how-do-i-calculate-the-azimuth-angle-to-north-between-two-wgs84-coordinates
#
#
sub getAngleBetweenPoints {
   my ($X1,$Y1,$X2,$Y2) = @_;
   my ($dx,$dy,$result) = (($X2 - $X1),($Y2 - $Y1),0);

   if ($dx > 0) {
      $result = (pi*0.5) - atan2($dy,$dx);
   } 
   elsif ($dx < 0 ) {
      $result = (pi*1.5) - atan2($dy,$dx);
   }
   elsif ($dy > 0) {
      $result = 0;
   }
   elsif ($dy < 0) {
      $result = pi; 
   }
   else {
      $result = -999;
   }
   $result = rad2deg($result) if ( $result != -999 );
   return $result;
}

Output:

Bearing using 'Clarke 1880':  90.1029398115509
Bearing using 'Australian National':  90.1029387649392
Difference between 'Clarke 1880' and 'Australian National':  1.0466117288388e-06