Converter as coordenadas X, Y em lat / long usando pyproj e Proj.4 retorna as coordenadas incorretas

10

Estou escrevendo um script python que lê vários arquivos XML contendo as coordenadas xey e os combina em um único arquivo csv. Latitude e Longitude são campos obrigatórios no csv, mas estou tendo dificuldades para converter as coordenadas x, y no plano do estado norte de Ohio usFt para 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 os métodos acima retornam o mesmo resultado, no entanto, esse período longo está em algum lugar na Baía de Hudson. Quando plogo as coordenadas no ArcMap, o comprimento correto do lat é: -81.142311,41.688205.

* Observe que todos os lat longs são fornecidos long, lat, pois esta é a ordem que o Proj usa

Alguém sabe por que eu estaria recebendo as coordenadas erradas do Proj.4 e pyproj?

Brian
fonte

Respostas:

8

Eu obtenho os mesmos resultados que o @geographika quando executo gdaltransforme a ferramenta 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

A reversão das coordenadas xey do seu ponto, no entanto, fornece o resultado que você está vendo no ArcMap:

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

Portanto, você precisará fazer uma verificação visual para garantir que as coordenadas x e y sejam o caminho certo. É um problema que tive no passado quando você obtém dois resultados semelhantes o suficiente para atribuir a erro de arredondamento ou algo assim.

MerseyViking
fonte
19

PyProj assume que suas coordenadas estão em metros. Eu acho que algo relacionado a pés / metros é a causa do problema.

Chamando uma instância da classe Proj com os argumentos lon, lat converterá lon / lat (em graus) em x / y coordenadas de projeção de mapa nativas (em metros)

Se a palavra-chave opcional 'preserve_units' for True, as unidades nas coordenadas de projeção do mapa não serão forçadas a serem metros.

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

Suas coordenadas iniciais estão em pés? Quando você carrega os dados no ArcMap, quais unidades o mapa usa?

Isso aproxima as coordenadas um pouco:

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)

Um problema semelhante pode ser encontrado aqui .

geographika
fonte
Muito obrigado. O argumento preserve_units definitivamente fez o truque, mas as coordenadas ainda estão incorretas. @MerseyViking Esta resposta me deu as coordenadas corretas. Eu gostaria de poder marcar as duas respostas como a resposta, porque ambas ajudaram.
27511 Brian
Bem, se as pessoas votam mais na resposta da @ geographika do que na minha, tudo vai dar certo :) Que bom que tudo deu certo.
MerseyViking
uma vez que a ligação é interrompida, pode ser útil para mostrar que você pode escrever:p1 = Proj( init="epsg:3734", preserve_units=True )
BenjaminGolder
4

Na verdade, eu estava tentando fazer a mesma coisa, exceto com a grade de avião do estado sul de OH e me deparei com sua pergunta. Eu estava obtendo resultados errados com o 3735, agora obtenho resultados corretos com o 3729. Espero que, se você mudar de 3734 para 3728, obterá os resultados corretos.

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

Eu usei o seu lat fornecido, longo e estou fora por menos de um pé.

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

p2 (-81.142311,41.688205)

(2339326.6558868014, 739401.4226131936)

vs 2339327.3, 739400.91

engenheiro de minas
fonte