incorrect conversion from ESRI WKT to Proj.4 for NAD83 UTM

[zhiyuli]

First I would like to say this lib helps me a lot. I am using it in conjunction with PyProj lib in several projects already.

Issue:
I have a ESRI WKT string, and in this particular case I know its corresponding EPSG code is 26916 because the ESRI WKT string represents NAD83 UTM Zone 16.

I was trying to re-project a point from above projection into EPSG:3857 using PyProj for a web mapping app, but the result looked incorrect. The following code shows the different results when PyProj objs are init'd by proj.4 string produced from PyCRS and EPSG code.

import pycrs
import pyproj

esri_prj_str_nad83_utm_zone16 = 'PROJCS["NAD_1983_UTM_Zone_16N",GEOGCS["GCS_North_American_1983",DATUM["D_NORTH_AMERICAN_1983",SPHEROID["GRS_1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-87],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["Meter",1]]'
fromcrs = pycrs.parser.from_esri_wkt(esri_prj_str_nad83_utm_zone16)
proj4_str = fromcrs.to_proj4()
prj_obj_from_proj4 = pyproj.Proj(proj4_str)

prj_obj_26916_from_epsg = pyproj.Proj("+init=EPSG:26916")
prj_obj_3857_from_epsg = pyproj.Proj("+init=EPSG:3857")

(-9748715.328601453, 3929480.5837150216)

print pyproj.transform(prj_obj_from_proj4, prj_obj_3857_from_epsg, 446495, 3681410)

(-9748750.584671738, 3931248.690782473)

print pyproj.transform(prj_obj_26916_from_epsg, prj_obj_3857_from_epsg, 446495, 3681410)

Thanks

[micahcochran]

I'm impressed that it is that close. This comes down to the definition of the UTM projection has changed in PROJ.4.

So PROJ.4 (which pyproj is built atop) takes this line of code:
prj_obj_26916_from_epsg = pyproj.Proj("+init=EPSG:26916")
It has a lookup in epsg and creates finds this proj4 string:

# NAD83 / UTM zone 16N
<26916> +proj=utm +zone=16 +datum=NAD83 +units=m +no_defs  <>

The proj4_str that you show above from pycrs results in the proj4_string:
+proj=tmerc +datum=NAD83 +ellps=GRS80 +a=6378137.0 +f=298.257222101 +pm=0.0 +lat_0=0.0 +lon_0=-87.0 +x_0=500000.0 +y_0=0.0 +units=m +axis=enu +no_defs

I replaced +proj=tmerc with etmerc, I get the first result that show. This is because PROJ.4 version 4.9.3 changed UTM projection algorithms from tmerc to etmerc (OSGeo/PROJ#316). At some point, PROJ.4 may change the definition of tmerc to be the same as etmerc, but I wouldn't hold my breath.

[karimbahgat]

Glad to hear you find it useful @zhiyuli. I'm just going to have to concur with @micahcochran, and props for keeping up with the latest on proj4.

So if I understand correctly, PyCRS uses the old way of referring to utm (tmerc)? How recent was this change and do you recommend that we change to the newer etmerc or stick with the old?

[micahcochran]

@karimbahgat I found out about this when I was creating the GIGS testing framework for Proj4.

That is correct. etmerc and tmerc are algorithms for Transverse Mercator projections.

Proj4 version 4.9.3 was released 2016-08-22, which made UTM use the etmerc version. The issue was raised in 2015-10-06. Their change was based on a report by the US National Geospatial Inteligence Agency on 2014-03-25. etmerc algorithm was added to Proj4 release 4.8.0 (source), which was released 2012-03-07.

The WKT itself doesn't not specify the method/algorithm used to reproject data. The WKT only specifies that the coordinates stored in Transverse Mercator projection. tmerc and etmerc both give valid results.

I don't know if there is enough information in the WKT to say that this is a +proj=utm. Another problem is that other spatial reference systems use the tmerc algorithm. I'm not sure if it is fair to those countries geospatial authorities to do an algorithm shift. Even the utm algorithm switch has already changed that algorithm around the world.

I haven't really thought through the right coarse of action.

[karimbahgat]

@micahcochran, thanks for the details. So since this is essentially a very recent legacy/incompatability issue, there is no clear solution of switching to the new etmerc method or sticking with the old tmerc, as this might lead to unexpected changes for projection definitions written before or not agreeing with the recent USGS recommendation.

Perhaps the best solution in this scenario is to let the user set a function arg or library constant deciding whether to use the old or new algo, defaulting to the old..

To be clear though, algorithms aside, setting the proj4 +init=utm is simply a shorthand to allow you to specify more easily w the +zone arg, but is functionally equivalent to using +etmerc with all the additional args (as was produced by pyproj)? Is this all that differentiates the utm projection from ordinary trans merc? This makes me wonder too, i cant remember taking into account utm zones when i wrote pycrs...

[micahcochran]

@karimbahgat Everything you have stated is correct.

What you propose seems like a good solution.

[karimbahgat]

Assigning this to myself to add the extra arg to choose between old or new.