Python – Transform Coordinates from Custom Projection to Lat/Lon Using PyProj

affine transformationcoordinate systempyprojpython

I am trying to convert a pair of coordinates specified in US survey feet on a custom transverse mercator projection to lat/long using PyProj. I must be doing something wrong because I cannot get the correct decimal degree values:

import pyproj

# Coordinates in ft and false easting/northing in ft
x1,y1 = [1383566.534, 1627687.047]
sourceProjection =  "+proj=tmerc +lat_0=42 +lon_0=-115.5 +k=0.9996 +x_0=1640416 +y_0=328083 +datum=NAD83 +units=us-ft +no_defs"

x2,y2 = pyproj.transform(pyproj.Proj(sourceProjection), pyproj.Proj(init='epsg:4326'), x1, y1)

print "Original:    {0}, {1}".format(x1, y1)
print "Transformed: {0}, {1}".format(x2, y2)

The output I get is:

Original:    1383566.534, 1627687.047
Transformed: -119.384930994, 53.6300677515

but the correct result is -116.503170 45.562201

What am I doing wrong please?

Best Answer

If I convert everything to metres I get a lot closer - but maybe I've used the us-foot instead of the imperial foot or mixed both. Anyway, I did this:

Convert your point to metres:

>>> x1m = x1 * 1200.0/3937.0 
>>> y1m = y1 * 1200.0/3937.0

Convert the projection offsets to metres:

>>> 328083*(1200.0/3937.0)
99999.89839979679
>>> 1649416*(1200.0/3937.0)
502743.00228600454

Create projection string with those offsets - remove "units":

>>> pm =  "+proj=tmerc +lat_0=42 +lon_0=-115.5 +k=0.9996 +x_0=502743 +y_0=100000 +datum=NAD83  +no_defs"

Transform:

>>> pyproj.transform(pyproj.Proj(pm), pyproj.Proj("+init=epsg:4326"), x1m,y1m)
(-116.53831300493289, 45.56188657919288)

SO close! I did also just try using imperial feet (conversion factor 0.3048) and its identical to about four decimal places, so its not that).

I don't understand why pyproj seems to ignore the units (try your original projection but with +units=m, the output is the same), but I'm not sure even with +units=us-ft that you can use feet in the x_0 type of parameters.

Related Question