I am trying to convert NAD83(2011) / Louisiana South (ftUS) coordinates on survey maps to lat/lon using Proj4 in Python.
Could someone point me in the right direction?
I have verified the lat/lon on earthpoint.us site. Their conversion seems to work. This is my python code which contains the Proj4 string.
===============
from pyproj import Proj
#EPSG:6479 from epsg-registry.org
southla = '+proj=lcc +lat_1=30.7 +lat_2=29.3 +lat_0=28.5 +lon_0=-91.33333333333333 +x_0=3280833.3333 +y_0=0 +ellps=GRS80 +datum=NAD83 +towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs'
p=Proj(southla)
lonlat = p(3037416.890, 709211.549,inverse=True)
print "%s, %s" % (lonlat[1],lonlat[0])
Answer
As as there is no support EPSG:6479 (NAD83(2011) / Louisiana South (ftUS)) in PROJ4 (look at the comment of mkennedy) I will illustrate the problem with EPSG:3452, NAD83 / Louisiana South (ftUS) because the problem is the same, the units of the projection.
Solution with Pyproj
import pyproj
southla = pyproj.Proj('+proj=lcc +lat_1=30.7 +lat_2=29.3 +lat_0=28.5 +lon_0=-91.33333333333333 +x_0=3280833.3333 +y_0=0 +ellps=GRS80 +datum=NAD83 +towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs')
# or southla = pyproj.Proj(init='epsg:3452')
wgs84 = pyproj.Proj('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')
longlat = pyproj.transform(southla,wgs84,3037416.890, 709211.549)
print longlat
(-93.986162555568697, 34.865093872789771)
Result
Solution with GDAL
from osgeo import osr
southla = osr.SpatialReference()
southla.ImportFromEPSG(3452)
wgs84 = osr.SpatialReference()
wgs84.ImportFromEPSG(4326)
x = 3037416.890
y = 709211.549
transformation = osr.CoordinateTransformation(southla,wgs84)
result = transformation.TransformPoint(x, y)
print result
(-92.105819299234454, 30.447920720765676, 0.0)
Result
Why the difference?
Because PyProj assumes that your coordinates are in meters and the unit of EPSG:3452 (or EPSG:6479) is "US survey foot" (units=us-ft) -> look at Converting X,Y coordinates to lat/long using pyproj and Proj.4 returns the wrong coordinates
preserve_units=True
fix the problem (pyproj silently changes '+units=' parameter)
foot_proj = pyproj.Proj(init="epsg:3452", preserve_units=True)
longlat = pyproj.transform(foot_proj,wgs84,3037416.890, 709211.549)>>> longlat
print longlat
(-92.105819299234469, 30.447920720765683)
Result
No comments:
Post a Comment