Conversion des coordonnées projetées en lat / lon en utilisant Python?

51

Je suppose que c'est une question fondamentale mais je n'arrive pas à trouver ou à reconnaître la solution.

Ce site retourne

Point:
X: -11705274.6374
Y: 4826473.6922

lorsque vous recherchez avec la première valeur de clé 000090 comme exemple. Je suppose que ceci est une référence spatiale et je comprends un peu ce que c'est.

Je cherche des instructions ou des exemples sur la façon de convertir cela en latitude et en longitude à l'aide de Python.

Vincent
la source

Réponses:

100

Le moyen le plus simple de transformer les coordonnées en Python est pyproj , c’est-à-dire l’ interface Python vers la bibliothèque PROJ.4 . En réalité:

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

résultats -105.150271116 39.7278572773

Antonio Falciano
la source
4
Oui. Pyproj jusqu'au bout.
sgillies
Cela fonctionne pour moi
lenhhoxung
36

Par défaut, le site auquel vous avez accédé utilise le système de référence spatiale EPSG 3857 (WGS84 Web Mercator). J'ai trouvé cette information ici .

Vous pouvez soit spécifier un autre système de référence spatiale en entrant le fichier EPSG souhaité dans le formulaire sous, Spatial Referencesoit convertir les coordonnées renvoyées avec Python.

Par exemple, vous pouvez utiliser les liaisons GDAL Python pour convertir ce point du système de coordonnées projeté (EPSG 3857) en un système de coordonnées géographiques (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()

Ceci retourne pour votre point les coordonnées de -105.150271116 39.7278572773.

ustroetz
la source
8

afalciano a la bonne réponse mais souhaite inclure une variante d'utilisation de pyproj.

Il fait besoin que vous connaissez la chaîne de proj4 et est un petit peu plus rapide.

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
la source
2
Vous n'avez pas besoin de la chaîne proj4, remplacez la deuxième ligne par p = pyproj.Proj(init='epsg:3857')et le résultat est identique.
alphabetasoup
1
Le résultat est le même, mais la dernière fois que j'ai vérifié, c'était un peu plus rapide.
Marcel Wilson
6

La sortie n'est pas un système de référence spatial / de coordonnées , c'est une paire de coordonnées. Vous devez savoir quelle est la référence spatiale pour reprojeter les coordonnées.

Cependant, ce n'est pas nécessaire dans ce cas. Il suffit de transmettre une référence spatiale de sortie appropriée au service et celui-ci renverra les coordonnées en lon / lat.

Voici la page avec les coordonnées de sortie au format Lon / Lat utilisant le système de référence géographique WGS-84 ( EPSG 4326 ).

utilisateur2856
la source
4

Essayé le code suggéré par Marcel Wilson et c'est plus rapide:

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
utilisateur1174460
la source
1

J'ai trouvé ce post en cherchant des moyens de le faire au sein de QGIS. Comme décrit ici , la méthode utilisée ressemble à ceci:

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
la source
2
Remarque, il y a un changement radical d'API dans QGIS 3, donc si v3.xvous en avez besoin, vous devrez utiliserxform = QgsCoordinateTransform(crsSrc, crsDest, QgsProject.instance())
Jonny
0

S'il vous plaît noter que la transformfonction de pyprojaccepte également arrays, ce qui est très utile quand il s'agit de dataframes

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
la source
-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
la source