Convertendo coordenadas projetadas para lat / lon usando Python?

51

Suponho que essa seja uma pergunta básica, mas não consigo encontrar ou reconhecer a solução.

Este site retorna

Point:
X: -11705274.6374
Y: 4826473.6922

ao pesquisar com o primeiro valor da chave 000090 como exemplo. Eu acho que essa é uma referência espacial e eu meio que entendo o que é isso.

Estou procurando instruções ou exemplos de como converter isso para Latitude e longitude usando Python.

Vincent
fonte

Respostas:

100

A maneira mais simples de transformar coordenadas no Python é pyproj , ou seja, a interface do Python para a biblioteca PROJ.4 . De fato:

from pyproj import Proj, transform

inProj = Proj(init='epsg:3857')
outProj = Proj(init='epsg:4326')
x1,y1 = -11705274.6374,4826473.6922
x2,y2 = transform(inProj,outProj,x1,y1)
print x2,y2

retorna -105.150271116 39.7278572773

Antonio Falciano
fonte
4
Sim. Pyproj todo o caminho.
sgillies
Funciona para mim
lenhhoxung
36

Por padrão, o site ao qual você vinculou usa o Sistema de Referência Espacial EPSG 3857 (WGS84 Web Mercator). Encontrei esta informação aqui .

Você pode especificar outro Sistema de Referência Espacial inserindo o EPSG desejado no formulário em Spatial Referenceou pode converter as coordenadas retornadas com Python.

Por exemplo, você pode usar as ligações GDAL Python para converter esse ponto do sistema de coordenadas projetadas (EPSG 3857) em um sistema de coordenadas geográficas (EPSG 4326).

import ogr, osr

pointX = -11705274.6374 
pointY = 4826473.6922

# Spatial Reference System
inputEPSG = 3857
outputEPSG = 4326

# create a geometry from coordinates
point = ogr.Geometry(ogr.wkbPoint)
point.AddPoint(pointX, pointY)

# create coordinate transformation
inSpatialRef = osr.SpatialReference()
inSpatialRef.ImportFromEPSG(inputEPSG)

outSpatialRef = osr.SpatialReference()
outSpatialRef.ImportFromEPSG(outputEPSG)

coordTransform = osr.CoordinateTransformation(inSpatialRef, outSpatialRef)

# transform point
point.Transform(coordTransform)

# print point in EPSG 4326
print point.GetX(), point.GetY()

Isso retorna para o seu ponto as coordenadas de -105.150271116 39.7278572773.

ustroetz
fonte
8

O afalciano tem a resposta certa, mas queria incluir um uso variante do pyproj.

Ele faz necessário que você sabe a corda proj4 e é um pouco mais rápido.

import pyproj
p = pyproj.Proj("+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378137 +b=6378137 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")
lon, lat = p(x, y, inverse=True)
print lat, lon
Marcel Wilson
fonte
2
Você não precisa da string proj4, substitua a segunda linha p = pyproj.Proj(init='epsg:3857')e o resultado é o mesmo.
alphabetasoup
11
O resultado é o mesmo, mas a última vez que verifiquei foi um pouco mais rápido.
Marcel Wilson
6

A saída não é um sistema de referência espacial / coordenada , é um par de coordenadas. Você precisa saber qual é a referência espacial para reprojetar as coordenadas.

No entanto, isso não é necessário neste caso. Basta passar uma referência espacial de saída apropriada para o serviço e ele retornará as coordenadas em Lon / Lat.

Aqui está a página com coordenadas de saída no formato Lon / Lat usando o sistema de referência geográfica geográfica WGS-84 ( EPSG 4326 ).

user2856
fonte
4

Tentei o código sugerido por Marcel Wilson e é mais rápido:

from pyproj import Proj, transform
import time
import pyproj


# Test 1 - 0.0006158 s
start=time.time()
inProj = Proj(init='epsg:3857')
outProj = Proj(init='epsg:4326')
x1,y1 = -11705274.6374,4826473.6922
x2,y2 = transform(inProj,outProj,x1,y1)
end=time.time()
print(y2,x2)
print('%.7f' % (end-start))

# Test 2 - 0.0000517 s --- factor 11,9
start=time.time()
x,y = -11705274.6374,4826473.6922
p = pyproj.Proj("+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378137 +b=6378137 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")
lon, lat = p(x, y, inverse=True)
end=time.time()
print(lat, lon)
print('%.7f' % (end-start))
-----------------

39.72785727727918 -105.15027111593008
0.0006158
39.72785727727918 -105.15027111593008
0.0000517
user1174460
fonte
1

Encontrei este post ao procurar maneiras de fazer isso no QGIS. Conforme descrito aqui , o método usado tem a seguinte aparência:

def convertProjection(self,x,y,from_crs,to_crs):
    crsSrc = QgsCoordinateReferenceSystem(from_crs)
    crsDest = QgsCoordinateReferenceSystem(to_crs)
    xform = QgsCoordinateTransform(crsSrc, crsDest)
    pt = xform.transform(QgsPoint(x,y))
    return pt.x, pt.y

# Remove the "EPSG:" part
from_crs = 3857
to_crs = 4326
x = -11705274.6374    
y = 4826473.6922
lon, lat = self.convertProjection(x,y,from_crs, to_crs)
Toivo Säwén
fonte
2
Observe que há uma alteração na API do QGIS 3, por isso, se estiver usando, v3.xvocê precisará usarxform = QgsCoordinateTransform(crsSrc, crsDest, QgsProject.instance())
Jonny
0

Observe que a transformfunção de accept pyprojtambém arrays, o que é bastante útil quando se trata de quadros de dados

import pandas as pd
from pyproj import Proj, transform

df = pd.DataFrame({'x': [-11705274.6374]*100, 
                   'y': [4826473.6922]*100})
inProj, outProj = Proj(init='epsg:3857'), Proj(init='epsg:4326')
df['x2'], df['y2'] = transform(inProj, outProj, df['x'].tolist(), df['y'].tolist())
J. Doe
fonte
-1
import cv2
import numpy as np
def onMouse(event, x, y, flags, param):
    if event == cv2.EVENT_LBUTTONDOWN:
       # draw circle here (etc...)
       print('x = %f, y = %f'%((x*2/100),(y*2/100)))
a=float(input("enter length of original dimension in cm"))
b=float(input("enter width of original dimension in cm"))
print("origional image coordinates")
im=cv2.imread('mask1.jpg',0)
re_im=cv2.resize(im,(round(a*100/2),round(b*100/2)))
cv2.imshow('image',re_im)
cv2.setMouseCallback('image', onMouse)
cv2.waitKey(0)
cv2.destroyAllWindows()
Subhash Verma
fonte