|
stevieb has asked for the wisdom of the Perl Monks concerning the following question:
I'm prototyping some code that I'll eventually port over to C for an Arduino project I've been working on (much quicker to prototype in Perl first) and am hoping for some feedback from the mathematicians.
So the code is designed to calculate a direction in compass degrees between two sets of GPS coordinates, then convert that degree into a compass direction. The following example is between Calgary, AB, Canada to Toronto, ON, Canada. From what I've found on the Internet, this direction is approximately 94 degrees, which the code outputs correctly.
I completely comprehend and had no problem at all writing the direction() function, but had to do a lot of reading and guessing with the trig stuff in bearing() as it's much beyond my mathematical understanding, and am wondering if it appears reasonably sane.
use warnings;
use strict;
use Math::Trig;
use constant PI => 3.14159265358979;
# calgary
my $lat1 = 51.32055555555556;
my $lon1 = -114.7297222222222221;
# toronto
my $lat2 = 43.65321111111111;
my $lon2 = -79.38321111111111;
my $bearing = bearing($lat1, $lon1, $lat2, $lon2);
my $direction = direction($bearing);
print "$bearing, $direction\n";
sub bearing {
my ($lat1, $lon1, $lat2, $lon2) = map { ($_ * PI) / 180 } @_;
my $d_lon = $lon2 - $lon1;
my $y = sin($d_lon) * cos($lat2);
my $x = cos($lat1) * sin($lat2)
- sin($lat1)
* cos($lat2) * cos($d_lon);
my $rad = atan2($y, $x);
my $deg = $rad * (180 / PI);
return $deg < 0 ? $deg += 360 : $deg;
}
sub direction {
my ($deg) = @_;
my @directions = qw(
N NNE NE ENE E ESE SE SSE S SSW SW WSW W WNW NW NNW N
);
my $calc = (($deg % 360) / 22.5) + .5;
return $directions[$calc];
}
__END__
94.0048844699519, E
|
Re: Calculate bearing between GPS coordinates
by RonW (Parson) on May 31, 2017 at 21:50 UTC
|
Your code looks reasonable and should work fine when translated to C.
Side note (just FYI): You are calculating the true heading, not the compass heading. A compass, whether magnetic or inertial (such as a gyroscope), will give you a different reading.
In the case of a magnetic compass, magnetic north and south are not aligned with the Earth's rotational north and south. Also, the magnetic north and south vary over time. The correction required depends on your location relative to the current magnetic north and south.
An inertial compass can be calibrated to the true heading, but accuracy decreases the further west or east of where the compass was calibrated. I don't remember how to calculate the correction.
(I learned about this back as a teen when I was a cadet in the Civil Air Patrol.)
| [reply] |
|
|
Thanks RonW!
I was a bit vague as I didn't want to be pedantic about details. What you're referring to is "magnetic declination" and I am aware of that, and understand it's not being considered here.
My main purpose for this little contraption (it's going to end up looking something like this) is for when I'm on my multi-kilometre hikes in the mountains, I turn the device on, press a button that writes to EEPROM the current lat/lon, then turn it off and throw it in my backpack.
In the extraordinary event that my land orienteering fails me with my map/compass due to exhaustion or lack of visibility because I'm in a large stand of 300' cedar trees and I get lost, I pull the unit back out, turn it back on, and press another button to "go back to the last save point". It'll then display the degrees and direction and distance I need to travel to get back to home base, which I'll use along with my compass if need be.
For these cases, if I'm only hiking say 30-40 km, declination corrections not being preset shouldn't be much of an issue. If I find that it is during trial runs, I'll calculate that into the software and go from there :)
I've already got duplicates of all the hardware I need (Arduino Trinket Pro, OLED screen, batteries (I'm going to test several), and the GPS sensor itself. I'm also adding in for convenience after the initial rig is completed and tested, micro USB charging port for the battery (so I can connect it to my portable solar rig to charge), an FTDI interface for the microproc so I can easily connect via serial for troubleshooting new code, and an SD card breakout so I'll be able to actually save my hiking paths if I ever so desire.
| [reply] |
|
|
For these cases, if I'm only hiking say 30-40 km, declination corrections not being preset shouldn't be much of an issue. If I find that it is during trial runs, I'll calculate that into the software and go from there :)
While i suspect that in most cases you will want to take the local deviation into account, in most cases you should also be able to get a general idea of it from your maps or just from memory.
But in specific cases that may not be enough as local geological factors such as nearby magnetic materials can cause extra and sometimes large changes. (see https://www.ngdc.noaa.gov/geomag/EMM/img/dec_WMM_EMM.jpg via https://www.ngdc.noaa.gov/geomag/EMM/index.html. There are places in the USA where the deviation is 15 degrees.
Being somewhat of a fan of geological tidbits i kinda remembered running into a site that had a database. While it is not what i remember https://www.ngdc.noaa.gov/geomag/emag2.html seems to have a world grid database available that you may find of interest to subset and somehow contain on your deviceEdit:wrong data, that was not declination, but still an interesting dataset.
| [reply] |
|
|
| [reply] |
|
|
|
|
This is largely academic due to the short distances involved, but... If you follow a constant compass bearing, you're not on a great circle path, you're following a rhumb line. For your Calgary-to-Toronto example, the great circle starts off at a bearing of 94 deg, but you need to adjust that along the way, so when you get to the destination, you're headed at 120.5 deg. The rhumb line is a constant 107.8 degrees the entire way. The rhumb bearing is calculated as follows:
my $bearing = atan2($lon2 - $lon1, aths($lat2) - aths($lat1));
sub aths { # atanh(sin(x))
my ($x) = @_;
my $y = sin($x);
return 0.5*log((1 + $y) / (1 - $y));
}
| [reply]
[d/l] |
|
|
So where are you going to hike?
In Canada?
| [reply] |
|
|
Re: Calculate bearing between GPS coordinates
by Anonymous Monk on May 31, 2017 at 20:23 UTC
|
Since you're using Math::Trig anyway, why don't you use its great_circle_direction() and rad2deg()? | [reply] |
|
|
I didn't even know about great_circle_direction() :)
As far as rad2deg() and deg2rad(), I had them in, but wanted to learn as many of the actual mathematical equations as I possibly could for those functions so I understood what was happening.
| [reply]
[d/l]
[select] |
Re: Calculate bearing between GPS coordinates
by pryrt (Abbot) on Jun 01, 2017 at 22:01 UTC
|
Your formulas have already been confirmed. But as for the meaning behind the trig: "Ask Dr. Math" has a derivation of that bearing forumla at http://mathforum.org/library/drmath/view/55417.html in the "Date: 06/19/2002 at 08:45:23" message. After reading it a few times, the steps I see:
- Given the points on the sphere NorthPole, Location1=(LON1, LAT1), and Location2=(LON2, LAT2)
- In general, a (LON,LAT) pair turns into (x,y,z) = ( cos(LAT)*cos(LON), cos(LAT)*sin(LON), sin(LAT)*1 ) on a sphere
→ we're wrongly assuming earth is a sphere, for the sake of calculation: it's close enough that you won't accidentally end up "unemployed, in Greenland"</Vizzini>)
- Rotate the earth (or, easier, your mathematical frame of reference) so that Location1=(0, LAT1), and Location2=(LON2-LON1, LAT2)=(DLON, LAT2)
→ this allows you to drop out a couple sines and cosines by knowing that cos(0)=1 and sin(0)=0
- Your coordinates become
N = <NORTH> = ( 0 , 0 , 1 )
A = <LOC#1> = ( cos(LAT1), 0, sin(LAT1) ) → this is greatly simplified by LON1 = 0
B = <LOC#2> = ( cos(LAT2)*cos(DLON), cos(LAT2)*sin(DLON), sin(LAT2) )
- Create two normal vectors for two planes (one plane thru N and A; the other thru A and B) by taking the cross products p=NxA and q=BxA (the direction of q keeps the bearing oriented properly)
p = NxA = (0, cos(LAT1), 0)
q = BxA = (sin(LAT1)*cos(LAT2)*sin(DLON),
cos(LAT1)*sin(LAT2)-sin(LAT1)*cos(LAT2)*cos(DLON),
-cos(LAT1)*cos(LAT2)*sin(DLON))
- The angle between these planes is your bearing. Normally that angle would be from cos(BEARING) = p DOT q / (|p|*|q|). But because NxA is all in the y-direction of your rotated earth, you can simplify to the arctangent function, already used in your code.
| [reply]
[d/l]
[select] |
Re: Calculate bearing between GPS coordinates
by stevieb (Canon) on Jun 03, 2017 at 21:29 UTC
|
So this is Arduino related and not really about Perl, but all the code I've been hacking on prototyping in Perl and porting to C has gotten me a lot of progress. I suppose it is Perl related to some degree, as if it wasn't for me knowing Perl and all the experience I've gained about writing C/XS, I wouldn't be here.
I have the GPS connected to my Trinket Pro on battery, with a GPS unit and a 0.96", 128x32 pixel OLED screen. I also have a button to switch between the "home" display contents, and the "return" display contents, and a button that writes to EEPROM that saves the coordinates of the current location.
The idea is that I can press the "save" button before going on a long hike, and turn the unit off. Then, if I ever get lost, I turn it back on, press the "return screen" button, and it guides me back to my starting location. Here are the screen layouts (in all cases I obscure some of the numbers as to not give out my home location :) Also, all coordinates have been shaved to 5 decimal points which is one metre. We don't need to be more precise than that, so it allowed me to shave some screen area) All of the N|W locations are generated dynamically based on whether the headings are negative or positive:
Home display. Top line is the current date and time in UTC. Second line is fixed-width north heading along with current heading (direction). Third line is fixed-width east heading, followed by the number of satellites we're currently locked onto. Fourth line is altitude, followed by current speed in KM/h.
17.06.03 21:17:16 UTC
N:051.11111 H:288
W:114.11111 s:6
A:1090 S:0.00
Return screen. First line is fixed-width north point of the destination we need to return to, followed by the distance in Km. Second line is the fixed-width west point we're returning to, followed by the direction in degrees we need to travel. Third line is the fixed-width north point that we're currently located followed by our current heading in degrees, and fourth line is our fixed-width west point followed by our current speed of travel in Km/h.
N:050.95355 K:38.687
W:114.58056 D:244.96
N:051.11111 H:288.07
W:114.11111 S:0.00
I have manually saved the coordinates of a location I'm driving to tomorrow (Bragg Creek, Alberta, Canada), and am going to give it a test run on that driving distance. Later this week, I will go for a hike about 8km and test the actual intended functionality. For now, I'm just going to do code cleanup, fix up a couple of tiny things, and if it all works after my tests this week, I'll be contemplating how I want to package it, and solder everything together.
I have not converted the awesome magnetic declination code into C yet, so the direction has not been calibrated. I'll hack together that code after the unit is packaged and field ready, before my first real life tests of it. | [reply]
[d/l]
[select] |
|
|
