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
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.
No comments:
Post a Comment