5

I am converting a point from wgs84 into a state plane coord system (epsg: 32051, which is defined in US survey feet) using pyproj, but the transform isn't changing the z values into feet:

point_4326 = (-81.23127933,37.70957554,668.105)
proj1 = pyproj.Proj(init='epsg:4326')
proj2 = pyproj.Proj(init='epsg:32051')
print pyproj.transform(proj1, proj2, *point_44326)

$ (589192.2346128735, 78765.53566868477, 668.105)

I have also tried specifying the exact proj4 string for the new coord system (which includes the USft scaling factor), but get the exact same result:

proj2 = pyproj.Proj('+proj=lcc +lat_1=37.48333333333333 +lat_2=38.88333333333333 +lat_0=37 +lon_0=-81 +x_0=609601.2192024384 +y_0=0 +ellps=clrk66 +datum=NAD27 +to_meter=0.3048006096012192 +no_defs')
print pyproj.transform(proj1, proj2, *point_4326) 

$ (589192.2346128735, 78765.53566868477, 668.105)

Any ideas how I can get pyproj to convert the elevation values to the correct units, or if there are similar tools / libraries which can do this?

PolyGeo
  • 65,136
  • 29
  • 109
  • 338
jeremyeastwood
  • 614
  • 5
  • 15

1 Answers1

5

Look at Proj4 String for NAD83(2011) / Louisiana South (ftUS)

preserve_units=True fix the problem (pyproj silently changes '+units=' parameter)

import pyproj
proj2 = pyproj.Proj(init='epsg:32051',preserve_units=True)
proj1 = pyproj.Proj(init='epsg:4326')
print pyproj.transform(proj1, proj2, *point_4326)
(1933041.523059067, 258416.59493967678, 2191.941154166668)

Control with GDAL/OSR

from osgeo import osr
wgs84 = osr.SpatialReference()
wgs84.ImportFromEPSG(4326)
proj2= osr.SpatialReference()
proj2.ImportFromEPSG(32051)
transformation = osr.CoordinateTransformation(wgs84,proj2)
result = transformation.TransformPoint(*point_4326)
print result
(1933041.5230590631, 258416.59493967678, 2191.941154166668)

And 2192 feets = 668.122 meters

Community
  • 1
gene
  • 54,868
  • 3
  • 110
  • 187