12

I need a function that maps gps positions to x/y values like this:

getXYpos(GeoPoint relativeNullPoint, GeoPoint p){
   deltaLatitude=p.latitude-relativeNullPoint.latitude;
   deltaLongitude=p.longitude-relativeNullPoint.longitude;
   ...
   resultX=latitude (or west to east) distance in meters from p to relativeNullPoint
   resultY=longitude (or south to north) distance in meters from p to relativeNullPoint
}

i have seen some implementations of "distance of two geoPoints" but they all just calculate the air-line distance. i think the deltaLongitude can be transformed into meters directly but the deltaLatitude depends in the Longitude. does anyone know how this problem can be solved?

Deestan
  • 16,738
  • 4
  • 32
  • 48
Simon
  • 13,173
  • 14
  • 66
  • 90

3 Answers3

22

To start with, I think you have your latitude and longitude reversed. Longitude measures X, and latitude measures Y.

The latitude is easy to turn into a north-south distance. We know that 360 degrees is a full circle around the earth through the poles, and that distance is 40008000 meters. As long as you don't need to account for the errors due to the earth being not perfectly spherical, the formula is deltaLatitude * 40008000 / 360.

The tricky part is converting longitude to X, as you suspected. Since it depends on the latitude you need to decide which latitude you're going to use - you could choose the latitude of your origin, the latitude of your destination, or some arbitrary point in between. The circumference at the equator (latitude 0) is 40075160 meters. The circumference of a circle at a given latitude will be proportional to the cosine, so the formula will be deltaLongitude * 40075160 * cos(latitude) / 360.

Edit: Your comment indicates you had some trouble with the longitude formula; you might have used degrees instead of radians in the call to cos, that's a common rookie mistake. To make sure there's no ambiguity, here's working code in Python.

def asRadians(degrees):
    return degrees * pi / 180

def getXYpos(relativeNullPoint, p):
    """ Calculates X and Y distances in meters.
    """
    deltaLatitude = p.latitude - relativeNullPoint.latitude
    deltaLongitude = p.longitude - relativeNullPoint.longitude
    latitudeCircumference = 40075160 * cos(asRadians(relativeNullPoint.latitude))
    resultX = deltaLongitude * latitudeCircumference / 360
    resultY = deltaLatitude * 40008000 / 360
    return resultX, resultY

I chose to use the relativeNullPoint latitude for the X calculation. This has the benefit that if you convert multiple points with the same longitude, they'll have the same X; north-south lines will be vertical.

Edit again: I should have pointed out that this is a very simple formula and you should know its limitations. Obviously the earth is not flat, so any attempt to map it to XY coordinates will involve some compromises. The formula I derived above works best when the area you're converting is small enough to consider flat, and where the slight curvature and non-parallelism of north-south lines can be ignored. There's a whole science to map projections; if you want to see some possibilities a good place to start would be Wikipedia. This specific projection is known as the Equirectangular projection, with some added scaling.

Mark Ransom
  • 299,747
  • 42
  • 398
  • 622
  • ok i tried this, and the the latitude calculation is very accurate :) but the longitude-meter values are always the factor 1.3 to big for the latitude im testing. i think its not a good idea to just divide by this factor because it will probably differ for different latitudes, maybe the earth circumfence has to be calculated relative to the current latitude? – Simon Jun 12 '10 at 10:31
  • I just checked some of my sources (like the link you send me http://geography.about.com/library/faq/blqzcircumference.htm ) and tried to calculate the earth cirumfence myself: c=2*PI*r=2*PI*6.378.137m=40.075.017m so i think the 40075160m value used everywhere is wrong.. i got r from http://en.wikipedia.org/wiki/World_Geodetic_System#A_new_World_Geodetic_System:_WGS_84 so it is probably correct. same problem for the circumfence through the poles: there i get 39.940.653 which is an extreme difference to the 40.008.000 value.. – Simon Sep 07 '10 at 23:58
  • ok the earth isn't a perfect sphere;) i think this will explain the deviant values, especially the one for the pole-circumference . but how is the circumference of the geoid calculated? – Simon Sep 08 '10 at 00:07
  • @Sponge, I don't know how the circumferences were calculated - I just trust my sources. You're free to use different values if that works for your application. You never specified the accuracy needs of your application - the difference in the two circumference through the poles figures is less than 0.2%. – Mark Ransom Sep 08 '10 at 00:59
  • What do you mean the earth is not flat ? – kubante Apr 27 '21 at 14:20
2

Harvesine function is what you need. Check it out at moveable-types There are Distance, Bearing, Midpoint and other stuff Javascript implementation working really good.

UPDATE

I've found a Java implementation of Harvesine function in another stackoverflow question

Community
  • 1
  • 1
user347594
  • 1,256
  • 7
  • 11
1

There are libraries on jstott.me.uk for PHP, Java and Javascript which do this, e.g.

var lld1 = new LatLng(40.718119, -73.995667); // New York
document.write("New York Lat/Long: " + lld1.toString() + "<br />");
var lld2 = new LatLng(51.499981, -0.125313);  // London
document.write("London Lat/Long: " + lld2.toString() + "<br />");
var d = lld1.distance(lld2);
document.write("Surface Distance between New York and London: " + d + "km");
skaffman
  • 398,947
  • 96
  • 818
  • 769