GML-Export von MultiPoint Geometry (OGR mit PyWPS)

[fPol]

Hallo Leute,

ich bastle gerade an einem PyWPS-Prozess, der ein TIN (triangulated irregular network) als GML
übergeben bekommt, die Stützpunkte des TIN extrahiert und in ausgedünnter Form (als GML) zurückliefert.
Die Ausdünnung der Punkte klappt wunderbar und wenn ich das Ergebnis zB. als WKT (well known text) zurückliefere,
kann ich dies auch in einem GIS (zB. Quantum GIS) darstellen. Nun versuche ich aber die Ergebnispunkte als GML zurückzuliefern,
sodass Quantum GIS diese direkt als neuen Layer hinzufügen kann. Dabei erscheint folgende Fehlermeldung:


Ich erstelle erst die MultiPoint-Geometry, füge ihr alle Punkte hinzu und exportiere dann zu GML.

Code: Alles auswählen

#create new point geoemtry from the resulting points
        pointGeom = ogr.Geometry(ogr.wkbMultiPoint)
        newPoint = ogr.Geometry(ogr.wkbPoint)
        for point in resultPoints:
            newPoint.AddPoint(point[0], point[1], point[2])
            pointGeom.AddGeometry(newPoint)
  
        gml = pointGeom.ExportToGML()

Der Fehler besagt ja "filename too long", es sieht so aus, als versuche OGR aus dem ganzen GML-File den Namen zu erzeugen.
Hat jemand vielleicht eine Idee, welcher Zwischenschritt mir fehlt, oder wo der Fehler sonst liegen könnnte?

Anmerkung: Das Zurückgeben des GML erfolgt über einen ComplexOutput, eingelesen wird über eine ComplexInput.
Wenn ich das GML nur zum Server sende und unverändert zurückgebe, wird es ordentlich in Quantum GIS dargestellt.
Der Fehler muss also irgendwo beim Erzeugen des neuen GML liegen.
Den vollständigen Code könnt ihr hier finden: http://pastebin.com/4L2yfzwz

Ich freue mich auf eure Antworten!

[fPol]

So, ich denke ich hab herausgefunden was der Fehler ist:

Die Funktion ExportToGML liefert wirklich nur die GML-Repräsentation der Geometrie, also in diesem Fall die Punkte
der MultiPoint-Geometrie, aber kein fertiges GML-File. Den Header usw. muss ich jetzt selber zusammenbasteln und
dann das Ergebnis des Exports hineinschreiben.

[BlackJack]

@fPol: Wo sieht man denn das es sich um ein `ComplexOutput` handelt? Ich sehe im Quelltext nur `LiteralOutput` als Name. (Ich habe keine Ahnung von der API).

Kann es sein, dass an der Stelle wo Du das XML übergibst eigentlich ein Dateiname erwartet wird? So sieht das jedenfalls aus.

Der Quelltext enthält einige sehr unschöne Sachen. Ein paar Methoden sollten eigentlich Funktionen sein, weil sie mit der Klasse überhaupt nichts zu tun haben. Auf der anderen Seite gibt es da schon eine Klasse, dann sollte es keinen Grund geben ``global`` zu verwenden. Den sollte es ja schon bei Funktionen schon nicht geben weil man da besser eine Klasse schreiben sollte wenn Argumente und Rückgabewerte als Sprachmittel nicht ausreichen um ``global`` zu vermeiden.

Funktionen/Methoden können Rückgabewerte haben. Das wird im Quelltext nie wirklich verwendet — stattdessen werden leere Listen übergeben die dann gefüllt werden. Das ist verwirrend weil der Leser so nicht weiss ob da grundsätzlich leere Listen rein kommen, und wenn nein was da denn schon drin ist oder sein darf. Jeder Aufrufer muss eine Liste erstellen und übergeben und kann da auch Fehler machen, wenn er eine falsche/nicht leere Liste übergibt. Das liesse sich alles vermeiden.

``return`` ohne Wert am Ende einer Funktion oder Methode zu schreiben ist sinnfrei. Fehler sollte man über Ausnahmen „nach oben” Melden und nicht durch spezielle Rückgabewerte.

An einer Stelle wird ein Wert in eine Zeichenkettendarstellung umgewandelt die dann über Zeichenkettenmethoden in Einzelteile/-werte zerlegt wird. Das ist gruselig und fehleranfällig. Der Datentyp hat sicher eine API mit der man auf einem *vernünftigen* und dafür vorgesehenen Weg an die Werte heran kommt. Der Code sieht sogar so aus als wenn es sich ursprünglich um ein Tupel handelt.

Die Namen sind teilweise schlecht gewählt. Man sollte `layer` nicht mit `lyr` abkürzen. Was bringt das? Ausser Verwirrung beim Leser. Wenn man statt ``for feat in lyr:`` besser ``for feature in layer:`` schreibt, kann man sich den Kommentar sparen der dem Leser erklärt wofür `feat` und `lyr` stehen. Andere Beispiele sind Namen die nicht zum Wert passen. Wenn man `points` sieht, erwartet man irgendwie geartete Punkt-Objekte und keine Zahl. Die Anzahl der Punkte sollte `point_count` heissen. Auch hier könnte man dann den zusätzlichen erklärenden Kommentar weglassen.

[fPol]

Danke für deine Hinweise BlackJack! Ich gebe zu, der Code ist alles andere als schön und sauber; ich werde da noch einiges überarbeiten.
Bezügl. des ComplexOutput: hab grad gesehen, dass ich bei Pastebin eine alte Version hochgeladen hab, die noch den LiteralInput nutzt.
Hier die aktuelle Version:
http://pastebin.com/jUzLCdpG

[BlackJack]

@fPol: Wenn ein Argument grundsätzlich eine Zeichenkette mit dem Namen des Moduls enthalten soll, dann kann man diesen dynamisch ermitteln. Das ist sicherer als ihn per Hand in den Quelltext zu schreiben.

Bei `self.Output2` stimmt die Beschreibung nicht mit dem Wert überein der am Ende dort an das Programm übergeben wird!? Das wäre an der Stelle wohl auch deutlich aufgefallen wenn es nicht einen so generischen, nummerierten Namen hätte, sondern einen der beschreibt was dort übergeben wird. Denn eine Zeile wie ``self.timeOutput.setValue(lowDists)`` würde schon irgendwo Fragen aufwerfen die sich bei ``self.Output2.setValue(lowDists)`` nicht stellen. Da `lowDists` nicht mehr gefüllt wird, denke ich an der Stelle sollte tatsächlich die verbrauchte Zeit übergeben werden.

Die Liste `allPoints` in `execute()` wird für gar nichts verwendet.

Man sollte auch nicht alle möglichen Namen am Anfang einer Funktion an Werte binden die erst viel später im Verlauf der Funktion benötigt werden. Das macht die Zusammenhänge schwerer verständlich weil man unnötig zurückschauen muss wo etwas definiert wird, und es ist schwieriger einzelne voneinander unabhängige Teile nachträglich in eigene Funktionen auszulagern wenn eine Funktion zu umfangreich wird.

Bei `fillCells()` wird das `points11`-Argument überhaupt nicht verwendet. Statt die Namen mit Nummern zu versehen könnte man hier die Listen ihrerseits wieder in eine Liste oder vielleicht sogar eine verschachtelte Liste stecken. Bei den Tests werden zwei Teilausdrücke immer und immer wieder berechnet die noch nicht einmal von der Schleife abhängen. Ausserdem werden Punkte die genau auf den Mittel liegen nicht berücksichtigt, soll das so sein?

Wenn man den Umstand ausnutzt, dass der Datentyp `bool` von `int` erbt und `True` und `False` den Zahlenwerten 1 und 0 entspricht und man Wahrheitswerte deshalb auch als Index in eine Liste verwenden kann,, lässt sich `fillCells()` *deutlich* kompakter und übersichtlicher ausdrücken.

Das ganze könnte dann ungefähr so aussehen: http://pastebin.com/WddebtBA

[fPol]

@BlackJack
Wow, danke, dass du dir soviel Arbeit gemacht hast das umzuschreiben!

Ich werde mir das gleich mal in Ruhe angucken. Du hast Recht: Punkte die genau in der Mitte, auf den Kanten und nah an
Kantengrenzen liegen, werden nicht berücksichtigt, bzw. wird die Entfernungsberechnung nur in ihrer jeweiligen Zelle
durchgeführt, wodurch Punkte erhalten bleiben, die direkt neben einem Punkt aus der Nachbarzelle liegen.
Das wollte ich noch durch einen prozentualen Overlap der Zellen übereinander korrigieren.

[BlackJack]

@fPol: Warum wird das denn *überhaupt* auf die Zellen aufgeteilt? Am Ende wird ja doch wieder alles zusammen geworfen. Der hauptsächliche Effekt ist doch, dass Punkte in benachbarten Zellen erhalten bleiben, welche den Mindestabstand unterschreiten. Aber warum wird das gemacht?

Edit: Die Funktion `filter_close_points()` etwas effizienter und auch kürzer:

Code: Alles auswählen

def filter_close_points(minimal_distance, points):
    # 
    # TODO: This first line is only neccessary if the points have to be looked
    #       at in the order they are given to this function.  If going through
    #       them in reverse would be okay we could spare this step.
    # 
    points = list(reversed(points))
    while points:
        point = points.pop()
        yield point
        points = [p for p in points if distance(point, p) >= minimal_distance]

[fPol]

Die Aufteilung in Zellen wird gemacht, um die Anzahl an überflüssigen Distanzberechnungen zu reduzieren.

[BlackJack]

@fPol: Die letzte Funktion war leider fehlerhaft. Die hat zu viele Punkte durchgelassen. Hier die korrigierte Version:

Code: Alles auswählen

def filter_close_points(minimal_distance, points):
    # 
    # TODO: The reversing is only neccessary if the points have to be looked
    #       at in the order they are given to this function.  If going through
    #       them in reverse would be okay we could spare this step.
    # 
    points = list(reversed(points))
    while points:
        point = points.pop()
        point_count = len(points)
        points = [p for p in points if distance(point, p) >= minimal_distance]
        if len(points) == point_count:
            yield point

Zur Aufteilung in Zellen: So wie es jetzt gemacht wird produziert das falsche Ergebnisse weil Abstände zwischen den Punkten in verschiedenen Zellen nicht berücksichtigt werden. Auch 5% Überlappung der Zellen bringt nicht automatisch korrekte Ergebnisse weil die nötige Überlappung ja nicht fest ist, sondern von der minimalen Distanz abhängt. Und auch dann kann es immer noch passieren dass man Punkte aus dem Überlappungsbereich im Ergebnis hat, die trotzdem in einer anderen Zelle Punkte haben die zu Nahe dran liegen, so dass sie eigentlich gar nicht in das Ergebnis gehören.

Ich würde diese Zellengeschichte ja weg lassen und die Funktion von oben verwenden. Die hat im schlechtesten Fall schon mal die Hälfte weniger Distanzberechnungen. Im allgemeinen noch weniger weil in jedem Durchlauf nicht nur der aktuelle Punkt verworfen wird, wenn es zu nahe Punkte gibt, sondern diese Punkte auch gleich mit, denn die kommen offensichtlich ja auch nicht in Frage.

Sollte das nicht reichen würde ich nichts selber programmieren sondern eine entsprechende Bibliothek verwenden die für solche Fälle, also „räumliche Anfragen” effiziente Algorithmen und Datenstrukturen bietet.

Der Distanzfunktion fehlt übrigens ein Wurzelziehen wenn man den Benutzer die minimale Distanz nicht als Quadrat der eigentlich gewünschten Distanz eingeben lassen will. Soll das so?

Vereinfacht, ohne die Zellen wäre ich dann bei: http://pastebin.com/e8L9Ezrq