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?

narmaps
  • 664
  • 5
  • 16
  • possible duplicate of https://gis.stackexchange.com/questions/213034/converting-elevations-into-correct-units-with-pyproj – AndreJ Oct 28 '18 at 06:28

1 Answers1

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.

Spacedman
  • 63,755
  • 5
  • 81
  • 115
  • 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. – narmaps Oct 18 '18 at 23:10
  • Mabe related: https://gis.stackexchange.com/questions/213034/converting-elevations-into-correct-units-with-pyproj – AndreJ Oct 28 '18 at 06:27