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
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
No comments:
Post a Comment