Ist es möglich, den EPSG-Wert mithilfe der OGR-Python-API von einer OSR-SpatialReference-Klasse abzurufen?

21

Beim Lesen eines Layers aus einer OGR-PostGIS-Verbindung kann ich die räumliche Referenz des Layers abrufen. Ist es jedoch möglich, den EPSG-Wert abzurufen? Gibt es dazu Unterlagen?

Beispielsweise:

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

Kehrt zurück:

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

Wie erhalte ich den EPSG-Wert für die Projektion? Z.B:

srs.GetEPSG()
print srs
27700

Ich habe es versucht srs.GetAttrValue('AUTHORITY'), aber das kommt gerade zurück 'EPSG'.

python postgis ogr Tomas
quelle
I've tried srs.GetAttrValue('AUTHORITY'), but this just returns 'EPSG'welches ist richtig. EPSG ist die Behörde
nmtoken

Antworten:

30

Es ist etwas vergraben, aber es gibt einen zweiten Parameter für GetAttrValue (), der den Wert an dieser Ordnungszahl zurückgibt. Also kann ich tun:

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'

Nach einigem Herumspielen habe ich festgestellt, dass Sie den Wert für jeden Parameter mithilfe einer Pipe |als Pfadtrennzeichen abrufen können:

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

Was beim Auffinden des geografischen Koordinatensystems eines projizierten CS hilfreich sein kann:

In [13]: srs.GetAttrValue("PROJCS|GEOGCS|AUTHORITY", 1)
Out[13]: '4277'
MerseyViking
quelle
1
Danke, das ist großartig. Ich werde es umsetzen. Ich hatte keine Zeit mehr, um weiter herumzuspielen - die schnelle Anwendungsentwicklung wurde wieder durch das Fehlen von GDAL / OGR-Dokumentation aufgehalten!
Tomas
Ich habe die GetAttrValue-Funktion mit den Argumenten "AUTHORITY" und "1" ausprobiert und festgestellt, dass der EPSG-Code nicht immer zurückgegeben wird, da der EPSG-Code nicht immer in der WKT enthalten ist. Ich bin immer noch ein bisschen verwirrt, warum das so ist. Ich fand die folgende Lösung gut für meine Bedürfnisse: gis.stackexchange.com/questions/7608/…
Burrow
5

Hier ist ein Codeausschnitt, der für mich funktioniert hat:

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
tomkralidis
quelle
0

SpatialReference.GetAuthorityCode()Nimmt Noneals Parameter, der einen Berechtigungsknoten auf dem Stammelement findet (dh je nach Bedarf projiziert / geografisch). Gleiches gilt für 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'
rcoup
quelle