GIS, UNIGIS und andere Kleinigkeiten

Seit dem Update auf ArcGIS 9.1 stehen uns nun auch Lizenzen des Geostatistical Analyst zur Verfügung. Letzte Woche hat sich dann auch eine erste Anwendungsmöglichkeit ergeben: die Interpolation der Temperaturmessungen aus dem NAQUA-Messstellennetz. Die Ergebnisse der Interpolation waren nicht zur weiteren Verwendung als Inputdaten in einem Modell gedacht, sondern als Visualisierung. Deshalb -- und auch weil ich nicht zuviel Zeit investieren konnte -- habe ich ein verkürztes (wahrscheinlich unwissenschaftliches) Verfahren gewählt. Nach einiger Recherche zum Thema Interpolation mit Google (die interessantesten Links folgen am Ende), bin ich zum Schluss gekommen, dass es die 'richtige' Interpolationsmethode gar nicht gibt und dass es zu aufwendig wäre, jede der Methoden ausführlich zu konfigurieren und zu überprüfen. Daher habe ich einfach die drei Methoden, die am häufigsten mit positiven Resultaten erwähnt wurden (Kriging, Inverse Distance Weighting und Radial Basis Function), mit den Defaultwerten übernommen.

Das Interface des Geostatistical Analyst ist gut gelungen. Es gibt ein Menu mit diversen Explorationstools (Histogramm, Semivariogramm etc.), die auch in anderen Zusammenhängen sehr von Nutzen sein werden. Daneben gibt es dann den Wizard, mit dem die Hauptarbeit erledigt wird. Damit kann man bequem in mehreren Schritten die Interpolationsmethode und deren Parameter auswählen. Gleichzeitig kann man damit auch noch die Fehlerberechnung durchführen und sich so über die Güte des Modells informieren. Hat man alles eingestellt, wird die Interpolation berechnet und als 'Geostatistical Analyst Output Layer' dargestellt. Dies ist kein normaler Layer und hat den Vorteil, dass man via 'Method Properties' die gewählten Parameter jederzeit wieder abändern kann. Die geänderten Parameter werden vom Layer dann umgehend dargestellt. Ideal zum Rumpröbeln und Optimieren. Hat man dann einmal ein fertiges Ergebnis erzielt, kann man den Layer in ein GRID (Raster) oder eine Feature Class (Vektor) exportieren und weiterverwenden. Diese speziellen Layer werden übrigens auch angezeigt, wenn der Analyst gar nicht geladen ist. Nur eine Veränderung der Parameter ist dann nicht mehr möglich.

Die drei enstandenen Layer habe ich dann rein optisch und subjektiv beurteilt und den 'besten' (Kriging) ausgewählt. Das Ergebnis ist optisch ansprechend. Ob es wissenschaftlichen Kriterien standhält, ist natürlich eher zu bezweifeln. Dennoch: der Geostatistical Analyst ist ein ausserordentlich mächtiges und gut designtes Tool. Für einen korrekten und sinnvollen Gebrauch ist aber eine umfassende Einarbeitung sicher notwendig.

Einige Links zum Thema Interpolation:

• A Comparison of Spatial Interpolation Techniques in Temperature Estimation
• An Evaluation of Spatial Interpolation Methods on Air Temperature in Phoenix, AZ
• A Guide to Predicting Spatial Distribution of Weed Emergence Using Geographic Information Systems (GIS)
• Mapping of temperature variables in Israel: a comparison of different interpolation methods
• The generation of monthly gridded datasets for a range of climatic variables over the United Kingdom
• Interpolation Techniques for Climate Variables
• Spatial Interpolation of Monthly Mean Climate Data for China

Vorletzte Woche haben sich ein paar Interessierte aus dem Süddeutschen und Schweizer Raum zu einer UNIGIS-Lerngruppe zusammengefunden und sich in Basel ein erstes Mal getroffen. Der Begriff 'Lerngruppe' tönte nach einem recht ambitiösen Programm. Dem war aber zum Glück gar nicht so. Wir haben über unsere ersten UNIGIS-Wochen diskutiert und sind auch zusammen die Aufgaben durchgegangen und haben dabei noch ein paar unklare Punkte klären können, z.B. wie sind andere den Aufsatz über GIS und CAD angegangen? Und natürlich haben wir noch über jede Menge anderer Themen aus dem UNIGIS-Universum gesprochen (Summer School, AGIT etc.). Obwohl wir nur zu dritt waren, hat sich dieses Treffen -- für mich jedenfalls -- gelohnt. Zwar hätte ich die paar offenen Punkte auch über die Jahrgangs-Mailingliste abklären können, aber nach ein paar Wochen 'einsamer' Arbeit vor dem eigenen Monitor war ein reales Treffen genau die richtige Abwechslung. Ich denke, wir werden das sicherlich in Zukunft noch weiterführen (ca. 1x pro Modul).

Mit Version 9.1 von ArcGIS gibt es ja die Möglichkeit, Geoprocessing-Scripts mit Python zu schreiben. Dabei greift man mit Python (mit dem win32com-Modul) auf das Geoprocessing-Objekt zu, das -- wie alle ArcObjects -- ein COM-Objekt ist. Ich habe mir nun die Frage gestellt, ob man nun mit Python nun auch auf die anderen ArcObjects zugreifen kann. Mit Programmiersprachen wie VB, C++ funktioniert das ja problemlos. Die Antwort ist ein Nein mit einem kleinen Aber...

Mit Hilfe des win32com-Moduls kann Python ohne Probleme und sehr komfortabel auf COM-Objekte zugreifen. Der Haken an der Sache ist einfach, dass mit dieser Methode nur auf COM-Objekte zugegriffen werden kann, die das IDispatch-Interface implementieren. Ist dies nicht der Fall, klappt der Zugriff nicht. Der Zugriff über das IDispatch-Interface kann variabel gestaltet werden und entweder als 'early binding' oder als 'late binding' erfolgen. Tatsache ist, dass das Geoprocessing-Objekt von ArcObjects dieses IDispatch-Interface implementiert. Damit ist es möglich mit Python auf dieses Objekt zuzugreifen:

gp = win32com.client.Dispatch("esriGeoprocessing.GpDispatch.1")

Es gibt aber nur wenig andere ArcObjects, die dieses Interface ebenfalls implementieren. Eine Übersicht über diese Objekte kann man sich verschaffen, wenn man in den Objektmodellen ganz einfach nach 'IDispatch' sucht. Es sind dies nicht genug, um damit eine einigermassen komplexe Applikation entwickeln zu können. In Python gibt es neben dem IDispatch-Interface noch andere Zugriffsmöglichkeiten auf COM-Objekte sogenannte 'native Interfaces'. Python kann aber damit nicht auf beliebige Objekte zugreifen, sondern nur auf solche, für die eine bestimmte Erweiterung d.h. ein Modul existiert. Für ArcObjects war eine solche Erweiterung einmal in Entwicklung (PyArcObjects), aber sie ist im Moment auf Eis gelegt.

Kurz gesagt: der beliebige Zugriff auf ArcObjects via COM gelingt z.Zt. mit Programmiersprachen wie VBA, VB, C++, VB.Net, C# etc. Ein möglicher Umweg wäre via .NET. Es gibt für Python mindestens zwei Implementationen (IronPython, Python for .NET), die eine enge Verzahnung mit der .NET-Technologie anstreben. Da die ArcObjects ja auch via .NET angesprochen werden können, bestünde hier vielleicht eine Möglichkeit. Problematisch ist allerdings die Tatsache, dass in ArcGIS 9.1 nur Python 2.1 installiert ist.

Ohne Probleme möglich ist die Erzeugung eines COM-Servers in Python. Win32com bietet hier alle notwendigen Werkzeuge (inkl. Systemregistrierung) um eigene COM-Objekte in Python entwickeln zu können. Dies ist sicherlich eine valable Möglichkeit, um das Hindernis der fehlenden Entwicklungsumgebung auf unserem GIS-Server zu umgehen.

Weitere Dokumente zum Thema Python und COM:

• Kapitel 12 ("Advanced Python and COM") aus dem O'Reily-Buch "Python Programming on Win32"
• ArcGIS Desktop Developer Guide: http://edndoc.esri.com/arcobjects/9.1/ => ArcGIS Desktop => Developer Guide

Um das Python-Scripting zu ermöglichen, wird mit der ArcGIS-Installation gleichzeitig auch eine Python-Umgebung installiert. In ArcGIS 9.1 wird die Python-Version 2.1 installiert. Die hat zwar schon einige Jahre auf dem Buckel (aktuell ist im Moment v2.4.2), das sollte aber zu keinen Problemen führen. Teil dieser Python-Umgebung ist die Pythonwin-IDE, mit der man die Scripts editieren und debuggen kann. Natürlich ist eine solche IDE sehr nützlich, nur sind mir bei dieser IDE einige Umstimmigkeiten aufgefallen:

Pythonwin hat -- zumindest auf unserem GIS-Server -- Probleme mit Umlauten (und wahrscheinlich auch mit anderen Sonderzeichen). Es spielt keine Rolle, ob die Umlaute in Strings oder auch nur in Kommentaren auftauchen, Pythonwin verschluckt sich daran. Das Programm stürzt zwar nicht ab, aber es funktioniert auch nicht mehr korrekt, so dass ein Arbeiten damit sehr erschwert ist.

Ab und zu hatte ich den Eindruck, dass Pythonwin sehr träge reagiert, v.a. beim Scrollen im Vollbildmodus. Hier ist aber nicht gesagt, dass das ein Pythonwin-Problem ist, vielleicht liegt das ja auch an unserer GIS-Umgebung.

Was mich an Pythonwin auch sehr enttäuscht, ist die Command Completion. Es existiert zwar eine solche, die zeigt einem bei der Eingabe aber nur diejenigen Objekte und Funktionen an, die im Script stehen. So werden mir zum Beispiel jeweils nicht alle möglichen Methoden oder Properties des Geoprocessing-Objekts angezeigt, sondern nur diejenigen, die im Script stehen. Zudem habe ich nicht herausgefunden, wie die Command Completion mit der Tastatur bedient werden kann. Weder die Enter- noch die Space- noch die KeyRight-Taste haben eine Wirkung darauf. Es funktioniert nur, wenn man mit der Maus auf den Vorschlag der Command Completion klickt. Erst dann wird dieser übernommen. [Update: Es ist natürlich die Tabulator-Taste, die den gewählten Befehl vervollständigt. Ich war wohl zu sehr an den VBA-Editor gewöhnt!] Kein Vergleich zu der um Klassen besseren Command Completion in der VBA-Umgebung. All diese Einschränkungen haben mich dazu gebracht, meine Scripts vorerst nicht mit Pythonwin zu schreiben, sondern mit einem simplen Texteditor (SciTE), den ich auch sonst häufig benutze. Der hat zwar nur Syntax Highlighting und keine Command Completion, dafür reagiert er sehr schnell und zuverlässig.

In dem vom Model Builder exportierten Python-Script taucht die folgende Anweisung auf:

gp.SetProduct("ArcInfo")

Auf den ersten Blick scheint diese Anweisung klar zu sein. Da im Script das Tool zum Importieren von Interchange-Files gebraucht wird (Import from Interchange File) und dieses nur mit einer ArcInfo-Lizenz benutzt werden kann, ging ich davon aus, dass die SetProduct-Anweisung eine solche ArcInfo-Lizenz -- sofern verfügbar -- anfordert. Der entsprechende Eintrag in der ArcGIS Desktop-Hilfe bestätigte diese Ansicht. Als ich aber das Script unter ArcView laufen liess, stürzte es trotzdem wegen der fehlenden ArcInfo-Lizenz ab. Eine ausführliche Suche in den ESRI-Foren ergab dann, dass in Tat und Wahrheit diese Funktion die Lizenzart nicht verändert. Sie teilt dem Python-Script einzig und allein mit, welche Lizenz gerade verwendet wird. Im Prinzip ist diese Funktion also nutzlos. ESRI empfiehlt aber, das benötigte Lizenz-Level doch mittels SetProduct dem Script mitzuteilen. Folgende Threads in den ESRI-Foren befassen sich auch mit diesem Thema:
Failure of SetProduct to get elevated permissions
Debugging Python scripts
Eine Möglichkeit, dieses Problem zu umgehen, wäre eventuell Folgendes: In der ArcToolbox werden diejenigen Scripts, die beim jeweils eingestellten Lizenz-Level nicht zur Verfügung stehen, entweder gar nicht oder dann besonders gekennzeichnet dargestellt. Ideal wäre es nun, wenn man dies auch für eigene Scripts implementieren könnte. Dann käme auch niemand auf die Idee, das Script mit der falschen Lizenz auszuführen. Bis jetzt habe ich aber noch nicht rausgefunden, ob und wie das machbar ist.

Die Geographic Information Technology Training Alliance (GITTA) -- ein Projekt des Swiss Virtual Campus (SVC) -- hat die von ihr entwickelten Lernmodule unter einer Creative Commons License veröffentlicht. Die bisher verfügbaren Module decken die Themen Einführung in GIS, Datenbanken, Datenerfassung, räumliches Modellieren, räumliche Analyse, Datenpräsentation und Fallstudien ab. Einige Module sind in deutsch geschrieben, einige in englisch, ein paar wenige sind in beiden Sprachen verfasst. Einige der Modultitel erinnern mich an Modultitel des UNIGIS-Studiums. Ich werde mir daher die GITTA-Module bei Gelegenheit mal näher anschauen. Da ich das erste UNIGIS-Modul bald beendet haben werde, kann ich da sicher schon einen ersten Vergleich anstellen.

Ein Wort zur Performance der Python-Scripts. Um hier zu wirklich aussagekräftigen Resultaten zu kommen, fehlte mir die Zeit. Zudem hätte ich mein Script noch als VBA-Programm implementieren sollen, um einen Vergleich zu machen. Die Performance ist sicherlich einigermassen zufriedenstellend. Es macht keinen Unterschied, ob das Script in ArcMap oder in ArcCatalog ausgeführt wird. Beide Male beträgt die Laufzeit ca. zwei Minuten (importiert wird u.a. immerhin ein 30 MB grosses e00-File). Deutlich langsamer läuft das Script ab, wenn ich es in der nackten DOS-Kommandozeile ausführe. Dann beträgt die Laufzeit ca. fünf Minuten. Dieser Unterschied ist zwar markant aber nicht weiter von Bedeutung, da ich Scripts meistens sowieso nur in ArcMap oder ArcCatalog ausführen werde. Insgesamt werde ich aber trotzdem das Gefühl nicht los, dass die Performance doch noch ein bisschen besser sein könnte. Darauf weist auch der Technical Article Nr. 30420 von ESRI hin.

Ganz im Gegensatz zu den USA sind in Europa die Geodaten der nationalen Landesvermessungbehörden nicht frei erhätlich sondern müssen zu einem z.T. erheblichen Preis erworben werden. Das ist meiner Meinung nach unverständlich, da diese Daten ja bereits mit Hilfe von Steuergeldern erhoben wurden. Ich sehe nicht ein, weshalb ich dafür noch einmal bezahlen muss. In den USA ist die freie Verfügbarkeit der nationalen Geodaten eine Selbstverständlichkeit. Zwar verliert der Staat vorübergehend Einnahmen aus den Lizenzgebühren für die Daten, dafür sind aber diese Geodaten die Grundlage für eine ausserordentlich breite Palette an neuen Nutzungsmöglichkeiten nicht nur im privaten Bereich sondern vor allem auch in wirtschaftlicher Hinsicht, wovon der Staat letztendlich wieder profitieren würde. Aus diesen und anderen Gründen läuft zur Zeit eine Online-Petition gegen die INSPIRE-Initiative der EU. Die INSPIRE-Initiative will im Prinzip eine bessere Vernetzung und Zusammenarbeit zwischen den nationalen Geoinformations-Behörden fördern, was sehr zu begrüssen ist. Allerdings ist der letzte Entwurf dazu so verwässert worden, dass er die freie Verfügbarkeit von Geodaten eher verhindert als fördert. Obwohl die Schweiz der EU nicht angehört, habe ich die Online-Petition unterzeichnet, da ich deren Forderungen zustimme und sich die Schweiz einer solchen EU-Initiative ja sowieso nicht verschliessen könnte. Jo Walsh - eine der InitiantInnen der Petition - erläutert ihre Argumente in einem lesenswerten Artikel auf Directionsmag. Alle Infos zur Petition sind auf petition.publicgeodata.org zu erfahren. Also: wer für die freie Verfügbarkeit von Geodaten ist, sollte die Petition unbedingt unterzeichnen!

Neben den im ersten Teil erwähnten Problemen, sind mir noch weitere interessante Punkte im Bezug auf Python-Scripting in ArcGIS aufgefallen. Bei ArcGIS musste ja man schon seit jeher aufpassen, dass man sich mit dem einen Programm nicht den Workspace sperrt ("lock"), so dass man mit einem anderen dort nichts mehr machen kann. Häufig tritt das bei gleichzeitigem Arbeiten mit ArcMap und ArcCatalog auf. Dieses Problem besteht natürlich auch, wenn man mit Scripts arbeitet. Eine häufige Fehlerquelle ist z.B. wenn man das Script in der Pythonwin-IDE schreibt und dort auch das Debugging macht und gleichzeitig noch ArcCatalog oder ArcMap geöffnet hat. Auch hier gilt demnach: wenn möglich immer nur eine Applikation zurselben Zeit geöffnet haben. Die von der Scripting-Umgebung in diesem Konfliktfall zurückgegebenen Fehlermeldungen lassen zudem auch nicht immer auf die tatsächliche Problematik schliessen.

Im einem früheren Post habe ich meine ersten Erfahrungen mit dem neuen Model Builder von ESRI beschrieben. Das dort erstellte Modell habe ich in ein Python-Script exportiert und mir das mal genauer angeschaut.

Die Struktur des Scripts ist sehr simpel, da ja auch mein zugrundeliegendes Modell nicht besonders kompliziert ist. Zuerst werden die benötigten Python-Module importiert:

import sys, string, os, win32com.client

dann das Geoprocessing-Objekt initialisiert:

gp = win32com.client.Dispatch("esriGeoprocessing.GpDispatch.1")

dann die Lizenz-Art und der Workspace gesetzt:

gp.SetProduct("ArcInfo")

dann die notwendigen Toolboxes importiert:

gp.AddToolbox("C:/Program Files/ArcGIS/ArcToolbox/Toolboxes/Conversion Tools.tbx")
gp.AddToolbox("C:/Program Files/ArcGIS/ArcToolbox/Toolboxes/Data Management Tools.tbx")
gp.AddToolbox("C:/ArcGIS/arcexe9x/Toolboxes/Coverage Tools.tbx")

dann die Eingabe-Parameter übernommen:

Segments_e00-File = sys.argv[1]
Catchments_e00-File = sys.argv[2]
BIL-File = sys.argv[3]
PGDB-Ordner = sys.argv[4]
PGDB-Name = sys.argv[5]

und danach die einzelnen Geoprocessing-Schritte ausgeführt. Soweit so gut.
Allerdings ist mir aufgefallen, dass der Export aus dem Modell in das Python-Script nicht ganz befriedigend ist. Im Script fand ich ziemlich viele Variablen, die absolute Pfade enthielten:

cat5237_polygon = "H:\\Data\\mars\\ccm\\catchments\\polygon"
BIL_projiziert = "H:\\Data\\mars\\ccm\\landcover"

Der Sinn solcher Scripts soll es ja sein, dass man sie immer wieder in verschiedensten Projekten benutzen kann und da sind fest eingebaute ("hard coded") absolute Pfade nicht erwünscht. Ich habe anschliessend versucht, mit Änderungen im Modell diese Variablen zum verschwinden zu bringen, was allerdings nicht möglich war.

Im Model Builder gibt es das Tool "Select Data". Wenn bei einem Tool der Output ein Coverage oder ein Ordner ist, dann muss im Model Builder das "Select Data"-Tool benutzt werden, um mit den im Coverage enthaltenen Feature Classes (arc, region, route) weiterarbeiten zu können. In der ArcGIS Desktop-Hilfe wird ausdrücklich erwähnt, dass dieses Tool nur im Model Builder verwendet werden soll und sonst nicht. Im exportierten Python-Script taucht allerdings dieses "Select Data"-Tool trotzdem auf. Es wird zwar ohne Fehler ausgeführt, ist aber vollkommen nutzlos. In den Python-Scripts kann man direkt auf Coverage Feature Classes zugreifen. Hat es z.B. im Coverage "Catchment" eine Region mit dem Namen seaoutlet, kann man auf diese ganz einfach zugreifen mit:

Catchment\region.seaoutlet

Dasselbe gilt natürlich auch für Geodatabases, Ordner etc.

In Zukunft werde ich Geoprocessing-Tasks wahrscheinlich ohne den Model Builder direkt als Python-Script implementieren. Die Möglichkeiten und die Kontrolle über den Code sind hier doch deutlich grösser als im Model Builder. Daher habe ich das exportierte Script auch manuell umgeschrieben, so dass die obigen Problemstellen nicht mehr enthalten sind. Das überarbeitete Script findet sich hier.