18

I've tried the following, input: lat/lon data then I'll calculate a box around it by, let's say 50 m, so +/- 50 m on easting/northing value.

Now I reconvert it to lat/lon and with a script:

http://robotics.ai.uiuc.edu/~hyoon24/LatLongUTMconversion.py I get a result that just can't be, lon before is around 7, afterwards around 2. 

zone, easting, northing = LLtoUTM(23, location.get_lat(), location.get_lon()) 

topUTM = northing + error
bottomUTM = northing - error
leftUTM = easting - error
rightUTM = easting + error
left, top = UTMtoLL(23, leftUTM, topUTM, zone)

Is the error in my code, or might be the script flawed?

So I've tried to use pyproj, just lat/lon to utm to lat/lon to see what happens

>>> p = pyproj.Proj(proj='utm', zone=32, ellps='WGS84')
>>> p
<pyproj.Proj object at 0x7ff9b8487dd0>
>>> x,y = p(47.9941214, 7.8509671)
>>> print x,y
5159550.36822 1114087.43925
>>> print p(x,y,inverse=True)
(47.971558538495991, 7.8546573140162605)

And here it's not as extremely far off as with the script from above, but it still seems strongly enough incorrect as to not be able to use it. How come? What can I do to get more exact results?

EDIT:

I ran test() and it all tests passed.

in epsg file there is no such thing. The closest I've found was this:

<32632> +proj=utm +zone=32 +ellps=WGS84 +datum=WGS84 +units=m +no_defs <>

no tmerc. Also What would I need to pass the towgs84 as parameters? The ones above?

Kara
  • 6,115
  • 16
  • 50
  • 57
luh
  • 181
  • 1
  • 1
  • 3

4 Answers4

62

I've created a small UTM conversion library for Python last week and uploaded it to the Python Package Index: http://pypi.python.org/pypi/utm

I have compared it to using pyproj and it is faster and more accurate. Given your sample data, this is the result:

>>> import utm

>>> u = utm.from_latlon(47.9941214, 7.8509671)
>>> print u
(414278, 5316285, 32, 'T')

>>> print utm.to_latlon(*u)
(47.994157948891505, 7.850963967574302)

UPDATE: Richards answer below describes the real solution for this issue.

Community
  • 1
  • 1
TBieniek
  • 4,858
  • 2
  • 24
  • 29
  • Thanks for the shout-out, @TBieniek! – Richard Sep 02 '15 at 19:47
  • Thanks a lot for this implementation, absolutely valuable! –  Nov 16 '17 at 10:19
  • 1
    Is utm.to_latlon broadcastable on a numpy array? – Ian Campbell Moore Mar 28 '18 at 19:41
  • 5
    Thanks for creating a whole freakin Python package and uploading it on pip to answer a StackOverflow question. Truly above and beyond – Jack Jul 31 '19 at 20:17
  • Thanks! So much easier to use (and install) than pyproj! – Simmovation Jul 28 '20 at 00:51
  • Hi, I've tried your pypi package, it does seem much easier to use. I have a conversion error however, utm.to_latlon(246000, 9855098, 37, 'S') returns (88.35340991983618, -159.09274162642123), am I doing something wrong in the input? – Olli Oct 07 '21 at 16:17
  • Careful that the conversion has flaws. I've run over a dataset in zone 31 V and there's 100m+ error... – Raf Feb 20 '23 at 08:02
34

The error is in your code.

First off, the PyProj issue listed in one of the other answers is real. You should check your epsg file and make sure it includes the line 

<2392> +proj=tmerc +lat_0=0 +lon_0=24 +k=1.000000 +x_0=2500000 +y_0=0 +ellps=intl +towgs84=-90.7,-106.1,-119.2,4.09,0.218,-1.05,1.37 +units=m +no_defs no_defs <>

Note the towgs84 parameter.

Your problem with PyProj stems from mis-using the projection command.

If we take 47.9941214N, 7.8509671E and convert to UTM we get Zone 32, 414278 Easting, 5316286 Northing.

You perform the following PyProj operations:

p = pyproj.Proj(proj='utm', zone=32, ellps='WGS84')
>>> x,y = p(47.9941214, 7.8509671)
>>> print x,y
5159550.36822 1114087.43925
>>> print p(x,y,inverse=True)
(47.971558538495991, 7.8546573140162605)

But, if we consult the PyProj documentation, we see the following:

Calling a Proj class instance with the arguments lon, lat will convert lon/lat (in degrees) to x/y native map projection coordinates (in meters).

Let's try running the OP's PyProj operations again, but switch the order of the lon/lat arguments:

p = pyproj.Proj(proj='utm', zone=32, ellps='WGS84')
>>> x,y = p(7.8509671, 47.9941214)
>>> print x,y
414278.16731 5316285.59492
>>> print p(x,y,inverse=True)
(7.850967099999812, 47.994121399999784)

The operation inverts itself (pretty much) perfectly!

To answer the first part of your question, if you look in http://robotics.ai.uiuc.edu/~hyoon24/LatLongUTMconversion.py at the definition of UTMtoLL, you find the following:

UTMtoLL(ReferenceEllipsoid, northing, easting, zone)

Yet you use UTMtoLL(23, leftUTM, topUTM, zone) where leftUTM is an Easting and topUTM is a Northing.

Therefore, in the case of both your first script and PyProj, you've used the wrong order of arguments.

It's a good reminder to always double- (or triple-) check your work before suggesting that someone else's is wrong. That said, Python's documentation is not the greatest and PyProj's documentation in this instance is cryptic at best. A nice web-based explanation of this command and accompanied by examples of its usage would have probably prevented angst on your part.

Richard
  • 56,349
  • 34
  • 180
  • 251
  • 1
    Is there any way to automatically detect UTM zone from latitude longitude. – BetterCallMe Jul 15 '21 at 06:48
  • Is there better doc on this since 2013 ? (GDAL seems to fall over on roundtrip lat-lon -> utm -> lat-lon, see this [plot](https://github.com/corteva/rioxarray/issues/396) ) – denis Sep 21 '21 at 18:29
  • @Lostman: `math.floor(((lon + 180) % 360) / 6) + 1` should work. See https://stackoverflow.com/questions/9186496/determining-utm-zone-to-convert-from-longitude-latitude – Eric Duminil Apr 29 '22 at 17:10
4

I have no problem with pyproj, try the following code

from pyproj import Proj

Lat = 52.063098675
Lon = -114.132980348 #Calgary

ZoneNo = "11" #Manually input, or calcuated from Lat Lon
myProj = Proj("+proj=utm +zone="+ZoneNo+",\
+north +ellps=WGS84 +datum=WGS84 +units=m +no_defs") #north for north hemisphere
UTMx, UTMy = myProj(Lon, Lat)

########################################

#UTM ==> Lat Lon:
ZoneNo = "11" #Manually input or from other sources
myProj = Proj("+proj=utm +zone="+\
ZoneNo+", +north +ellps=WGS84 +datum=WGS84 +units=m +no_defs")
Lon2, Lat2 = myProj(UTMx, UTMy,inverse=True)

print Lat2
print Lon2
Ting On Chan
  • 121
  • 2
2

Your issue with pyProj sounds just like the one described here:

http://code.google.com/p/pyproj/issues/detail?id=3

which is resolved:

solved! in epsg file there must be

<2392> +proj=tmerc +lat_0=0 +lon_0=24 +k=1.000000 +x_0=2500000 +y_0=0 +ellps=intl +towgs84=-90.7,-106.1,-119.2,4.09,0.218,-1.05,1.37 +units=m +no_defs no_defs <>

note the towgs84 parameter!

Check that thread out if you want to continue to use pyproj.

Also, does the test() function of the module work? Have you tried any of the scripts that come with it in the test directory?

agf
  • 171,228
  • 44
  • 289
  • 238
  • Those hoping that Ubuntu 13.04 Raring would include an up-to-date PyProj would find themselves mistaken. Install from source! – Richard Sep 04 '13 at 17:46