Is it possible to get the EPSG value from an OSR SpatialReference class using the OGR Python API?

Question

When reading a layer from a OGR PostGIS connection I can get the SpatialReference of the layer, but is it possible to get the EPSG value? Is there any documentation on this?

For example:

lyr = conn.GetLayerByName(tbl) # Where conn is OGR PG connection
srs = ly.GetSpatialRef()
print srs

Returns:

PROJCS["OSGB 1936 / British National Grid",
GEOGCS["OSGB 1936",
 DATUM["OSGB_1936",
 SPHEROID["Airy 1830",6377563.396,299.3249646,
 AUTHORITY["EPSG","7001"]],
 AUTHORITY["EPSG","6277"]],
 PRIMEM["Greenwich",0,
 AUTHORITY["EPSG","8901"]],
 UNIT["degree",0.01745329251994328,
 AUTHORITY["EPSG","9122"]],
 AUTHORITY["EPSG","4277"]],
UNIT["metre",1,
 AUTHORITY["EPSG","9001"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",49],
PARAMETER["central_meridian",-2],
PARAMETER["scale_factor",0.9996012717],
PARAMETER["false_easting",400000],
PARAMETER["false_northing",-100000],
AUTHORITY["EPSG","27700"],
AXIS["Easting",EAST],
AXIS["Northing",NORTH]]

So how do I get the EPSG value for the projection? E.g.:

srs.GetEPSG()
print srs
27700

I've tried srs.GetAttrValue('AUTHORITY'), but this just returns 'EPSG'.

• nmtoken: I've tried srs.GetAttrValue('AUTHORITY'), but this just returns 'EPSG' which is correct. EPSG is the authority

Answer 1

It's a bit buried, but there is a second parameter to GetAttrValue() which returns the value at that ordinal. So I can do:

In [1]: import osgeo.osr as osr
In [2]: srs = osr.SpatialReference()
In [3]: srs.SetFromUserInput("EPSG:27700")
Out[3]: 0
In [4]: print srs
PROJCS["OSGB 1936 / British National Grid",
 GEOGCS["OSGB 1936",
 DATUM["OSGB_1936",
 SPHEROID["Airy 1830",6377563.396,299.3249646,
 AUTHORITY["EPSG","7001"]],
 TOWGS84[375,-111,431,0,0,0,0],
 AUTHORITY["EPSG","6277"]],
 PRIMEM["Greenwich",0,
 AUTHORITY["EPSG","8901"]],
 UNIT["degree",0.0174532925199433,
 AUTHORITY["EPSG","9122"]],
 AUTHORITY["EPSG","4277"]],
 PROJECTION["Transverse_Mercator"],
 PARAMETER["latitude_of_origin",49],
 PARAMETER["central_meridian",-2],
 PARAMETER["scale_factor",0.9996012717],
 PARAMETER["false_easting",400000],
 PARAMETER["false_northing",-100000],
 UNIT["metre",1,
 AUTHORITY["EPSG","9001"]],
 AXIS["Easting",EAST],
 AXIS["Northing",NORTH],
 AUTHORITY["EPSG","27700"]]
In [5]: srs.GetAttrValue("AUTHORITY", 0)
Out[5]: 'EPSG'
In [6]: srs.GetAttrValue("AUTHORITY", 1)
Out[6]: '27700'

After a bit of playing about, I've found you can get the value for any parameter using a pipe | as a path separator:

In [12]: srs.GetAttrValue("PRIMEM|AUTHORITY", 1)
Out[12]: '8901'

Which may be of use in finding the geographic coordinate system of a projected CS:

In [13]: srs.GetAttrValue("PROJCS|GEOGCS|AUTHORITY", 1)
Out[13]: '4277'

• Tomas: Thanks, that's great. I will implement it. I had run out of time for further 'playing about' - rapid application development being held up by lack of GDAL/OGR documentation once again!
• Burrow: I tried the GetAttrValue function with the "AUTHORITY" and "1" arguments and noticed it does not always return the EPSG code because the EPSG code is not always included in the WKT. I'm still a little fuzzy about why this is the case. I found the following solution to work well for my needs: gis.stackexchange.com/questions/7608/…

Answer 2

SpatialReference.GetAuthorityCode() takes None as a parameter, that finds an authority node on the root element (ie. projected/geographic as appropriate). Same applies to GetAuthorityName():

In [1]: import osgeo.osr as osr
In [2]: srs = osr.SpatialReference()
In [3]: srs.SetFromUserInput("EPSG:27700")
Out[3]: 0
In [4]: srs.GetAuthorityCode(None)
Out[4]: '27700'
In [5]: srs.GetAuthorityName(None)
Out[5]: 'EPSG'

• nickves: None means that it will search in the root by itself to find the authority code (which in this context is epsg)

Answer 3

Here's a code snippet that has worked for me:

def wkt2epsg(wkt, epsg='/usr/local/share/proj/epsg', forceProj4=False):
''' Transform a WKT string to an EPSG code
Arguments
---------
wkt: WKT definition
epsg: the proj.4 epsg file (defaults to '/usr/local/share/proj/epsg')
forceProj4: whether to perform brute force proj4 epsg file check (last resort)
Returns: EPSG code
'''
code = None
p_in = osr.SpatialReference()
s = p_in.ImportFromWkt(wkt)
if s == 5: # invalid WKT
 return None
if p_in.IsLocal() == 1: # this is a local definition
 return p_in.ExportToWkt()
if p_in.IsGeographic() == 1: # this is a geographic srs
 cstype = 'GEOGCS'
else: # this is a projected srs
 cstype = 'PROJCS'
an = p_in.GetAuthorityName(cstype)
ac = p_in.GetAuthorityCode(cstype)
if an is not None and ac is not None: # return the EPSG code
 return '%s:%s' % \
 (p_in.GetAuthorityName(cstype), p_in.GetAuthorityCode(cstype))
else: # try brute force approach by grokking proj epsg definition file
 p_out = p_in.ExportToProj4()
 if p_out:
 if forceProj4 is True:
 return p_out
 f = open(epsg)
 for line in f:
 if line.find(p_out) != -1:
 m = re.search('<(\\d+)>', line)
 if m:
 code = m.group(1)
 break
 if code: # match
 return 'EPSG:%s' % code
 else: # no match
 return None
 else:
 return None