Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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?

share|improve this question

3 Answers 3

up vote 6 down vote accepted

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.

share|improve this answer

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

share|improve this answer

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.

share|improve this answer
    
Thank you very much. The preserve_units argument definitely did the trick, but the coordinates are still incorrect. @MerseyViking This answer gave me the correct coordinates. I wish I could mark both answers as the answer because they both helped. –  Brian May 25 '11 at 14:46
    
Well if people upvote @geographika's answer more than mine, it'll all pan out :) Glad it all worked tho. –  MerseyViking May 25 '11 at 15:19
    
since the link is broken, it might be helpful to show that you can write: p1 = Proj( init="epsg:3734", preserve_units=True ) –  BenjaminGolder Jan 20 at 5:21

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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