Shapefile (Vector) Handling

This section introduces geospatial analysis of shapefiles with gdal, ogr, and osr. For interactive reading and executing code blocks and find geo01-shp.ipynb, or install Python and JupyterLab locally.

Laden Sie eine bestehende Shapefile¶

OSGeo’s ogr module is an excellent source for handling shapefiles. After importing the libraries, the correct driver for handling shapefiles ("ESRI Shapefile") with Python needs to be called. Next, with the driver object (ogr.GetDriverByName("SHAPEFILE")), we can open (instantiate) a shapefile (object with shp_driver.Open("SHAPEFILE")), which contains layer information. It is precisely the layer information (i.e., references to shapefile attributes) that we want to work with. Therefore we have to instantiate a shapefile layer object with shp_dataset.GetLayer().

Erstellen Sie eine neue Formdatei¶

ogr ermöglicht auch eine neue Punkt-, Zeilen- oder Polygonformdatei zu erstellen. Der folgende Codeblock definiert eine Funktion zur Erstellung einer Formdatei, wobei das optionale Keyword-Argument overwrite verwendet wird, um zu kontrollieren, ob eine bestehende Formdatei mit dem gleichen Namen überschrieben werden soll (Standard: True). Der Befehl shp_driver.CreateDataSource(SHP-FILE-DIR) erstellt eine neue Formdatei und der Rest der Funktion fügt der Formdatei eine Schicht hinzu, wenn die optionalen Keyword-Argumente layer_name und layer_type bereitgestellt werden. Beide optionalen Schlüsselwörter müssen strings sein, wobei layer_name für die neue Schicht jeder Name sein kann. layer_type muss entweder "point", "line" oder "polygon" sein, um einen Punkt, (Poly)line oder Polygon formfile zu erstellen. Die Funktion verwendet den geometry_dict dictionary, um der Schicht den richtigen ogr.SHP-TYPE zuzuordnen. Weitere Optionen zur Erweiterung der create_shp(...)-Funktion finden Sie unter pcjericks’ Github page.

The create_shp function is also provided with flusstools (flusstools.geotools.create_shp()) and aids to create a new shapefile (make sure to get the directory right):

Shapefiles können auch in QGIS erstellt und gezeichnet werden und die folgenden Figuren führen durch die Prozedur der Erstellung einer Polygon-Formdatei. Wir brauchen nicht mehr die resultierende Formdatei in diesem Abschnitt, sondern später, damit sie mit Rasterdatensätzen interagiert.

The first step to making a shapefile with QGIS is obviously to run QGIS and create a new project. The following example uses water depth and flow velocity raster data as background information to delineate the so-called morphological unit of slackwater. Both the water depth and flow velocity rasters are part of the River Architect sample datasets (precisely located at RiverArchitect/SampleData/01_Conditions/2100_sample/). After downloading the sample data, they can be imported in QGIS by dragging the files from the Browser panel into the Layers panel. Then:

Holen Sie sich und setzen Sie Shapefile-Projektionen¶

Die in den .prj-Dateien einer Formdatei verwendete Terminologie entspricht den Definitionen im Abschnitt geospatial data. In Python erhalten Sie Informationen über das Koordinatensystem unter shp_layer.GetSpatialRef() der ogr Bibliothek:

Diese GEOGCS-Definition der obigen Formdatei entspricht Esri’s well-known Text (WKT). Da das Formdateiformat von Esri entwickelt wurde, muss esris WKT (esriwkt) Format in .prj Dateien verwendet werden. Das Open Geospatial Consortium (OGC) verwendet einen anderen bekannten Text in ihren EPSG:XXXXDefinitionen (z.B. erhältlich unter spatialreference.org).

Um das Koordinatensystem einer Formdatei neu definieren oder neu definieren zu können, können wir spatialreference.org durch Pythons Standard urllibbibliothek verwenden.

Diese Funktion kann verwendet werden, um eine neue Projektionsdatei zu erstellen:

Eine Offline-Alternative zur Generierung einer .prj-Datei ist die osr-Bibliothek, die zusammen mit gdal:

Transformieren (Re-Projekt) a Shapefile¶

Eine Re-Projektion kann beispielsweise erforderlich sein, um eine Formdatei in EPSG:4326 (z.B. mit QGIS erstellt) in einer Web-GIS-Anwendung (z.B. offene Straßenkarten) zu verwenden, die typischerweise EPSG:3857 verwendet. Um eine andere Projektion auf geometrische Objekte einer Formdatei anzuwenden, reicht es nicht aus, die .prj-Datei einfach neu zu schreiben. Das folgende Beispiel zeigt daher die Neuprojektion der countries.shp formfile aus der Natural Earth quick start kit. Der folgende Workflow führt die Reprojektion durch:

• Die zu transformierende Formdatei befindet sich im Unterverzeichnis geodata/shapefiles/countries.shp und öffnet sie wie oben beschrieben im Python-Skript.

• Lesen und identifizieren Sie das räumliche Referenzsystem in der Eingabeformdatei:

  • Erstellen Sie ein räumliches Referenzobjekt mit in_sr = osr.SpatialReference(str(shapefile.GetSpatialRef())).

  • Erkennen Sie das räumliche Referenzsystem im EPSG-Format mit AutoIdentifyEPSG().

  • Zuordnen des EPSG-formatierten räumlichen Referenzsystems an das räumliche Referenzobjekt der Eingabeformdatei (ImportFromEPSG(int(in_sr.GetAuthorityCode(None)))).

• Erstellen Sie die Ausgaberaumreferenz mit out_sr = osr.SpatialReference() und wenden Sie den Ziel-EPSG-Code an (out_sr.ImportFromEPSG(3857)).

• Erstellen Sie ein Koordinatentransformationsobjekt (coord_trans = osr.CoordinateTransformation(in_sr, out_sr)), das später eine Neuprojektierung von Geometrieobjekten ermöglicht.

• Erstellen Sie die Ausgabeformdatei, die einer Kopie der Eingabeformdatei entspricht (verwenden Sie die oben definierte create_shp Funktion mit layer_name="basemap" und layer_type="line").

• Kopieren Sie die Feldnamen und Typen der Eingabeformdatei:

  • Lesen Sie die Attributschicht aus den Layer Definitionen der Eingabedatei mit in_lyr_def = in_shp_lyr.GetLayerDefn()

  • Iterate durch die Felddefinitionen und an die Ausgabeformfile-Schicht anhängen (out_shp_lyr)

• Iterate durch die Geometriemerkmale in der Eingabeformdatei:

  • Verwenden Sie die Layer-Definitionen der neuen (Output-)Formdatei (out_shp_lyr_def = out_shp_lyr.GetLayerDefn()) um transformierte Geometrieobjekte später anzuhängen.

  • Definieren Sie eine Iterationsvariable in_feature als Instanz von in_shp_lyr.GetNextFeature.

  • In einer while-Schleife wandeln Sie jede Geometrie (geometry = in_feature.GetGeometryRef()) in der Eingabeformdatei um, wandeln Sie die geometry (geometry.Transform(coord_trans)) in eine ogr.Feature() mit der SetGeometry(geometry)-Methode, kopieren Sie Felddefinitionen (nästed for-loop) um und fügen Sie die neue Funktion an die Ausgabeformfile-Schicht (out_shp_lyr_def.CreateFeature(out_feature).

  • Am Ende des while-loop suchen Sie das nächste Feature in den Attributen der Eingabeformdatei mit in_feature = in_shp_lyr.GetNextFeature()

• Entsperren Sie alle Schichten und Formdateien, indem Sie die Objekte mit None überschreiben (Nichts wird in jede Datei geschrieben, solange diese Variablen existieren!).

• Senden Sie die neue Projektion EPSG:3857 mit der oben definierten get_wkt Funktion.

In the case of success, the code sequence in_sr.AutoIdentifyEPSG() returns 0 (i.e., it successfully recognized the EPSG number), but unfortunately, many EPSG numbers are not known to AutoIdentifyEPSG(). In the case that AutoIdentifyEPSG() did not function properly, the method does not return 0, but 7, for example. A workaround for the limited functionality of srs.AutoIdentifyEPSG() is srs.FindMatches. srs.FindMatches returns a matching srs_match from a larger database, which is somewhat nested; for instance, use:

matches = srs.FindMatches()

Dann sieht matches so aus: [(osgeo.osr.SpatialReference, INT)]. So sieht ein kompletter Workaround für srs.AutoIdentifyEPSG() (oder in_sr.AutoIdentifyEPSG() im obigen Codeblock) so aus:

Fügen Sie Fields und Point Features zu einer Shapefile hinzu¶

Ein Shapefile-Feature kann ein Punkt, eine Zeile oder ein Polygon sein, das Feld-Attribute hat (z.B. "id"=1, um zu beschreiben, dass es sich um Polygon Nummer 1 handelt oder einem id Block 1 zugeordnet ist). Feldattribute können mehr als nur ein IDentifier sein und z.B. die Polygon- oder Stadtmarken enthalten, wie im oben dargestellten Beispiel, das cities.shp (shp_driver.Open("geodata/shapefiles/cities.shp")) illustriert.

Um ** einen Punkt formfile** zu erstellen, können wir die oben vorgestellte create_shp-Funktion (oder flusstools.geotools.create_shp(shp_file_dir="")) verwenden und die Projektion mit der get_wkt()-Funktion (auch oben vorgestellt) einstellen. Der folgende Codeblock zeigt die Nutzung beider Funktionen, um einen river zu erstellen. schp Punkt formfile, die drei Punkte an drei Flüssen in Mitteleuropa enthält. Der Codeblock verdeutlicht auch die Erstellung eines Feldes in der Attributtabelle und die Erstellung von drei Punktmerkmalen. Hier ist, wie der Code funktioniert:

• Die Shapefile befindet sich in der Variable rivers_pts. Beachten Sie, dass die layer_type bereits die Art der Geometrien bestimmt, die in der Formdatei verwendet werden können. So wird z.B. das Hinzufügen einer Zeile oder eines Polygons zu einer ogr.wkbPointSchicht zu einer ERROR 1-Nachricht führen.

• Die basemap (Schicht) wird der Variablen lyr = river_pts.GetLayer() zugeschrieben.

• Ein string-Feld wird hinzugefügt und der Attributtabelle angepasst:

  • sofort ein neues Feld mit field_gname = ogr.FieldDefn("FIELD-NAME", ogr.OFTString) (Recall: der Feldname darf nicht mehr als 10 Zeichen haben!)

  • das neue Feld mit lyr.CreateField(field_gname)

  • andere Feldtypen als OFTString können sein:OFTInteger,OFTReal,OFTDate,OFTTime,OFTDateTime, OFTBinary,@@@@@,OFTIntegerList,OFTRealList, oder OFTStringList.

• Fügen Sie drei Punkte hinzu, die unter pt_names = {RIVER-NAME: (x-coordinate, y-coordinate)} in einer Schleife über die Wörterbuchschlüssel gespeichert sind:

  • für jeden neuen Punkt, erstellen Sie eine Funktion als Kind der Schichtdefinitionen mit feature = ogr.Feature(lyr.GetLayerDefn())

  • den Wert des Feldnamens für jeden Punkt mit feature.SetField(FIELD-NAME, FIELD-VALUE)

  • eine Reihe des neuen Punktes im WKT-Format mit wkt = "POINT(X-COORDINATE Y-COORDINATE)"

  • den WKT-formatierten Punkt in eine Punktgeometrie mit point = ogr.CreateGeometryFromWkt(wkt)

  • den neuen Punkt als die Geometrie der neuen Funktion mit feature.SetGeometry(point)

  • die neue Funktion mit lyr.CreateFeature(feature)

• Entsperren Sie die Formdatei, indem Sie die Variable lyr und river_pts mit None überschreiben.

Die daraus resultierende rivers.shp formfile kann in QGIS zusammen mit einer Digitales Oberflächenmodell (DOM) aus dem Natural Earth quick start kit] importiert werden.

Multiline (Polyline) Shapefile¶

Ähnlich wie bei dem Verfahren zum Erstellen und Hinzufügen von Punkten zu einer neuen Punktformdatei kann eine (multi) Linie (oder Polyline) zu einer Shapefile hinzugefügt werden. Die create_shp-Funktion erstellt eine Multi-Line-Formdatei, wenn der Schichttyp als "line" (flusstools.geotools.create_shp(shp_file_dir="", layer_type="line")) definiert ist. Das Koordinatensystem wird mit der oben definierten get_wkt()-Funktion erstellt.

Der folgende Codeblock verwendet die Koordinaten der Städte entlang des Rheins, die in einem Dictionary namens station_names gespeichert sind. Die Stadtnamen werden nicht verwendet und nur die Koordinaten werden mit line.AddPoint(X, Y) angehängt. Wie zuvor wird ein Feld geschaffen, um dem Fluss einen Namen zu geben. Die eigentliche Zeilenfunktion wird wieder als Kind der Schicht mit line_feature = ogr.Feature(lyr.GetLayerDefn()) erstellt. Dieser Codeblock führt eine Linie hervor, die etwa dem Rhein zwischen Frankreich und Deutschland folgt.

Die daraus resultierende rhine_proxy.shp formfile kann in QGIS zusammen mit einer Digitales Oberflächenmodell (DOM) importiert werden und die Städte weisen formfile aus der Natural Earth quick start kit.

Polygon Shapefile¶

Polygone sind Oberflächenpatches, die punktweise, linienweise oder aus einer "Multipolygon"WKB-Definition erstellt werden können. Bei der Erstellung von Polygonen aus Punkten oder Linien wollen wir eine Oberfläche schaffen, und deshalb wird der entsprechende Geometrietyp wkbLinearRing zum Aufbau von Polygonen aus beiden Punkten oder Linien (anstatt wkbPoint bzw. wkbLine) genannt. Der folgende Codeblock bietet ein Beispiel für den Aufbau einer Polygon-Formdatei, die das Hydrauliklabor der Universität Stuttgart abgrenzt. Der Unterschied zwischen dem obigen Beispiel zum Erstellen einer Linienformdatei sind:

• Die Projektion ist EPSG:4326.

• Die Punktkoordinaten werden innerhalb eines ogr.wkbLinearRing Objekt Schritt für Schritt anstatt in einer Schleife über diktionäre Einträge erzeugt.

• Datei, Variable und Feldnamen.

Gestaltdatei von JSON erstellen¶

Die Belastung von Geometriedaten aus in-line definierten Variablen ist in der Praxis umständlich, wo Geospatialdaten häufig auf öffentlichen Plattformen (z.B. Landnutzung oder Abdeckung) zur Verfügung gestellt werden. Das folgende Beispiel verwendet eine JSON-Datei, die mit Kartendienstdaten des Baden-Württemberg State Institute for the Environment, Survey and Nature Conservation (LUBW) erstellt wurde, in der Polygonknoten im WKT-Polygongeometrieformat gespeichert werden ("MultiPolygon (((node1_x node1_y, nodej_x, nodej_y, ... ...)))"):

• Die JSON-Datei (herunterladen hq100-dreisam.json und speichern Sie sie an ein Unterverzeichnis namens /geodata/json/) wird mit Pandas gelesen und die Formdatei wird wie zuvor mit der create_shp-Funktion erstellt.

• Die Projektion ist EPSG:25832.

• Zwei Felder werden in Form von

  • "tbg_name" (der ursprüngliche Stringname der Polygone in den LUBW-Daten) und

  • "area" (ein echtes Nummernfeld, in dem der Polygonbereich unter m² unter polygon.GetArea() berechnet wird).

• Die Polygongeometrien stammen aus den WKT-formatierten Definitionen im Feld "wkt_geom" des pandas Dataframe-Objekts dreisam_inundation.

Auch GeoJSON Daten können verwendet werden, um eine ogr.Geometry mit ogr.createFromGeoJson(FILENAME) zu erstellen:

Geometrische Attribute berechnen¶

Der oben genannte Codeblock verdeutlicht die Nutzung von polygon.GetArea(), um den Polygonbereich in m$^2$ zu berechnen. Die ogr Bibliothek bietet viele weitere Funktionen, um geometrische Attribute von Features zu berechnen und hier eine Zusammenfassung:

• Mehrere Polygone
  wkt_... = ...
polygon_a = ogr.CreateGeometryFromWkt(wkt_1)
polygon_b = ogr.CreateGeometryFromWkt(wkt_2)
polygon_union = polygon_a.Union(polygon_b)

• Intersect zwei Polygone
  polygon_intersection = polygon_a.Intersection(polygon_b)

• Envelope (minimale und maximale Ausdehnung) eines Polygons
  env = polygon.GetEnvelope()
print("minX: %d, minY: %d, maxX: %d, maxY: %d" % (env[0],env[2],env[1],env[3])

• Convex Rumpf (Ummantelungsfläche) mehrerer Geometrien (Punkte, Linien, Polygone)
  all_polygons = ogr.Geometry(ogr.wkbGeometryCollection)
for feature in POLYGON-SOURCE-LAYER: all_polygons.AddGeometry(feature.GetGeometryRef())
convexhull = all_polygons.ConvexHull()
  Save convexhull to shapefile (use the create_shp function as shown in the above examples or read more at pcjericks’ Github pages)
  Tip: To create a tight hull (e.g., of a point cloud), look for concavehull functions.

• Länge (eine Zeile)
  wkt = "LINESTRING (415128.5 5320979.5, 415128.6 5320974.5, 415129.75 5320974.7)"
line = ogr.CreateGeometryFromWkt(wkt)
print("Length = %d" % line.Length())

• Bereich (von einem Polygon): polygon.GetArea() (siehe oben Beispiel)

• Beispiel für die Berechnung Zentrenkoordinaten von polygons.

Export in andere Formate¶

Die obigen Beispiele behandeln nur .shp-Dateien, aber andere Formate können nützlich sein (z.B. Web-Anwendungen erstellen oder auf Google Earth exportieren). Zu diesem Zweck illustrieren die folgenden Absätze die Erstellung von GeoJSON- und KML-Dateien. Mehrere andere Konvertierungen können durchgeführt werden, nicht nur zwischen Dateiformaten, sondern auch zwischen Merkmalstypen. So können beispielsweise Polygone aus Punktwolken erzeugt werden (u.a. mit der oben genannten ConvexHull-Methode). Interessierte Studenten können mehr über Umwandlungen in Michael Diener’s Python Geospatial Analysis Cookbook.

GeoJSON¶

GeoJSON-Dateien können wie zuvor einfach erstellt werden, auch ohne einen Treiber zu aktivieren:

Für eine robustere Bearbeitung und Definition einer Projektion aktivieren Sie den Treiber ogr.GetDriverByName("GeoJSON"). So funktioniert die Erstellung und Manipulation von GeoJSON-Dateien ähnlich wie die oben gezeigten Shapefile-Handler.

KML (Google Earth)¶

Um Punkt-, Zeilen- oder Polygon-Features in Google Earth anzuzeigen, können Funktionen in Googles KML (Keyhole Markup Language), ähnlich der Erstellung einer GeoJSON-Datei, und mit der Funktion geometry.ExportToKML:

Darüber hinaus können ähnlich wie GeoJSON-Dateien und Formdateien KML-Dateien robuster generiert werden (z.B. mit einer definierten Projektion). Sie müssen nur den KML-Treiber (kml_driver = ogr.GetDriverByName("KML")) einleiten und eine KML-Datenquelle definieren (kml_file = kml_driver.CreateDataSource(FILENAME.KML)).