Distance Calculation – Calculating Latitude/Longitude X Miles from Point

algorithmdistancehaversinespherical-geometry

I am wanting to find a latitude and longitude point given a bearing, a distance, and a starting latitude and longitude.

This appears to be the opposite of this question (Distance between lat/long points).

I have already looked into the haversine formula and think it's approximation of the world is probably close enough.

I am assuming that I need to solve the haversine formula for my unknown lat/long, is this correct? Are there any good websites that talk about this sort of thing? It seems like it would be common, but my googling has only turned up questions similar to the one above.

What I am really looking for is just a formula for this. I'd like to give it a starting lat/lng, a bearing, and a distance (miles or kilometers) and I would like to get out of it a lat/lng pair that represent where one would have ended up had they traveled along that route.

Best Answer

I'd be curious how results from this formula compare with Esri's pe.dll.

(citation).

A point {lat,lon} is a distance d out on the tc radial from point 1 if:

 lat=asin(sin(lat1)*cos(d)+cos(lat1)*sin(d)*cos(tc))
 IF (cos(lat)=0)
    lon=lon1      // endpoint a pole
 ELSE
    lon=mod(lon1-asin(sin(tc)*sin(d)/cos(lat))+pi,2*pi)-pi
 ENDIF

This algorithm is limited to distances such that dlon < pi/2, i.e those that extend around less than one quarter of the circumference of the earth in longitude. A completely general, but more complicated algorithm is necessary if greater distances are allowed:

 lat =asin(sin(lat1)*cos(d)+cos(lat1)*sin(d)*cos(tc))
 dlon=atan2(sin(tc)*sin(d)*cos(lat1),cos(d)-sin(lat1)*sin(lat))
 lon=mod( lon1-dlon +pi,2*pi )-pi

Here's an html page for testing.

Related Solutions

Trilateration Latitude Longitude – Using 3 Points and Distances for Trilateration

After some looking around at Wikipedia and the same question/answer at StackOverflow, I figured I would take a stab at it, and try to fill in the gaps.

First off, Not sure where you got the output, but it appears to be wrong. I plotted the points in ArcMap, buffered them to the distances specified, ran intersect to on the buffers, and then captured the vertex of intersection to get the solutions. Your proposed output is the point in green. I calculated the value in the callout box, which is about 3 meters of what ArcMap gave for solution derived from the intersect.

alt text

The math on the wikipedia page isn't too bad, just need to covert your geodetic coordinates to the cartesian ECEF, which can be found here. the a/x +h terms can be replaced by the authalic sphere radius, if you aren't using an ellipsoid.

Probably easiest just give you some well(?) documented code, so here it is in python

import math
import numpy

#assuming elevation = 0
earthR = 6371
LatA = 37.418436
LonA = -121.963477
DistA = 0.265710701754
LatB = 37.417243
LonB = -121.961889
DistB = 0.234592423446
LatC = 37.418692
LonC = -121.960194
DistC = 0.0548954278262

#using authalic sphere
#if using an ellipsoid this step is slightly different
#Convert geodetic Lat/Long to ECEF xyz
#   1. Convert Lat/Long to radians
#   2. Convert Lat/Long(radians) to ECEF
xA = earthR *(math.cos(math.radians(LatA)) * math.cos(math.radians(LonA)))
yA = earthR *(math.cos(math.radians(LatA)) * math.sin(math.radians(LonA)))
zA = earthR *(math.sin(math.radians(LatA)))

xB = earthR *(math.cos(math.radians(LatB)) * math.cos(math.radians(LonB)))
yB = earthR *(math.cos(math.radians(LatB)) * math.sin(math.radians(LonB)))
zB = earthR *(math.sin(math.radians(LatB)))

xC = earthR *(math.cos(math.radians(LatC)) * math.cos(math.radians(LonC)))
yC = earthR *(math.cos(math.radians(LatC)) * math.sin(math.radians(LonC)))
zC = earthR *(math.sin(math.radians(LatC)))

P1 = numpy.array([xA, yA, zA])
P2 = numpy.array([xB, yB, zB])
P3 = numpy.array([xC, yC, zC])

#from wikipedia
#transform to get circle 1 at origin
#transform to get circle 2 on x axis
ex = (P2 - P1)/(numpy.linalg.norm(P2 - P1))
i = numpy.dot(ex, P3 - P1)
ey = (P3 - P1 - i*ex)/(numpy.linalg.norm(P3 - P1 - i*ex))
ez = numpy.cross(ex,ey)
d = numpy.linalg.norm(P2 - P1)
j = numpy.dot(ey, P3 - P1)

#from wikipedia
#plug and chug using above values
x = (pow(DistA,2) - pow(DistB,2) + pow(d,2))/(2*d)
y = ((pow(DistA,2) - pow(DistC,2) + pow(i,2) + pow(j,2))/(2*j)) - ((i/j)*x)

# only one case shown here
z = numpy.sqrt(pow(DistA,2) - pow(x,2) - pow(y,2))

#triPt is an array with ECEF x,y,z of trilateration point
triPt = P1 + x*ex + y*ey + z*ez

#convert back to lat/long from ECEF
#convert to degrees
lat = math.degrees(math.asin(triPt[2] / earthR))
lon = math.degrees(math.atan2(triPt[1],triPt[0]))

print lat, lon

Coordinate System Algorithm – Offset Latitude/Longitude by Some Amount of Meters

If your displacements aren't too great (less than a few kilometers) and you're not right at the poles, use the quick and dirty estimate that 111,111 meters (111.111 km) in the y direction is 1 degree (of latitude) and 111,111 * cos(latitude) meters in the x direction is 1 degree (of longitude).

Best Answer

Related Solutions

Trilateration Latitude Longitude – Using 3 Points and Distances for Trilateration

Coordinate System Algorithm – Offset Latitude/Longitude by Some Amount of Meters

Related Question