convert latitude and longitude to x and y grid system using python

9
votes

I have file with latitude and longitude values, that i want to convert the x and y in km I want to measure the distance from each point.

for instance I make the first points of latitude and longitude(which are 51.58, -124.6 respectfully)

to (0,0) in my x and y system so than basically i want to find out what the other points are and their location from the origin so i want to find what 51.56(lat) -123.64(long) is in (x,y) in km and so on for the rest of the file.

I want to do this all in python, is there some sort code ?

I found sites online , for instance

http://www.whoi.edu/marine/ndsf/cgi-bin/NDSFutility.cgi?form=0&from=LatLon&to=XY

does exactly want i want to do, I just don't know how they do it.

5
pythonlatitude-longitudecoordinate-systems
You do realize that the earth is a sphere, so there is no general solution to your problem (in other words, the world cannot be mapped to a 2D surface without some approximations / distortions). You can find the formulas you need at movable-type.co.uk/scripts/latlong.html – Floris
If you want to try to solve the problem yourself, you should start coding, and come back when you have a specific, narrow question. If you just want a ready-made answer, this is probably not the site to be looking on. – Mark Reed
@Floris technically an oblate spheroid, but your point stands :) – dabhaid
The site linked is made in HTML and js you could translate it to python? – Noelkd
@Floris in this case i am assuming that this certain location that corresponds with my data is flat, just for this conversion of lat, long to x and y – learner

5 Answers

10
votes

The following gets you pretty close (answer in km). If you need to be better than this, you have to work harder at the math - for example by following some of the links given.

import math
dx = (lon1-lon2)*40000*math.cos((lat1+lat2)*math.pi/360)/360
dy = (lat1-lat2)*40000/360

Variable names should be pretty obvious. This gives you 

dx = 66.299 km (your link gives 66.577)
dy = 2.222 km (link gives 2.225)

Once you pick coordinates (for example, lon1, lat1) as your origin, it should be easy to see how to compute all the other XY coordinates.

Note - the factor 40,000 is the circumference of the earth in km (measured across the poles). This gets you close. If you look at the source of the link you provided (you have to dig around a bit to find the javascript which is in a separate file) you find that they use a more complex equation:

function METERS_DEGLON(x)
{  
   with (Math)
   {
      var d2r=DEG_TO_RADIANS(x);
      return((111415.13 * cos(d2r))- (94.55 * cos(3.0*d2r)) + (0.12 * cos(5.0*d2r)));
   }
}

function METERS_DEGLAT(x)
{
   with (Math)
   {
      var d2r=DEG_TO_RADIANS(x);
      return(111132.09 - (566.05 * cos(2.0*d2r))+ (1.20 * cos(4.0*d2r)) - (0.002 * cos(6.0*d2r)));
   }
}

It looks to me like they are actually taking account of the fact that the earth is not exactly a sphere... but even so when you are making the assumption you can treat a bit of the earth as a plane you are going to have some errors. I'm sure with their formulas the errors are smaller...

4
votes

UTM projections are in meters. So you could use something like the utm lib at this link:

https://pypi.python.org/pypi/utm

Googling python lat lon to UTM will point to several options.

UTM zones are 6 degrees of longitude wide and start from 0 at the prime meridian. The origin of each UTM zone is on the equator (x-axis) with the y-axis at the western most degree of longitude. This makes the grid positive to the north and east. You could calculate your distance from these results. Values are most accurate in the middle of the UTM zone.

You also should know what datum your original lat lon values are based on and use the same datum in your conversion.

2
votes

You can get the distance between GPS points using the Great Circle Distance formula. Latitude and longitude are in an geodectic coordinate system, so you can't just convert to a flat 2D grid and use euclidean distances. You can convert sufficiently close points to an approximate grid by taking an arbitrary point like your (X,Y), setting it to the origin (like you've done) and then using great circle distance together with bearing to plot the points relative to each other on the plane, but it's an approximation.

2
votes

if you were to use a 3D system, these functions will do:

def arc_to_deg(arc):
    """convert spherical arc length [m] to great circle distance [deg]"""
    return float(arc)/6371/1000 * 180/math.pi

def deg_to_arc(deg):
    """convert great circle distance [deg] to spherical arc length [m]"""
    return float(deg)*6371*1000 * math.pi/180

def latlon_to_xyz(lat,lon):
    """Convert angluar to cartesian coordiantes

    latitude is the 90deg - zenith angle in range [-90;90]
    lonitude is the azimuthal angle in range [-180;180] 
    """
    r = 6371 # https://en.wikipedia.org/wiki/Earth_radius
    theta = math.pi/2 - math.radians(lat) 
    phi = math.radians(lon)
    x = r * math.sin(theta) * math.cos(phi) # bronstein (3.381a)
    y = r * math.sin(theta) * math.sin(phi)
    z = r * math.cos(theta)
    return [x,y,z]

def xyz_to_latlon (x,y,z):
    """Convert cartesian to angular lat/lon coordiantes"""
    r = math.sqrt(x**2 + y**2 + z**2)
    theta = math.asin(z/r) # https://stackoverflow.com/a/1185413/4933053
    phi = math.atan2(y,x)
    lat = math.degrees(theta)
    lon = math.degrees(phi)
    return [lat,lon]
0
votes

You can use UTM:

pip install utm

Here is an example:

>>> import utm
>>> utm.from_latlon(51.2, 7.5)
(395201.3103811303, 5673135.241182375, 32, 'U')

The return has the form (EASTING, NORTHING, ZONE_NUMBER, ZONE_LETTER).

Notes

It works with NumPy arrays too:

>>> utm.from_latlon(np.array([51.2, 49.0]), np.array([7.5, 8.4]))
(array([395201.31038113, 456114.59586214]),
 array([5673135.24118237, 5427629.20426126]),
 32,
 'U')

And in reverse:

>>> utm.to_latlon(340000, 5710000, 32, 'U')
(51.51852098408468, 6.693872395145327)