Python - OGR: Transform coordinates from meter to decimal degrees

I have a point in a shapefile. The shapefile is in a WGS_1984_Web_Mercator Projected Coordinate System. I want to get the coordinates of that point in decimal degrees with GDAL/OGR in Python.

I am able to get the coordinates in meters with the code below.

I tried to transform the shapefile from a Projected Coordinate System to a Geographic Coordinate System without success. I assume it is not possible to remove a Projected Coordinate System with OGR. Is that right?

What other ways are there within OGR/Python to get the coordinates in decimal degrees?

import ogr

driver = ogr.GetDriverByName('ESRI Shapefile')

shp = driver.Open('testpoint.shp', 0)
lyr = shp.GetLayer()

feat = lyr.GetNextFeature()
geom = feat.GetGeometryRef()

x = geom.GetX() #Get X coordinates
y = geom.GetY() #Get Y cooridnates

print x,y
• I'm not familiar enough with OGR to know if there's something out of the box, but you can get the WGS84 lat/long values from the Pseudo-Mercator XYs. Warning: there's quite a bit of math involved if you want to roll your own conversion. – Mintx May 20 '13 at 20:44
• Yes, I am aware of that. There is a short explanation here how to do the math: gis.stackexchange.com/a/30437/15607 But I would rather do it without calculating the coordinates manually. – ustroetz May 20 '13 at 20:48
• Looks like that link has the answer with ogr2ogrand -s_srs EPSG:3857 with -t_srs EPSG:4326 – Mintx May 20 '13 at 20:53
• Yeah, I saw that too. However I like to do it in Python and not in the command line. Do you know how to use ogr2ogr within Python? – ustroetz May 20 '13 at 20:55
• You could make a system call with os.system() – Mintx May 20 '13 at 20:58

In order to get the coordinates in decimal degrees, the data needs to be reprojected to WGS84.

import ogr, osr

driver = ogr.GetDriverByName('ESRI Shapefile')

shp = driver.Open('testpoint.shp', 0)
lyr = shp.GetLayer()

feat = lyr.GetNextFeature()
geom = feat.GetGeometryRef()

# Transform from Web Mercator to WGS84
sourceSR = lyr.GetSpatialRef()
targetSR = osr.SpatialReference()
targetSR.ImportFromEPSG(4326) # WGS84
coordTrans = osr.CoordinateTransformation(sourceSR,targetSR)
geom.Transform(coordTrans)

x = geom.GetX()
y = geom.GetY()

print x,y
• using this code, i only get 0.0 from all coordinates, should this code convert X_Coord and Y_Coord fields of the shape file ? – Elgert Jun 20 '18 at 10:43