Artikel getaggt mit ‘SAGA GIS’

Januar 20th, 2012

Geoprocessing mit Open Source GIS Tools und Python

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

8 Leute mögen diesen Artikel.

Dezember 23rd, 2010

Uphill/Downhill-Berechnung mit SAGA GIS

Zielsetzung

In diesem Artikel möchte ich zeigen, wie man ausgehend von einem Weg und einem Digitalen Geländemodell jeweils die Bergseite und die Talseite des Weges mit Hilfe von GIS bestimmt. Dies bezeichnet man als Uphill/Downhill-Berechnung. In der Forstwirtschaft ist es beispielsweise interessant, ob das geschlagene Holz eines Bestandes bergauf oder bergab auf den Weg gerückt (transportiert) wird. Bergab wird die Reichweite der Rückung bei relativ gleichbleibenden Rückekosten und somit auch die Erschließungswirkung eines Weges mit ca. 300m angenommen. Bergauf jedoch kann man für vergleichbare Rückekosten nur bis maximal 150m kalkulieren. Daher ist es interessant, die Berg – und Talseite eines Weges zu bestimmen. Die Methoden wurden dem Dokument „Predicting wood skidding direction on Steep Terrain by DEM and Forest Road Network Extension“ entnommen und auf die Software  SAGA GIS in der Version 2.0.6 übertragen.

weiterlesen »

6 Leute mögen diesen Artikel.

März 25th, 2010

Regenwald in Google Earth (8) – Vektorisierung

Von Raster- zu Vektordaten

Bislang hatten wir es mit Rasterdaten zu tun. Die Klimadaten wurden entsprechend gefiltert, klassifiziert und verschnitten. Das Ergebnis könnten wir als Rasterdatensatz (bspw. GeoTIFF) exportieren und auf Google Earth als Bild-Overlay einfügen. Die Farbgebung ist dann allerdings festgelegt und auch mit der Transparenz von Bildern haperts nach meinen Erfahrungen noch bei Google Earth.

Deshalb werden wir unsere Ergebnisdaten vom Rasterformat ins Vektorformat konvertieren und anschließend ins KML-Format überführen.

weiterlesen »

1 andere Person mag diesen Artikel.

März 25th, 2010

Regenwald in Google Earth (7) – Analyse II


Globale Landbeckung

Jetzt können wir die Fläche mit den klimatischen Voraussetzungen für tropischen Regenwald mit der tatsächlichen Fläche des immergrünen Waldes abgleichen. Im ersten Teil dieses Artikels haben wir uns bereits Daten über die weltweite Landbedeckung des GLCF besorgt. Diese müssen nur noch in SAGA GIS importiert und weiterverarbeitet werden.

weiterlesen »

Be the first to like.

März 14th, 2010

Regenwald in Google Earth (6) – Analyse I

Verschneidung der Klimadaten

In diesem Schritt werden wir nun die Daten der weltweiten Jahresdurchschnittstemperatur und des weltweiten gemittelten Jahresniederschlags miteinander kombinieren, denn ein Parameter allein hat ja noch keine Aussagekraft bezüglich des potentiellen Wachstums tropischen Regenwaldes. Beispielsweise ist der Temperaturbereich in der Sahara prinzipiell im tropischen Bereich – die Feuchtigkeit fehlt hier jedoch. Analog dazu gibt es auch in Alaska genügend Niederschlag, nur ist es dort natürlich ein bisschen zu frisch für den tropischen Regenwald.

Also: Nur dort, wo sich beide Bereiche räumlich decken sind die klimatischen Voraussetzungen für tropischen Regenwald gegeben.

weiterlesen »

Be the first to like.

März 13th, 2010

Regenwald in Google Earth (5) – Klassifizierung

Klassifizierung von Temperatur und Niederschlag

Mit den Daten sind wir nun soweit fertig. Wir haben Temperatur und Niederschläge als Grids vorliegen. Was jetzt noch fehlt ist eine sinnvolle Klassifizierung der Daten. Die ursprüngliche Fragestellung war: Wo auf der Welt sind die klimatischen Voraussetzungen für tropischen Regenwald gegeben?

weiterlesen »

Be the first to like.

März 13th, 2010

Regenwald in Google Earth (4) – Niederschlag

Niederschlag

Die Jahresdurchschnittstemperatur haben wir bereits im vorhergehenden Teil berechnet. Nun soll es um den weltweiten kumulativen Jahresniederschlag gehen.

weiterlesen »

Be the first to like.

März 3rd, 2010

Regenwald in Google Earth (2) – SAGA GIS

SAGA GIS

SAGA GIS ist ein Open Source GIS mit dem Schwerpunkt in der Rasteranalyse. Mit diesem Werkzeug werden die Daten im Folgenden bearbeitet. Rasterdatenimport kann immer zunächst über die simple Variante mit Hilfe des GDAL-Treibers versucht werden:

weiterlesen »

Be the first to like.

März 2nd, 2010

Regenwald in Google Earth (1) – Daten

Klimadaten

Vor einiger Zeit wollte ich die Regenwaldflächen der Erde mit den klimatischen Voraussetzungen für tropisches Klima vergleichen. Also erstmal Treibstoff fürs GIS holen: Niederschlagsdaten können vom GPCC bezogen werden. Der Deutsche Wetterdienst (DWD) bietet auf seinen Seiten einen Visualizer/Downloaddienst vor allem für Niederschlagsdaten an. Hier kann für jedes Monat und Jahr der gewünschte Mittelwert (Temperatur / Niederschlag) in vielen Formaten ausgegeben werden. Wir lassen uns das Ganze im Format ESRI ASCII ausgeben.

weiterlesen »

1 andere Person mag diesen Artikel.