21

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

  • I've tried srs.GetAttrValue('AUTHORITY'), but this just returns 'EPSG' which is correct. EPSG is the authority – nmtoken Jun 17 at 14:36
29

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'
  • 1
    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! – Tomas Feb 14 '12 at 14:14
  • 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/… – Burrow Jun 28 '17 at 20:42
5

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
0

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.GetAuthorityCode(None)
Out[5]: 'EPSG'

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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