Converting projected coordinates to lat/lon using Python

Question

Point:
X: -11705274.6374
Y: 4826473.6922

when you search with the first key value of 000090, as an example. I assume that this is a spatial reference.

I am looking for instructions, or examples, of how to convert this to latitude and longitude using Python.

That site returns: "Access to the application you were trying to use has been blocked in accordance with County policy. Please contact your organization's Service Desk if you believe this is in error." — Daniel F, Jul 06 '22 at 04:43

score 159 · Accepted Answer · edited May 13 '23 at 12:12

159

The simplest way to transform coordinates in Python is pyproj, i.e. the Python interface to PROJ.4 library.

Deprecated solution

The solution originally posted in this answer is shown below. This method has been discouraged since version 2.0, and some parts now result in deprecation warnings.

from pyproj import Proj, transform
inProj = Proj(init='epsg:3857') # deprecated since 2.4
outProj = Proj(init='epsg:4326') # deprecated since 2.4
inProj = Proj('epsg:3857')

outProj = Proj('epsg:4326') 
x1,y1 = -11705274.6374,4826473.6922
x2,y2 = transform(inProj,outProj,x1,y1) # deprecated some time later than 2.4
print x2,y2 # -105.150271116 39.7278572773

edited May 13 '23 at 12:12

Vince

20,017
15
45
64

answered Nov 29 '13 at 10:04

Antonio Falciano

14,333
2
36
66

pyproj 2.4 gives a FutureWarning about deprecated Proj() initialization with the init= syntax. The updated syntax is identical but without the init=. Like this:
inProj = Proj('epsg:3857') and outProj = Proj('epsg:4326')
– Marc Compere Dec 19 '19 at 21:35
1

I removed the 'init=' based on the edit, but my output of transform(inProj,outProj,x1,y1) gave me (latitude, longitude) instead of (longitude, latitude). Anybody knows why ? – gdelab Feb 18 '20 at 15:18
2

@gdelab See https://github.com/pyproj4/pyproj/issues/538#issuecomment-585734088 – Antonio Falciano Feb 22 '20 at 17:17
1

This might be helpful to reference: https://pyproj4.github.io/pyproj/stable/gotchas.html#upgrading-to-pyproj-2-from-pyproj-1 – snowman2 Jan 26 '21 at 01:20
1

This will also be good to reference: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6 – snowman2 Jan 26 '21 at 01:21
1

Now transform is also deprecated: https://pyproj4.github.io/pyproj/stable/gotchas.html#upgrading-to-pyproj-2-from-pyproj-1 – zabop Jan 13 '22 at 16:10
As zabob said, now, recommended usage is: transformer = Transformer.from_crs("EPSG:4326", "EPSG:3857") x2, y2 = transformer.transform(x1, y1) – Alexey Mar 23 '23 at 16:20

score 41 · Answer 2 · edited Jan 28 '22 at 12:55

By default the site you linked to uses the Spatial Reference System EPSG 3857 (WGS84 Web Mercator). I found this information here.

You can either specify another Spatial Reference System by entering the desired EPSG into the form under Spatial Reference or you can convert the returned coordinates with Python.

For instance you can use the GDAL Python bindings to convert this point from the projected coordinate system (EPSG 3857) to a geographic coordinate system (EPSG 4326).

from osgeo import ogr
from osgeo import 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()

This returns for your point the coordinates of -105.150271116 39.7278572773.

@Thomas it is returning 39.7278572773 -105.150271116, why is it flipping coordinates? — Ayaz49, Mar 27 '22 at 10:23
The only comment is about the print, just add brackets to print without errors (for python 3 +) — Rafa, Apr 01 '23 at 16:15

score 20 · Answer 3 · answered Feb 10 '20 at 02:28

When using pyproj, note the differences from various releases in how it is used to transform data. Here are a few examples using new/old capabilities based on the question:

Using pyproj >= 2.2.0

import pyproj
print(pyproj.__version__)  # 2.4.1
print(pyproj.proj_version_str) # 6.2.1

proj = pyproj.Transformer.from_crs(3857, 4326, always_xy=True)

x1, y1 = (-11705274.6374, 4826473.6922)
x2, y2 = proj.transform(x1, y1)
print((x2, y2))  # (-105.15027111593008, 39.72785727727918)

Using pyproj <= 1.9.6

import pyproj
print(pyproj.__version__)  # 1.9.6
print(pyproj.proj_version_str) # 4.9.3

inProj = pyproj.Proj(init='epsg:3857')
outProj = pyproj.Proj(init='epsg:4326')

x1, y1 = (-11705274.6374, 4826473.6922)
x2, y2 = pyproj.transform(inProj, outProj, x1, y1)
print((x2, y2))  # (-105.15027111593008, 39.72785727727918)

There are a few considerations with the different versions of PROJ/pyproj:

transform for older versions always return the same axis order of "x, y" or "longitude, latitude", whereas PROJ 6+ is the order as defined by EPSG, unless an option like always_xy=True is specified
Proj is limited to converting between geographic and projection coordinates within one datum, whereas the newer Transformer takes into account datum shifts, and is recommended for newer pyproj use

score 13 · Answer 4 · answered Apr 15 '16 at 16:26

13

afalciano has the right answer but wanted to include a variant usage of pyproj.

It does require you know the proj4 string and is a tiny bit faster.

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

answered Apr 15 '16 at 16:26

Marcel Wilson

329
3
6

3

You don't need the proj4 string, substitute the second line for p = pyproj.Proj(init='epsg:3857') and the result is the same. – alphabetasoup Sep 15 '16 at 00:08
1

The result is the same but last I checked this was a tiny bit faster. – Marcel Wilson Sep 17 '16 at 17:03
4

Even though it was not the intention of this answer, it is particularly useful when you have a totally custom projection, which is not listed in the EPSG-repository. – Andreas L. Jan 13 '20 at 13:51

user2856 · Answer 5 · 2013-11-28T04:21:11.270

The output is not a spatial/coordinate reference system, it's a pair of coordinates. You need to know what the spatial reference is to reproject the coordinates.

However, that's not required in this case. Just pass an appropriate output spatial reference to the service and it will return the coordinates in Lon/Lat.

Here is the page with output coordinates in Lon/Lat format using the WGS-84 geographic spatial reference system (EPSG 4326).

score 7 · Answer 6 · answered Feb 20 '19 at 10:24

Tried the code suggested by Marcel Wilson and it is faster:

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

score 5 · Answer 7 · answered Dec 09 '16 at 10:27

5

I found this post when looking for ways of doing this within QGIS. As described here, the method used looks like so:

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)

answered Dec 09 '16 at 10:27

Toivo Säwén

325
2
9

3

Note, there's a breaking API change in QGIS 3, so if using v3.x you'll need to use xform = QgsCoordinateTransform(crsSrc, crsDest, QgsProject.instance()) – Jonny Jan 17 '19 at 15:38
1

Also see PyQgis Cookbook >> Projection Support >> CRS Transformation, https://docs.qgis.org/3.16/en/docs/pyqgis_developer_cookbook/crs.html#crs-transformation, and how to transform an entire layer here on Stack: https://gis.stackexchange.com/questions/109235/how-to-make-crs-transformation-of-a-qgsvectorlayer – matt wilkie Jan 30 '21 at 18:48

score 5 · Answer 8 · answered Jul 01 '19 at 07:33

Please note that the transform function of pyproj accepts also arrays, which is quite useful when it comes to 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())

score 1 · Answer 9 · answered Jul 06 '22 at 04:52

This is a shorter version of ustroetz's answer since the ogr.wkbPoint geometry is not required for a point.

from osgeo import osr
pointX = -11705274.6374 
pointY = 4826473.6922
pointZ = 0.0
Spatial Reference System
inputEPSG = 3857
outputEPSG = 4326
create coordinate transformation
inSpatialRef = osr.SpatialReference()
inSpatialRef.ImportFromEPSG(inputEPSG)
outSpatialRef = osr.SpatialReference()
outSpatialRef.ImportFromEPSG(outputEPSG)
coordTransform = osr.CoordinateTransformation(inSpatialRef, outSpatialRef)
transform point, height z is optional
x, y, z = coordTransform.TransformPoint(pointX, pointY, pointZ)
print point in EPSG 4326
print(x, y)

score -2 · Answer 10 · answered Jan 16 '22 at 15:32

convert esri extent bbox to lat long

use esri javascipt api module ( client side projection, fast, no need ajax, no need internet, projection done in your browser, extremely fast)

first include css, js

link(rel='stylesheet',href='https://js.arcgis.com/4.17/esri/themes/light/main.css') script(src='https://js.arcgis.com/4.17/')

This is my working example:

proj4js ( need ajax, internet to read projection string from website, due to website does not allow jsonp, not allow cors, you have to setup your own proxy to bypass CORS problem, not easy for normal user )

But if you already known from what to what, then you can just copy past projection string, no internet needed, all projectin done in your browser.

If you don't know project from what to what, projection wkid is dynamic generated in the runtime, then like I said previously, you must need ajax, need internet to read projection string from website on the fly, in the run time. It is not easy, because you have to bypass CORS problem.

Here is my working example:

use Arcgis.com hosted web service ( need internet, very slow, not reliable, if ESRI server is crowded, you might get it very slow. Every ajax request can only conver one pair of lat lng)

Here is how to project on website

you can get xmin, ymin from here

Downvoted because the question is tagged python and this answer uses JavaScript exclusively. Also the screenshots... please don't do this on Stack Exchange. — Daniel F, Jul 06 '22 at 05:03

Converting projected coordinates to lat/lon using Python

10 Answers10

Recommended solution since pyproj 2

Deprecated solution

inProj = Proj(init='epsg:3857') # deprecated since 2.4

outProj = Proj(init='epsg:4326') # deprecated since 2.4

Spatial Reference System

create a geometry from coordinates

create coordinate transformation

transform point

print point in EPSG 4326

Using pyproj >= 2.2.0

Using pyproj <= 1.9.6

Spatial Reference System

create coordinate transformation

transform point, height z is optional

print point in EPSG 4326

Linked