Aufgabenstellung
Momentan soll ich eine ganze Menge an Geodaten für einen wissenschaftlichen Auftraggeber verarbeiten. Dabei geht es um die Potentialabschätzung der Wuchsleistung von schnellwachsenden Baumarten. Ich freute mich riesig auf den professionellen Einsatz von Open Source GIS Tools unter etwas anspruchsvolleren Bedingungen – denn in letzter Zeit hatte ich beruflich wenig mit dem OS GIS Universum zu tun.
Große Datensätze
Anspruchsvoll? Oh ja – es handelt sich um eine bayernweite Analyse und da kommt schnell einiges zusammen…
Hier eine kleine Übersicht:
- 4 Sachdatenbanken (MS Access)
- 9 vektorielle Geodatensätze in insgesamt über 750 einzelnen ESRI Shapefiles (ca. 6,8 GB)
- 3 Rasterdatensätze als ESRI ASCII bzw. Binary GRID (ca. 4,6 GB)
Nicht übel. Diese Menge von rund 11 GB sollte also mit Open Source GIS Tools verarbeitet werden. Das war im Auftrag natürlich nicht gefordert, da dieser softwareunabhängig formuliert war – aber für mich war es klar, dass ich FOSSGIS-Tools einsetze.
Ehrlich gesagt habe ich bislang mit FOSSGIS-Werkzeugen keine aufwendigen GIS-Analysen durchgeführt, sondern eher ausprobiert, was damit alles geht. Nun also eine vollständige Analyse. Dass das spannend wird, dachte ich mir schon…
Welche Werkzeuge?
Nun gut – welches Werkzeug solls denn nun sein? Da die Daten – vor allem die Shapefiles – ziemlich zerhackstückelt über ganz Bayern daher kamen (Schönen Gruß ans LfU und das LVG), mussten diese erstmal zu bayernweiten Datensätzen zusammengeführt werden. Das war vor Angebotsabgabe leider nicht bekannt – aber gut, ich nahm die Herausforderung an.
QGIS – fTools
Bei QGIS gibts eine Erweiterung, welche auch Shapefiles zusammenführen kann. Die sog. fTools können dies und auch noch viel mehr (z.B. Intersection, Dissolve, Multipart to single part etc.). Dabei basiert fTools auf QGIS und GDAL/OGR, um die Geoprocessing-Funktionen anbieten zu können. Außerdem ist es in Python geschrieben.
Die Variante der fTools in QGIS kam aufgrund der doch langen Laufzeiten nicht in Frage. Die GUI von QGIS frisst da anscheinend doch einiges an Performance — außerdem muss fTools dann als Plugin von QGIS mit dem Mutterprozess hin und her kommunizieren. Das wollt ich mir ersparen.
Man muss aber sagen, dass die fTools für kleinere Batch-Arbeiten gut geignet sind. Das Pythonskript, das ich mit der Vorlage der fTools geschrieben habe liest sich einfach. Das QGIS-Objektmodell folgt dem üblichen „Feature – Geometry – Attributtable“-Schema und ist somit leicht nachzuvollziehen.
Bis die QGIS-Pythonbindings dann aber mal ansprechbar sind kann ne Zeit vergehen. Gary Sherman hat dazu eine ausführliche Anleitung hier. Dazu mal an anderer Stelle mehr.
GDAL/OGR – kennste eine – kennste alle?
GDAL/OGR ist ja das ETL (Extract-Transform-Load) Werkzeug schlechthin. Sagt man. Ist auch so sag ich jetzt, aber noch vor ein paar Wochen hatte ich noch nicht viel Ahnung, wie man die Bibliothek einsetzt.
Ich hatte GDAL/OGR schon immer unbedacht benutzt, wenn ich z.B. mit QGIS oder GRASS arbeitete, da viele Desktop-GIS diese Bibliothek unter der Haube haben – allerdings eben immer „nur“ als Datenprovider für Desktop GIS.
Nun war klar, dass ich die Python-Bindings von GDAL/OGR genauer anschauen muss. Also erstmal recherchieren und viel ausprobieren (die Python-Doku für GDAL/OGR ist derzeit noch sehr spärlich) und irgendwann klappts dann aber auch. Das grundsätzliche Modell bei Geoprocessing-Bibliotheken ist ja immer relativ ähnlich: es gibt ein Feature-Objekt mit Sachdaten und ner Geometrie und man hat ein paar GIS-Methoden (welche im FOSSGIS-Universum ja eigentlich fast immer auf GEOS und somit der JTS basieren) zur Verfügung (Intersection, Buffer, etc.). Wie soll es auch mordsmäßig verschieden sein?
Genau so ists auch bei den Python-Bindings von GDAL/OGR. Die Extract und Load-Komponente kommt von GDAL/OGR selbst. Die geometrischen Verschneidungsoperatoren werden bei GDAL/OGR mit GEOS bewerkstelligt, während Transformationen mit PROJ4 behandelt werden.
Wer – wie ich – eher aus der ESRI-Ecke mit dem ArcGIS Geoprocessor kommt, wird sich bei GDAL/OGR in Python auch recht schnell zurecht finden. Die FDO API in C# war da schon ne andere (weil abstraktere) Nuss.
Interessante Links für GDAL/OGR mit Python:
Letztlich konnten die Daten dann mit einem Pythonskript und der GDAL/OGR-Bibliothek zu bayernweiten Datensätzen zusammengeführt werden. Und wie schnell das geht! Da ist GDAL/OGR wirklich in ihrem Element. Wahnsinn, was die Bibliothek wegschreibt:
- rund 100 Shapefiles mit insgesamt 2,1 Mio Features in 20 Minuten (800MB)
- rund 100 Shapefiles mit insgesamt 3,1 Mio Features in 35 Minuten (4,2 GB)
Diese Zahlen überzeugen.
PostGIS
Auch PostGIS fällt nicht aus der Reihe der FOSSGIS-Tools und basiert auf GEOS und der PROJ4-Bibliothek. Die Kombination von SQL als Datenabfragesprache mit geometrischen Funktionen (Spatial SQL) fand ich immer spannend. Nachdem die Shapefiles also zusammengeführt waren, wollte ich einfach alle Vektordaten in eine PostGIS-DB importieren und die Aufbereitung in PostGIS machen. Auch dass mit PL/Python meine Lieblingssprache Python in einer Datenbank als prozedurale Sprache zur Verfügung steht find ich klasse (ich denke hier an die Lesbarkeit von PL/SQL Code vs. Python…).
Allerdings kam ich nicht weit. PostGIS hat definitiv als professionelle Geodatenhaltung seine Stärken und einige Anfragen an Daten kann es auch schnell beantworten, aber PostGIS ist nicht „post GIS“ (im Sinne von „nach GIS“) wie es im Vorwort des Buchs „PostGIS in Action“ hies. Man braucht bei solchen Datenmengen anscheinend schon noch ein extra Analysewerkzeug. Vielleicht liegts auch an meiner Rechenpower oder meinen PostGIS/SQL-Kenntnissen, aber GIS-Analysen mit dieser Datenmenge klappt nach meiner Erfahrung nicht in PostGIS. Nach Laufzeiten von 2 Tagen habe ich die aufwendigen Verschneidungen (bayernweit) sein gelassen. Schade eigentlich. Ich hätte irgendwie mehr Performance von PostGIS erwartet.
SAGA GIS
Raster? Nimm SAGA GIS. Eben – für Rasterdaten ist SAGA GIS spitze. Performant, gut bedienbar, sogar per Python oder Kommandozeile ansprechbar für Batchabläufe. Bei Vektordaten hat SAGA aber so seine Macken. Außerdem werden nicht viele Möglichkeiten in der Vektoranalyse angeboten. Diejenigen, die dabei sind basieren interessanterweise jedoch nicht auf GEOS, sondern auf CGAL.
Aber gut, SAGA GIS ist eben vor allem ein Raster-GIS. In Ordnung.
GRASS GIS – hat so nen Bart?
Allerdings. GRASS GIS hat schon einige Jahre auf dem Buckel (30+?) – wird aber laufend fortentwickelt. Seit geraumer Zeit sieht auch die GUI etwas ansprechender aus. Von den Analysemöglichkeiten her ist GRASS GIS im FOSSGIS-Universum meiner Meinung nach einsame spitze.
Nimm GRASS GIS – wenn Du nicht ewig Zeit hast. Im Zusammenhang mit Scripting ist das Werkzeug ein verdammt performantes Analyse-Tool. Natürlich sollte man an dieser Stelle die etwas höhere Lernkurve bei GRASS ggü. anderen Werkzeugen berücksichtigen. Aber die Zeitinvestition lohnt sich!
Ich bin in den o.g. Ausführungen immer auf die Geometriebibliothek der Tools eingegangen. Da fast alle OS GIS Werkzeuge letztlich die geometrischen Funktionen an die GEOS-Bibliothek weiterreichen, ist auch die Performance bei großen Datenmengen vergleichbar.
Damit ist es relativ egal, ob man GDAL/OGR oder die fTools von QGIS oder PostGIS verwendet. Alle basieren auf GEOS. GEOS ist ja ein C++ Port der JTS (Java Topology Suite). Und diese ist wiederum in anderen FOSSGIS-Tools verbaut wie gvSIG, OpenJUMP, uDIG oder Kosmo. GIS-Analysen mit derartigen Datenmengen in Java-Tools durchzuführen hatte bei mir leider regelmäßig einen „Out of heap space“ Fehler zur Folge. Deshalb schied auch die JTS aus.
Die GRASS GIS Overlay-Operatoren sind jedoch als C-Algorithmen anscheinend richtig speicheroptimiert geschrieben. Trotzdem bleibt ein kleiner Wermutstropfen: Auch in GRASS GIS konnte ich die Datenmengen nicht bayernweit importieren bzw. verschneiden. Mit einem 64bit System und mehr RAM wäre das vielleicht schon gegangen. So aber blieb mir nur ein „Tiled Processing“, in dem ich die bayernweiten Datensätze in 16 Kacheln aufgeteilt habe. Mit GRASS GIS hats dann im Batch-Modus per Pythonbindings allerdings recht flott (rund 1,5 Std. für einen bayernweiten Verschneidungsdurchlauf) geklappt. Vielleicht mach ich demnächst auch nochmal nen Vergleich mit einem Tiled Processing der Daten in PostGIS, aber ich habe wenig Hoffnung, dass dies ähnlich performant ist wie GRASS GIS.
Und ne kleine GRASS Python Utility Klasse ist auch noch herausgekommen:
# -*- coding: latin-1 -*-
# Helper for GRASS GIS Python Scripting
# Johannes Sommer, 18.01.2012
class utility():
import core as grass
''' Utility-Class for GRASS GIS Python Scripting.
Copyright by Johannes Sommer 2012.
- doAddColumn
- doCalcColumnValue
- doRenameColumn
- doIntersection
- doDifference
- getMapsets
- getVectorLayers
- getRasterLayers
- getVectorLayersFromMapset
- getRasterLayersFromMapset
- getColumnList
- doFilterLayerList
Usage:
import grass.utils as utils
gUtil = util.utility()
gUtil.dodifference('forest', lakes')
'''
def doAddColumn(self, _lyrA, _colName, _colType):
''' Adds a column to the input vector layer with the given type.
GRASS datatypes depend on underlying database, but all
support:
- varchar (length)
- double precision
- integer
- date
Usage:
doAddColumn('lakes', 'width', 'double precision') '''
retVal = grass.run_command('v.db.addcol', map=_lyrA,
columns='%s %s' % (_colName, _colType) )
return retVal
def doCalcColumnValue(self, _lyrA, _colName, _newVal):
''' Updates a column of the input vector layer with the passed value.
@_newVal can be a certain value or a calculation of existing columns.
Usage:
doCalcColumnValue('lakes', 'width', 'length/2') '''
retVal = grass.run_command('v.db.update', map=_lyrA,
column=_colName, value=_newVal)
return retVal
def doRenameColumn(self, _lyrA, _colNameA, _colNameB):
''' Renames a column of the input vector layer.
Usage:
doRenameColumn('lakes', 'width', 'width_old') '''
retVal = grass.run_command('v.db.renamecol', map=_lyrA,
column='%s,%s' %(_colNameA,_colNameB))
return retVal
def doIntersection(self, _lyrA, _lyrB, OVERWRITE_OUTPUT==False):
''' Performs a geometric intersection ('AND'-Operator) of two vector layers
and writes the output to a layer called 'layerNameA_X_layerNameB'.
@OVERWRITE_OUTPUT can be set to 'Y' or 'N'
Usage:
doIntersection('lakes', 'forest') '''
outputLyr = _lyrA.split('@')[0] + '_X_' + _lyrB.split('@')[0]
optArg = ''
if OVERWRITE_OUTPUT == True:
optArg = '--overwrite'
retVal = grass.run_command('v.overlay %s' % optArg, ainput=_lyrA,
binput=_lyrB, output=outputLyr,
operator='and')
return retVal
def doDifference(self, _lyrA, _lyrB, OVERWRITE_OUTPUT==False):
''' Performs a geometric difference ('NOT'-Operator) of two vector layers
and writes the output to a layer called 'layerNameA_DIFF_layerNameB'.
@OVERWRITE_OUTPUT can be set to 'Y' or 'N'
Usage:
doDifference('lakes', 'forest') '''
outputLyr = _lyrA.split('@')[0] + '_DIFF_' + _lyrB.split('@')[0]
optArg = ''
if OVERWRITE_OUTPUT == True':
optArg = '--overwrite'
retVal = grass.run_command('v.overlay %s' % optArg, ainput=_lyrA,
binput=_lyrB, output=outputLyr,
operator='not')
return retVal
def getMapsets(self):
''' Returns all GRASS mapsets in current GRASS location as list.'''
return grass.mapsets()
def getVectorLayers(self):
''' Returns all vector layers in current GRASS location as list.'''
return grass.list_strings('vect')
def getRasterLayers(self):
''' Returns all raster layers in current GRASS location as list.'''
return grass.list_strings('rast')
def getVectorLayersFromMapset(self, _mapset):
''' Returns all vector layers in given mapset as list.'''
vLyrList = getVectorLayers()
# Filter List by MAPSET
vLyrListFiltered = []
for item in vLyrList:
if _mapset in item:
vLyrListFiltered.append(item)
return vLyrListFiltered
def getRasterLayersFromMapset(self, _mapset):
''' Returns all raster layers in given mapset as list.'''
rLyrList = getRasterLayers()
# Filter List by MAPSET
rLyrListFiltered = []
for item in rLyrList:
if _mapset in item:
rLyrListFiltered.append(item)
return rLyrListFiltered
def getColumnList(self, _vLyr):
''' Returns all column names of a vector layer as list.'''
# Import vector wrapper from %GISBASE%\etc\python\script
import vector as vector
out_colList = []
out_colList = vector.vector_columns(_vLyr, getDict = False)
return out_colList
def doFilterLayerList(self, _lyrList, _lyrkey):
''' Filters input layer list per _lyrkey and returns filtered list.'''
out_lyrList = []
for item in _lyrList:
if _lyrkey in item:
out_lyrList.append(item)
return out_lyrList