proj4 not respecting +units?

[JamesFysh]

Hi,

I was somewhat surprised to find that I was getting different coordinates back when transforming from WGS84 to NAD83 StatePlane Missouri Central (i.e. EPSG:102697), than what I was seeing in QGIS and my customers GIS.

The long and the short of it is that the EPSG:102697 uses +units=us-ft, however proj4 appears to ignore this and all projected coordinates appear to be given in meter-units. I'm not sure if this is expected behavior, and why it differs to that exhibited by GDAL?

Following is an IPython notebook than demonstrates the issue:

import pyproj
from osgeo import osr

# Given a particular coordinate in Springfield, Missouri
#
# -93.28, 37.2
LONGITUDE = -93.28
LATITUDE = 37.2

# and taking a "sensible" projected coordinate system for that region: EPSG:102697
PROJ4_INIT = "+proj=tmerc +lat_0=35.83333333333334 +lon_0=-92.5 +k=0.999933333 +x_0=500000 +y_0=0 +datum=NAD83 +units=us-ft +no_defs"

# .. Note the use of +units=us-ft in the above

def meters_to_feet(value):
    return value / 0.30480060960122

# In Proj4
p1 = pyproj.Proj("+init=epsg:4326")
p2 = pyproj.Proj(PROJ4_INIT)

print("PyPROJ: {}, {} -> {}".format(LONGITUDE, LATITUDE, pyproj.transform(p1, p2, LONGITUDE, LATITUDE)))
x, y = pyproj.transform(p1, p2, LONGITUDE, LATITUDE)
print("Converting to us-ft: {}, {}".format(meters_to_feet(x), meters_to_feet(y)))
print("Backwards: {}, {} -> {}".format(x, y, pyproj.transform(p2, p1, x, y)))

# Observe: The result appears to come back in meters (not us-ft), and a transform from 
# projected coordinates -> geographic expects units in meters also.

PyPROJ: -93.28, 37.2 -> (430757.2001325557, 151931.97713024777)
Converting to us-ft: 1413242.5807682227, 498463.49496815325
Backwards: 430757.2001325557, 151931.97713024777 -> (-93.28000000000004, 37.19999999999998)

# In GDAL
p1 = osr.SpatialReference()
p1.ImportFromEPSG(4326)
p2 = osr.SpatialReference()
p2.ImportFromProj4(PROJ4_INIT)
project = osr.CoordinateTransformation(p1, p2)
reverse = osr.CoordinateTransformation(p2, p1)

print("GDAL: {}, {} -> {}".format(LONGITUDE, LATITUDE, project.TransformPoint(LONGITUDE, LATITUDE)))
x, y, _ = project.TransformPoint(LONGITUDE, LATITUDE)
print("Already appears in us-ft: {}, {}".format(x, y))
print("Backwards: {}, {} -> {}".format(x, y, reverse.TransformPoint(x, y)))

# Observe: The result appears to come back in us-ft, as expected, and a transform from
# projected coordinates -> geographic expects units in us-ft also.

GDAL: -93.28, 37.2 -> (1413242.5807682204, 498463.49496815726, 0.0)
Already appears in us-ft: 1413242.5807682204, 498463.49496815726
Backwards: 1413242.5807682204, 498463.49496815726 -> (-93.28000000000004, 37.19999999999998, 0.0)

# If we then revise the PROJ4_INIT string and re-run the GDAL test, we see the same (meter) units as per Proj.4
PROJ4_INIT = "+proj=tmerc +lat_0=35.83333333333334 +lon_0=-92.5 +k=0.999933333 +x_0=500000 +y_0=0 +datum=NAD83 +no_defs"

p1 = osr.SpatialReference()
p1.ImportFromEPSG(4326)
p2 = osr.SpatialReference()
p2.ImportFromProj4(PROJ4_INIT)
project = osr.CoordinateTransformation(p1, p2)
reverse = osr.CoordinateTransformation(p2, p1)

print("GDAL: {}, {} -> {}".format(LONGITUDE, LATITUDE, project.TransformPoint(LONGITUDE, LATITUDE)))
x, y, _ = project.TransformPoint(LONGITUDE, LATITUDE)
print("Converting to us-ft: {}, {}".format(meters_to_feet(x), meters_to_feet(y)))
print("Backwards: {}, {} -> {}".format(x, y, reverse.TransformPoint(x, y)))

GDAL: -93.28, 37.2 -> (430757.2001325534, 151931.9771302485, 0.0)
Converting to us-ft: 1413242.5807682152, 498463.49496815563
Backwards: 430757.2001325534, 151931.9771302485 -> (-93.28000000000004, 37.19999999999998, 0.0)

[rouault]

This looks like an issue with pyproj at first sight, since GDAL uses the proj C API and doesn't do itself the units to foot conversion but delegates that to proj. You'd bette reporting to https://github.com/jswhit/pyproj . Closing under the assumption that's a pyproj issue. Please reopen if you don't think it is (and have more elements to "prove" it)

[JamesFysh]

You are indeed correct, and I should have done the following simple test first...

jfysh@ampere:~$ cs2cs +init=epsg:4326 +to +proj=tmerc +lat_0=35.83333333333334 +lon_0=-92.5 +k=0.999933333 +x_0=500000 +y_0=0 +datum=NAD83 +units=us-ft +no_defs
-93.28 37.2
1413242.58      498463.49 0.00

Will raise against pyproj

Thanks!

EDIT: Already a known-issue for pyproj: pyproj4/pyproj#67