{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "(chpt-shp)=\n", "# Shapefile (Vector) Handling\n", "\n", "This section introduces geospatial analysis of shapefiles with gdal, ogr, and osr. For interactive reading and executing code blocks [![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/hydro-informatics/jupyter-python-course/main?filepath=jupyter) and find *geo01-shp.ipynb*, or install {ref}`Python ` and {ref}`JupyterLab ` locally.\n", "\n", "```{admonition} Requirements\n", "* Make sure to understand {ref}`shapefiles ` and {ref}`vector data `\n", "* Accomplish the {ref}`QGIS tutorial `\n", "```\n", "\n", "```{admonition} Use flusstools\n", ":class: tip\n", "The core functions featured in this section are also implemented in {{ ft_url }}. To use those functions, make sure flusstools is installed and import it as follows: `from flusstools import geotools`. Some of the functions shown in this tutorial can then be used with `geotools.function_name()`.\n", "```\n", "\n", "```{admonition} Watch this section and the Python tutorials in video formats\n", ":class: tip, dropdown\n", "\n", "

Watch this section as a video on the @Hydro-Morphodynamics channel on YouTube.

\n", "```\n", "\n", "\n", "## Load an Existing Shapefile\n", "\n", "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()`." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "from osgeo import ogr\n", "shp_driver = ogr.GetDriverByName(\"ESRI Shapefile\")\n", "shp_dataset = shp_driver.Open(\"geodata/shapefiles/cities.shp\")\n", "shp_layer = shp_dataset.GetLayer()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "```{tip}\n", "To get a full list of supported `ogr` drivers (e.g., for `DXF`, `ESRIJSON`, `GPS`, `PDF`, `SQLite`, `XLSX`, and many more), [download the script `get_ogr_drivers.py`](https://raw.githubusercontent.com/hydro-informatics/material-py-codes/main/geo/get_ogr_drivers.py) from hydro-informatics on Github.\n", "```\n", "\n", "```{admonition} Do not use **from gdal import ogr**\n", ":class: important\n", "\n", "Importing `ogr` with `from gdal import ogr` is deprecated since GDAL v3. Therefore, `ogr` must be imported from `osgeo`.\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "(create-shp)=\n", "## Create a New Shapefile \n", "\n", "`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`).\n", "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 *string*s, 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 [*pcjerick*'s Github pages](https://pcjericks.github.io/py-gdalogr-cookbook/geometry.html)." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "from osgeo import ogr\n", "import os\n", "\n", "\n", "def create_shp(shp_file_dir, overwrite=True, *args, **kwargs):\n", " \"\"\"\n", " :param shp_file_dir: STR of the (relative shapefile directory (ends on \".shp\")\n", " :param overwrite: [optional] BOOL - if True, existing files are overwritten\n", " :kwarg layer_name: [optional] STR of the layer_name - if None: no layer will be created (max. 13 chars)\n", " :kwarg layer_type: [optional] STR (\"point, \"line\", or \"polygon\") of the layer_name - if None: no layer will be created\n", " :output: returns an ogr shapefile layer\n", " \"\"\"\n", " shp_driver = ogr.GetDriverByName(\"ESRI Shapefile\")\n", "\n", " # check if output file exists if yes delete it\n", " if os.path.exists(shp_file_dir) and overwrite:\n", " shp_driver.DeleteDataSource(shp_file_dir)\n", "\n", " # create and return new shapefile object\n", " new_shp = shp_driver.CreateDataSource(shp_file_dir)\n", "\n", " # create layer if layer_name and layer_type are provided\n", " if kwargs.get(\"layer_name\") and kwargs.get(\"layer_type\"):\n", " # create dictionary of ogr.SHP-TYPES\n", " geometry_dict = {\"point\": ogr.wkbPoint,\n", " \"line\": ogr.wkbMultiLineString,\n", " \"polygon\": ogr.wkbMultiPolygon}\n", " # create layer\n", " try:\n", " new_shp.CreateLayer(str(kwargs.get(\"layer_name\")),\n", " geom_type=geometry_dict[str(kwargs.get(\"layer_type\").lower())])\n", " except KeyError:\n", " print(\"Error: Invalid layer_type provided (must be 'point', 'line', or 'polygon').\")\n", " except TypeError:\n", " print(\"Error: layer_name and layer_type must be string.\")\n", " except AttributeError:\n", " print(\"Error: Cannot access layer - opened in other program?\")\n", " return new_shp" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The `create_shp` function is also provided with {{ ft_url }} ([`flusstools.geotools.create_shp()`](https://flusstools.readthedocs.io/en/latest/geotools.html#module-flusstools.geotools.shp_mgmt)) and aids to create a new shapefile (make sure to get the directory right):" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "a_new_shp_file = create_shp(r\"\" + os.getcwd() + \"/geodata/shapefiles/new_polygons.shp\", layer_name=\"basemap\", layer_type=\"polygon\")\n", "\n", "# release data source\n", "a_new_shape_file = None" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "```{important}\n", "A **shapefile name** may **not** have **more than 13 characters** and a **field name** may **not** have **more than 10 characters** (read more in [Esri's shapefile docs](http://resources.arcgis.com/en/help/main/10.1/index.html#//005600000013000000)).\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Shapefiles can also be created and drawn in QGIS and the following figures guide through the procedure of creating a polygon shapefile. We will not need the resulting shapefile in this section anymore, but later for making it interact with raster datasets.\n", "\n", "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*](https://www.sciencedirect.com/science/article/pii/S0169555X14000099). Both the water depth and flow velocity rasters are part of the [*River Architect* sample datasets](https://github.com/RiverArchitect/SampleData/archive/master.zip) (precisely located at [`RiverArchitect/SampleData/01_Conditions/2100_sample/`](https://github.com/RiverArchitect/SampleData/tree/master/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:\n", "\n", "```{figure} ../img/qgis-create-shp.png\n", ":alt: qgis create new shapefile\n", ":name: qgis-create-shp-py\n", "\n", "QGIS: Create Layer > New Shapefile Layer...\n", "```\n", "\n", "```{figure} ../img/qgis-new-shp.png\n", ":alt: qgis define new shapefile\n", ":name: qgis-new-shp-py\n", "\n", "Define New Shapefile Layer\n", "```\n", "\n", "```{figure} ../img/qgis-toggle-editing.png\n", ":alt: qgis shapefile toggle-editing\n", ":name: qgis-toggle-editing-py\n", "\n", "Enable shapefile editing\n", "```\n", "\n", "```{figure} ../img/qgis-draw-polygon.png\n", ":alt: qgis draw polygon shapefile\n", ":name: qgis-draw-polygon-py\n", "\n", "Draw polygons in a shapefile in QGIS.\n", "```\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "(prj-shp)=\n", "## Get and Set Shapefile Projections \n", "\n", "The terminology used in the `.prj` files of a shapefile corresponds to the definitions in the {ref}`geospatial data ` section. In Python, information about the coordinate system is available through `shp_layer.GetSpatialRef()` of the `ogr` library:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "GEOGCS[\"WGS 84\",\n", " DATUM[\"WGS_1984\",\n", " SPHEROID[\"WGS 84\",6378137,298.257223563,\n", " AUTHORITY[\"EPSG\",\"7030\"]],\n", " AUTHORITY[\"EPSG\",\"6326\"]],\n", " PRIMEM[\"Greenwich\",0,\n", " AUTHORITY[\"EPSG\",\"8901\"]],\n", " UNIT[\"degree\",0.0174532925199433,\n", " AUTHORITY[\"EPSG\",\"9122\"]],\n", " AXIS[\"Latitude\",NORTH],\n", " AXIS[\"Longitude\",EAST],\n", " AUTHORITY[\"EPSG\",\"4326\"]]\n" ] } ], "source": [ "shp_srs = shp_layer.GetSpatialRef()\n", "print(shp_srs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This `GEOGCS` definition of the above shapefile corresponds to Esri's *well-known* text (WKT). Since the shapefile format was developed by Esri, Esri's WKT (**esriwkt**) format must be used in `.prj` files. The *Open Geospatial Consortium* (*OGC*) uses a different well-known text in their `EPSG:XXXX` definitions (e.g., available at [spatialreference.org](http://www.spatialreference.org)). " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "```java\n", "GEOGCS[\"WGS 84\",\n", " DATUM[\"WGS_1984\", SPHEROID[\"WGS84\", 6378137, 298.257223563, AUTHORITY[\"EPSG\", \"7030\"]], AUTHORITY[\"EPSG\",\"6326\"]],\n", " PRIMEM[\"Greenwich\", 0, AUTHORITY[\"EPSG\", \"8901\"]],\n", " UNIT[\"degree\",0.01745329251994328, AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]]\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To redefine or newly define the coordinate system of a shapefile, we can use [spatialreference.org](http://www.spatialreference.org) through Python's default `urllib` library.\n", "\n", "```{note}\n", "Running the following code block requires an internet connection.\n", "```" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "import urllib\n", "\n", "# function to get spatialreferences with epsg code\n", "def get_esriwkt(epsg): \n", " # usage get_epsg_code(4326)\n", " try:\n", " with urllib.request.urlopen(\"http://spatialreference.org/ref/epsg/{0}/esriwkt/\".format(epsg)) as response:\n", " return str(response.read()).strip(\"b\").strip(\"'\")\n", " except Exception:\n", " pass\n", " try:\n", " with urllib.request.urlopen(\"http://spatialreference.org/ref/sr-org/epsg{0}-wgs84-web-mercator-auxiliary-sphere/esriwkt/\".format(epsg)) as response:\n", " return str(response.read()).strip(\"b\").strip(\"'\")\n", " # sr-org codes are available at \"https://spatialreference.org/ref/sr-org/{0}/esriwkt/\".format(epsg)\n", " # for example EPSG:3857 = SR-ORG:6864 -> https://spatialreference.org/ref/sr-org/6864/esriwkt/ = EPSG:3857\n", " except Exception:\n", " print(\"ERROR: Could not find epsg code on spatialreference.org. Returning default WKT(epsg=4326).\")\n", " return 'GEOGCS[\"GCS_WGS_1984\",DATUM[\"D_WGS_1984\",SPHEROID[\"WGS_1984\",6378137,298.257223563]],PRIMEM[\"Greenwich\",0],UNIT[\"Degree\",0.017453292519943295],UNIT[\"Meter\",1]]'" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This function can be used to create a new projection file:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Wrote projection file : GEOGCS[\"GCS_WGS_1984\",DATUM[\"D_WGS_1984\",SPHEROID[\"WGS_1984\",6378137,298.257223563]],PRIMEM[\"Greenwich\",0],UNIT[\"Degree\",0.017453292519943295]]\n" ] } ], "source": [ "# open the hypy-area shapefile\n", "shp_file = \"hypy-area\"\n", "\n", "# create new .prj file for the shapefile (.shp and .prj must have the same name)\n", "with open(\"geodata/shapefiles/{0}.prj\".format(shp_file), \"w\") as prj:\n", " epsg_code = get_esriwkt(4326)\n", " prj.write(epsg_code)\n", " print(\"Wrote projection file : \" + epsg_code)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "An offline alternative for generating a `.prj` file is the `osr` library that comes along with `gdal`:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "from osgeo import osr\n", "\n", "def get_wkt(epsg, wkt_format=\"esriwkt\"):\n", " default = 'GEOGCS[\"GCS_WGS_1984\",DATUM[\"D_WGS_1984\",SPHEROID[\"WGS_1984\",6378137,298.257223563]],PRIMEM[\"Greenwich\",0],UNIT[\"Degree\",0.017453292519943295],UNIT[\"Meter\",1]]'\n", " spatial_ref = osr.SpatialReference()\n", " try:\n", " spatial_ref.ImportFromEPSG(epsg)\n", " except TypeError:\n", " print(\"ERROR: epsg must be integer. Returning default WKT(epsg=4326).\")\n", " return default\n", " except Exception:\n", " print(\"ERROR: epsg number does not exist. Returning default WKT(epsg=4326).\")\n", " return default\n", " if wkt_format==\"esriwkt\":\n", " spatial_ref.MorphToESRI()\n", " # return a nicely formatted WKT string (alternatives: ExportToPCI(), ExportToUSGS(), or ExportToXML())\n", " return spatial_ref.ExportToPrettyWkt()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "(reproject-shp)=\n", "## Transform (Re-project) a Shapefile\n", "\n", "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](http://naciscdn.org/naturalearth/packages/Natural_Earth_quick_start.zip). The following workflow performs the reprojection:\n", "\n", "* Save shapefile to transform is located in the subdirectory `geodata/shapefiles/countries.shp` and open it in the Python script as described above.\n", "* Read and identify the spatial reference system used in the input shapefile:\n", " - Create a spatial reference object with `in_sr = osr.SpatialReference(str(shapefile.GetSpatialRef()))`.\n", " - Detect the spatial reference system in EPSG format with `AutoIdentifyEPSG()`.\n", " - Assign the EPSG-formatted spatial reference system to the spatial reference object of the input shapefile (`ImportFromEPSG(int(in_sr.GetAuthorityCode(None)))`).\n", "* Create the output spatial reference with `out_sr = osr.SpatialReference()` and apply the target EPSG code (`out_sr.ImportFromEPSG(3857)`).\n", "* Create a coordinate transformation object (`coord_trans = osr.CoordinateTransformation(in_sr, out_sr)`) that enables re-projecting geometry objects later.\n", "* 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\"`).\n", "* Copy the field names and types of the input shapefile:\n", " - Read the attribute layer from the input file's layer definitions with `in_lyr_def = in_shp_lyr.GetLayerDefn()`\n", " - Iterate through the field definitions and append them to the output shapefile layer (`out_shp_lyr`)\n", "* Iterate through the geometry features in the input shapefile:\n", " - Use the new (output) shapefile's layer definitions (`out_shp_lyr_def = out_shp_lyr.GetLayerDefn()`) to append transformed geometry objects later.\n", " - Define an iteration variable `in_feature` as an instance of `in_shp_lyr.GetNextFeature`.\n", " - 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)`).\n", " - 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()`\n", "* Unlock (release) all layers and shapefiles by overwriting the objects with `None` (nothing is written to any file as long as these variables exist!).\n", "* Assign the new projection EPSG:3857 using the above-defined `get_wkt` function." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "from osgeo import ogr\n", "from osgeo import osr\n", "\n", "shp_driver = ogr.GetDriverByName(\"ESRI Shapefile\")\n", "\n", "# open input shapefile and layer\n", "in_shp = shp_driver.Open(r\"\" + os.path.abspath(\"\") + \"/geodata/shapefiles/countries.shp\")\n", "in_shp_lyr = in_shp.GetLayer()\n", "\n", "# get input SpatialReference\n", "in_sr = osr.SpatialReference(str(in_shp_lyr.GetSpatialRef()))\n", "# auto-detect epsg\n", "in_sr.AutoIdentifyEPSG()\n", "# assign input SpatialReference\n", "in_sr.ImportFromEPSG(int(in_sr.GetAuthorityCode(None)))\n", "\n", "# create SpatialReference for new shapefile\n", "out_sr = osr.SpatialReference()\n", "out_sr.ImportFromEPSG(3857)\n", "\n", "# create a CoordinateTransformation object\n", "coord_trans = osr.CoordinateTransformation(in_sr, out_sr)\n", "\n", "# create output shapefile and get layer\n", "out_shp = create_shp(r\"\" + os.path.abspath(\"\") + \"/geodata/shapefiles/countries-web.shp\", layer_name=\"basemap\", layer_type=\"line\")\n", "out_shp_lyr = out_shp.GetLayer()\n", "\n", "# look up layer (features) definitions in input shapefile\n", "in_lyr_def = in_shp_lyr.GetLayerDefn()\n", "# copy field names of input layer attribute table to output layer\n", "for i in range(0, in_lyr_def.GetFieldCount()):\n", " out_shp_lyr.CreateField(in_lyr_def.GetFieldDefn(i))\n", "\n", "# instantiate feature definitions object for output layer (currently empty)\n", "out_shp_lyr_def = out_shp_lyr.GetLayerDefn()\n", "\n", "# iteratively append all input features in new projection\n", "in_feature = in_shp_lyr.GetNextFeature()\n", "while in_feature:\n", " # get the input geometry\n", " geometry = in_feature.GetGeometryRef()\n", " # re-project (transform) geometry to new system\n", " geometry.Transform(coord_trans)\n", " # create new output feature\n", " out_feature = ogr.Feature(out_shp_lyr_def)\n", " # assign in-geometry to output feature and copy field values\n", " out_feature.SetGeometry(geometry)\n", " for i in range(0, out_shp_lyr_def.GetFieldCount()):\n", " out_feature.SetField(out_shp_lyr_def.GetFieldDefn(i).GetNameRef(), in_feature.GetField(i))\n", " # add the feature to the shapefile\n", " out_shp_lyr.CreateFeature(out_feature)\n", " # prepare next iteration\n", " in_feature = in_shp_lyr.GetNextFeature()\n", "\n", "# release shapefiles and layers\n", "in_shp = None\n", "in_shp_lyr = None\n", "out_shp = None\n", "out_shp_lyr = None\n", "\n", "# create .prj file for output shapefile for web application references\n", "with open(r\"\" + os.path.abspath('') + \"/geodata/shapefiles/countries-web.prj\", \"w+\") as prj:\n", " prj.write(get_wkt(3857))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "```{admonition} Challenge\n", "Re-write the above code block in a `re_project(shp_file, target_epsg)` function.\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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 propperly, 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:
\n", "\n", "```python\n", "matches = srs.FindMatches()\n", "```\n", "\n", "Then, `matches` looks like this: `[(osgeo.osr.SpatialReference, INT)]`. Therefore, a complete workaround for `srs.AutoIdentifyEPSG()` (or `in_sr.AutoIdentifyEPSG()` in the code block above) looks like this:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "# set epsg and create spatial reference object\n", "epsg = 3857\n", "srs = osr.SpatialReference()\n", "srs.ImportFromEPSG(epsg)\n", "\n", "# identify spatial reference\n", "auto_detect = srs.AutoIdentifyEPSG()\n", "if auto_detect is not 0:\n", " srs = srs.FindMatches()[0][0] # Find matches returns list of tuple of SpatialReferences\n", " srs.AutoIdentifyEPSG() # Re-perform auto-identification" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "(add-field)=\n", "## Add Fields and Point Features to a Shapefile \n", "\n", "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 *ID*entifier 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\")`). \n", "\n", "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:\n", "\n", "* 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.\n", "* The `basemap` (layer) is attributed to the variable `lyr = river_pts.GetLayer()`.\n", "* A *string*-type field is added and appended to the attribute table:\n", " - instantiate a new field with `field_gname = ogr.FieldDefn(\"FIELD-NAME\", ogr.OFTString)` (recall: the field name may not have more than 10 characters!)\n", " - append the new field to the shapefile with `lyr.CreateField(field_gname)` \n", " - other field types than `OFTString` can be: `OFTInteger`, `OFTReal`, `OFTDate`, `OFTTime`, `OFTDateTime`, `OFTBinary`, `OFTIntegerList`, `OFTRealList`, or `OFTStringList`.\n", "* Add three points stored in `pt_names = {RIVER-NAME: (x-coordinate, y-coordinate)}` in a loop over the dictionary keys:\n", " - for every new point, create a feature as a child of the layer defintions with `feature = ogr.Feature(lyr.GetLayerDefn())`\n", " - set the value of the field name for every point with `feature.SetField(FIELD-NAME, FIELD-VALUE)`\n", " - create a string of the new point in WKT format with `wkt = \"POINT(X-COORDINATE Y-COORDINATE)\"`\n", " - convert the WKT-formatted point into a point-type geometry with `point = ogr.CreateGeometryFromWkt(wkt)`\n", " - set the new point as the new feature's geometry with `feature.SetGeometry(point)`\n", " - append the new feature to the layer with `lyr.CreateFeature(feature)`\n", "* Unlock (release) the shapefile by overwriting the `lyr` and `river_pts` variable with `None`.\n", "\n", "```{important}\n", "The operations are not written to the shapefile if the `lyr` and `river_pts` objects are not overwritten with `None`.\n", "```" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "shp_dir = r\"\" + os.path.abspath('') + \"/geodata/shapefiles/rivers.shp\"\n", "river_pts = create_shp(shp_dir, layer_name=\"basemap\", layer_type=\"point\")\n", "\n", "# create .prj file for the shapefile for web application references\n", "with open(shp_dir.split(\".shp\")[0] + \".prj\", \"w+\") as prj:\n", " prj.write(get_wkt(3857))\n", "\n", "# get basemap layer\n", "lyr = river_pts.GetLayer()\n", "\n", "# add string field \"rivername\"\n", "field_gname = ogr.FieldDefn(\"rivername\", ogr.OFTString)\n", "lyr.CreateField(field_gname)\n", "\n", "# names and coordinates of central EU rivers in EPSG:3857 WG84 / Pseudo-Mercator\n", "pt_names = {\"Aare\": (916136.03, 6038687.72),\n", " \"Ain\": (623554.12, 5829154.69),\n", " \"Inn\": (1494878.95, 6183793.83)}\n", "\n", "# add the three rivers as points to the basemap layer\n", "for n in pt_names.keys():\n", " # create Feature as child of the layer\n", " pt_feature = ogr.Feature(lyr.GetLayerDefn())\n", " # define value n (river) in the rivername field\n", " pt_feature.SetField(\"rivername\", n)\n", " # use WKT format to add a point geometry to the Feature\n", " wkt = \"POINT(%f %f)\" % (float(pt_names[n][0]), float(pt_names[n][1]))\n", " point = ogr.CreateGeometryFromWkt(wkt)\n", " pt_feature.SetGeometry(point)\n", " # append the new feature to the basement layer\n", " lyr.CreateFeature(pt_feature)\n", " \n", "# release files\n", "lyr = None\n", "river_pts = None" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The resulting `rivers.shp` shapefile can be imported in QGIS along with a DEM from the [Natural Earth quick start kit](http://naciscdn.org/naturalearth/packages/Natural_Earth_quick_start.zip).\n", "\n", "```{figure} ../img/qgis-rivers.png\n", ":alt: river map shapefile\n", ":name: qgis-rivers-py\n", "\n", "The river names bound to points in a point shapefile shown in QGIS.\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "(line-create)=\n", "## Multiline (Polyline) Shapefile \n", "\n", "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.\n", "\n", "```{tip}\n", "The term *multi-line* is used in *OGC* and `ogr`, while *polyline* is used in Esri GIS environments.\n", "```\n", "\n", "The following code block uses the coordinates of cities along the Rhine River, stored in a *dictionary* named `station_names`. The city names are not used and only the coordinates are appended with `line.AddPoint(X, Y)`. As before, a field is created to give the river a name. The actual line feature is again created as a child of the layer with `line_feature = ogr.Feature(lyr.GetLayerDefn())`. Running this code block produces a line that approximately follows the Rhine River between France and Germany." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "shp_dir = r\"\" + os.path.abspath(\"\") + \"/geodata/shapefiles/rhine_proxy.shp\"\n", "rhine_line = create_shp(shp_dir, layer_name=\"basemap\", layer_type=\"line\")\n", "\n", "# create .prj file for the shapefile for web application references\n", "with open(shp_dir.split(\".shp\")[0] + \".prj\", \"w+\") as prj:\n", " prj.write(get_wkt(3857))\n", "\n", "# get basemap layer\n", "lyr = rhine_line.GetLayer()\n", "\n", "# coordinates for EPSG:3857 WG84 / Pseudo-Mercator\n", "station_names = {\"Basel\": (844361.68, 6035047.42),\n", " \"Kembs\": (835724.27, 6056449.76),\n", " \"Breisach\": (842565.32, 6111140.43),\n", " \"Rhinau\": (857547.04, 6158569.58),\n", " \"Strasbourg\": (868439.31, 6203189.68)}\n", "\n", "# create line object and add points from station names\n", "line = ogr.Geometry(ogr.wkbLineString)\n", "for stn in station_names.values():\n", " line.AddPoint(stn[0], stn[1])\n", "\n", "# create field named \"river\"\n", "field_name = ogr.FieldDefn(\"river\", ogr.OFTString)\n", "lyr.CreateField(field_name)\n", "\n", "# create feature, geometry, and field entry\n", "line_feature = ogr.Feature(lyr.GetLayerDefn())\n", "line_feature.SetGeometry(line)\n", "line_feature.SetField(\"river\", \"Rhine\")\n", "\n", "# add feature to layer\n", "lyr.CreateFeature(line_feature)\n", "\n", "lyr = None\n", "rhine_line = None" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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](http://naciscdn.org/naturalearth/packages/Natural_Earth_quick_start.zip).\n", "\n", "```{figure} ../img/qgis-rhine.png\n", ":alt: rhine line shapefile\n", ":name: qgis-rhine-py\n", "\n", "The line approximating the Rhine river in QGIS.\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "(poly-create)=\n", "## Polygon Shapefile\n", "\n", "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:\n", "\n", "* The projection is `EPSG:4326`.\n", "* The point coordinates are generated within an `ogr.wkbLinearRing` object step-by-step rather than in a loop over *dictionary* entries.\n", "* File, variable, and field names." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "shp_dir = r\"\" + os.path.abspath(\"\") + \"/geodata/shapefiles/va4wasserbau.shp\"\n", "va_geo = create_shp(shp_dir, layer_name=\"basemap\", layer_type=\"polygon\")\n", "\n", "# create .prj file for the shapefile for GIS map applications\n", "with open(shp_dir.split(\".shp\")[0] + \".prj\", \"w+\") as prj:\n", " prj.write(get_wkt(4326))\n", "\n", "# get basemap layer\n", "lyr = va_geo.GetLayer()\n", "\n", "# create polygon points\n", "pts = ogr.Geometry(ogr.wkbLinearRing)\n", "pts.AddPoint(9.103686, 48.744251)\n", "pts.AddPoint(9.104689, 48.744198)\n", "pts.AddPoint(9.104667, 48.743960)\n", "pts.AddPoint(9.103557, 48.744009)\n", "\n", "# create polygon geometry\n", "poly = ogr.Geometry(ogr.wkbPolygon)\n", "# build polygon geometry from points\n", "poly.AddGeometry(pts)\n", "\n", "# add field to classify building type\n", "field = ogr.FieldDefn(\"building\", ogr.OFTString)\n", "lyr.CreateField(field)\n", "\n", "poly_feature_defn = lyr.GetLayerDefn()\n", "poly_feature = ogr.Feature(poly_feature_defn)\n", "poly_feature.SetGeometry(poly)\n", "poly_feature.SetField(\"building\", \"Versuchsanstalt\")\n", "\n", "lyr.CreateFeature(poly_feature)\n", "\n", "lyr = None\n", "va_geo = None" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "(json-create)=\n", "## Build Shapefile from JSON\n", "\n", "Loading geometry data from in-line defined variables is cumbersome in practice, where geospatial data are often provided on public platforms (e.g., land use or cover). The following example uses a JSON file generated with map service data from the [Baden-Württemberg State Institute for the Environment, Survey and Nature Conservation (LUBW)](https://udo.lubw.baden-wuerttemberg.de/), where polygon nodes are stored in WKB polygon geometry format (`\"MultiPolygon (((node1_x node1_y, nodej_x, nodej_y, ... ...)))\"`):\n", "\n", "* The JSON file ([download hq100-dreisam.json](https://raw.githubusercontent.com/hydro-informatics/jupyter-python-course/main/geodata/json/hq100-dreisam.json) and save it to a subdirectory called `/geodata/json/`)is read with {ref}`pandas` and the shapefile is created, as before, with the `create_shp` function.\n", "* The projection is `EPSG:25832`.\n", "* Two fields are added in the form of\n", " - `\"tbg_name\"` (the original string name of the polygons in the LUBW data), and\n", " - `\"area\"` (a real number field, in which the polygon area is calculated in m2 using `polygon.GetArea()`).\n", "* The polygon geometries are derived from the WKB-formatted definitions in the `\"wkb_geom\"` field of the pandas dataframe object `dreisam_inundation`." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "# get data from json file\n", "dreisam_inundation = pd.read_json(r\"\" + os.path.abspath(\"\") + \"/geodata/json/hq100-dreisam.json\")\n", "\n", "# create shapefile\n", "shp_dir = r\"\" + os.path.abspath('') + \"/geodata/shapefiles/dreisam_hq100.shp\"\n", "dreisam_hq100 = create_shp(shp_dir, layer_name=\"basemap\", layer_type=\"polygon\")\n", "\n", "# create .prj file for the shapefile for GIS map applications\n", "with open(shp_dir.split(\".shp\")[0] + \".prj\", \"w+\") as prj:\n", " prj.write(get_wkt(25832))\n", "\n", "# get basemap layer\n", "lyr = dreisam_hq100.GetLayer()\n", "\n", "# add string field \"tbg_name\"\n", "lyr.CreateField(ogr.FieldDefn(\"tbg_name\", ogr.OFTString))\n", "\n", "# add string field \"area\"\n", "lyr.CreateField(ogr.FieldDefn(\"area\", ogr.OFTReal))\n", "\n", "for wkt, tbg in zip(dreisam_inundation[\"wkt_geom\"], dreisam_inundation[\"TBG_NAME\"]):\n", " # create Feature as child of the layer\n", " feature = ogr.Feature(lyr.GetLayerDefn())\n", " # assign tbg_name\n", " feature.SetField(\"tbg_name\", tbg)\n", " # use WKT format to add a polygon geometry to the Feature\n", " polygon = ogr.CreateGeometryFromWkt(wkt)\n", " # define default value of 0 to the area field\n", " feature.SetField(\"area\", polygon.GetArea())\n", "\n", " feature.SetGeometry(polygon)\n", " # append the new feature to the basement layer\n", " lyr.CreateFeature(feature)\n", "\n", "lyr = None\n", "dreisam_hq100 = None" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "```{tip}\n", "Open the new `dreisam_hq100.shp` in QGIS and explore the attribute table.\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Also *GeoJSON* data can be used to create an `ogr.Geometry` with `ogr.createFromGeoJson(FILENAME)`:" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "X=1013452, Y=6231540 (EPSG:3857)\n" ] } ], "source": [ "from osgeo import ogr\n", "geojson_data = \"\"\"{\"type\":\"Point\",\"coordinates\":[1013452.282805,6231540.674235]}\"\"\"\n", "point = ogr.CreateGeometryFromJson(geojson_data)\n", "print(\"X=%d, Y=%d (EPSG:3857)\" % (point.GetX(), point.GetY()))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "(calc)=\n", "## Calculate Geometric Attributes\n", "\n", "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: \n", "\n", "* Unify multiple polygons
\n", " `wkt_... = ...`
\n", " `polygon_a = ogr.CreateGeometryFromWkt(wkt_1)`
\n", " `polygon_b = ogr.CreateGeometryFromWkt(wkt_2)`
\n", " `polygon_union = polygon_a.Union(polygon_b)`\n", "* Intersect two polygons
\n", " `polygon_intersection = polygon_a.Intersection(polygon_b)`\n", "* Envelope (minimum and maximum extents) of a polygon
\n", " `env = polygon.GetEnvelope()`
\n", " `print(\"minX: %d, minY: %d, maxX: %d, maxY: %d\" % (env[0],env[2],env[1],env[3])`\n", "* Convex hull (envelope surface) of multiple geometries (points, lines, polygons)
\n", " `all_polygons = ogr.Geometry(ogr.wkbGeometryCollection)`
\n", " `for feature in POLYGON-SOURCE-LAYER: all_polygons.AddGeometry(feature.GetGeometryRef())`
\n", " `convexhull = all_polygons.ConvexHull()`
\n", " Save `convexhull` to shapefile (use the `create_shp` function as shown in the above examples or read more at [pcjerick's Github pages](https://pcjericks.github.io/py-gdalogr-cookbook/vector_layers.html#save-the-convex-hull-of-all-geometry-from-an-input-layer-to-an-output-layer))
\n", " Tip: To create a tight hull (e.g., of a point cloud), look for `concavehull` functions.\n", "* Length (of a line)
\n", " `wkt = \"LINESTRING (415128.5 5320979.5, 415128.6 5320974.5, 415129.75 5320974.7)\"`
\n", " `line = ogr.CreateGeometryFromWkt(wkt)`
\n", " `print(\"Length = %d\" % line.Length())`\n", "* Area (of a polygon): `polygon.GetArea()` (see above example)\n", "* Example to calculate [centroid coordinates of polygons](https://pcjericks.github.io/py-gdalogr-cookbook/geometry.html#quarter-polygon-and-create-centroids).\n", "\n", "```{important}\n", "The units of the geometric attribute (e.g., m$^2$, U.S. feet, or others) are calculated based on the definitions in the `.prj file ` (recall the definition of projections in WKT format in the {ref}`prj` section).\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "(export)=\n", "## Export to Other Formats\n", "\n", "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*](https://github.com/mdiener21/python-geospatial-analysis-cookbook)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### GeoJSON\n", "GeoJSON files can be easily created as before, even without activating a driver:" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "triangle = ogr.Geometry(ogr.wkbLinearRing)\n", "triangle.AddPoint(-11717151.498691, 2356192.894805)\n", "triangle.AddPoint(-11717120.446149, 2355586.175893)\n", "triangle.AddPoint(-11719392.059083, 2354012.050842)\n", "\n", "polygon = ogr.Geometry(ogr.wkbPolygon)\n", "polygon.AddGeometry(triangle)\n", "\n", "with open(r\"\" + os.path.abspath('') + \"/geodata/geojson/pitillal-triangle.geojson\", \"w+\") as gjson:\n", " gjson.write(polygon.ExportToJson())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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. " ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "gjson_driver = ogr.GetDriverByName(\"GeoJSON\")\n", "\n", "# make spatial reference\n", "sr = osr.SpatialReference()\n", "sr.ImportFromEPSG(3857)\n", "\n", "# create GeoJSON file\n", "gjson = gjson_driver.CreateDataSource(\"pitillal-full.geojson\")\n", "gjson_lyr = gjson.CreateLayer(\"pitillal-full.geojson\", geom_type=ogr.wkbPolygon, srs=sr)\n", "\n", "# get layer feature definitions\n", "feature_def = gjson_lyr.GetLayerDefn()\n", "# create new feature\n", "new_feature = ogr.Feature(feature_def)\n", "# assign the triangle from the above code block\n", "new_feature.SetGeometry(polygon)\n", "# add new feature to Layer\n", "gjson_lyr.CreateFeature(new_feature)\n", "\n", "# close links to data sources\n", "gjson = None\n", "gjson_lyr = None" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### KML (Google Earth)\n", "To display point, line, or polygon features in Google Earth, features can be plugged into Google's [KML](https://developers.google.com/kml/documentation/kml_tut) (Keyhole Markup Language), similar to the creation of a GeoJSON file, and with the function `geometry.ExportToKML`:" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "triangle = ogr.Geometry(ogr.wkbLinearRing)\n", "triangle.AddPoint(-11717151.498691, 2356192.894805)\n", "triangle.AddPoint(-11717120.446149, 2355586.175893)\n", "triangle.AddPoint(-11719392.059083, 2354012.050842)\n", "\n", "polygon = ogr.Geometry(ogr.wkbPolygon)\n", "polygon.AddGeometry(triangle)\n", "\n", "#geojson = poly.ExportToJson()\n", "with open(r\"\" + os.path.abspath(\"\") + \"/geodata/pitillal-triangle.kml\", \"w+\") as gjson:\n", " gjson.write(polygon.ExportToKML())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Moreover, similar to GeoJSON files and shapefiles, KML files can be generated more robustly (with a defined projection, for example). All you need to do is initiating the KML driver (`kml_driver = ogr.GetDriverByName(\"KML\")`) and define a KML data source (`kml_file = kml_driver.CreateDataSource(FILENAME.KML)`)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "```{admonition} Exercise\n", "Familiarize with shapefile handling in the {ref}`geospatial ecohydraulics ` exercise.\n", "```" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.12" } }, "nbformat": 4, "nbformat_minor": 4 }