6

This question already has an answer here:

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?

marked as duplicate by user2856, Kersten, whyzar, nmtoken, Simbamangu Jan 9 '18 at 13:25

This question has been asked before and already has an answer. If those answers do not fully address your question, please ask a new question.

  • 1
    Use osr.SpatialReference - gis.stackexchange.com/questions/60371/… – user2856 Jan 9 '18 at 4:34
  • 1
    @Luke why was this marked as duplicate? These aren't duplicate questions at all. One finds the EPSG value, one finds the name of the Spatial Reference. Two completely different use cases. Two different workflows. I would mark it as duplicate if the other question was related to EPSG, but its not....so these two questions aren't really duplicates are they. I actually don't see another post on stackexchange discussing how to programatically extract an EPSG value with the gdal python bindings. – onakua Jan 17 '18 at 16:41
  • @onakua both questions are essentially "How do I get some attribute from a raster coordinate system" and all the answers to this question are applicable to the other question. You are welcome to raise this further by asking a question in meta or even nominate the question for reopening. – user2856 Jan 17 '18 at 23:22
  • 2
    Wait I'm confused none of the other answers in the other question extract an EPSG value :). GDAL/OGR has very poor documentation. Parsing different elements out of a SRS object is valuable, there are many elements stored in each SRS object. Having documented use-cases for extracting individual elements is valuable, especially for a site geared towards providing answers and information to an audience who isn't as familiar with GDAL/OGR and may not fully understand the nuances of working with coordinate systems. You are splitting hairs for the purpose of splitting hairs. – onakua Jan 17 '18 at 23:48
  • @onakua please don't use comments to try to argue about this. As I said, you are welcome to raise this in meta. – user2856 Jan 18 '18 at 2:33
9

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

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

  • It is a workaround, and perhaps that is the intent of GDAL developers, lets see if any direct answers come up, but thanks! – icypy Jan 9 '18 at 5:06
  • 3
    Just to add to this: I recognized that in some cases it might be necessary to call proj.AutoIdentifyEPSG() before calling proj.GetAttrValue('AUTHORITY',1). – s6hebern Mar 26 '18 at 13:28
7

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.

  • +1 but note this only works with GDAL 2.0+ – user2856 Jan 9 '18 at 8:10
5

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'})
  • I totally agree but the question specifically asks for GDAL. – bugmenot123 Jan 9 '18 at 10:01
  • That is good to know about rasterio, but yes the question was about gdal. And now I realize, that ogr is the main module for spatial reference and not gdal, it is weird. – icypy Jan 17 '18 at 21:44

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