177

I'm looking for an algorithm which when given a latitude and longitude pair and a vector translation in meters in Cartesian coordinates (x,y) would give me a new coordinate.

Sort of like a reverse Haversine. I could also work with a distance and a heading transformation, but this would probably be slower and not as accurate. Ideally, the algorithm should be fast as I'm working on an embedded system. Accuracy is not critical, within 10 meters would be good.

2
  • So you'd be fine modeling the earth as a sphere?
    – underdark
    Commented Oct 26, 2010 at 22:54
  • 1
    Yeah, that would be fine as I'm expecting <1km offsets.
    – Thomas O
    Commented Oct 26, 2010 at 22:58

9 Answers 9

167

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).

12
  • 6
    @Thomas: Actually, you can be very close to the poles. I checked against a UTM calculation using equal x- and y-displacements of 1400 m (so the total displacement is 2 km). The results are good to 8.6 meters or better. The worst latitude (for this direction and amount of displacement) is 81 degrees: the approximation actually gets more accurate as you move north and its error stays below 10 meters until you get beyond 89.6 degrees!
    – whuber
    Commented Oct 27, 2010 at 21:36
  • 95
    Incidentally, these magic numbers of 111,111 are easy to remember by knowing some history: the French originally defined the meter so that 10^7 meters would be the distance along the Paris meridian from the equator to the north pole. Thus, 10^7 / 90 = 111,111.1 meters equals one degree of latitude to within the capabilities of French surveyors two centuries ago.
    – whuber
    Commented Oct 27, 2010 at 21:38
  • 6
    So with the formula if I wanted to move +100m in the y direction from say 10.0 N, 10.0 E, would I just add 100 / 111111? If moving in the x direction +100m, would it be 100÷(111,111×(cos 10))? Just making sure I've got this right.
    – Thomas O
    Commented Oct 27, 2010 at 21:47
  • 8
    @Thomas Yes, that's right. Note how the second formula expands the apparent x-displacement (by virtue of dividing by a number less than 1) as it should, because a degree of longitude gets smaller as you move towards the poles from the equator. The only potential hitch is to make sure you and your software platform agree on what "cos" means: it had better interpret cos(10) as the cosine of 10 degrees, not 10 radians! (If not, 10 degrees = 10 * pi / 180 radians illustrates the simple conversion.) At this point, the code offered by @haakon_d should make complete sense to you.
    – whuber
    Commented Oct 27, 2010 at 21:54
  • 13
    Somebody attempted to edit this answer to replace "meters" by "km." They probably were reading the comma "," in the European sense of a decimal point. I follow the American convention (which I believe is the convention of international publications, too) of using a comma to separate long digit strings into groups of three and a decimal point "." instead of the comma. (This usage is clearly shown in previous comments.) To avoid any ambiguity I have edited the answer to show clearly what the comma and point mean.
    – whuber
    Commented Feb 11, 2013 at 16:20
86

As Liedman says in his answer Williams’s aviation formulas are an invaluable source, and to keep the accuracy within 10 meters for displacements up to 1 km you’ll probably need to use the more complex of these.

But if you’re willing to accept errors above 10m for points offset more than approx 200m you may use a simplified flat earth calculation. I think the errors still will be less than 50m for offsets up to 1km.

 //Position, decimal degrees
 lat = 51.0
 lon = 0.0

 //Earth’s radius, sphere
 R=6378137

 //offsets in meters
 dn = 100
 de = 100

 //Coordinate offsets in radians
 dLat = dn/R
 dLon = de/(R*Cos(Pi*lat/180))

 //OffsetPosition, decimal degrees
 latO = lat + dLat * 180/Pi
 lonO = lon + dLon * 180/Pi 

This should return:

 latO = 51,00089832
 lonO = 0,001427437
8
  • 13
    I just want to point out that this is identical to the answer I provided except you have replaced my value of 111,111 meters per degree by 111,319.5. Your value is slightly better at high latitudes but slightly worse in the lower latitudes (from 0 to about 40 degrees). Either value meets the stated requirements of accuracy.
    – whuber
    Commented Oct 27, 2010 at 21:42
  • 3
    +1 for providing code. Note that it's more accurate than you suspect (the error is typically less than 5 m over 2000 m).
    – whuber
    Commented Oct 27, 2010 at 21:55
  • 1
    I wondered if I should add a remark in my answer that this is an identical solution to yours except for the value of R, but left it out due to brevity. When it comes to precision you’re right as long as you don’t add any rotational errors to the system. Using offsets measured in a local projected coordinate system the rotational errors may grow quite large.
    – haakon_d
    Commented Oct 28, 2010 at 8:33
  • 1
    That's an excellent point: we have implicitly assumed the x-displacement is at least close to true east-west and the y-displacement is close to north-south. If not, they have to be converted into equivalent E-W and N-S displacements (not just "eastings" and "northings") before computing their lat-lon equivalents.
    – whuber
    Commented Nov 2, 2010 at 21:35
  • The d distance parameter of the Aviation Formulary equations is in radians, e.g. (distance/radius-of-earth). Commented Nov 15, 2012 at 21:10
27

I find that Aviation Formulary, here is great for these types of formulas and algorithms. For your problem, check out the "lat/long given radial and distance":here

Note that this algorithm might be a bit too complex for your use, if you want to keep use of trigonometry functions low, etc.

1
  • Thanks for this - looks ideal. Though I can't figure out if the distance is in meters or some other measurement.
    – Thomas O
    Commented Oct 27, 2010 at 11:25
4

I created a simple custom map on Google Maps that illustrates the estimation algorithm mentioned by the accepted answer (1/111111 == one meter). Feel free to see and play with it here:

https://drive.google.com/open?id=1XWlZ8BM00PIZ4qk43DieoJjcXjK4z7xe&usp=sharing

enter image description here

To give more context - this map shows the coordinates 0,0 and then shows two more pins which are 1 meter north and 1 meter east of 0,0.

3

It might make sense to project the point first. You could make something like this pseudo-code:

falt_coordinate = latlon_to_utm(original_koordinate)
new_flat_coordinate = flat_coordinate + (x,y)
result_coordinate = utm_to_latlon(new_flat_coordinate)

where (x,y) is the desired offset.

You don't need to use utm, any flat coordinate system, that makes sense in your area will do.

What software are you working with?

2
  • 1
    This would not work if you are at the edge of an UTM zone and if the move would make it cross to another UTM zone
    – Goldorak84
    Commented Jan 15, 2021 at 14:38
  • wow how can convert the last point in degree distance relative to first point ? Commented Jul 14, 2022 at 15:16
3

Here is the Python code for @whuber's answer

from math import cos, radians

def meters_to_lat_lon_displacement(m, origin_latitude):
    lat = m / 111111
    lon = m / (111111 * cos(radians(origin_latitude)))
    return lat, lon

Here is the R code:

deg2rad = function(deg) {(deg * pi) / (180)}

meters_to_lat_lon_displacement = function(m, origin_latitude){
  lat = m / 111111
  lon = m / (111111 * cos((deg2rad(origin_latitude))))
  return(list(lat=lat,lon=lon))
}
1
  • wow how can convert the last point in degree distance relative to first point ? Commented Jul 14, 2022 at 15:16
0

Here is Swift version used on kodisha answer on SO:

extension CLLocationCoordinate2D {
  
  func earthRadius() -> CLLocationDistance {
    let earthRadiusInMetersAtSeaLevel = 6378137.0
    let earthRadiusInMetersAtPole = 6356752.314
    
    let r1 = earthRadiusInMetersAtSeaLevel
    let r2 = earthRadiusInMetersAtPole
    let beta = latitude
    
    let earthRadiuseAtGivenLatitude = (
      ( pow(pow(r1, 2) * cos(beta), 2) + pow(pow(r2, 2) * sin(beta), 2) ) /
        ( pow(r1 * cos(beta), 2) + pow(r2 * sin(beta), 2) )
    )
    .squareRoot()
    
    return earthRadiuseAtGivenLatitude
  }
  
  func locationByAdding(
    distance: CLLocationDistance,
    bearing: CLLocationDegrees
  ) -> CLLocationCoordinate2D {
    let latitude = self.latitude
    let longitude = self.longitude
    
    let earthRadiusInMeters = self.earthRadius()
    let brng = bearing.degreesToRadians
    var lat = latitude.degreesToRadians
    var lon = longitude.degreesToRadians
    
    lat = asin(
      sin(lat) * cos(distance / earthRadiusInMeters) +
        cos(lat) * sin(distance / earthRadiusInMeters) * cos(brng)
    )
    lon += atan2(
      sin(brng) * sin(distance / earthRadiusInMeters) * cos(lat),
      cos(distance / earthRadiusInMeters) - sin(lat) * sin(lat)
    )
    
    let newCoordinate = CLLocationCoordinate2D(
      latitude: lat.radiansToDegrees,
      longitude: lon.radiansToDegrees
    )
    
    return newCoordinate
  }
}

extension FloatingPoint {
  var degreesToRadians: Self { self * .pi / 180 }
  var radiansToDegrees: Self { self * 180 / .pi }
}
1
  • wow how can convert the last point in degree distance relative to first point ? Commented Jul 14, 2022 at 15:16
0

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).

This is just an estimate and is not good enough as it starts changing as we go towards the poles (not just at the poles).

The more accurate formula is as below

# Meters per degree at given latitude
meters_per_degree_lat = (111132.92 - 559.82 * np.cos(2 * lat0_rad) +
                             1.175 * np.cos(4 * lat0_rad) - 0.0023 * np.cos(6 * lat0_rad))

meters_per_degree_lon = (111412.84 * np.cos(lat0_rad) -
                            93.5 * np.cos(3 * lat0_rad) + 0.118 * np.cos(5 * lat0_rad))

Reference

Length of Degree @ https://en.wikipedia.org/wiki/Geographic_coordinate_system

whubers more detailed answer @ https://gis.stackexchange.com/a/75535/61719

-3

Vincenty's direct formula should do the job.

1
  • 4
    Welcome to GIS SE! As a new user please take the tour. A good answer should be self-contained, and include all the steps and details required. Please edit your answer to add explanation as to how to use your suggestion, and why.
    – Midavalo
    Commented Jul 27, 2020 at 1:54

Not the answer you're looking for? Browse other questions tagged or ask your own question.