OGR in Python

[frabron]

Hallo,

nun will ich auch mal einen Schnippsel posten. Er ist wohl hauptsächlich für Geografen und Geoinformatiker interessant, die die GDAL Bibliothek benutzen und sich über die Python Bindings freuen. GDAL ist eine im GIS Bereich weit verbreitete Bibliothek zur Manipulation von Geodaten (Raster- & Vektordaten), die in C/C++ geschrieben ist und mittels SWIG eine Anbindung an Python hat. Ich brauche für ein Projekt ein Skript, dass mir Vektordaten von einem Format in ein anderes konvertiert, und dabei auch das räumliche Referenzsystem ändert (das Koordinatensystem sozusagen). Dazu benutzt man den OGR Teil von GDAL.

Wer sich einmal ein Python - GDAL Skript (wie z.B. ogr2ogr.py) angesehen hat, weiss, dass diese z.T. schwer zu durchschauen sind, denn sie sind bewusst am C Code angelehnt, damit man direkt vergleiche kann, wie die Mappings zwischen Python OGR und der darunterliegenden GDAL/OGR Bibliothek sind. Optisch ist das aber keine Augenweide

Deshalb jetzt mein Code, der zwar weniger Funktionsumfang hat, aber hoffentlich besser zu verstehen ist

http://paste.pocoo.org/show/336317/
zeigt den eigentlichen Code, wobei die convert_datasource Funktion der interessante Teil ist

http://paste.pocoo.org/show/336318/
ist ein Skript, dass ich zum Testen verwende und als Beispiel dient.

Vielleicht ist das ja hilfreich für den einen oder anderen. Über Verbesserungshinweise freue ich mich natürlich auch

Frank

[BlackJack]

@frabron: Ungeordnete Gedanken und Fragen:

Ist `check_resource()` nötig? Was passiert wenn man nicht existierende Ressourcen benutzt? Bekommt man dann eine Ausnahme von der zugrunde liegenden Bibliothek.

Kann man `is_available()` nicht auch durch `driver_by_name()` ausdrücken ohne über alle Treiber zu iterieren? Also einfach den Treiber holen und wenn es eine Ausnahme gibt, ist er halt nicht verfügbar. Die Liste `drivers` in der Funktion wird nicht verwendet.

`is_available()` teilt sich viel Quelltext mit `list_all()`. Ich tendiere oft dazu lieber einen Iterator zu schreiben -- wenn man eine Liste braucht, kann man ja einfach `list()` auf den Iterator anwenden. Dann könnte man `is_availabe()` auch auf ein `iter_drivers()` anwenden. Den Index würde ich da nicht mitliefern -- den kann man sich bei Bedarf auch mit `enumerate()` generieren. Ungetestet:

Code: Alles auswählen

def is_available(name):
    return any(d.GetName() == name for d in iter_drivers())

def iter_drivers():
    for index in xrange(ogr.GetDriverCount()):
        yield ogr.GetDriver(index)

Was gibt `layer.GetNextFeature()` nach dem letzten Feature zurück? Ist das etwas festes, dokumentiertes, wie zum Beispiel `None`? Dann kann man das nämlich in einer ``for``-Schleife ausdrücken:

Code: Alles auswählen

def iter_features(layer):
    for feature in iter(layer.GetNextFeature, None):
        yield feature
    layer.ResetReading()

Die ``for``-Schleife in `get_geometry_type()` finde ich etwas verwirrend. Das verschleiert ein wenig die Tatsache das man eigentlich nur das erste Feature haben möchte und es wird IMHO auch nicht wirklich deutlich genug, dass die Funktion `None` zurückgibt, sollte es keine Features geben. (Ist das Absicht?) Folgendes fände ich da etwas klarer:

Code: Alles auswählen

def get_geometry_type(layer):
    try:
        feature = iterate_features(layer).next()
    except StopIteration:
        return None
    else:
        return feature.GetGeometryRef().GetGeometryType()

[lunar]

Man könnte auch "next()" verwenden:

Code: Alles auswählen

def get_geometry_type(layer):
    feature = next(iterate_features(layer), None)
    if feature is not None:
         return feature.GetGeometriyRef().GetGeometryType()
    return None

[BlackJack]

@lunar: Wenn man ein Python verwendet das `next()` schon als Funktion hat. Ich bin da ja immer noch etwas konservativ.

[frabron]

Hallo Blackjack,

danke für deine Anregungen:

Ist `check_resource()` nötig? Was passiert wenn man nicht existierende Ressourcen benutzt? Bekommt man dann eine Ausnahme von der zugrunde liegenden Bibliothek.

Es passiert erst mal nix, driver.Open() gibt `None` zurück, wenn der Pfad falsch ist oder wenn sonst was schief läuft. Dennoch kann man sich fragen, ob `check_resource()` nötig ist, ist ja nur ein Einzeiler.

Kann man `is_available()` nicht auch durch `driver_by_name()` ausdrücken ohne über alle Treiber zu iterieren? ...

Schicke Idee! Auf die Idee mit `any()` wäre ich nie gekommen, habe ich auch so noch nie benutzt.

Was gibt `layer.GetNextFeature()` nach dem letzten Feature zurück? ...

Ich bin mir nicht sicher. Ich schau mal, was da zurückkommt. Leider gibt es keine Dokumentation zum Python-Binding, man muss sich mit der C oder C++ Dokumentation helfen.

Die ``for``-Schleife in `get_geometry_type()` finde ich etwas verwirrend. ...

Code: Alles auswählen

# ...
feature.GetGeometryRef().GetGeometryType()

Man muss wohl ein bisschen aufpassen mit den Referenzen auf Objekte, siehe hier:

Python crashes if you use an object after deleting an object it has a relationship with
...

Aber ich schau mir die Funktion mal an, denn du hast recht, dass die nicht besonders eindeutig ist was den Rückgabewert angeht. Ich mach mal die neue Version fertig und poste sie dann noch mal hier.

vielen Dank für deine Tipps bisher

[frabron]

Lunar: ich hab wohl wie Blackjack noch kein next()

[frabron]

Schade, der hier

Code: Alles auswählen

    for feature in iter(layer.GetNextFeature, None):
        yield feature
    layer.ResetReading()

geht leider nicht. Da geht wohl was verloren

for feature in iter(layer.GetNextFeature, None):
File "C:\Python26\lib\site-packages\osgeo\ogr.py", line 2110, in __cmp__
return self.Equal(other)
File "C:\Python26\lib\site-packages\osgeo\ogr.py", line 1505, in Equal
return _ogr.Feature_Equal(*args)
RuntimeError: Pointer 'hOtherFeat' is NULL in 'OGR_F_Equal'.

[frabron]

So, hier die neue Version. Es ist nicht mehr viel passiert, ich habe lediglich noch eine Funktion zum Erzeugen von SRS hinzugefügt (Zeile 46ff), da die Funktionalität mehrfach benötigt wird. Und natürlich, soweit möglich, die Verbesserungstipps von Blackjack eingebaut. Da der Code unter 2.5 auf meinem Server laufen muss, konnte ich Lunars Tipp leider nicht berücksichtigen, zumal bei iterate_feature() die Konstruktion mit while der einzig lauffähige Code war, den ich erzeugen konnte. Versuche mit iter() lieferten leider beim iterieren über die Features einen

RuntimeError: Pointer 'hOtherFeat' is NULL in 'OGR_F_Equal'

Kann ich wohl nix machen, das geht in die unbekannte Welt der Pointer ...