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));

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.

## All 3 Replies

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 += 360 if ( \$heading < 0 );

\$XH = 1.0;
\$YH = 7.2;
\$XT = 2.0;
\$YT = 8.2;
\$angle = atan2((\$YT - \$YH) , (\$XT - \$XH));
\$heading += 360 if ( \$heading < 0 );

\$XH = 756979.0;
\$YH = 7951269.2;
\$XT = 765484.9;
\$YT = 7951229.6;
\$angle = atan2((\$YT - \$YH) , (\$XT - \$XH));
\$heading += 360 if ( \$heading < 0 );

Output:

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

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

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\':  \$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\':  \$ans2\n";
print "Difference between \'\$names\' and \'\$names\':  ".(\$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``````
Be a part of the DaniWeb community

We're a friendly, industry-focused community of developers, IT pros, digital marketers, and technology enthusiasts meeting, learning, and sharing knowledge.