7

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
6
  • 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
    Commented May 20, 2013 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
    Commented May 20, 2013 at 20:48
  • Looks like that link has the answer with ogr2ogrand -s_srs EPSG:3857 with -t_srs EPSG:4326
    – Mintx
    Commented May 20, 2013 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
    Commented May 20, 2013 at 20:55
  • You could make a system call with os.system()
    – Mintx
    Commented May 20, 2013 at 20:58

1 Answer 1

12

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
1
  • 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
    Commented Jun 20, 2018 at 10:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.