Inaccuracy of the UTM to/from lat, lon conversion

[carlodef]

When converting coordinates from UTM to lat, lon then back to UTM with utm.to_latlon and utm.from_latlon you don't end up on the initial point. The error is about 20 to 30 centimeters, depending on the position of the initial point.

import utm
import pyproj
import numpy as np

def utm_latlon_utm_error_with_utm_package(x, y, number, letter):
    lat, lon = utm.to_latlon(x, y, number, letter)
    x1, y1 = utm.from_latlon(lat, lon, force_zone_number=number)[:2]
    return np.linalg.norm([x1 - x, y1 - y])

def utm_latlon_utm_error_with_pyproj_package(x, y, number, letter):
    proj = '{:02d}{}'.format(number, letter)
    lon, lat = pyproj.transform(pyproj.Proj('+proj=utm +zone={}'.format(proj)),
                                pyproj.Proj('+proj=latlong'), x, y)
    x1, y1 = pyproj.transform(pyproj.Proj('+proj=latlong'),
                              pyproj.Proj('+proj=utm +zone={}'.format(proj)),
                              lon, lat)
    return np.linalg.norm([x1 - x, y1 - y])

def lon_lat_to_utm(lon, lat):
    n = utm.latlon_to_zone_number(lat, lon)
    l = utm.latitude_to_zone_letter(lat)
    proj = '{:02d}{}'.format(n, l)
    x, y = pyproj.transform(pyproj.Proj('+proj=latlong'),
                            pyproj.Proj('+proj=utm +zone={}'.format(proj)),
                            lon, lat)
    return x, y, n, l

lat, lon = 2.339404, 48.857767
x, y, n, l = lon_lat_to_utm(lon, lat)
print('Distance after UTM --> (lon, lat) --> UTM conversion')
print('\twith the utm package: {:.6f} m'.format(utm_latlon_utm_error_with_utm_package(x, y, n, l)))
print('\twith the pyproj package: {:.6f} m'.format(utm_latlon_utm_error_with_pyproj_package(x, y, n, l)))

Running this Python snippet gives:

Distance after UTM --> (lon, lat) --> UTM conversion
	with the utm package: 0.280308 m
	with the pyproj package: 0.000000 m

Any idea on where this might come from?

[agressin]

I also see such inaccurate transformation when comparing utm package with the wheel known osr package https://pcjericks.github.io/py-gdalogr-cookbook/projection.html,
On my sample data, I observe shift up to 60cm in X (eastern) direction, Y (northern) seam good.

[dzanaga]

I also encountered the same issue:

>>> import utm

>>> utm.from_latlon(*utm.to_latlon(193185, 1833395, 29, 'Q'))
(193184.4488629346, 1833394.9978964194, 29, 'Q')

With geopandas instead:

>>> import geopandas as gpd

>>> gpd.GeoSeries(Point(193185, 1833395), crs={'init': 'epsg:32629'}).to_crs(epsg=4326).to_crs(epsg=32629)
0    POINT (193185.0000000001 1833395)
dtype: object

[tjwilli58]

If you can fix this issue without requiring pyproj, please do. I can install this module under msys2, but I can't install pyproj under msys2. (At least until someone with more experience with me builds it!)

[tjwilli58]

Similar results with geopy:
(python >>> removed)

import utm
import geopy
from geopy import distance
lat, lon = 2.339404, 48.857767
pt1 = geopy.Point(lat,lon)
pt1
Point(2.339404, 48.857767, 0.0)
u = utm.from_latlon(lat, lon)
u
(261764.74612451534, 258757.73393756224, 39, 'N')
pt2 = geopy.Point(utm.to_latlon(*u))
pt2
Point(2.3394040038379655, 48.8577644806325, 0.0)
distance.distance(pt1,pt2).meters
0.2802228495994366

[tjwilli58]

Apparently PROJ default datum is GRS-80, so

distance.distance(pt1,pt2,ellipsoid='WGS-84').meters
0.2802228495994366
distance.distance(pt1,pt2,ellipsoid='GRS-80').meters
0.28022284959482346
distance.distance(pt1,pt2,ellipsoid='GRS-80').meters - distance.distance(pt1,pt2,ellipsoid='WGS-84').meters
-4.613143200771219e-12

[tjwilli58]

Looking the the PROJ documentation. It mentions that the algorithms are given in [Karney2013] and the addenda points to a link of implementations in Python and other languages. Is utm.to_latlon() using this? Maybe geographiclib on PyPi could be incorporated?

[chlxt]

I fixed a bug, see pull #49

this bug may cause in-accuracy mentioned above.

[chlxt]

import utm
import pyproj
import numpy as np

def utm_latlon_utm_error_with_utm_package(x, y, number, letter):
    lat, lon = utm.to_latlon(x, y, number, letter)
    x1, y1 = utm.from_latlon(lat, lon, force_zone_number=number)[:2]
    return np.linalg.norm([x1 - x, y1 - y])

def utm_latlon_utm_error_with_pyproj_package(x, y, number, letter):
    proj = '{:02d}{}'.format(number, letter)
    lon, lat = pyproj.transform(pyproj.Proj('+proj=utm +zone={}'.format(proj)),
                                pyproj.Proj('+proj=latlong'), x, y)
    x1, y1 = pyproj.transform(pyproj.Proj('+proj=latlong'),
                              pyproj.Proj('+proj=utm +zone={}'.format(proj)),
                              lon, lat)
    return np.linalg.norm([x1 - x, y1 - y])

def lon_lat_to_utm(lon, lat):
    n = utm.latlon_to_zone_number(lat, lon)
    l = utm.latitude_to_zone_letter(lat)
    proj = '{:02d}{}'.format(n, l)
    x, y = pyproj.transform(pyproj.Proj('+proj=latlong'),
                            pyproj.Proj('+proj=utm +zone={}'.format(proj)),
                            lon, lat)
    return x, y, n, l

lat, lon = 2.339404, 48.857767
x, y, n, l = lon_lat_to_utm(lon, lat)
print('Distance after UTM --> (lon, lat) --> UTM conversion')
print('\twith the utm package: {:.6f} m'.format(utm_latlon_utm_error_with_utm_package(x, y, n, l)))
print('\twith the pyproj package: {:.6f} m'.format(utm_latlon_utm_error_with_pyproj_package(x, y, n, l)))

Now has the following result:

Distance after UTM --> (lon, lat) --> UTM conversion
        with the utm package: 0.001527 m
        with the pyproj package: 0.000000 m

I think this error(~1e-3) is reasonable when using this kind of approximate algorithm.

[chlxt]

and

utm.from_latlon(*utm.to_latlon(193185, 1833395, 29, 'Q'))

now gives

(193185.00006130006, 1833395.009825141, 29, 'Q')

[dzanaga]

@chlxt nice! is it available already in the latest version?

[chlxt]

I don't think so, currently It hasn't be approved to merged to the master branch Turbo87/utm, but you can try my fork of utm

[tjwilli58]

I think this hasn't been merged because the build failed for some older versions of Python for your pull request
I not familiar with the process, but maybe a comment on your pull request is in order?

[tjwilli58]

Looks like PRs #46 and #47 have the same issue.