Shapefile (Vector) Handling

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

Laden Sie ein vorhandenes Shapefil¶

OSGeo’s ogr module handles shapefiles. After importing the library, obtain the "ESRI Shapefile" driver with ogr.GetDriverByName("ESRI Shapefile"). The driver can open a shapefile dataset with shp_driver.Open("path/to/file.shp"); obtain its layer with shp_dataset.GetLayer().

Erstellen Sie ein neues Shapefile¶

ogr also enables creating a new point, line, or polygon shapefile. The following code block defines a function for creating a shapefile, where the optional keyword argument overwrite is used to control whether an existing shapefile with the same name should be overwritten (default: True). The command shp_driver.CreateDataSource(SHP-FILE-DIR) creates a new shapefile and the rest of the function adds a layer to the shapefile if the optional keyword arguments layer_name and layer_type are provided. Both optional keywords must be strings, where layer_name can be any name for the new layer. layer_type must be either "point", "line", or "polygon" to create a point, (poly)line, or polygon shapefile, respectively. The function uses the geometry_dict dictionary to assign the correct ogr.SHP-TYPE to the layer. There are more options for extending the create_shp(...) function listed on pcjericks’ Github pages.

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 das Erstellen eines Polygon-Shapefiles. Wir werden die resultierende Shapefile in diesem Abschnitt nicht mehr benötigen, sondern später, um sie mit Rasterdatensätzen interagieren zu lassen.

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:

Erhalten und Set Shapefile Projections¶

The terminology used in the .prj files of a shapefile corresponds to the definitions in the geospatial data section. In Python, information about the coordinate system is available through shp_layer.GetSpatialRef() of the ogr library:

Diese GEOGCS Definition der obigen Shapefile entspricht Esris bekanntem Text (WKT). Da das Shapefile-Format von Esri entwickelt wurde, muss Esris WKT (esriwkt) Format in .prj Dateien verwendet werden. Das Open Geospatial Consortium (OGC) verwendet in seinen EPSG:XXXX-Definitionen einen anderen bekannten Text (z. B. verfügbar unter spatialreference.org].

Um das Koordinatensystem einer Shape-Datei neu zu definieren oder neu zu definieren, können wir spatialreference.org über Pythons Standardbibliothek urllib verwenden.

Diese Funktion kann zum Erstellen einer neuen Projektionsdatei verwendet werden:

Eine Offline-Alternative zum Generieren einer .prj-Datei ist die osr-Bibliothek, die zusammen mit gdal kommt:

Transformation (Re-Projekt) eines Shapefiles¶

A re-projection may be needed, for instance, to use a shapefile in EPSG:4326 (e.g., created with QGIS) in a web-GIS application (e.g., open street maps) that typically uses EPSG:3857. To apply a different projection to geometric objects of a shapefile, it is not enough to simply rewrite the .prj file. Therefore, the following example shows the re-projection of the countries.shp shapefile from the Natural Earth quick start kit. The following workflow performs the reprojection:

• The shapefile to transform is located in the subdirectory geodata/shapefiles/countries.shp and open it in the Python script as described above.

• Lesen und Identifizieren des in der Eingabe-Formdatei verwendeten räumlichen Bezugssystems:

  • Create a spatial reference object with in_sr = osr.SpatialReference(str(shapefile.GetSpatialRef())).

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

  • Weisen Sie das EPSG-formatierte räumliche Bezugssystem dem räumlichen Bezugsobjekt der Eingabe-Shapefile zu (ImportFromEPSG(int(in_sr.GetAuthorityCode(None)))).

• Create the output spatial reference with out_sr = osr.SpatialReference() and apply the target EPSG code (out_sr.ImportFromEPSG(3857)).

• Erstellen Sie ein Koordinatentransformationsobjekt (coord_trans = osr.CoordinateTransformation(in_sr, out_sr)), mit dem Geometrieobjekte später neu projiziert werden können.

• Create the output shapefile, which will correspond to a copy of the input shapefile (use the above-defined create_shp function with layer_name="basemap" and layer_type="line").

• Kopieren Sie die Feldnamen und Typen der Eingabe-Shapefile:

  • Read the attribute layer from the input file’s layer definitions with in_lyr_def = in_shp_lyr.GetLayerDefn()

  • Iterieren Sie die Felddefinitionen und fügen Sie sie an die Ausgabe-Shapefile-Schicht an (out_shp_lyr)

• Iterieren Sie durch die Geometriemerkmale in der Eingabe-Shapefile:

  • Verwenden Sie die Layerdefinitionen der neuen Shapefile (out_shp_lyr_def = out_shp_lyr.GetLayerDefn()), um später transformierte Geometrieobjekte anzuhängen.

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

  • In a while loop, instantiate every geometry (geometry = in_feature.GetGeometryRef()) in the input shapefile, transform the geometry (geometry.Transform(coord_trans)), convert it to an ogr.Feature() with the SetGeometry(geometry) method, copy field definitions (nested for-loop), and append the new feature to the output shapefile layer (out_shp_lyr_def.CreateFeature(out_feature)).

  • At the end of the while-loop, look for the next feature in the input shapefile’s attributes with in_feature = in_shp_lyr.GetNextFeature()

• Unlock (release) all layers and shapefiles by overwriting the objects with None (nothing is written to any file as long as these variables exist!).

• Assign the new projection EPSG:3857 using the above-defined get_wkt function.

Im Erfolgsfall gibt die Codesequenz in_sr.AutoIdentifyEPSG() 0 zurück (d.h. sie erkennt erfolgreich die EPSG-Nummer), aber leider sind viele EPSG-Nummern AutoIdentifyEPSG() nicht bekannt. Falls AutoIdentifyEPSG() nicht richtig funktioniert hat, gibt die Methode nicht 0 zurück, sondern beispielsweise 7. Ein Workaround für die eingeschränkte Funktionalität von srs.AutoIdentifyEPSG() ist srs.FindMatches. srs.FindMatches gibt ein matching srs_match aus einer größeren Datenbank zurück, die etwas verschachtelt ist; zum Beispiel:

matches = srs.FindMatches()

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

Hinzufügen von Feldern und Punktfunktionen zu einem Shapefile¶

A shapefile feature can be a point, a line, or a polygon, which has field attributes (e.g., "id"=1 to describe that this is polygon number 1 or associated to an id block 1). Field attributes can be more than just an IDentifier and include, for example, the polygon area or city labels as in the above-shown example illustrating cities.shp (shp_driver.Open("geodata/shapefiles/cities.shp")).

To create a point shapefile, we can use the above-introduced create_shp function (or flusstools.geotools.create_shp(shp_file_dir="")) and set the projection with the get_wkt() function (also above introduced). The following code block shows the usage of both functions to create a river.shp point shapefile that contains three points located at three rivers in central Europe. The code block also illustrates the creation of a field in the attribute table and the creation of three point features. Here is how the code works:

• The shapefile is located in the rivers_pts variable. Note that the layer_type already determines the type of geometries that can be used in the shapefile. For instance, adding a line or polygon feature to an ogr.wkbPoint layer will result in an ERROR 1 message.

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

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

  • instantiate a new field with field_gname = ogr.FieldDefn("FIELD-NAME", ogr.OFTString) (recall: the field name may not have more than 10 characters!)

  • append the new field to the shapefile with lyr.CreateField(field_gname)

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

• Add three points stored in pt_names = {RIVER-NAME: (x-coordinate, y-coordinate)} in a loop over the dictionary keys:

  • for every new point, create a feature as a child of the layer definitions with feature = ogr.Feature(lyr.GetLayerDefn())

  • set the value of the field name for every point with feature.SetField(FIELD-NAME, FIELD-VALUE)

  • create a string of the new point in WKT format with wkt = "POINT(X-COORDINATE Y-COORDINATE)"

  • convert the WKT-formatted point into a point-type geometry with point = ogr.CreateGeometryFromWkt(wkt)

  • set the new point as the new feature’s geometry with feature.SetGeometry(point)

  • append the new feature to the layer with lyr.CreateFeature(feature)

• Entsperren (Freigeben) der Shapefile durch Überschreiben der Variablen lyr und river_pts mit None.

The resulting rivers.shp shapefile can be imported in QGIS along with a DEM from the Natural Earth quick start kit.

Multiline (Polyline) Shapefil¶

Similar to the procedure for creating and adding points to a new point shapefile, a (multi) line (or polyline) can be added to a shapefile. The create_shp function creates a multi-line shapefile when the layer type is defined as "line" (flusstools.geotools.create_shp(shp_file_dir="", layer_type="line")). The coordinate system is created with the above-defined get_wkt() function.

Der folgende Codeblock verwendet die Koordinaten der Städte entlang des Rheins, gespeichert in einem Wörterbuch namens station_names. Die Stadtnamen werden nicht verwendet und nur die Koordinaten werden mit line.AddPoint(X, Y) angehängt. Wie zuvor wird ein Feld angelegt, um dem Fluss einen Namen zu geben. Die eigentliche Linienfunktion wird wieder als Kind der Schicht mit line_feature = ogr.Feature(lyr.GetLayerDefn()) erstellt. Durch den Betrieb dieses Codeblocks entsteht eine Linie, die ungefähr dem Rhein zwischen Frankreich und Deutschland folgt.

The resulting rhine_proxy.shp shapefile can be imported in QGIS along with a DEM and the cities point shapefile from the Natural Earth quick start kit.

Polygon Shapefil¶

Polygons are surface patches that can be created point-by-point, line-by-line, or from a "Multipolygon" WKB definition. When creating polygons from points or lines, we want to create a surface and this is why the corresponding geometry type is called wkbLinearRing for building polygons from both points or lines (rather than wkbPoint or wkbLine, respectively). The following code block features an example for building a polygon shapefile delineating the hydraulic laboratory of the University of Stuttgart. The difference between the above example for creating a line shapefile are:

• Die Projektion ist EPSG:4326.

• The point coordinates are generated within an ogr.wkbLinearRing object step-by-step rather than in a loop over dictionary entries.

• Datei-, Variablen- und Feldnamen.

Build Shapefile von JSON¶

Das Laden von Geometriedaten aus inline definierten Variablen ist in der Praxis umständlich, da Geodaten häufig auf öffentlichen Plattformen bereitgestellt werden (z. B. Landnutzung oder Abdeckung). Im folgenden Beispiel wird eine JSON-Datei verwendet, die mit Kartendienstdaten des Baden-Württemberg State Institute for the Environment, Survey and Nature Conservation (LUBW) generiert wurde, wobei Polygonknoten im WKT-Polygongeometrieformat ("MultiPolygon (((node1_x node1_y, nodej_x, nodej_y, ... ...)))") gespeichert sind:

• Die JSON-Datei (hq100-dreisam.json herunterladen und in einem Unterverzeichnis namens /geodata/json/ speichern) wird mit Pandas gelesen und die Shape-Datei 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 reelles Zahlenfeld, in dem der Polygonbereich in m² unter Verwendung von polygon.GetArea() berechnet wird).

• The polygon geometries are derived from the WKT-formatted definitions in the "wkt_geom" field of the pandas dataframe object dreisam_inundation.

Also GeoJSON data can be used to create an ogr.Geometry with ogr.createFromGeoJson(FILENAME):

Geometrische Attribute berechnen¶

The above code block illustrates the usage of polygon.GetArea() to calculate the polygon area in m$^2$. The ogr library provides many more functions to calculate geometric attributes of features and here is a summary:

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

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

• Envelope (minimum and maximum extents) of a polygon
  env = polygon.GetEnvelope()
print("minX: %d, minY: %d, maxX: %d, maxY: %d" % (env[0],env[2],env[1],env[3])

• Konvexer Rumpf (Hüllflä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 (einer Linie)
  wkt = "LINESTRING (415128.5 5320979.5, 415128.6 5320974.5, 415129.75 5320974.7)"
line = ogr.CreateGeometryFromWkt(wkt)
print("Length = %d" % line.Length())

• Area (of a polygon): polygon.GetArea() (see above example)

• Beispiel für die Berechnung von Zentrumskoordinaten von Polygons].

Export in andere Formate¶

The above examples deal with .shp files only, but other formats can be useful (e.g., to create web applications or export to Google Earth). To this end, the following paragraphs illustrate the creation of GeoJSON and KML files. Several other conversions can be performed, not only between file formats but also between feature types. For instance, polygons can be created from point clouds (among others with the ConvexHull method mentioned above). Interested students can learn more about conversions in Michael Diener’s Python Geospatial Analysis Cookbook.

GeoJSON¶

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

For more robust file handling and defining a projection, activate the driver ogr.GetDriverByName("GeoJSON"). Thus, the creation and manipulation of GeoJSON files work similarly to the shapefile handlers shown above.

KML (Google Earth)¶

Um Punkt-, Linien- oder Polygonfunktionen in Google Earth anzuzeigen, können Funktionen in Googles KML (Keyhole Markup Language) eingefügt werden, ähnlich der Erstellung einer GeoJSON-Datei und mit der Funktion geometry.ExportToKML:

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