2

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?

2

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.

  • I think you nailed it. I tweaked the false easting and northing to the actual values shown in QGIS which are in metres, similar to the values you used. I read somewhere that proj4 definitions should always be in SI units (metres) and then your solution produces a sufficiently close result to be considered correct. – user2723146 Oct 18 '18 at 23:10
  • Mabe related: gis.stackexchange.com/questions/213034/… – AndreJ Oct 28 '18 at 6:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.