[GIS] Converting lat,lon to x,y coordinates with pyproj

convertcoordinate systempyproj

I am trying to convert lat,lon to x,y coordinates in France, using the pyproj library and what I am doing is obviously wrong.

Here is some of my code:

import pyproj
import math

PROJ='+proj=utm +zone=31, +north +ellps=WGS84 +datum=WGS84 +units=m +no_defs'

def LatLon_To_XY(Lat,Lon):
    p1=pyproj.Proj(PROJ,preserve_units=True)
    (x,y)=p1(Lat,Lon)
    return(x,y)


def XY_To_LatLon(x,y):
    p1=pyproj.Proj(PROJ,preserve_units=True)
    (lat,lon)=p1(x,y,inverse=True)
    return(lat,lon)

def distance(x1,y1,x2,y2):
    d1=x1-x2
    d2=y1-y2
    out=math.sqrt(d1*d1+d2*d2)
    return(out)

Problem is when I am validating this computing the distance between Paris and Toulouse, I get 816km where the distance is supposed to be around 580km.

I guess I have wrong parameters for the projection definition.

Here are the values displayed:

  • For Paris lat,lon (48.856614, 2.3522219) and x,y: 6253348.594037977 374425.11288885673

  • For Toulouse lat,lon (43.604652, 1.444209) and x,y: 5453729.038474649 210744.18988320036

Best Answer

Nope, there's nothing wrong with your calculations; you are looking at projection distortions!

I´d suggest to use the built in Geod class in pyproj for geodetic computations. e.g.

import pyproj
import math

P = pyproj.Proj(proj='utm', zone=31, ellps='WGS84', preserve_units=True)
G = pyproj.Geod(ellps='WGS84')

def LatLon_To_XY(Lat,Lon):
    return=P(Lat,Lon)    

def XY_To_LatLon(x,y):
    return P(x,y,inverse=True)    

def distance(Lat1, Lon1, Lat2, Lon2):
    return G.inv(Lon1, Lat1, Lon2, Lat2)[2]

or a framework that implements geodetic algebra, e.g. vincenty, GeoPy...

Related Question