Sunday, 18 February 2018

gdal - 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

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

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...