Czy można uzyskać wartość EPSG z klasy OSR SpatialReference przy użyciu interfejsu API Python OGR?

21

Czytając warstwę z połączenia OGR PostGIS, mogę uzyskać parametr SpatialReference warstwy, ale czy można uzyskać wartość EPSG? Czy jest na to jakaś dokumentacja?

Na przykład:

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

Zwroty:

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

Jak więc uzyskać wartość EPSG dla projekcji? Na przykład:

srs.GetEPSG()
print srs
27700

Próbowałem srs.GetAttrValue('AUTHORITY'), ale to po prostu powraca 'EPSG'.

python postgis ogr Tomas
źródło
I've tried srs.GetAttrValue('AUTHORITY'), but this just returns 'EPSG'który jest poprawny. EPSG jest autorytetem
nmtoken

Odpowiedzi:

30

Jest nieco zakopany, ale istnieje drugi parametr GetAttrValue (), który zwraca wartość na tej porządkowej. Więc mogę zrobić:

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'

Po krótkiej zabawie odkryłem, że możesz uzyskać wartość dowolnego parametru za pomocą potoku |jako separatora ścieżek:

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

Które mogą być przydatne do znalezienia geograficznego układu współrzędnych rzutowanego CS:

In [13]: srs.GetAttrValue("PROJCS|GEOGCS|AUTHORITY", 1)
Out[13]: '4277'
MerseyViking
źródło
1
Dzięki, to świetnie. Zaimplementuję to. Zabrakło mi czasu na dalszą zabawę - szybkie tworzenie aplikacji zostało wstrzymane przez brak dokumentacji GDAL / OGR!
Tomas
Próbowałem funkcji GetAttrValue z argumentami „AUTHORITY” i „1” i zauważyłem, że nie zawsze zwraca kod EPSG, ponieważ kod EPSG nie zawsze jest zawarty w WKT. Nadal jestem trochę niejasny, dlaczego tak jest. Znalazłem następujące rozwiązanie, które działa dobrze dla moich potrzeb: gis.stackexchange.com/questions/7608/…
Burrow
5

Oto fragment kodu, który zadziałał dla mnie:

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
źródło
0

SpatialReference.GetAuthorityCode()przyjmuje Nonejako parametr znajdujący węzeł uprawnień na elemencie głównym (tj. rzutowany / geograficzny, zależnie od przypadku). To samo dotyczy 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
źródło