Tuesday 24 October 2017

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


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

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...