Saturday 1 July 2017

GDAL/Python: How do I get coordinate system name from SpatialReference?


In Python, using GDAL, I've extracted a raster's projection as a WKT string as follows:



wkt = dataset.GetProjection()
# wkt is 'PROJCS["GDA_1994_Transverse_Mercator",GEOGCS["GDA_1994",DATUM["GDA_1994",SPHEROID["GRS_1980",6378137,298.2572221010002],TOWGS84[0,0,0,0,0,0,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",117],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["Meter",1]]'

Using the WKT string, I can create a SpatialReference instance as follows:


src = osr.SpatialReference()
src.ImportFromWkt(wkt)

This is easy-peasy. I can extract several parameters of the projection from src such as the UTM zone, etc., easily enough. But I can't figure out how to extract something like the name of the projection, i.e. "GDA_1994_Transverse_Mercator". This has surely got to be possible, but the Python API documentation may as well be nonexistent for all the use it is.


How do I extract the names of the projection and the geographic coordinate system?



Answer




See the OGR Projections tutorial and the OGRSpatialReference class. In particular, the GetAttrValue method.


Here's a worked example.


from osgeo import gdal,osr
ds=gdal.Open(r'SOMERASTER.TIF')
prj=ds.GetProjection()
print prj

srs=osr.SpatialReference(wkt=prj)
if srs.IsProjected:
print srs.GetAttrValue('projcs')

print srs.GetAttrValue('geogcs')

For my raster this prints:


PROJCS["WGS 84 / UTM zone 55N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",147],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","32655"]]
'WGS 84 / UTM zone 55N'
'WGS 84'

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