La conversión de coordenadas X, Y a lat / long usando pyproj y Proj.4 devuelve las coordenadas incorrectas

10

Estoy escribiendo un script de Python que lee múltiples archivos XML que contienen coordenadas x e y y los combina a todos en un solo archivo csv. Latitud y longitud son campos obligatorios en el csv, pero tengo dificultades para convertir las coordenadas x, y en el plano del estado de Ohio North usFt a 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)

Ambos métodos anteriores devuelven el mismo resultado, sin embargo, este último tiempo está en algún lugar de la Bahía de Hudson. Cuando trazo las coordenadas en ArcMap, el lat lat correcto es: -81.142311,41.688205.

* Observe que todos los lat largos se proporcionan largos, lat ya que este es el orden que usa Proj

¿Alguien sabe por qué obtendría las coordenadas incorrectas de Proj.4 y pyproj?

Brian
fuente

Respuestas:

8

Obtengo los mismos resultados que @geographika cuando ejecuto gdaltransformy la herramienta proj.4 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

Sin embargo, invertir las coordenadas x e y de su punto da el resultado que está viendo en ArcMap:

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

Por lo tanto, deberá realizar una verificación visual para asegurarse de que sus coordenadas x e y sean correctas. Es un problema que tuve en el pasado cuando obtienes dos resultados que son lo suficientemente similares como para atribuirlo a un error de redondeo o algo así.

MerseyViking
fuente
19

PyProj supone que sus coordenadas están en metros. Supongo que algo relacionado con pies / metros es la causa del problema.

Al llamar a una instancia de clase Proj con los argumentos lon, lat convertirá lon / lat (en grados) a coordenadas de proyección del mapa nativo x / y (en metros)

Si la palabra clave opcional 'preserve_units' es True, las unidades en las coordenadas de proyección del mapa no están obligadas a ser metros.

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

¿Tus coordenadas iniciales están en pies? Cuando carga los datos en ArcMap, ¿qué unidades usa el mapa?

Esto hace que las coordenadas estén un poco más cerca:

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)

Un problema similar se puede encontrar aquí .

geographika
fuente
Muchas gracias. El argumento preserve_units definitivamente hizo el truco, pero las coordenadas aún son incorrectas. @MerseyViking Esta respuesta me dio las coordenadas correctas. Desearía poder marcar ambas respuestas como la respuesta porque ambas ayudaron.
Brian
Bueno, si la gente vota más la respuesta de @ geographika que la mía, todo saldrá bien :) Me alegro de que todo haya funcionado.
MerseyViking
Dado que el enlace está roto, puede ser útil mostrar que puede escribir:p1 = Proj( init="epsg:3734", preserve_units=True )
BenjaminGolder
4

En realidad estaba tratando de hacer lo mismo, excepto con la cuadrícula del avión del estado del sur de OH y me encontré con su pregunta. Estaba obteniendo resultados incorrectos con 3735, ahora obtengo resultados correctos con 3729. Espero que si cambias de 3734 a 3728, obtendrás los resultados correctos.

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)

Utilicé su lat proporcionado, largo y estoy fuera por menos de un pie.

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

p2 (-81.142311,41.688205)

(2339326.6558868014, 739401.4226131936)

vs 2339327.3, 739400.91

ingeniero de minas
fuente