Artikel getaggt mit ‘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

9 Leute mögen diesen Artikel.

Mai 29th, 2011

Feature Data Object API und DotSpatial in .NET (3)

Die Agenda beim Thema FDO in .NET sah folgende Punkte vor:

  1. Get the list of installed providers
  2. Create a connection manager
  3. Create a connection
  4. Set the connection properties
  5. Open a connection
  6. Get the connection state
  7. Fetch Data

Die Punkte 1-6 wurden bereits im letzten Artikel abgehandelt. Wir werden uns nun dem 7. Punkt widmen. Und der hat es in sich. Denn wenn man die FDO API zum Laden von Daten verwendet, bekommt man FDO Objekte. Ziel soll es aber sein, die FDO Objekte in einem DotSpatial Viewer darzustellen.

Dazu muss man das Ergebnis der FDO API zu einem DotSpatial FeatureSet zusammenbauen.

7. Fetch Data

Zunächst müssen wir ein Command erzeugen, das wir später dann ausführen werden. (Eine Liste der verfügbaren Commands ist im Namespace OSGeo.FDO.Commands.CommandType zu finden. Außerdem sind die Commands auch immer abhängig vom Datenprovider. Die Dokumentation der FDO hilft hier weiter.)

if (connState == OSGeo.FDO.Connections.ConnectionState.ConnectionState_Open)
{
	ISelect sel = (ISelect)conn.CreateCommand(OSGeo.FDO.Commands.CommandType.CommandType_Select);
    sel.SetFeatureClassName("world_0");
	IFeatureReader reader = sel.Execute();
	FgfGeometryFactory geoFac = new FgfGeometryFactory();

Wir erzeugen also ein Select-Command, mit dem man die Daten einer Datenquelle abfragen kann. Dann weisen wir noch dem Select-Command den Namen der FeatureClass zu, die wir abfragen wollen. Im Falle des Shapefiles ist es der blanke Dateiname ohne die Endung „.shp“. Da wir keinen Filter eingebaut haben, werden im nächsten Befehl alle Daten per FeatureReader zurückgegeben.
Außerdem instanziieren wir eine GeometryFactory, die wir später noch benötigen werden.

	//FDO Get Feature Class Defintion (Attribut schema)
	DataTable dt = new DataTable();
	ClassDefinition cDef = reader.GetClassDefinition();
	foreach (PropertyDefinition pDef in cDef.Properties)
	{
		if (pDef.PropertyType == PropertyType.PropertyType_DataProperty)
		{
			DataPropertyDefinition dDef = pDef as DataPropertyDefinition;
			// Get .NET Datatype because of FDO's Enum of datatypes
			Type datatype = GetNETDataType(dDef);
			if (!dt.Columns.Contains(dDef.Name))
				dt.Columns.Add(dDef.Name, datatype);
		}
	}

Im oben aufgezeigten Code-Block wird für das Attributschema der zu erzeugenden FeatureClass zuerst eine DataTable erzeugt und mit den entsprechenden Columns und den zugehörigen Datentypen versehen. Da die FDO API eigene Enumerations für die Datentypen einer Class (FeatureClass oder Tabelle) zurückgibt, für die Definition eines Attributschemas der DataTable jedoch .NET Datentypen gebraucht werden, habe ich eine Hilfsfunktion (GetNETDataType) geschrieben, die eben für ein Element aus dem Enum von FDO einen .NET Datentyp zurückgibt. Den Code dieser Funktion findet man ganz unten auf dieser Seite.

Weiter gehts mit der Instanziierung eines FeatureSets. Ein FeatureSet entspricht in der Terminologie von DotSpatial einer Sammlung von Features. Das Äquivalent der FDO API oder auch von ESRI wäre hier eine FeatureClass. Da wir bereits eine DataTable mit dem passenden Attributschema des Shapefiles aus der FDO API haben, können wir die Struktur der DataTable mit einer Methode dem FeatureSet übergeben.

	DotSpatial.Data.FeatureSet fs = new DotSpatial.Data.FeatureSet();
	fs.CopyTableSchema(dt);

Jetzt holen wir uns die Attributdaten und die Geometry per FDO-FeatureReader aus dem Datensatz. Für jede Datenzeile erzeugen wir ein neues DotSpatial-Feature, das später dann ins FeatureSet wandert. Anschließend gehe ich in einer foreach-Schleife die FDO PropertyDefinition durch, was im Endeffekt einem Attributschema entspricht. Dabei kann man mit der FDO API unterscheiden, ob es eine DataProperty oder eine GeometryProperty ist.

Dann kommt eine etwas unschöne switch-case-Anweisung auf die ich nicht Stolz bin, aber es geht meines Wissens mit der FDO API nicht anders. Man muss beim FeatureReader der FDO wissen, welcher Datentyp dem Attribut zu Grunde liegt, um es mit der entsprechenden Methode (z.B.: reader.GetString()) korrekt auslesen zu können. Am Besten versteht man das im Code. Sind alle Attribute durchlaufen und die entsprechenden Attributwerte ausgelesen, wird die DataRow an die DataRow des DotSpatial-Features gehängt.

	//FDO Get Data & Geometry
	while (reader.ReadNext())
	{
		DataRow dr = dt.NewRow();
		DotSpatial.Data.Feature feature = new DotSpatial.Data.Feature();

		foreach (PropertyDefinition pDef in cDef.Properties)
		{
			if (pDef.PropertyType == PropertyType.PropertyType_DataProperty)
			{
				DataPropertyDefinition dDef = pDef as DataPropertyDefinition;
				switch (dDef.DataType.ToString())
				{
					case "DataType_String":
						if (!reader.IsNull(dDef.Name))
							dr[dDef.Name] = reader.GetString(dDef.Name);
						break;
					case "DataType_Int16":
						if (!reader.IsNull(dDef.Name))
							dr[dDef.Name] = reader.GetInt16(dDef.Name);
						break;
					case "DataType_Int32":
						if (!reader.IsNull(dDef.Name))
							dr[dDef.Name] = reader.GetInt32(dDef.Name);
						break;
					case "DataType_Int64":
						if (!reader.IsNull(dDef.Name))
							dr[dDef.Name] = reader.GetInt64(dDef.Name);
						break;
					case "DataType_DateTime":
						if (!reader.IsNull(dDef.Name))
							dr[dDef.Name] = reader.GetDateTime(dDef.Name);
						break;
					case "DataType_Single":
						if (!reader.IsNull(dDef.Name))
							dr[dDef.Name] = reader.GetSingle(dDef.Name);
						break;
					case "DataType_Double":
						if (!reader.IsNull(dDef.Name))
							dr[dDef.Name] = reader.GetDouble(dDef.Name);
						break;
					case "DataType_Decimal":
						if (!reader.IsNull(dDef.Name))
							dr[dDef.Name] = reader.GetDouble(dDef.Name);
						break;
					case "DataType_Boolean":
						if (!reader.IsNull(dDef.Name))
							dr[dDef.Name] = reader.GetBoolean(dDef.Name);
						break;
					default:
						if (!reader.IsNull(dDef.Name))
							dr[dDef.Name] = reader.GetByte(dDef.Name);
						break;
				}
				feature.DataRow = dr;
			}

Beim Auslesen der Geometrie wirds auch nochmal spannend. Hier wird per FDO API die Geometrie binär ausgelesen und in eine FDO-Geometrie mit Hilfe der anfangs erzeugten GeometryFactory verwandelt. Diese Bytes sind aber eben nicht standardmäßig vom Typ Well-Known-Binary (WKB). Das muss man erst noch einmal aus einer FDO-Geometrie erzeugen. WKB ist deswegen so wichtig, weil es der kleinste gemeinsame Nenner beim Austausch der Geometrie zwischen der FDO API und DotSpatial ist. Hab ich also eine WKB-Geometrie kann ich sie mit dem enstprechenden WKBReader von DotSpatial in eine DotSpatial-Geometrie verwandeln und dem Feature zuweisen. Dann bin ich mit den Attributdaten und der Geometrie durch und muss das Feature dem FeatureSet hinzufügen.

Achtung: Wichtig ist hier noch, dass man nicht die Methode fs.AddFeature() verwendet. Diese akzeptiert auch Features – unterschlägt jedoch die Attributdaten des übergebenen Features.

		if (pDef.PropertyType == PropertyType.PropertyType_GeometricProperty)
		{
			//Get Raw Binary Data and convert it to WKB-DotSpatial
			Byte[] bytes = reader.GetGeometry("Geometry");
			IGeometry geom = geoFac.CreateGeometryFromFgf(bytes);
			Byte[] wkb = geoFac.GetWkb(geom);

			DotSpatial.Data.WKBReader wkbReader = new DotSpatial.Data.WKBReader();
			DotSpatial.Topology.IGeometry dSgeom = wkbReader.Read(wkb);
			feature.BasicGeometry = dSgeom;
		}
	}
	fs.Features.Add(feature);
}
//TODO: Add DotSpatial FeatureSet to map
//for example:
//fs.Name = "test";
//this.map1.Layers.Add(fs);
//this.map1.ZoomToExtents();

Hier noch die Übersetzermethode für die FDO-Datentypen in .NET-Datentypen.

/// <summary>
/// Method for getting the .NET DataType from FDO DataPropertyDefintion
/// </summary>
/// <param name="dDef">FDO DataPropertyDefinition</param>
/// <returns>.NET DataType</returns>
public Type GetNETDataType(OSGeo.FDO.Schema.DataPropertyDefinition dDef)
{
    switch (dDef.DataType.ToString())
    {
	case "DataType_String":
	    return System.Type.GetType("System.String");
	case "DataType_Int16":
	    return System.Type.GetType("System.Int16");
	case "DataType_Int32":
	    return System.Type.GetType("System.Int32");
	case "DataType_Int64":
	    return System.Type.GetType("System.Int64");
	case "DataType_DateTime":
	    return System.Type.GetType("System.DateTime");
	case "DataType_Single":
	    return System.Type.GetType("System.Single");
	case "DataType_Double":
	    return System.Type.GetType("System.Double");
	case "DataType_Decimal":
	    return System.Type.GetType("System.Decimal");
	case "DataType_Boolean":
	    return System.Type.GetType("System.Boolean");
	default:
	    return System.Type.GetType("System.Byte");
    }
}

Soweit die Verwendung der FDO API in .NET für DotSpatial. Ich habe das Ganze anhand eines Shapefiles durchgespielt, was natürlich nicht notwendig wäre, weil DotSpatial ESRI Shapefiles bereits liest. Interessant wird die FDO API aber dann bei OGC Webservices, GDAL/OGR und vor allem bei Oracle Spatial und der ArcSDE. Und um die Verwendung der FDO API bei der ArcSDE soll es im nachfolgenden Artikel gehen..

2 Leute mögen diesen Artikel.

Mai 29th, 2011

Feature Data Object API und DotSpatial in .NET (2)

Ich möchte zunächst einmal eine Verbindung zu einem ESRI Shapefile mit der FDO API erstellen. Die Verbindung zu einer ESRI ArcSDE mit der FDO werde ich später aufzeigen.
Die Dokumentation von FDO gibt unter der Rubrik „Getting Started“ folgende Hinweise:

Write the Code to Connect to a Provider

Do the following:

  1. Get the list of installed providers
  2. Create a connection manager
  3. Create a connection
  4. Get the connection state
  5. Get the connection properties
  6. Get values for the connection properties
  7. Set the connection properties
  8. Open a connection
  9. Open a pending connection

Diese Auflistung werde ich etwas verkürzen. Die Stufen 5-6 kann man nämlich auch mit Hilfe der Dokumentation lösen. Interessant ist jedoch, dass man auch per Code die jeweiligen Verbindungseinstellungen und Möglichkeiten abrufen könnte.

Meine Agenda sieht folgendes vor:

  1. Get the list of installed providers
  2. Create a connection manager
  3. Create a connection
  4. Set the connection properties
  5. Open a connection
  6. Get the connection state
  7. Fetch Data

Und das wollen wir nun mit ein paar Code-Schnippsel füllen.

1. Get the list of installed providers

Mit den folgenden Zeilen kann man sich alle registrierten Provider anzeigen lassen.

//FDO Provider Registry
IProviderRegistry FDORegistry = (IProviderRegistry)FeatureAccessManager.GetProviderRegistry();
ProviderCollection providers = FDORegistry.GetProviders();
string s = "Registered FDO Providers: n";
foreach (Provider p in providers)
{
	s += "# " + p.DisplayName + "n";
}
MessageBox.Show(s, "Registered FDO providers", MessageBoxButtons.OK, MessageBoxIcon.Information);

Hat alles geklappt, sollte man folgende MessageBox sehen:

Achtung: Dabei werden die Provider ausgegeben, wie sie in der Datei providers.xml aufgelistet sind. Das bedeutet jedoch nicht, dass alle Provider und die Assemblies dafür auch wirklich vorhanden sind. Die Provider-Registry liest also nur das XML aus, ohne zu prüfen, ob die referenzierten DLLs auch existieren.

 

2. Create a connection manager / 3. Create a connection / 4. Set the connection properties / 5. Open a connection

Wir brauchen also erst einmal einen ConnectionManager. Mit diesem erzeugen wir eine Verbindung mit der Angabe des entsprechenden Datenproviders als string, z.B. "OSGeo.SHP.3.6".

//FDO Connection - ShapeFile example
IConnectionManager connManager = FeatureAccessManager.GetConnectionManager();
IConnection conn = connManager.CreateConnection("OSGeo.SHP.3.6");
string shapeFile = @"D:\GISDATA\world_adm0.shp";
conn.ConnectionInfo.ConnectionProperties.SetProperty("DefaultFileLocation", shapeFile);
OSGeo.FDO.Connections.ConnectionState connState = conn.Open();

Mit der Wahl des Datenproviders ändern sich die Möglichkeiten, die man bei den Verbindungseigenschaften (ConnectionProperties) setzen kann. Eine Liste der Eigenschaften, die je nach Provider gesetzt werden müssen bzw. können findet man in der FDO Doku – The Open Source Providers ->Connection API. Bei ESRI Shapefile muss man mindestens "DefaultFileLocation" oder "TemporaryFileLocation" und den Pfad zum Shapefile angeben (s.o.).

In der letzten Zeile des obigen Code-Blocks merkt man spätestens beim Aufruf der Methode Open(), ob die benötigten Datenprovider auch wirklich an der Stelle sind, wie sie in der Datei providers.xml beschrieben wurden. Erst dann werden die DLLs samt Abhängigkeiten nämlich zur Laufzeit geladen.

6. Get the connection state

Den Status der Verbindung sollte man abrufen, um eine Fehlerquelle auszuschließen. Das geht auch recht einfach. Anschließend werden wir Daten aus der Datenquelle holen.

//FDO Check ConnectionState
if (connState == OSGeo.FDO.Connections.ConnectionState.ConnectionState_Open)
{
	//Fetch Data
}

Wie man nun aus der FDO-Verbindung Daten herausbekommt (Punkt 7. Fetch Data) bespreche ich im dritten Teil

2 Leute mögen diesen Artikel.

Mai 29th, 2011

Feature Data Object API und DotSpatial in .NET (1)

DotSpatial

DotSpatial ist eine sehr gut dokumentierte, quelloffene GIS – Bibliothek für .NET 4. In diesem Projekt sind einige Open Source basierte GIS-Bibliotheken zusammen geflossen – unter anderem MapWindow6, ProjNET oder GPS.NET3. Auch Entwickler von SharpMap sind bei dem groß angelegten Projekt dabei. Gute Tutorials zu DotSpatial für den Einstieg findet man hier. SharpMap ist zwar ganz toll – es existiert ja auch schon länger für .NET. Aber DotSpatial ist wirklich der Hammer!

Project Vision: DotSpatial aims to provide a free, open source, consistent and dependable set of libraries for the .NET, Silverlight and Mono platforms, enabling developers to easily incorporate spatial data, analysis, and mapping into their applications thereby unleashing the massive potential of GIS in solutions for organizations and communities of all types in a nonrestrictive way.

DotSpatial Architektur

Mit DotSpatial kann man – ohne viel Programmierkenntnisse zu besitzen – ein kleines GIS zusammenbasteln. Mit der entsprechenden Entwicklungsumgebung (MS Visual Studio Express oder SharpDevelop) kommt man schon mit Drag&Drop so weit, dass man die meisten Funktionen von DotSpatial nutzen kann. Bei DotSpatial werden die entsprechenden GIS-Controls mitgeliefert – so wie man es vielleicht von anderen kommerziellen GIS-Bibliotheken kennt. Man merkt schon: Ich bin begeistert..

Im aktuellen Release von DotSpatial kann die Bibliothek mehrere Geodatenformate (sowohl Vektor- als auch Rasterformate) lesen. Schreiben kann es auf jeden Fall ESRI Shapefiles. Seit kurzem gibt es auch einen experimentellen lesenden Datenprovider für PostGIS und SpatiaLite. Schöne Sache – aber es gibt ja noch mehr Geodatenbanken bzw. Middleware auf diesen.  Ich wollte also lesenden Zugriff auf die ESRI ArcSDE in einem eigenen GIS-Viewer realisieren.

Nach etwas Recherche kam ich auf die Feature Data Object API.

Feature Data Object API

Die Feature Data Object API ist ein OSGeo-Projekt, das ursprünglich von Autodesk entwickelt wurde. Es stellt ein abstraktes Objektmodell für Geodaten dar, das viele Geodatenformate lesen und schreiben kann.

FDO Data Access Technology is an API for manipulating, defining and analyzing geospatial information regardless of where it is stored. FDO uses a provider-based model for supporting a variety of geospatial data sources, where each provider typically supports a particular data format or data store. FDO (“Feature Data Object”) is free, open source software licensed under the LGPL.

Feature Data Objects ArchitekturUnter anderem liest diese Bibliothek ESRI ArcSDE, OGR/GDAL, OGC Web Services, MySQL und mit einem proprietären Provider auch Oracle Spatial. Eine vollständige Liste und die Möglichkeiten der Provider findet man hier.

Interessanterweise setzt auch safe die FDO API in ihrer FME-Software ein. Auf jeden Fall ist die Bibliothek Open Source und sowohl in C++ als auch in .NET einsetzbar. Eine GUI für die API ist die FDO Toolbox.

So wunderbar die Dokumentation von DotSpatial für den Anfang auch ist – bei FDO fuer .NET muss man schon länger suchen und einiges ausprobieren. Zwar gibt es einiges an Doku auf der Seite der FDO. Allerdings fehlten mir praktische Beispiele für den schnellen Einstieg. Hilfreich war aber dieses Tutorial:

„Taking Geospatial Data Access to the next level with the FDO API“

Voraussetzungen für die Verwendung der FDO API in .NET

Im Ausgabeverzeichnis (bin\Debug oder bin\Release) des Visual Studio-Projekts müssen folgende Assemblies von FDO liegen:

  • Fdo.dll
  • FDOCommon.dll
  • FDOGeometry.dll
  • FDOMessages.dll
  • FDOSpatial.dll
  • OSGeo.FDO.dll
  • OSGeo.FDO.Common.dll
  • OSGeo.FDO.Geometry.dll
  • Xalan-C_1_11.dll
  • XalanMessages_1_11.dll
  • xerces-c_3_1.dll

Direkt arbeitet man in .NET nur mit den Assemblies im Namespace OSGeo, die dem entsprechend auch im Projekt referenziert werden müssen. Die anderen oben aufgeführten DLLs werden von diesen Assemblies jedoch auch im selben Verzeichnis benötigt. Darüber hinaus existieren noch mehr Assemblies der Provider für die FDO, die jedoch nicht im Ausgabeverzeichnis sein müssen. In der Datei providers.xml werden die Datenprovider und der Pfad zu diesen DLLs und den weiteren abhängigen DLLs angegeben. Dieses XML sollte den Bedürfnissen (v.a. beim Tag LibraryPath) entsprechend angepasst werden und auch im Ausgabepfad der Anwendung gespeichert sein.

Eine mögliche Konfiguration ist dabei, dass die Core-Assemblies (s.o.) und die Datei providers.xml im Ausgabeverzeichnis sind und die eigentlichen Datenprovider mit deren Abhängigkeiten in einem Unterverzeichnis des Ausgabepfads. Im Beispiel oben habe ich einen Unterordner FDO_3.6.0 angelegt.

Außerdem sollte man gleich alle Namespaces von FDO per using-Direktive ins eigene Projekt einbinden:

//FDO
using OSGeo.FDO;
using OSGeo.FDO.Connections;
using OSGeo.FDO.Common;
using OSGeo.FDO.Schema;
using OSGeo.FDO.Commands;
using OSGeo.FDO.Geometry;
using OSGeo.FDO.ClientServices;
using OSGeo.FDO.Commands.DataStore;
using OSGeo.FDO.Commands.Schema;
using OSGeo.FDO.Expression;
using OSGeo.FDO.Commands.Feature;
using OSGeo.FDO.Filter;

Weiter gehts im nächsten Teil..

6 Leute mögen diesen Artikel.

Februar 16th, 2011

G#.Viewer – Laden und Anzeigen eines Shapefiles (3)

Im Dritten Teil dieses Artikel über den G#.Viewer möchte ich auf das Laden und Anzeigen eines Shapefiles mit Sharpmap in C# eingehen. Nun werden wir das erste Mal C# Code nicht nur über den Designer der IDE erstellen, sondern auch selbst Hand anlegen.

Zunächst einmal definiere ich einige Felder am Anfang der Klasse, die im Folgenden verwendet werden:

public partial class frmMain : Form
{
	// Map
	private SharpMap.Map map;
	private bool isMapOpen = false;
	private SharpMap.Layers.VectorLayer vLayer;
	public int width;
	public int height;

	//TreeView
	private TreeNode rootNode;

	public frmMain()
	{
		InitializeComponent();

		//Größe der Karte = Größe des TabPages tabMap
		width = tabMap.Width;
		height = tabMap.Height;
	}

weiterlesen »

Be the first to like.

Februar 16th, 2011

G#.Viewer – GUI – Layerübersicht und OpenFileDialog (2)

Erweiterung der GUI

Im letzten Artikel habe ich einen ersten Entwurf der GUI für den G#.Viewer mit der IDE SharpDevelop angefertigt. Diesen Entwurf möchte ich nun etwas erweitern. Was noch gefehlt hat war eine Übersicht der Layer, die in den Viewer geladen werden. In dieser Übersicht sollen die Layer verschiebbar sein – also die Reihenfolge geändert werden können, in der sie im Kartenbild überlagert werden und die Sichtbarkeit verändert werden können. Außerdem braucht man die eine Layerübersicht wenn man Funktionen auf einzelne Layer beziehen möchte – beispielsweise das Zoomen auf eine bestimmte Ebene oder etwa auch das Entfernen einer Ebene aus dem Viewer.

Deshalb holen wir uns jetzt noch ein TreeView-Control hinzu. Mit Hilfe eines zusätzlichen SplitContainers bauen wir den TreeView in unseren GIS-Viewer ein. Anschließend müssen die Eigenschaften bearbeitet werden, so dass aus dem TreeView auch eine Layerübersicht werden kann. Die weniger interessanten Eigenschaften des TreeViews habe ich ausgeblendet. Ansonsten sollten die fett markierten Eigenschaften wie folgt aussehen:

weiterlesen »

Be the first to like.

Dezember 29th, 2010

G#.Viewer – GUI (1)

G#.Viewer – ein eigener GIS-Datenviewer in C#

Mit der Open Source .NET-Bibliothek SharpMap kann man eigene einfache GIS-Anwendungen mit C# relativ schnell entwickeln. In den nun folgenden Artikeln werde ich die Entwicklung eines einfachen GIS-Viewers beschreiben. Das Projekt heisst „G#.Viewer“ und soll das können: Geodaten darstellen und abfragen.

Als Entwicklungsumgebung (IDE) empfehle ich das Werkzeug SharpDevelop. Auch dies ist Open Source und ist eine profunde Alternative zu Microsoft Visual Studio C# Express Edition oder zur Vollversion Microsoft Visual Studio. In diesen Artikeln werde ich die portable Version 3.2 von SharpDevelop nutzen.

Los geht’s

Zunächst sollten wir uns mit der aktuellen Version der SharpMap-Bibliothek eindecken. Unter „Downloads“ auf der Seite von SharpMap wird man schnell fündig. Am besten läd man sich den aktuellen Quellcode (ich benutze im folgenden das ChangeSet 81947). Es werden zwar auch bereits kompilierte DLLs zum Download angeboten – dabei fehlen jedoch einige nette Sachen wie der Namespace SharpMap.UI, in dem alle nützlichen GUI-Kartenwerkzeuge für Windows Forms zu finden sind (nämlich das MapImage-Control).

Deshalb ist das Motto: selber kompilieren! Ist ja auch mit SharpDevelop ohne weiteres möglich..

weiterlesen »

1 andere Person mag 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.

Dezember 11th, 2010

SharpMap – Open Source GIS mit C#

Ich wollte schon immer eigene Applikationen schreiben. Dabei habe ich mir an verschiedenen Programmiersprachen die Zähne ausgebissen und habe es immer wieder verworfen. Einen guten Einstieg in die Programmierung fand ich dann durch Python: es war leicht zu erlernen und man konnte schnell etwas interessantes entwickeln.

Nachdem ich mit Python und Geoprocessing etwas Erfahrung gesammelt hatte und ein Praktikum bei ESRI Deutschland absolvierte, bei dem ich unter anderem ein bisschen in die Sprache C# schnuppern durfte, wollte ich C# weiter vertiefen.  Vor allem aber wollte ich eigentlich C# mit GIS verbinden. Python spielt sich ja meist auf der Konsole ab. Was ich wollte, war eine Benutzeroberfläche! Einen eigenen GIS-Datenviewer in C# – das war mein Ziel. Nur gab es da ein Problem: .NET-Bibliotheken oder Engines für GIS-Applikationen waren nur kommerziell verfügbar. Ich wollte aber doch nur spielen und lernen – also was tun?

Nach einer längeren Recherche fand ich dann SharpMap. SharpMap ist eine Mapping-Bibilothek für C#. Sie steht unter der LGPL-Lizenz und ist damit für kommerzielle, closed-source Produkte ohne Schwierigkeiten verwendbar – was manch einer als einen Nachteil sehen möchte. (Meine Meinung dazu ist eher pragmatisch: mit der LGPL kann man sowohl Open Source bleiben, als auch Closed Source. Das ist eine Möglichkeit mehr.)

Mit SharpMap ist es relativ leicht, eigene GIS-Applikationen (oder zumindest GIS-Datenviewer) in C# zu erstellen.  Und genau das werde ich in den nächsten Beiträgen etwas vorstellen.

2 Leute mögen diesen Artikel.

Tags: , , , ,
August 18th, 2010

ArcSDE – eine Multiuser-Geodatenbank? (4)

Fazit

Offensichtlich ist es wohl so, dass die Mechanismen einer SDE-Instanz auf der gesamten Datenbankebene verlaufen. Beispielsweise gibt es ein Datenbankschema SDE, in dem alle Verwaltungstabellen residieren.

Dabei müssten die Funktionalitäten doch auf Schemaebene getrennt werden. Es wäre super, wenn der Nutzer A ein eigenes Datenbankschema erhält und dort auch völlig getrennt von anderen Nutzern agieren könnte.

Ich habe versucht zwei verschiedene Schemas auf Datenbankebene anzulegen und auch mit den Rechten der Versionen (Private, Public und Protected) experimentiert – ohne Erfolg. Der Compress-Befehl wirkt sich einfach auf die gesamte Instanz aus und alle Add– und Delete-Tabellen werden im Schema SDE geführt nicht im jeweiligen Schema des Datenthemabesitzers.

In Sachdatenbanken kann man (nach meinem Kenntnisstand) völlig autark in seinem Schema Daten speichern, verändern, Prozeduren anlegen und Historisierungstabellen führen – warum nicht in einer mit ArcSDE erweiterten Datenbank?

Im Falle der ArcSDE lautet das Motto von ESRI leider wieder einmal: „Spatial is special“.

1 andere Person mag diesen Artikel.