Python – How to Transform Coordinates from Meter to Decimal Degrees using OGR

coordinate systemgdalogrpython

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

Best Answer

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
Related Question