Extracting EPSG from a raster using gdal bindings in Python [duplicate]

Question

My goal is to check that a shapefile and raster have the same projection and datum using EPSG number.

One can easily extract EPSG number from a shapefile by the following:

ds=ogr.Open(r'....satates48.shp')
utm_sr=ds.GetSpatialReference()
utm_sr.GetAttrValue('AUTHORITY',1)

This works great but how does one do the same for a raster file, you cannot, to my knowledge get a spatial reference of a raster.

One can get projection, or projecitonref (I am not sure if there is any difference), but you cannot do more, like:

fn_raster=gdal.open(r'raster_file.tif')
raster_sr=fn_raster.GetProjection().GetAttrValue('AUTHORITY',1)

This does not work. Perhaps there is an alternative way to do this?

Answer 1

I found the following workaround. I am unsure if it is the most efficient, but it does work for me.

from osgeo import gdal, osr

path = r"C:\temp\test2.tif"
d = gdal.Open(path)
proj = osr.SpatialReference(wkt=d.GetProjection())
print(proj.GetAttrValue('AUTHORITY',1))

EDIT: Slightly more condensed

Answer 2

You can also do this in one line using gdal.Info as so:

epsg = int(gdal.Info(input, format='json')['coordinateSystem']['wkt'].rsplit('"EPSG","', 1)[-1].split('"')[0])

This does essentially the same thing as the answer provided by Lennert. It reads the WKT representation of the file's spatial reference, then parses the string to extract the EPSG.

Note that this method is slightly more robust because gdal.Info is very flexible with the types of inputs it can take. input may be a datasource, valid filepath, or anything within gdal's extensive virtual file systems, which allows this method to accommodate a large number of potential use cases.

Answer 3

I would highly recommend using rasterio, which has a handy python API. For example:

>>> with rasterio.open('/path/to/your/geotiff.tif') as src:
        print (src.crs)

CRS({'init': u'epsg:32615'})