Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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
share|improve this question
    
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

1 Answer 1

up vote 10 down vote accepted

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
share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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