Has anyone done the 'Haversine' formula?

[Andre]

Thanks for all this wonderful codes, which I'm successfully using to calculate the distance between two geographic points!

But can you show me the other way around, meaning I have one point (latitude,longitude) and a specific distance (in km) around this point. How can I calculate the most-northern/southern latitude and most-western/eastern longitude matching this distance from the given point?

(I've seen a lot formulas on internet, but don't get it... latitude seems easy, but longitude... )

Thanks in advance!

[Paul]

Code: Select all

Define.d lat1,lon1,Distance,Bearing,Radius=6372.795477598
Macro DistFromPoint (lat1,lon1,Distance,Bearing)
	lat2.d=Degree(ASin(Sin(Radian(lat1))*Cos(Distance/Radius)+Cos(Radian(lat1))*Sin(Distance/Radius)*Cos(Radian(Bearing)) ))
	lon2.d=lon1+Degree(ATan2(Cos(Distance/Radius)-Sin(Radian(lat1))*Sin(Radian(lat2)),Sin(Radian(Bearing))*Sin(Distance/Radius)*Cos(Radian(lat1)) ))
EndMacro

DistFromPoint(49,-94,20,0) ;<-- Most northern point 20km from lat/lon
Debug lat2
Debug lon2

DistFromPoint(49,-94,20,180) ;<-- Most southern point 20km from lat/lon
DistFromPoint(49,-94,20,90) ;<-- Most eastern point 20km from lat/lon
DistFromPoint(49,-94,20,270) ;<-- Most western point 20km from lat/lon

[Andre]

Thank you very much, Paul!

I tried to encapsulate the complete function and variables a bit more, and extended the examples.
Hopefully all is done the right way and the results are true
(e.g I'm not sure, why the latitude changes a bit in the 3rd example..., or if the +/- of the resulting lat/long are always handled correctly...)

My own test inside my project will take a bit...

[Update: 9th Feb. 2019]

Code: Select all

Define.d Distance, Bearing

Structure geopoint
  Latitude.d
  Longitude.d
EndStructure
Define pos.geopoint

#MostNorthern = 0
#MostSouthern = 180
#MostEastern  = 90
#MostWestern  = 270

Procedure DistanceFromPoint(*pos.geopoint, Distance, Bearing)
  ; This function calculates the new geographic point in the 'Distance' (of xxx km) from the
  ; given current geographic point '@pos'.
  ; 'Bearing' determines the direction, where the new position should be calculated:
  ;   #MostNorthern and #MostSouthern will give a new latitude value
  ;   #MostEastern and #MostWestern will give a new longitude value
  Protected Radius.d = 6372.795477598  ;  a fixed value
  Protected OldLatitude.d = *pos\Latitude, OldLongitude.d = *pos\Longitude
  Protected NewLatitude.d, NewLongitude.d
  
  ; Now we use the given values (old Lat/Long position, Distance and Bearing/Direction) for calculating the new position:
  NewLatitude  = Degree(ASin(Sin(Radian(OldLatitude))*Cos(Distance/Radius)+Cos(Radian(OldLatitude))*Sin(Distance/Radius)*Cos(Radian(Bearing)) ))
  NewLongitude = OldLongitude + Degree(ATan2(Cos(Distance/Radius)-Sin(Radian(OldLatitude))*Sin(Radian(NewLatitude)),Sin(Radian(Bearing))*Sin(Distance/Radius)*Cos(Radian(OldLatitude))))
  
  ; Return the new position by saving it in the structured variable (which can be accessed in the main program):
  *pos\Latitude  = NewLatitude
  *pos\Longitude = NewLongitude
EndProcedure

; 1st example - going to the North:
pos\Latitude = 49
pos\Longitude = -94
Distance = 200
Debug "Old latitude/longitude: " + pos\Latitude + " / " + pos\Longitude
DistanceFromPoint(@pos, Distance, #MostNorthern) ;<-- Most northern point 200 km from lat/lon
Debug "..new latitude/longitude " + Distance + " km into North: " + pos\Latitude + " / " + pos\Longitude
Debug "----------------"

; 2nd example - going to the South:
pos\Latitude = 49
pos\Longitude = -94
Distance = 20
Debug "Old latitude/longitude: " + pos\Latitude + " / " + pos\Longitude
DistanceFromPoint(@pos, Distance, #MostSouthern) ;<-- Most southern point 20 km from lat/lon
Debug "..new latitude/longitude " + Distance + " km into South: " + pos\Latitude + " / " + pos\Longitude
Debug "----------------"

; 3rd extended example - going from Northpole (latitude = +90) to Southpole (latitude = -90) and calculate the new longitude for going 200 km to the East:
Debug "Here a more extended example: starting point moves from North to South (latitude decreases 10 degress everytime), but with constant longitude (20):"
For a = 85 To -85 Step -10
  pos\Latitude = a
  pos\Longitude = 20
  Distance = 200
  Debug "Old latitude/longitude: " + pos\Latitude + " / " + pos\Longitude
  DistanceFromPoint(@pos, Distance, #MostEastern) ;<-- Most eastern point 200 km from lat/lon
  Debug "..new latitude/longitude " + Distance + " km into East: " + pos\Latitude + " / " + pos\Longitude
Next
Debug "----------------"

; 4th example - going to the West:
pos\Latitude = 49
pos\Longitude = -94
Distance = 20
Debug "Old latitude/longitude: " + pos\Latitude + " / " + pos\Longitude
DistanceFromPoint(@pos, Distance, #MostWestern) ;<-- Most western point 20 km from lat/lon
Debug "..new latitude/longitude " + Distance + " km into West: " + pos\Latitude + " / " + pos\Longitude

[Josh]

... latitude seems easy, but longitude...

If you take a closer look, you will find that calculating the latitude is probably the bigger problem.

For example, if you want to drive 2000 km to the west, you have to follow the great circle to really reach the farthest point to the west. To follow the great circle means that on the northern hemisphere you have to go north first and then south again after half of your distance. In the southern hemisphere the opposite is true, at the equator it doesn't matter.

At GPS starting time I wrote a simple navigation program (shit, I'm old). I always calculated the globe to be in the middle of the circle. So all calculations were much easier. I haven't thought it over for your case right now, but maybe this can help.

[Little John]

A riddle:
You are somewhere on earth. You go 5 km south, then 5 km west and then 5 km north. Then you have arrived back at the place where you started.
You see a bear. What color is it?

[Josh]

You see a bear. What color is it?

I have never seen a polar bear so far north

[Paul]

You see a bear. What color is it?

I have never seen a polar bear so far north

What about Santa ?

[Michael Vogel]

Tried to do a quick check, don't think, the results are okay (maybe I am mixing bears and penguins here)...

Here's the output...
Old latitude/longitude: 49 / -94
New latitude/longitude 20 km into West: 48.9996754155 / -94.4177635481
30.568590932620445

...and here's the code:

Code: Select all

 :
pos\Latitude = 49
pos\Longitude = -94
new.geopoint=pos
Distance = 20
DistanceFromPoint(@new, Distance, #MostWestern) ;<-- Most western point 20 km from lat/lon
Debug "Old latitude/longitude: " + pos\Latitude + " / " + pos\Longitude
Debug "New latitude/longitude " + Distance + " km into West: " + new\Latitude + " / " + new\Longitude
Debug distVincenty2(pos\Latitude,pos\Longitude,new\Latitude,new\Longitude)

[Paul]

When Andre converted my snippet, he missed something...

Code: Select all

  NewLongitude = OldLongitude + Degree(ATan2(Cos(Distance/Radius)-Sin(Radian(OldLatitude))*Sin(Radian(NewLatitude)),
  Sin(Radian(Bearing))*Sin(Distance/Radius)*Cos(Radian(lat1))))
                                                        ^^^

He renamed all the 'lat1' to 'OldLongitude' except this one

[Andre]

Thanks Paul, for pointing this out

I've updated my code snippet above, including a slightly modified 3rd example.

But I still think too, that moving to the West (4th example) shouldn't change the latitude, when it doesn't for moving to the East (2nd example)!?

[Paul]

This is a good read...
https://www.movable-type.co.uk/scripts/latlong.html

Explains why final bearing differs from the initial heading by varying degrees according to distance and latitude.


(calculations would probably be much easier if everyone would just accept the 'flat earth' theory )

[Andre]

This is a good read...
https://www.movable-type.co.uk/scripts/latlong.html

...

(calculations would probably be much easier if everyone would just accept the 'flat earth' theory )

Thanks for the link, Paul.

But like other links too, this math section is too high for me....
So I'm someone for the 'flat earth'....

I will post again, when I can report more about testing the code in my project (will take a while).