# Calculating Lat/Lng X miles from point

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.

-
Are you looking for a tool that does this (like Esri's pe.dll) or a formula? – Kirk Kuykendall Feb 4 '11 at 19:06
Sorry I wasn't specific... I am looking for a formula. I'll update my question to be more specific. – Jason Whitehorn Feb 4 '11 at 19:08
Bunch of worked out math samples are here <a href="movable-type.co.uk/scripts/latlong.html">Calculate distance, bearing and more between Latitude/Longitude points</a> which includes the solution to "Destination point given distance and bearing from start point". – Stephen Quan Jan 24 '12 at 23:33
Closely related: gis.stackexchange.com/questions/2951/…. – whuber Jul 21 '12 at 17:07
here's the page that link to lat/long calculations [Lat/long calculations ](movable-type.co.uk/scripts/latlong.html) also this page Lat/long calculations there's a code + calculator – Abo gaseem Jan 16 '15 at 21:58

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.

-
Thank you for the quick reply. Let me digest some of this information, and I will get back with you. On the surface, though, this looks spot on. – Jason Whitehorn Feb 4 '11 at 20:21
I tried the direct case using pe.dll (actually libpe.so on solaris) after retrieving the distance and the forward azimuth from the html page for 32N,117W to 40N,82W. Using 32N,117W,d=3255.056515890041,azi=64.24498012065699, I got 40N,82W (actually 82.00000000064). – mkennedy Feb 4 '11 at 22:17
Awesome! Thank you very much for the link to the Aviation Formulary article by Ed Williams, I had not seen this before but it has thus far proven to be a great read. Just a note for anyone looking at this in the future, the inputs and outputs of this formula are ALL in radians, even the distance. – Jason Whitehorn Feb 5 '11 at 6:09
What is the unit of distance in this formula? – Karthik Murugan Apr 2 '15 at 6:35

If you were in a plane, then the point that is r meters away at a bearing of a degrees east of north is displaced by r * cos(a) in the north direction and r * sin(a) in the east direction. (These statements more or less define the sine and cosine.)

Although you are not in a plane--you're working on the surface of a curved ellipsoid that models the Earth's surface--any distance less than a few hundred kilometers covers such a small part of the surface that for most practical purposes it can be considered flat. The only remaining complication is that one degree of longitude does not cover the same distance as a degree of latitude. In a spherical Earth model, one degree of longitude is only cos(latitude) as long as a degree of latitude. (In an ellipsoidal model, this is still an excellent approximation, good to about 2.5 significant figures.)

Finally, one degree of latitude is approximately 10^7 / 90 = 111,111 meters. We now have all the information needed to convert meters to degrees:

The northwards displacement is r * cos(a) / 111111 degrees;

The eastwards displacement is r * sin(a) / cos(latitude) / 111111 degrees.

For example, at a latitude of -0.31399 degrees and a bearing of a = 30 degrees east of north, we can compute

``````cos(a) = cos(30 degrees) = cos(pi/6 radians) = Sqrt(3)/2 = 0.866025.
sin(a) = sin(30 degrees) = sin(pi/6 radians) = 1/2 = 0.5.
cos(latitude) = cos(-0.31399 degrees) = cos(-0.00548016 radian) = 0.999984984.
r = 100 meters.
east displacement = 100 * 0.5 / 0.999984984 / 111111 = 0.000450007 degree.
north displacement = 100 * 0.866025 / 111111 = 0.000779423 degree.
``````

Whence, starting at (-78.4437, -0.31399), the new location is at (-78.4437 + 0.00045, -0.31399 + 0.0007794) = (-78.4432, -0.313211).

A more accurate answer, in the modern ITRF00 reference system, is (-78.4433, -0.313207): this is 0.43 meters away from the approximate answer, indicating the approximation errs by 0.43% in this case. To achieve higher accuracy you must use either ellipsoidal distance formulas (which are far more complicated) or a high-fidelity conformal projection with zero divergence (so that the bearing is correct).

-
+1 for properly understanding the mathematical context (i.e., its locally plane) – Frank Conry Aug 10 '14 at 5:46

If you're interested in better accuracy, there's Vincenty. (Link is to the 'direct' form, which is exactly what you're after).

There are quite a few existing implementations, if you're after code.

Also, a question: You're not assuming the traveller maintains the same bearing for the entire trip, are you? If so, then these methods aren't answering the right question. (You'd be better off projecting to mercator, drawing a straight line, then un-projecting the result.)

-
Very good question, despite the wording in my question that may have indicated that I was calculating a destination for a traveler, I am not. Good point though. This was mainly so that I could calculate a bounding area (on a small order, say 50 miles)... I hope that makes sense. – Jason Whitehorn Feb 6 '11 at 3:45
gis.stackexchange.com/questions/3264/… had the same question (constructing bounding areas from a point & distance) – Dan S. Feb 7 '11 at 19:46

If you need a JavaScript solution consider these `functions` and this fiddle:

``````var gis = {
/**
* All coordinates expected EPSG:4326
* @param {Array} start Expected [lon, lat]
* @param {Array} end Expected [lon, lat]
* @return {number} Distance - meter.
*/
calculateDistance: function(start, end) {
var lat1 = parseFloat(start[1]),
lon1 = parseFloat(start[0]),
lat2 = parseFloat(end[1]),
lon2 = parseFloat(end[0]);

return gis.sphericalCosinus(lat1, lon1, lat2, lon2);
},

/**
* All coordinates expected EPSG:4326
* @param {number} lat1 Start Latitude
* @param {number} lon1 Start Longitude
* @param {number} lat2 End Latitude
* @param {number} lon2 End Longitude
* @return {number} Distance - meters.
*/
sphericalCosinus: function(lat1, lon1, lat2, lon2) {
var radius = 6371e3; // meters
var dLon = gis.toRad(lon2 - lon1),
distance = Math.acos(Math.sin(lat1) * Math.sin(lat2) +
Math.cos(lat1) * Math.cos(lat2) * Math.cos(dLon)) * radius;

return distance;
},

/**
* @param {Array} coord Expected [lon, lat] EPSG:4326
* @param {number} bearing Bearing in degrees
* @param {number} distance Distance in meters
* @return {Array} Lon-lat coordinate.
*/
createCoord: function(coord, bearing, distance) {
/** http://www.movable-type.co.uk/scripts/latlong.html
* φ is latitude, λ is longitude,
* θ is the bearing (clockwise from north),
* δ is the angular distance d/R;
* d being the distance travelled, R the earth’s radius*
**/
var

var φ2 = Math.asin(Math.sin(φ1)*Math.cos(δ) +
Math.cos(φ1)*Math.sin(δ)*Math.cos(θ));

var λ2 = λ1 + Math.atan2(Math.sin(θ) * Math.sin(δ)*Math.cos(φ1),
Math.cos(δ)-Math.sin(φ1)*Math.sin(φ2));
// normalise to -180..+180°
λ2 = (λ2 + 3 * Math.PI) % (2 * Math.PI) - Math.PI;

return [gis.toDeg(λ2), gis.toDeg(φ2)];
},
/**
* All coordinates expected EPSG:4326
* @param {Array} start Expected [lon, lat]
* @param {Array} end Expected [lon, lat]
* @return {number} Bearing in degrees.
*/
getBearing: function(start, end){
var
dLong = endLong - startLong;

var dPhi = Math.log(Math.tan(endLat/2.0 + Math.PI/4.0) /
Math.tan(startLat/2.0 + Math.PI/4.0));

if (Math.abs(dLong) > Math.PI) {
dLong = (dLong > 0.0) ? -(2.0 * Math.PI - dLong) : (2.0 * Math.PI + dLong);
}

return (gis.toDeg(Math.atan2(dLong, dPhi)) + 360.0) % 360.0;
},
toDeg: function(n) { return n * 180 / Math.PI; },
toRad: function(n) { return n * Math.PI / 180; }
};
``````

So if you want to calculate a new coordinate, it can be like so:

``````var start = [15, 38.70250];
var end = [21.54967, 38.70250];
var total_distance = gis.calculateDistance(start, end); // meters
var percent = 10;
// this can be also meters
var distance = (percent / 100) * total_distance;
var bearing = gis.getBearing(start, end);
var new_coord = gis.createCoord(icon_coord, bearing, distance);
``````
-

Here is a Python solution:

``````    def displace(self, theta, distance):
"""
Displace a LatLng theta degrees counterclockwise and some
meters in that direction.
Notes:
http://www.movable-type.co.uk/scripts/latlong.html
0 DEGREES IS THE VERTICAL Y AXIS! IMPORTANT!
Args:
theta:    A number in degrees.
distance: A number in meters.
Returns:
A new LatLng.
"""
theta = np.float32(theta)

return np.divide(np.dot(theta, np.pi), np.float32(180.0))

def to_degrees(theta):
return np.divide(np.dot(theta, np.float32(180.0)), np.pi)

lat2 = np.arcsin( np.sin(lat1) * np.cos(delta) +
np.cos(lat1) * np.sin(delta) * np.cos(theta) )

lng2 = lng1 + np.arctan2( np.sin(theta) * np.sin(delta) * np.cos(lat1),
np.cos(delta) - np.sin(lat1) * np.sin(lat2))

lng2 = (lng2 + 3 * np.pi) % (2 * np.pi) - np.pi

return LatLng(to_degrees(lat2), to_degrees(lng2))
``````
-

I got this working in ObjectiveC. The key here is knowing that you need lat/lng points in radians and you need to convert them back to lat/lng after applying the equation. Also, know that distance and tc are in radians.

Here's the original equation:

`````` 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
``````

Here it is implemented in ObjC where radian is a radian measured counter clockwise from N (e.g. PI / 2 is W, PI is S, 2 PI / 3 is E) and distance is in kilometers.

``````+ (CLLocationCoordinate2D)displaceLatLng:(CLLocationCoordinate2D)coordinate2D withRadian:(double)radian
withDistance:(CGFloat)distance {
double lat1Radians = coordinate2D.latitude * (M_PI / 180);
double lon1Radians = coordinate2D.longitude * (M_PI / 180);
double distanceRadians = distance / 6371;
double lon;
if (cos(lat) == 0) {
} else {
(float) (2 * M_PI)) - M_PI;
}
double newLat = lat * (180 / M_PI);
double newLon = lon * (180 / M_PI);
return CLLocationCoordinate2DMake(newLat, newLon);
}
``````
-

/*
I use the approach described below in determining the next coordinate given the bearing and distance from previous coordinate. I have problem on accuracy with other approach I read from the internet. I use this in determining the area of land, which is a polygon, and plot that polygon in google earth. A land Title has bearings and distances written in this manner: "NorthOrSouth x degrees y minutes EastOrWest, z meters to point n;". So, starting from the given coordinates of the reference point, I first compute the distance per one degree latitude and one degree longitude using haversine formula because this varies depending on location. Then I determine the next coordinate from the trigonometry sine and cosine formula. Below is the javascript:
*/
var mapCenter = new google.maps.LatLng(referencePointLatitude, referencePointLongitude); //the ref point lat and lon must be given, usually a land mark (BLLM)
var latDiv = latDiv(mapCenter); //distance per one degree latitude in this location
var lngDiv = lngDiv(mapCenter); //distance per one degree longitude in this location
var LatLng2 = NextCoord(PrevCoord,NorthOrSouth,x,y,EastOrWest,z); //next coordinate given the bearing and distance from previous coordinate
var Lat2 = LatLng2.lat(); //next coord latitude in degrees
var Lng2 = LatLng2.lng(); //next coord longitude in degrees
var polygon=[p1,p2,...,pn-1,pn,p1]; //p1,p2,etc. are the coordinates of points of the polygon, i.e. the land Title. Be sure to close the polygon to the point of beginning p1
var area = Area(polygon); //area of the polygon in sq.m.
function NextCoord(PrevCoord,NorthOrSouth,x,y,EastOrWest,z) {
var angle = ( x + ( y / 60 ) ) * Math.PI / 180;
var a = 1;
var b = 1;
if (NorthOrSouth == 'South') { a = -1; }
if (EastOrWest == 'West') { b = -1; }
var nextLat = PrevCoord.lat() + ( a * z * Math.cos(angle) / latDiv );
var nextLng = PrevCoord.lng() + ( b * z * Math.sin(angle) / lngDiv );
var nextCoord = new google.maps.LatLng(nextLat, nextLng);
return nextCoord;
}
function latDiv(mapCenter) {
var p1 = new google.maps.LatLng(mapCenter.lat()-0.5, mapCenter.lng());
var p2 = new google.maps.LatLng(mapCenter.lat()+0.5, mapCenter.lng());
return dist(p1,p2);
}
function lngDiv(mapCenter) {
var p3 = new google.maps.LatLng(mapCenter.lat(), mapCenter.lng()-0.5);
var p4 = new google.maps.LatLng(mapCenter.lat(), mapCenter.lng()+0.5);
return dist(p3,p4);
}
function dist(pt1, pt2) {
var dLat = ( pt2.lat() - pt1.lat() ) * Math.PI / 180;
var dLng = ( pt2.lng() - pt1.lng() ) * Math.PI / 180;
var a = Math.sin(dLat/2) * Math.sin(dLat/2) +
Math.sin(dLng/2) * Math.sin(dLng/2);
var R = 6372800; //earth's radius
var distance = R * 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));
return distance;
}
function Area(polygon) {
var xPts=[];
for (i=1; i<polygon.length; i++) { xPts[i-1] = ( polygon[i].lat() - polygon[0].lat() ) * latDiv; }
var yPts=[];
for (i=1; i<polygon.length; i++) { yPts[i-1] = ( polygon[i].lng() - polygon[0].lng() ) * lngDiv; }
var area = 0;
j = polygon.length-2;
for (i=0; i<polygon.length-1; i++) {
area = area + ( xPts[j] + xPts[i] ) * ( yPts[j] - yPts[i] );
j = i;
}
area = Math.abs(area/2);
return area;
}

-
The Cartesian formula for polygon area you attempt to apply here is not applicable to distances and angles computed on a curved surface (such as a spheroid). This formula makes an additional error by using latitude and longitude as if they were Cartesian coordinates. The only circumstances under which its use could be considered would be exactly those (for very small polygons) where the haversine formula is unnecessary anyway. Overall it appears this code works much too hard for no gain. – whuber Jun 28 '14 at 14:27

## protected by Community♦Oct 23 '15 at 12:10

Thank you for your interest in this question. Because it has attracted low-quality or spam answers that had to be removed, posting an answer now requires 10 reputation on this site.