Calculating Latitude/Longitude X miles from point?
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.
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).
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),
lat1 = gis.toRad(lat1),
lat2 = gis.toRad(lat2),
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
radius = 6371e3, // meters
δ = Number(distance) / radius, // angular distance in radians
θ = gis.toRad(Number(bearing));
φ1 = gis.toRad(coord[1]),
λ1 = gis.toRad(coord[0]);
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
startLat = gis.toRad(start[1]),
startLong = gis.toRad(start[0]),
endLat = gis.toRad(end[1]),
endLong = gis.toRad(end[0]),
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);