Converting X,Y coordinates to lat/long using pyproj and Proj.4 returns the wrong coordinates

Question

I'm writing a python script that reads multiple XML files containing x and y coordinates and combines them all into a single csv file. Latitude and Longitude are required fields in the csv, but I am having difficulty converting the x,y coordinates in Ohio North State Plane usFt to WGS84.

>>> p = Proj(r'+proj=lcc +lat_1=41.7 +lat_2=40.43333333333333 +lat_0=39.66666666666666 +lon_0=-82.5 +x_0=600000 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=us-ft +no_defs') #Nad83 State Plane Ohio North US Feet Proj object using parameters
>>> p(739400.91,2339327.3,inverse=True)
(-80.138057868777224, 60.278230707978487)
>>> p1 = Proj(init="epsg:3734") #Nad83 State Plane Ohio North US Feet Proj object using EPSG code
>>> p1(739400.91,2339327.3,inverse=True)
(-80.138057868777224, 60.278230707978487)

Both methods above return the same result, however this lat long is somewhere in Hudson Bay. When I plot the coordinates in ArcMap, the correct lat long is: -81.142311,41.688205.

*Notice all lat longs are provided long,lat as this is the order Proj uses

Does anyone know why I would be getting the wrong coordinates from Proj.4 and pyproj?

Answer 1

PyProj assumes that your coordinates are in meters. I'd guess something relating to feet/meters is the cause of the issue.

Calling a Proj class instance with the arguments lon, lat will convert lon/lat (in degrees) to x/y native map projection coordinates (in meters)

If the optional keyword 'preserve_units' is True, the units in map projection coordinates are not forced to be meters.

http://pyproj.googlecode.com/svn/trunk/docs/pyproj.Proj-class.html

Are your initial coordinates in feet? When you load the data into ArcMap what units does the map use?

This gets the coordinates a little closer:

p1 = Proj(init="epsg:3734")
#1 foot = 0.3048 meters
conv = 0.3048
print p1(739400.91 * conv,2339327.3 * conv,inverse=True)
(-87.3195533069909, 45.98605408134072)

A similar issue can be found here.

Answer 2

I get the same results as @geographika when I run gdaltransform and the proj.4 tool cs2cs:

$ gdaltransform -s_srs EPSG:3734 -t_srs EPSG:4326
739400.9 2339327.3             
-87.3195485720169 45.9860670658218 0

cs2cs +proj=lcc +lat_1=41.7 +lat_2=40.43333333333333 +lat_0=39.66666666666666 +lon_0=-82.5 +x_0=600000 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=us-ft +no_defs +to +proj=lonlat +datum=WGS84
739400.9 2339327.3
87d19'10.375"W 45d59'9.841"N 0.000

Reversing the x and y coordinates of your point however gives the result that you're seeing in ArcMap:

gdaltransform -s_srs EPSG:3734 -t_srs EPSG:4326
2339327.3 739400.9
-81.1423086719059 41.6882035384526 0

So you'll need to do a visual check to ensure you do have your x and y coordinates the right way round. It's a problem I've had in the past when you get two results that are similar enough you put it down to rounding error or something.

Answer 3

I was actually trying to do the same thing except with the OH south state plane grid and I came across your question. I was getting wrong results with 3735, now I get correct results with 3729. I expect if you change from 3734 to 3728, you will get the correct results.

EPSG:3728: NAD83(NSRS2007) / Ohio North (ftUS) EPSG:3729: NAD83(NSRS2007) / Ohio South (ftUS) EPSG:3734: NAD83 / Ohio North (ftUS) EPSG:3735: NAD83 / Ohio South (ftUS)

I used your provided lat, long and am off by less than one foot.

p2 = pyproj.Proj( init="epsg:3728", preserve_units=True )

p2(-81.142311,41.688205)

(2339326.6558868014, 739401.4226131936)

vs 2339327.3, 739400.91