Apakah mungkin untuk mendapatkan nilai EPSG dari kelas OSR SpatialReference menggunakan OGR Python API?

21

Saat membaca sebuah layer dari koneksi OGR PostGIS saya bisa mendapatkan SpatialReference dari layer, tetapi apakah mungkin untuk mendapatkan nilai EPSG? Apakah ada dokumentasi tentang ini?

Sebagai contoh:

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

Pengembalian:

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

Jadi, bagaimana cara mendapatkan nilai EPSG untuk proyeksi? Misalnya:

srs.GetEPSG()
print srs
27700

Saya sudah mencoba srs.GetAttrValue('AUTHORITY'), tetapi ini baru saja kembali 'EPSG'.

python postgis ogr Tomas
sumber
I've tried srs.GetAttrValue('AUTHORITY'), but this just returns 'EPSG'yang mana yang benar. EPSG adalah otoritas
nmtoken

Jawaban:

30

Agak terkubur, tetapi ada parameter kedua untuk GetAttrValue () yang mengembalikan nilai pada ordinal itu. Jadi saya bisa melakukan:

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'

Setelah sedikit bermain, saya menemukan Anda bisa mendapatkan nilai untuk parameter apa pun menggunakan pipa |sebagai pemisah jalur:

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

Yang mungkin berguna dalam menemukan sistem koordinat geografis dari CS yang diproyeksikan:

In [13]: srs.GetAttrValue("PROJCS|GEOGCS|AUTHORITY", 1)
Out[13]: '4277'
MerseyViking
sumber
1
Terima kasih, itu hebat. Saya akan menerapkannya. Saya kehabisan waktu untuk 'bermain-main' lebih lanjut - pengembangan aplikasi yang cepat terhambat oleh kurangnya dokumentasi GDAL / OGR sekali lagi!
Tomas
Saya mencoba fungsi GetAttrValue dengan argumen "AUTHORITY" dan "1" dan memperhatikan bahwa itu tidak selalu mengembalikan kode EPSG karena kode EPSG tidak selalu termasuk dalam WKT. Saya masih sedikit bingung tentang mengapa ini terjadi. Saya menemukan solusi berikut ini berfungsi dengan baik untuk kebutuhan saya: gis.stackexchange.com/questions/7608/…
Burrow
5

Berikut cuplikan kode yang berhasil bagi saya:

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
sumber
0

SpatialReference.GetAuthorityCode()mengambil Nonesebagai parameter, yang menemukan simpul otoritas pada elemen root (mis. diproyeksikan / geografis jika diperlukan). Hal yang sama berlaku untuk 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
sumber