diff --git a/05-raster-vector.ipynb b/05-raster-vector.ipynb deleted file mode 100644 index b055f3c5..00000000 --- a/05-raster-vector.ipynb +++ /dev/null @@ -1,1875 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Raster-vector interactions {#sec-raster-vector}\n", - "\n", - "## Prerequisites\n" - ], - "id": "7fd26792" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "#| echo: false\n", - "import pandas as pd\n", - "import matplotlib.pyplot as plt\n", - "pd.options.display.max_rows = 6\n", - "pd.options.display.max_columns = 6\n", - "pd.options.display.max_colwidth = 35\n", - "plt.rcParams['figure.figsize'] = (5, 5)" - ], - "id": "522ade54", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Let's import the required packages:\n" - ], - "id": "6f3a3895" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "import numpy as np\n", - "import shapely\n", - "import matplotlib.pyplot as plt\n", - "import geopandas as gpd\n", - "import rasterio\n", - "import rasterio.mask\n", - "import rasterstats\n", - "import rasterio.plot\n", - "import rasterio.features\n", - "import math\n", - "import os" - ], - "id": "1d5dd0d4", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "and load the sample data:\n" - ], - "id": "6b738b3d" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "src_srtm = rasterio.open('data/srtm.tif')\n", - "src_nlcd = rasterio.open('data/nlcd.tif')\n", - "src_grain = rasterio.open('data/grain.tif')\n", - "src_elev = rasterio.open('data/elev.tif')\n", - "src_dem = rasterio.open('data/dem.tif')\n", - "zion = gpd.read_file('data/zion.gpkg')\n", - "zion_points = gpd.read_file('data/zion_points.gpkg')\n", - "cycle_hire_osm = gpd.read_file('data/cycle_hire_osm.gpkg')\n", - "us_states = gpd.read_file('data/us_states.gpkg')\n", - "nz = gpd.read_file('data/nz.gpkg')\n", - "src_nz_elev = rasterio.open('data/nz_elev.tif')" - ], - "id": "919d20e5", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Introduction\n", - "\n", - "This Chapter focuses on interactions between raster and vector geographic data models, both introduced in @sec-spatial-class.\n", - "It includes four main techniques: \n", - "\n", - "* raster cropping and masking using vector objects (@sec-raster-cropping), \n", - "* extracting raster values using different types of vector data (Section @sec-raster-extraction), and \n", - "* raster-vector conversion (@sec-rasterization and @sec-spatial-vectorization).\n", - "\n", - "These concepts are demonstrated using data from in previous chapters, to understand their potential real-world applications.\n", - "\n", - "## Raster masking and cropping {#sec-raster-cropping}\n", - "\n", - "Many geographic data projects involve integrating data from many different sources, such as remote sensing images (rasters) and administrative boundaries (vectors). \n", - "Often the extent of input raster datasets is larger than the area of interest. \n", - "In this case raster *masking*, *cropping*, or both, are useful for unifying the spatial extent of input data (@fig-raster-crop (b) and (c), and the following two examples, illustrate the difference between masking and cropping). \n", - "Both operations reduce object memory use and associated computational resources for subsequent analysis steps, and may be a necessary preprocessing step before creating attractive maps involving raster data.\n", - "\n", - "We will use two layers to illustrate raster cropping:\n", - "\n", - "* The `srtm.tif` raster representing elevation (meters above sea level) in south-western Utah (a **rasterio** file connection named `src_srtm`) (see @fig-raster-crop (a))\n", - "* The `zion.gpkg` vector layer representing the Zion National Park boundaries (a `GeoDataFrame` named `zion`)\n", - "\n", - "Both target and cropping objects must have the same projection. \n", - "Since it is easier and more precise to reproject vector layers, compared to rasters, we use the following expression to reproject the vector layer `zion` into the CRS of the raster `src_srtm`:\n" - ], - "id": "f77993e9" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "zion = zion.to_crs(src_srtm.crs)" - ], - "id": "b22a06c0", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "To mask the image, i.e., convert all pixels which do not intersect with the `zion` polygon to \"No Data\", we use the [`rasterio.mask.mask`](https://rasterio.readthedocs.io/en/stable/api/rasterio.mask.html#rasterio.mask.mask) function, as follows:\n" - ], - "id": "1941d3b0" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "out_image_mask, out_transform_mask = rasterio.mask.mask(\n", - " src_srtm, \n", - " zion.geometry, \n", - " crop=False, \n", - " nodata=9999\n", - ")" - ], - "id": "056044a4", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Note that we need to choose and specify a \"No Data\" value, within the valid range according to the data type. \n", - "Since `srtm.tif` is of type `uint16` (how can we check?), we choose `9999` (a positive integer that is guaranteed not to occur in the raster). \n", - "Also note that **rasterio** does not directly support **geopandas** data structures, so we need to pass a \"collection\" of **shapely** geometries: a `GeoSeries` (see below) or a `list` of **shapely** geometries (see next example) both work.\n", - "\n", - "The result is the `out_image` array with the masked values: \n" - ], - "id": "aa8fddce" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "out_image_mask" - ], - "id": "1b63c2cb", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "and the new transformation matrix `out_transform`:\n" - ], - "id": "c8f8e963" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "out_transform_mask" - ], - "id": "6f6a363b", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Note that masking (without cropping!) does not modify the raster spatial configuration. \n", - "Therefore, the new transform is identical to the original:\n" - ], - "id": "982cbba7" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "src_srtm.transform" - ], - "id": "1e414d34", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Unfortunately, the `out_image` and `out_transform` object do not contain any information indicating that `9999` represents \"No Data\". \n", - "To associate the information with the raster, we must write it to file along with the corresponding metadata. \n", - "For example, to write the cropped raster to file, we need to modify the \"No Data\" setting in the metadata:\n" - ], - "id": "0249dff0" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "dst_kwargs = src_srtm.meta\n", - "dst_kwargs.update(nodata=9999)\n", - "dst_kwargs" - ], - "id": "5dc5d0cf", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Then we can write the cropped raster to file:\n" - ], - "id": "cb239291" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "new_dataset = rasterio.open('output/srtm_masked.tif', 'w', **dst_kwargs)\n", - "new_dataset.write(out_image_mask)\n", - "new_dataset.close()" - ], - "id": "c02ade45", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Now we can re-import the raster:\n" - ], - "id": "1e86778e" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "src_srtm_mask = rasterio.open('output/srtm_masked.tif')" - ], - "id": "e53c66b9", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The `.meta` property contains the `nodata` entry. \n", - "Now, any relevant operation (such as plotting) will take \"No Data\" into account:\n" - ], - "id": "690e1046" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "src_srtm_mask.meta" - ], - "id": "a4e1e955", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "*Cropping* means reducing the raster extent to the extent of the vector layer:\n", - "\n", - "* To just crop, *without* masking, we can derive the bounding box polygon of the vector layer, and then crop using that polygon (see further below), also combined with `crop=True`\n", - "* To crop *and* mask, we can use `rasterio.mask.mask` same shown as above for masking, just setting `crop=True` instead of the default `crop=False` (see further below)\n", - "\n", - "For the example of cropping only, here is how we can obtain the extent polygon of `zion`, as a `shapely` geometry object (@fig-zion-bbox):\n" - ], - "id": "f90ff908" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "#| label: fig-zion-bbox\n", - "#| fig-cap: Bounding box `'Polygon'` geometry of the `zion` layer\n", - "\n", - "bb = zion.unary_union.envelope\n", - "bb" - ], - "id": "fig-zion-bbox", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The extent can now be used for masking. \n", - "Here, we are also using the `all_touched=True` option so that pixels partially overlapping with the extent are included:\n" - ], - "id": "d52f328a" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "out_image_crop, out_transform_crop = rasterio.mask.mask(\n", - " src_srtm, \n", - " [bb], \n", - " crop=True, \n", - " all_touched=True, \n", - " nodata=9999\n", - ")" - ], - "id": "e61505dd", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Finally, here is how we crop and mask, using `rasterio.mask.mask` with `crop=True`. \n", - "When writing to file, it is now crucual to update the transform and dimensions, since they were modified as a result of cropping. \n", - "Also note that `out_image_mask_crop` is a three-dimensional array (even tough it has one band in this case), so the number of rows and columns are in `.shape[1]` and `.shape[2]` (rather than `.shape[0]` and `.shape[1]`), respectively:\n" - ], - "id": "7426d6b4" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "out_image_mask_crop, out_transform_mask_crop = rasterio.mask.mask(\n", - " src_srtm, \n", - " zion.geometry, \n", - " crop=True, \n", - " nodata=9999\n", - ")\n", - "dst_kwargs = src_srtm.meta\n", - "dst_kwargs.update({\n", - " 'nodata': 9999,\n", - " 'transform': out_transform_mask_crop,\n", - " 'width': out_image_mask_crop.shape[2],\n", - " 'height': out_image_mask_crop.shape[1]\n", - "})\n", - "new_dataset = rasterio.open('output/srtm_masked_cropped.tif', 'w', **dst_kwargs)\n", - "new_dataset.write(out_image_mask_crop)\n", - "new_dataset.close()\n", - "src_srtm_mask_crop = rasterio.open('output/srtm_masked_cropped.tif')" - ], - "id": "0005dd34", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "out_image_mask_crop.shape" - ], - "id": "7f1f9e46", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "@fig-raster-crop shows the original raster, and the cropped and masked results:\n" - ], - "id": "6eb62399" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "#| label: fig-raster-crop\n", - "#| fig-cap: Raster masking and cropping\n", - "#| layout-ncol: 2\n", - "#| fig-subcap:\n", - "#| - Original\n", - "#| - Masked\n", - "#| - Cropped\n", - "#| - Masked+Cropped\n", - "\n", - "# Original\n", - "fig, ax = plt.subplots(figsize=(3.5, 3.5))\n", - "rasterio.plot.show(src_srtm, ax=ax)\n", - "zion.plot(ax=ax, color='none', edgecolor='black');\n", - "\n", - "# Masked\n", - "fig, ax = plt.subplots(figsize=(3.5, 3.5))\n", - "rasterio.plot.show(src_srtm_mask, ax=ax)\n", - "zion.plot(ax=ax, color='none', edgecolor='black');\n", - "\n", - "# Cropped\n", - "fig, ax = plt.subplots(figsize=(3.5, 3.5))\n", - "rasterio.plot.show(out_image_crop, transform=out_transform_crop, ax=ax)\n", - "zion.plot(ax=ax, color='none', edgecolor='black');\n", - "\n", - "# Masked+Cropped\n", - "fig, ax = plt.subplots(figsize=(3.5, 3.5))\n", - "rasterio.plot.show(src_srtm_mask_crop, ax=ax)\n", - "zion.plot(ax=ax, color='none', edgecolor='black');" - ], - "id": "fig-raster-crop", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Raster extraction {#sec-raster-extraction}\n", - "\n", - "Raster extraction is the process of identifying and returning the values associated with a 'target' raster at specific locations, based on a (typically vector) geographic 'selector' object. \n", - "The reverse of raster extraction---assigning raster cell values based on vector objects---is rasterization, described in @sec-rasterization.\n", - "\n", - "In the following examples, we use a third-party package called `rasterstats`, which is specifically aimed at extracting raster values: \n", - "\n", - "* to *points* (@sec-extraction-to-points) or to *lines* (@sec-extraction-to-lines), via the `rasterstats.point_query` function, or \n", - "* to *polygons* (@sec-extraction-to-polygons), via the `rasterstats.zonal_stats` function.\n", - "\n", - "### Extraction to points {#sec-extraction-to-points}\n", - "\n", - "The simplest type of raster extraction is extracting the values of raster cells at specific points. \n", - "To demonstrate extraction to points, we will use `zion_points`, which contains a sample of 30 locations within the Zion National Park (@fig-zion-points): \n" - ], - "id": "9c393db0" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "#| label: fig-zion-points\n", - "#| fig-cap: '30 point locations within the Zion National Park, with elevation in the background'\n", - "\n", - "fig, ax = plt.subplots()\n", - "rasterio.plot.show(src_srtm, ax=ax)\n", - "zion_points.plot(ax=ax, color='black');" - ], - "id": "fig-zion-points", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The following expression extracts elevation values from `srtm.tif` according to `zion_points`, using `rasterstats.point_query`:\n" - ], - "id": "4f102560" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "result = rasterstats.point_query(\n", - " zion_points, \n", - " src_srtm.read(1), \n", - " nodata = src_srtm.nodata, \n", - " affine = src_srtm.transform,\n", - " interpolate='nearest'\n", - ")" - ], - "id": "f9f39edb", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The resulting object is a `list` of raster values, corresponding to `zion_points`. \n", - "For example, here are the elevations of the first five points:\n" - ], - "id": "da4a805d" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "result[:5]" - ], - "id": "c22658ef", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "To get a `GeoDataFrame` with the original points geometries (and other attributes, if any), as well as the extracted raster values, we can assign the extraction result into a new column:\n" - ], - "id": "fd21d413" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "zion_points['elev'] = result\n", - "zion_points" - ], - "id": "8d34af99", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Extraction to lines {#sec-extraction-to-lines}\n", - "\n", - "Raster extraction is also applicable with line selectors. \n", - "The typical line extraction algorithm is to extract one value for each raster cell touched by a line. \n", - "However, this particular approach is not recommended to obtain values along the transects, as it is hard to get the correct distance between each pair of extracted raster values.\n", - "\n", - "For line extraction, a better approach is to split the line into many points (at equal distances along the line) and then extract the values for these points using the \"extraction to points\" technique (@sec-extraction-to-points). \n", - "To demonstrate this, the code below creates (see @sec-vector-data for recap) `zion_transect`, a straight line going from northwest to southeast of the Zion National Park (@fig-zion-transect-shapely):\n" - ], - "id": "8a05ba94" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "#| label: fig-zion-transect-shapely\n", - "#| fig-cap: Transect through the Zion National Park\n", - "\n", - "coords = [[-113.2, 37.45], [-112.9, 37.2]]\n", - "zion_transect = shapely.LineString(coords)\n", - "zion_transect" - ], - "id": "fig-zion-transect-shapely", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Here is a printout demonstrating that this is a `'LineString'` geometry representing a straight line between two points:\n" - ], - "id": "62eb49ad" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "print(zion_transect)" - ], - "id": "04507a5d", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The line is also illustrated in the context of the raster in @fig-zion-transect.\n", - "\n", - "The utility of extracting heights from a linear selector is illustrated by imagining that you are planning a hike. \n", - "The method demonstrated below provides an 'elevation profile' of the route (the line does not need to be straight), useful for estimating how long it will take due to long climbs.\n", - "\n", - "First, we need to create a layer `zion_transect_pnt` consisting of points along our line (`zion_transect`), at specified intervals (e.g., `250`). \n", - "To do that, we need to transform the line into a projected CRS (so that we work with true distances, in $m$), such as UTM. \n", - "This requires going through a `GeoSeries`, as **shapely** geometries have no CRS definition nor concept of reprojection (see @sec-vector-layer-from-scratch):\n" - ], - "id": "cf2b2f10" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "zion_transect_utm = gpd.GeoSeries(zion_transect, crs=4326).to_crs(32612)\n", - "zion_transect_utm = zion_transect_utm.iloc[0]" - ], - "id": "8aac55f3", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The printout of the new geometry shows this is still a straight line between two points, only with coordinates in a projected CRS:\n" - ], - "id": "8d2e609e" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "print(zion_transect_utm)" - ], - "id": "e99c9ded", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Next, we need to calculate the distances, along the line, where points are going to be generated, using [`np.arange`](https://numpy.org/doc/stable/reference/generated/numpy.arange.html). \n", - "This is a numeric sequence starting at `0`, going up to line `.length`, in steps of `250` ($m$):\n" - ], - "id": "4401c407" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "distances = np.arange(0, zion_transect_utm.length, 250)\n", - "distances[:7] ## First 7 distance cutoff points" - ], - "id": "559efe50", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The distances cutoffs are used to sample (\"interpolate\") points along the line. \n", - "The **shapely** [`.interpolate`](https://shapely.readthedocs.io/en/stable/manual.html#object.interpolate) method is used to generate the points. \n", - "The points are then reprojected back to the geographic CRS of the raster (EPSG:`4326`):\n" - ], - "id": "d5df418c" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "zion_transect_pnt = [zion_transect_utm.interpolate(distance) for distance in distances]\n", - "zion_transect_pnt = gpd.GeoSeries(zion_transect_pnt, crs=32612)\n", - "zion_transect_pnt = zion_transect_pnt.to_crs(src_srtm.crs)\n", - "zion_transect_pnt" - ], - "id": "f9b5ee73", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Finally, we extract the elevation values for each point in our transect and combine the information with `zion_transect_pnt` (after \"promoting\" it to a `GeoDataFrame`, to accomodate extra attributes), using the point extraction method shown earlier (@sec-extraction-to-points). \n", - "We also attach the respective distance cutoff points `distances`:\n" - ], - "id": "865b9759" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "result = rasterstats.point_query(\n", - " zion_transect_pnt, \n", - " src_srtm.read(1), \n", - " nodata = src_srtm.nodata, \n", - " affine = src_srtm.transform,\n", - " interpolate='nearest'\n", - ")\n", - "zion_transect_pnt = gpd.GeoDataFrame(geometry=zion_transect_pnt)\n", - "zion_transect_pnt['dist'] = distances\n", - "zion_transect_pnt['elev'] = result\n", - "zion_transect_pnt" - ], - "id": "bf58c88b", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The information in `zion_transect_pnt`, namely the `'dist'` and `'elev'` attributes, can now be used to draw an elevation profile, as illustrated in @fig-zion-transect (b):\n" - ], - "id": "7af186d9" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "#| label: fig-zion-transect\n", - "#| fig-cap: Extracting a raster values profile to line\n", - "#| layout-ncol: 2\n", - "#| fig-subcap:\n", - "#| - Raster and a line transect\n", - "#| - Extracted elevation profile\n", - "\n", - "# Raster and a line transect\n", - "fig, ax = plt.subplots()\n", - "rasterio.plot.show(src_srtm, ax=ax)\n", - "gpd.GeoSeries(zion_transect).plot(ax=ax, color='black')\n", - "zion.plot(ax=ax, color='none', edgecolor='darkgrey');\n", - "\n", - "# Elevation profile\n", - "fig, ax = plt.subplots()\n", - "zion_transect_pnt.set_index('dist')['elev'].plot(ax=ax)\n", - "ax.set_xlabel('Distance (m)')\n", - "ax.set_ylabel('Elevation (m)');" - ], - "id": "fig-zion-transect", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Extraction to polygons\n", - "\n", - "The final type of geographic vector object for raster extraction is polygons. \n", - "Like lines, polygons tend to return many raster values per polygon. \n", - "For continuous rasters (@fig-raster-extract-to-polygon (a)), we typically want to generate summary statistics for raster values per polygon, for example to characterize a single region or to compare many regions. \n", - "The generation of raster summary statistics, by polygons, is demonstrated in the code below, which creates a list of summary statistics (in this case a list of length 1, since there is just one polygon), using `rasterstats.zonal_stats`:\n" - ], - "id": "cfa24ad9" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "result = rasterstats.zonal_stats(\n", - " zion, \n", - " src_srtm.read(1), \n", - " nodata = src_srtm.nodata, \n", - " affine = src_srtm.transform, \n", - " stats = ['mean', 'min', 'max']\n", - ")\n", - "result" - ], - "id": "d42a1890", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Transformation of the `list` to a `DataFrame` (e.g., to attach the derived attributes to the original polygon layer), is straightforward:\n" - ], - "id": "3091b6a3" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "pd.DataFrame(result)" - ], - "id": "c19a6179", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Because there is only one polygon in the example, a `DataFrame` with a single row is returned. \n", - "However, if `zion` was composed of more than one polygon, we would accordingly get more rows in the `DataFrame`.\n", - "The result provides useful summaries, for example that the maximum height in the park is around `2661` $m$ above see level. \n", - "\n", - "Note the `stats` argument, where we determine what type of statistics are calculated per polygon. \n", - "Possible values other than `'mean'`, `'min'`, `'max'` are:\n", - "\n", - "* `'count'`---The number of valid (i.e., excluding \"No Data\") pixels\n", - "* `'nodata'`---The number of pixels with 'No Data\"\n", - "* `'majority'`---The most frequently occurring value\n", - "* `'median'`---The median value\n", - "\n", - "See the [documentation](https://pythonhosted.org/rasterstats/manual.html#statistics) of `rasterstats.zonal_stats` for the complete list. \n", - "Additionally, the `rasterstats.zonal_stats` function accepts user-defined functions for calculating any custom statistics.\n", - "\n", - "To count occurrences of categorical raster values within polygons (@fig-raster-extract-to-polygon (b)), we can use masking (@sec-raster-cropping) combined with `np.unique`, as follows:\n" - ], - "id": "40cf29e8" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "out_image, out_transform = rasterio.mask.mask(\n", - " src_nlcd, \n", - " zion.geometry.to_crs(src_nlcd.crs), \n", - " crop=False, \n", - " nodata=9999\n", - ")\n", - "counts = np.unique(out_image, return_counts=True)" - ], - "id": "4f7fb532", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "According to the result, for example, pixel value `2` (\"Developed\" class) appears in `4205` pixels within the Zion polygon:\n" - ], - "id": "5f5a893b" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "counts" - ], - "id": "5e6e48ea", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The input data for raster to polygon extraction examples in this section are illustrated in @fig-raster-extract-to-polygon:\n" - ], - "id": "f9c04cb6" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "#| label: fig-raster-extract-to-polygon\n", - "#| fig-cap: Sample data used for continuous and categorical raster extraction to a polygon\n", - "#| layout-ncol: 2\n", - "#| fig-subcap:\n", - "#| - Continuous raster\n", - "#| - Categorical raster\n", - "\n", - "# Continuous raster\n", - "fig, ax = plt.subplots()\n", - "rasterio.plot.show(src_srtm, ax=ax)\n", - "zion.plot(ax=ax, color='none', edgecolor='black');\n", - "\n", - "# Categorical raster\n", - "fig, ax = plt.subplots()\n", - "rasterio.plot.show(src_nlcd, ax=ax, cmap='Set3')\n", - "zion.to_crs(src_nlcd.crs).plot(ax=ax, color='none', edgecolor='black');" - ], - "id": "fig-raster-extract-to-polygon", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Rasterization {#sec-rasterization}\n", - "\n", - "### Rasterizing points {#sec-rasterizing-points}\n", - "\n", - "Rasterization is the conversion of vector objects into their representation in raster objects. \n", - "Usually, the output raster is used for quantitative analysis (e.g., analysis of terrain) or modeling. \n", - "As we saw in @sec-spatial-class, the raster data model has some characteristics that make it conducive to certain methods. \n", - "Furthermore, the process of rasterization can help simplify datasets because the resulting values all have the same spatial resolution: rasterization can be seen as a special type of geographic data aggregation.\n", - "\n", - "The `rasterio` package contains the [`rasterio.features.rasterize`](https://rasterio.readthedocs.io/en/stable/api/rasterio.features.html#rasterio.features.rasterize) function for doing this work. \n", - "To make it happen, we need to have the \"template\" grid definition, i.e., the \"template\" raster defining the extent, resolution and CRS of the output, in the `out_shape` (the output dimensions) and `transform` (the transformation matrix) arguments of `rasterio.features.rasterize`. \n", - "In case we have an existing template raster, we simply need to query its `.shape` and `.transform`. \n", - "In case we create a custom template, e.g., covering the vector layer extent with specified resolution, there is some extra work to calculate both of these objects (see next example). \n", - "\n", - "Furthermore, the `rasterio.features.rasterize` function requires the input vector shapes in the form for a generator of `(geom,value)` tuples, where:\n", - "\n", - "* `geom` is the given geometry (**shapely** geometry object)\n", - "* `value` is the value to be \"burned\" into pixels coinciding with the geometry (`int` or `float`)\n", - "\n", - "Again, this will be made clear in the next example.\n", - "\n", - "The geographic resolution of the \"template\" raster has a major impact on the results: if it is too low (cell size is too large), the result may miss the full geographic variability of the vector data; if it is too high, computational times may be excessive. \n", - "There are no simple rules to follow when deciding an appropriate geographic resolution, which is heavily dependent on the intended use of the results. \n", - "Often the target resolution is imposed on the user, for example when the output of rasterization needs to be aligned to the existing raster.\n", - "\n", - "To demonstrate rasterization in action, we will use a template raster that has the same extent and CRS as the input vector data `cycle_hire_osm_projected` (a dataset on cycle hire points in London, illustrated in @fig-rasterize-points (a)) and a spatial resolution of 1000 $m$. \n", - "First, we take the vector layer and tansform it to a projected CRS:\n" - ], - "id": "f32b3f1c" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "cycle_hire_osm_projected = cycle_hire_osm.to_crs(27700)" - ], - "id": "d8d931c9", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Next, we need to calculate the `out_shape` and `transform` of our template raster. \n", - "To calculate the transform, we combine the top-left corner of the `cycle_hire_osm_projected` bounding box with the required resolution (e.g., 1000 $m$):\n" - ], - "id": "5f063752" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "bounds = cycle_hire_osm_projected.total_bounds\n", - "res = 1000\n", - "transform = rasterio.transform.from_origin(\n", - " west=bounds[0], \n", - " north=bounds[3], \n", - " xsize=res, \n", - " ysize=res\n", - ")\n", - "transform" - ], - "id": "eb03a032", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "To calculate the `out_shape`, we divide the x-axis and y-axis extent by the resolution, and take the ceiling of the results:\n" - ], - "id": "f9f1e755" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "rows = math.ceil((bounds[3] - bounds[1]) / res)\n", - "cols = math.ceil((bounds[2] - bounds[0]) / res)\n", - "shape = (rows, cols)\n", - "shape" - ], - "id": "3f08740b", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Now, we can rasterize. \n", - "Rasterization is a very flexible operation: the results depend not only on the nature of the template raster, but also on the type of input vector (e.g., points, polygons), the pixel \"activation\" method, and the function applied when there is more than one match.\n", - "\n", - "To illustrate this flexibility we will try three different approaches to rasterization (@fig-rasterize-points (b)-(d)). \n", - "First, we create a raster representing the presence or absence of cycle hire points (known as presence/absence rasters). \n", - "In this case, we transfer the value of `1` to all pixels where at least one point falls in. \n", - "To transform the point `GeoDataFrame` into a generator of `shapely` geometries and the (fixed) values, we use the following expression:\n" - ], - "id": "8669278a" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "g = ((g, 1) for g in cycle_hire_osm_projected.geometry.to_list())\n", - "g" - ], - "id": "29bc994b", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Then, the rasterizing expression is:\n" - ], - "id": "6dfd2714" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "ch_raster1 = rasterio.features.rasterize(\n", - " shapes=g, \n", - " out_shape=shape, \n", - " transform=transform\n", - ")" - ], - "id": "c68dea95", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The result `ch_raster1` is an `ndarray` with the burned values of `1` where the pixel coincides with at least one point, and `0` in \"unaffected\" pixels:\n" - ], - "id": "a36ce54b" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "ch_raster1" - ], - "id": "d6e3055b", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "To count the number of bike hire stations, we can use the fixed value of `1` combined with the `merge_alg=rasterio.enums.MergeAlg.add`, which means that multiple values burned into the same pixel are *summed*, rather than replaced keeping last (the default):\n" - ], - "id": "443c1dd3" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "g = ((g, 1) for g in cycle_hire_osm_projected.geometry.to_list())\n", - "ch_raster2 = rasterio.features.rasterize(\n", - " shapes=g,\n", - " out_shape=shape,\n", - " transform=transform,\n", - " merge_alg=rasterio.enums.MergeAlg.add\n", - ")" - ], - "id": "db504758", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Here is the resulting array of point counts `ch_raster2`:\n" - ], - "id": "2ffcbfa1" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "ch_raster2" - ], - "id": "da022a43", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The new output, `ch_raster2`, shows the number of cycle hire points in each grid cell. \n", - "The cycle hire locations have different numbers of bicycles described by the capacity variable, raising the question, what's the capacity in each grid cell? \n", - "To calculate that, we must sum the field (`'capacity'`) rather than the fixed values of `1`. \n", - "This requires using an expanded generator of geometries and values, where we (1) extract both geometries and attribute values, and (2) filter out \"No Data\" values, as follows:\n" - ], - "id": "ecf77a98" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "g = ((g, v) for g, v in cycle_hire_osm_projected[['geometry', 'capacity']] \\\n", - " .dropna(subset='capacity')\n", - " .to_numpy() \\\n", - " .tolist())\n", - "ch_raster3 = rasterio.features.rasterize(\n", - " shapes=g,\n", - " out_shape=shape,\n", - " transform=transform,\n", - " merge_alg=rasterio.enums.MergeAlg.add\n", - ")" - ], - "id": "8ec1f073", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Here is the resulting array of capacity sums `ch_raster3`:\n" - ], - "id": "a48ca47f" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "ch_raster3" - ], - "id": "719a8a96", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The input point layer `cycle_hire_osm_projected` and the three variants of rasterizing it `ch_raster1`, `ch_raster2`, and `ch_raster3` are shown in @fig-rasterize-points:\n" - ], - "id": "0ef19a6d" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "#| label: fig-rasterize-points\n", - "#| fig-cap: Three variants of point rasterization\n", - "#| layout-ncol: 2\n", - "#| fig-subcap:\n", - "#| - Input points\n", - "#| - Presence/Absence\n", - "#| - Point counts\n", - "#| - Summed attribute values\n", - "\n", - "# Input points\n", - "fig, ax = plt.subplots()\n", - "cycle_hire_osm_projected.plot(column='capacity', ax=ax);\n", - "\n", - "# Presence/Absence\n", - "fig, ax = plt.subplots()\n", - "rasterio.plot.show(ch_raster1, transform=transform, ax=ax);\n", - "\n", - "# Point counts\n", - "fig, ax = plt.subplots()\n", - "rasterio.plot.show(ch_raster2, transform=transform, ax=ax);\n", - "\n", - "# Summed attribute values\n", - "fig, ax = plt.subplots()\n", - "rasterio.plot.show(ch_raster3, transform=transform, ax=ax);" - ], - "id": "fig-rasterize-points", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Rasterizing lines and polygons\n", - "\n", - "Another dataset based on California's polygons and borders (created below) illustrates rasterization of lines. \n", - "There are three preliminary steps. First, we subset the California polygon:\n" - ], - "id": "9f49067b" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "california = us_states[us_states['NAME'] == 'California']\n", - "california" - ], - "id": "14ed5b76", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Second, we \"cast\" the polygon into a `'MultiLineString'` geometry, using the [`.boundary`](https://geopandas.org/en/stable/docs/reference/api/geopandas.GeoSeries.boundary.html) property that `GeoSeries` have: \n" - ], - "id": "82f262ba" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "california_borders = california.geometry.boundary\n", - "california_borders" - ], - "id": "31decabf", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Third, we create a template raster with a resolution of a `0.5` degree, using the same approach as in @sec-rasterizing-points:\n" - ], - "id": "ecc2da1e" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "bounds = california_borders.total_bounds\n", - "res = 0.5\n", - "transform = rasterio.transform.from_origin(\n", - " west=bounds[0], \n", - " north=bounds[3], \n", - " xsize=res, \n", - " ysize=res\n", - ")\n", - "rows = math.ceil((bounds[3] - bounds[1]) / res)\n", - "cols = math.ceil((bounds[2] - bounds[0]) / res)\n", - "shape = (rows, cols)\n", - "shape" - ], - "id": "76547b03", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Finally, we rasterize `california_borders` based on the calculated template's `shape` and `transform`. \n", - "When considering line or polygon rasterization, one useful additional argument is `all_touched`. \n", - "By default it is `False`, but when changed to `True`---all cells that are touched by a line or polygon border get a value. \n", - "Line rasterization with `all_touched=True` is demonstrated in the code below (@fig-rasterize-lines-polygons, left). \n", - "We are also using `fill=np.nan` to set \"background\" values as \"No Data\":\n" - ], - "id": "b5ad8982" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "california_raster1 = rasterio.features.rasterize(\n", - " ((g, 1) for g in california_borders.to_list()),\n", - " out_shape=shape,\n", - " transform=transform,\n", - " all_touched=True,\n", - " fill=np.nan\n", - ")" - ], - "id": "ff26576a", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Compare it to a polygon rasterization, with `all_touched=False` (the default), which selects only raster cells whose centroids are inside the selector polygon, as illustrated in @fig-rasterize-lines-polygons (right):\n" - ], - "id": "7156eae4" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "california_raster2 = rasterio.features.rasterize(\n", - " ((g, 1) for g in california.geometry.to_list()),\n", - " out_shape=shape,\n", - " transform=transform,\n", - " fill=np.nan\n", - ")" - ], - "id": "76e33b42", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "To illustrate which raster pixels are actually selected as part of rasterization, we also show them as points. \n", - "This also requires the following code section to calculate the points, which we explain in @sec-spatial-vectorization:\n" - ], - "id": "97da3262" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "ids = california_raster1.copy()\n", - "ids = np.arange(0, california_raster1.size) \\\n", - " .reshape(california_raster1.shape) \\\n", - " .astype(np.int32)\n", - "shapes = rasterio.features.shapes(ids, transform=transform)\n", - "pol = list(shapes)\n", - "pnt = [shapely.geometry.shape(i[0]).centroid for i in pol]\n", - "pnt = gpd.GeoSeries(pnt)" - ], - "id": "b5861134", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The input vector layer, the rasterization results, and the points `pnt`, are shown in @fig-rasterize-lines-polygons:\n" - ], - "id": "78cf7748" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "#| label: fig-rasterize-lines-polygons\n", - "#| fig-cap: Examples of line and polygon rasterization\n", - "#| layout-ncol: 2\n", - "#| fig-subcap:\n", - "#| - Line rasterization w/ `all_touched=True`\n", - "#| - Polygon rasterization w/ `all_touched=False`\n", - "\n", - "# Line rasterization\n", - "fig, ax = plt.subplots()\n", - "rasterio.plot.show(california_raster1, transform=transform, ax=ax, cmap='Set3')\n", - "gpd.GeoSeries(california_borders).plot(ax=ax, edgecolor='darkgrey', linewidth=1)\n", - "pnt.plot(ax=ax, color='black', markersize=1);\n", - "\n", - "# Polygon rasterization\n", - "fig, ax = plt.subplots()\n", - "rasterio.plot.show(california_raster2, transform=transform, ax=ax, cmap='Set3')\n", - "california.plot(ax=ax, color='none', edgecolor='darkgrey', linewidth=1)\n", - "pnt.plot(ax=ax, color='black', markersize=1);" - ], - "id": "fig-rasterize-lines-polygons", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Spatial vectorization {#sec-spatial-vectorization}\n", - "\n", - "Spatial vectorization is the counterpart of rasterization (@sec-rasterization), but in the opposite direction. \n", - "It involves converting spatially continuous raster data into spatially discrete vector data such as points, lines or polygons.\n", - "\n", - "There are three standard methods to convert a raster to a vector layer, which we cover next:\n", - "\n", - "* Raster to polygons (@sec-raster-to-polygons)\n", - "* Raster to points (@sec-raster-to-points)\n", - "* Raster to contours (@sec-raster-to-contours)\n", - "\n", - "The most straightforward form of vectorization is the first one, converting raster cells to polygons, where each pixel is represented by a rectangular polygon. The second method, raster to points, has the additional step of calculating polygon centroids. The third method, raster to contours, is somewhat unrelated. Let us demonstrate the three in the given order.\n", - "\n", - "### Raster to polygons {#sec-raster-to-polygons}\n", - "\n", - "The `rasterio.features.shapes` function can be used to access to the raster pixel as polygon geometries, as well as raster values. The returned object is a generator, which yields `geometry,value` pairs. The additional `transform` argument is used to yield true spatial coordinates of the polygons, which is usually what we want. \n", - "\n", - "For example, the following expression returns a generator named `shapes`, referring to the pixel polygons:\n" - ], - "id": "2f84e2a0" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "shapes = rasterio.features.shapes(\n", - " src_grain.read(), \n", - " transform=src_grain.transform\n", - ")\n", - "shapes" - ], - "id": "b0875b16", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "We can generate all shapes at once, into a `list` named `pol`, as follows:\n" - ], - "id": "22ce150e" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "pol = list(shapes)" - ], - "id": "502903ea", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Each element in `pol` is a `tuple` of length 2, containing:\n", - "\n", - "* The GeoJSON-like `dict` representing the polygon geometry\n", - "* The value of the pixel(s) which comprise the polygon\n", - "\n", - "For example:\n" - ], - "id": "7611e716" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "pol[0]" - ], - "id": "1d6400a0", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Note that each raster cell is converted into a polygon consisting of five coordinates, all of which are stored in memory (explaining why rasters are often fast compared with vectors!).\n", - "\n", - "To transform the `list` into a `GeoDataFrame`, we need few more steps of data reshaping:\n" - ], - "id": "ea863cd3" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "# Create 'GeoSeries' with the polygons\n", - "geom = [shapely.geometry.shape(i[0]) for i in pol]\n", - "geom = gpd.GeoSeries(geom, crs=src_grain.crs)\n", - "# Create 'Series' with the values\n", - "values = [i[1] for i in pol]\n", - "values = pd.Series(values)\n", - "# Combine the 'Series' and 'GeoSeries' into a 'DataFrame'\n", - "result = gpd.GeoDataFrame({'value': values, 'geometry': geom})\n", - "result" - ], - "id": "77dd4039", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The resulting polygon layer is shown in @fig-raster-to-polygons. As shown using the `edgecolor='black'` option, neighboring pixels sharing the same raster value are dissolved into larger polygons. The `rasterio.features.shapes` function does not offer a way to avoid this type of dissolving. One way to work around that is to convert an array with consecutive IDs, instead of the real values, to polygons, then extract the real values from the raster (similarly to the \"raster to points\" example, see below).\n" - ], - "id": "ad849c1a" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "#| label: fig-raster-to-polygons\n", - "#| fig-cap: '`grain.tif` converted to a polygon layer'\n", - "\n", - "result.plot(column='value', edgecolor='black', legend=True);" - ], - "id": "fig-raster-to-polygons", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Raster to points {#sec-raster-to-points}\n", - "\n", - "To transform raster to points, we can use `rasterio.features.shapes`, as in conversion to polygons, only with the addition of the `.centroid` method to go from polygons to their centroids. However, to avoid dissolving nearby pixels, we will actually convert a raster with consecutive IDs, then extract the \"true\" values by point (it is not strictly necessary in this example, since the values of `elev.tif` are all unique):\n" - ], - "id": "c74416a5" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "# Prepare IDs array\n", - "r = src_elev.read(1)\n", - "ids = r.copy()\n", - "ids = np.arange(0, r.size).reshape(r.shape).astype(np.int32)\n", - "ids" - ], - "id": "1e699c5d", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "# IDs raster to points\n", - "shapes = rasterio.features.shapes(ids, transform=src_elev.transform)\n", - "pol = list(shapes)\n", - "geom = [shapely.geometry.shape(i[0]).centroid for i in pol]\n", - "geom = gpd.GeoSeries(geom, crs=src_elev.crs)\n", - "result = gpd.GeoDataFrame(geometry=geom)" - ], - "id": "a6fda307", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "# Extract values to points\n", - "result['value'] = rasterstats.point_query(\n", - " result, \n", - " r, \n", - " nodata = src_elev.nodata, \n", - " affine = src_elev.transform,\n", - " interpolate='nearest'\n", - ")" - ], - "id": "0858e6a3", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The result is shown in @fig-raster-to-points:\n" - ], - "id": "9c1b6d94" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "#| label: fig-raster-to-points\n", - "#| fig-cap: Raster and point representation of `elev.tif`\n", - "#| layout-ncol: 2\n", - "#| fig-subcap:\n", - "#| - Input raster\n", - "#| - Points\n", - "\n", - "# Input raster\n", - "fig, ax = plt.subplots()\n", - "result.plot(column='value', legend=True, ax=ax)\n", - "rasterio.plot.show(src_elev, ax=ax);\n", - "\n", - "# Points\n", - "fig, ax = plt.subplots()\n", - "result.plot(column='value', legend=True, ax=ax)\n", - "rasterio.plot.show(src_elev, cmap='Greys', ax=ax);" - ], - "id": "fig-raster-to-points", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Raster to contours {#sec-raster-to-contours}\n", - "\n", - "Another common type of spatial vectorization is the creation of contour lines representing lines of continuous height or temperatures (isotherms) for example. We will use a real-world digital elevation model (DEM) because the artificial raster elev produces parallel lines (task for the reader: verify this and explain why this happens). Plotting contour lines is straightforward, using the `contour=True` option of `rasterio.plot.show` (@fig-raster-contours1):\n" - ], - "id": "d61226ed" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "#| label: fig-raster-contours1\n", - "#| fig-cap: Displaying raster contours\n", - "\n", - "fig, ax = plt.subplots()\n", - "rasterio.plot.show(src_dem, ax=ax)\n", - "rasterio.plot.show(\n", - " src_dem, \n", - " ax=ax, \n", - " contour=True, \n", - " levels=np.arange(0,1200,50), \n", - " colors='black'\n", - ");" - ], - "id": "fig-raster-contours1", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Unfortunately, `rasterio` does not provide any way of extracting the contour lines in the form of a vector layer, for uses other than plotting. There are two possible workarounds:\n", - "\n", - "1. Using `gdal_contour` on the [command line](https://gdal.org/programs/gdal_contour.html) (see below), or through its Python interface [`osgeo`](https://gis.stackexchange.com/questions/360431/how-can-i-create-contours-from-geotiff-and-python-gdal-rasterio-etc-into-sh)\n", - "1. Writing a custom function to export contour coordinates generated by, e.g., [`matplotlib`](https://www.tutorialspoint.com/how-to-get-coordinates-from-the-contour-in-matplotlib) or [`skimage`](https://gis.stackexchange.com/questions/268331/how-can-i-extract-contours-from-a-raster-with-python)\n", - "\n", - "We hereby demonstrate the first and easiest approach, using `gdal_contour`. Although we deviate from exclusively using the Python language, the benefit of `gdal_contour` is the proven algorithm, customized to spatial data, and with many relevant options. `gdal_contour` (along with other GDAL programs) should already be installed on your system since this is a dependency of `rasterio`. For example, generating 50 $m$ contours of the `dem.tif` file can be done as follows: \n" - ], - "id": "99e7f506" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "#| eval: false\n", - "os.system('gdal_contour -a elev data/dem.tif output/dem_contour.gpkg -i 50.0')" - ], - "id": "1c0cf0a4", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Note that we ran the `gdal_contour` command through `os.system`, in order to remain in the Python environment. You can also run the standalone command in the command line interface you are using, such as the Anaconda Prompt:\n", - "\n", - "```sh\n", - "gdal_contour -a elev data/dem.tif output/dem_contour.gpkg -i 50.0\n", - "```\n", - "\n", - "Like all GDAL programs, `gdal_contour` works with files. Here: \n", - "\n", - "* The input is the `data/dem.tif` file\n", - "* The result is exported to the `output/dem_contour.gpkg` file\n", - "\n", - "To illustrate the result, let's read the result back into the Python environment. Note that the layer contains an arrtibute named `elev` (as specified using `-a elev`) with the contour elevation values:\n" - ], - "id": "4e031c5a" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "contours = gpd.read_file('output/dem_contour.gpkg')\n", - "contours" - ], - "id": "e06a90bd", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Here is a plot of the contour layer in `dem_contour.gpkg` (@fig-raster-contours2):\n" - ], - "id": "bf13030c" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "#| label: fig-raster-contours2\n", - "#| fig-cap: Raster contours calculated with the `gdal_contour` program\n", - "\n", - "fig, ax = plt.subplots()\n", - "rasterio.plot.show(src_dem, ax=ax)\n", - "contours.plot(ax=ax, edgecolor='black');" - ], - "id": "fig-raster-contours2", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Distance to nearest geometry {#sec-distance-to-nearest-geometry}\n", - "\n", - "Calculating a raster of distances to the nearest geometry is an example of a \"global\" raster operation (@sec-global-operations-and-distances). \n", - "To demonstrate it, suppose that we need to calculate a raster representing the distance to the nearest coast in New Zealand. \n", - "This example also wraps many of the concepts introduced in this chapter and in previous chapter, such as raster aggregation, raster conversion to points, and rasterizing points.\n", - "\n", - "For the coastline, we will dissolve the New Zealand administrative division polygon layer and \"extract\" the boundary (which is a `\"MultiLineString\"` geometry):\n" - ], - "id": "44797bec" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "coastline = gpd.GeoSeries(nz.unary_union, crs=nz.crs) \\\n", - " .to_crs(src_nz_elev.crs) \\\n", - " .boundary\n", - "coastline" - ], - "id": "d381c236", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "For a \"template\" raster, we will aggregate the New Zealand DEM, in the `nz_elev.tif` file, to 5 times coarser resolution. \n", - "The code section below follows the aggeregation example in @sec-raster-agg-disagg, then replaces the original (aggregated) values with unique IDs, which is a preliminary step when converting to points, as explained in @sec-raster-to-points. \n", - "Finally, we also replace \"erase\" (i.e., replace with `np.nan`) IDs which were `np.nan` in the aggregated elevation raster, i.e., beyond the land area of New Zealand:\n" - ], - "id": "fdb0731d" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "factor = 0.2\n", - "# Reading aggregated array\n", - "r = src_nz_elev.read(1,\n", - " out_shape=(\n", - " int(src_nz_elev.height * factor),\n", - " int(src_nz_elev.width * factor)\n", - " ),\n", - " resampling=rasterio.enums.Resampling.average\n", - ")\n", - "# Updating the transform\n", - "new_transform = src_nz_elev.transform * src_nz_elev.transform.scale(\n", - " (src_nz_elev.width / r.shape[1]),\n", - " (src_nz_elev.height / r.shape[0])\n", - ")\n", - "# Generating unique IDs per cell\n", - "ids = r.copy()\n", - "ids = np.arange(0, r.size).reshape(r.shape).astype(np.float32)\n", - "# \"Erasing\" irrelevant IDs\n", - "ids[np.isnan(r)] = np.nan\n", - "ids" - ], - "id": "3cec443d", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The result is an array named `ids` with the IDs, and the corresponding `new_transform`, as plotted in @fig-raster-distances1:\n" - ], - "id": "f1b7e5d3" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "#| label: fig-raster-distances1\n", - "#| fig-cap: Template with cell IDs to calculate distance to nearest geometry\n", - "\n", - "fig, ax = plt.subplots()\n", - "rasterio.plot.show(ids, transform=new_transform, ax=ax)\n", - "gpd.GeoSeries(coastline).plot(ax=ax, edgecolor='black');" - ], - "id": "fig-raster-distances1", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "To calculate distances, we must convert each pixel to a vector (point) geometry. \n", - "We use the technique demonstrated in @sec-raster-to-points: \n" - ], - "id": "ce71b2b0" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "shapes = rasterio.features.shapes(ids, transform=new_transform)\n", - "pol = list(shapes)\n", - "pnt = [shapely.geometry.shape(i[0]).centroid for i in pol]" - ], - "id": "094ee0b3", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The result `pnt` is a `list` of `shapely` geometries, representing raster cell centroids (excluding `np.nan` pixels):\n" - ], - "id": "ae131e68" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "print(pnt[0])" - ], - "id": "377015cf", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Next we calculate the correspinding `list` of distances, using the `.distance` method from `shapely`:\n" - ], - "id": "2a6e5358" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "distances = [(i, i.distance(coastline)) for i in pnt]\n", - "distances[0]" - ], - "id": "1058e355", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Finally, we rasterize (see @sec-rasterizing-points) the distances into our raster template:\n" - ], - "id": "1b0096d4" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "image = rasterio.features.rasterize(\n", - " distances,\n", - " out_shape=ids.shape,\n", - " dtype=np.float_,\n", - " transform=new_transform,\n", - " fill=np.nan\n", - ")\n", - "image" - ], - "id": "f72bf7da", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The resulting raster of distances is shown in @fig-raster-distances2:\n" - ], - "id": "8d4a5b3b" - }, - { - "cell_type": "code", - "metadata": {}, - "source": [ - "#| label: fig-raster-distances2\n", - "#| fig-cap: Distance to nearest coastline in New Zealand\n", - "\n", - "fig, ax = plt.subplots()\n", - "rasterio.plot.show(image, transform=new_transform, ax=ax)\n", - "gpd.GeoSeries(coastline).plot(ax=ax, edgecolor='black');" - ], - "id": "fig-raster-distances2", - "execution_count": null, - "outputs": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Exercises\n" - ], - "id": "017dc18b" - } - ], - "metadata": { - "kernelspec": { - "name": "python3", - "language": "python", - "display_name": "Python 3" - } - }, - "nbformat": 4, - "nbformat_minor": 5 -} \ No newline at end of file diff --git a/05-raster-vector.qmd b/05-raster-vector.qmd index b47b4b05..b398d872 100644 --- a/05-raster-vector.qmd +++ b/05-raster-vector.qmd @@ -30,6 +30,7 @@ import rasterio.plot import rasterio.features import math import os +import osgeo ``` and load the sample data: @@ -728,17 +729,22 @@ There are three standard methods to convert a raster to a vector layer, which we * Raster to points (@sec-raster-to-points) * Raster to contours (@sec-raster-to-contours) -The most straightforward form of vectorization is the first one, converting raster cells to polygons, where each pixel is represented by a rectangular polygon. The second method, raster to points, has the additional step of calculating polygon centroids. The third method, raster to contours, is somewhat unrelated. Let us demonstrate the three in the given order. +The most straightforward form of vectorization is the first one, converting raster cells to polygons, where each pixel is represented by a rectangular polygon. +The second method, raster to points, has the additional step of calculating polygon centroids. +The third method, raster to contours, is somewhat unrelated. +Let us demonstrate the three in the given order. ### Raster to polygons {#sec-raster-to-polygons} -The `rasterio.features.shapes` function can be used to access to the raster pixel as polygon geometries, as well as raster values. The returned object is a generator, which yields `geometry,value` pairs. The additional `transform` argument is used to yield true spatial coordinates of the polygons, which is usually what we want. +The [`rasterio.features.shapes`](https://rasterio.readthedocs.io/en/stable/api/rasterio.features.html#rasterio.features.shapes) function can be used to access to the raster pixel as polygon geometries, as well as raster values. +The returned object is a generator, which yields `geometry,value` pairs. +The additional `transform` argument is used to yield true spatial coordinates of the polygons, which is usually what we want. For example, the following expression returns a generator named `shapes`, referring to the pixel polygons: ```{python} shapes = rasterio.features.shapes( - src_grain.read(), + rasterio.band(src_grain, 1), transform=src_grain.transform ) shapes @@ -763,21 +769,24 @@ pol[0] Note that each raster cell is converted into a polygon consisting of five coordinates, all of which are stored in memory (explaining why rasters are often fast compared with vectors!). -To transform the `list` into a `GeoDataFrame`, we need few more steps of data reshaping: +To transform this `list` into a `GeoDataFrame`, we need few more steps of data reshaping: ```{python} -# Create 'GeoSeries' with the polygons +# 'GeoSeries' with the polygons geom = [shapely.geometry.shape(i[0]) for i in pol] geom = gpd.GeoSeries(geom, crs=src_grain.crs) -# Create 'Series' with the values +# 'Series' with the values values = [i[1] for i in pol] values = pd.Series(values) -# Combine the 'Series' and 'GeoSeries' into a 'DataFrame' +# Combine into a 'GeoDataFrame' result = gpd.GeoDataFrame({'value': values, 'geometry': geom}) result ``` -The resulting polygon layer is shown in @fig-raster-to-polygons. As shown using the `edgecolor='black'` option, neighboring pixels sharing the same raster value are dissolved into larger polygons. The `rasterio.features.shapes` function does not offer a way to avoid this type of dissolving. One way to work around that is to convert an array with consecutive IDs, instead of the real values, to polygons, then extract the real values from the raster (similarly to the "raster to points" example, see below). +The resulting polygon layer is shown in @fig-raster-to-polygons. +As shown using the `edgecolor='black'` option, neighboring pixels sharing the same raster value are dissolved into larger polygons. +The `rasterio.features.shapes` function does not offer a way to avoid this type of dissolving. +One way to work around that is to convert an array with consecutive IDs, instead of the real values, to polygons, then extract the real values from the raster (see the example in @sec-raster-to-points). ```{python} #| label: fig-raster-to-polygons @@ -788,27 +797,32 @@ result.plot(column='value', edgecolor='black', legend=True); ### Raster to points {#sec-raster-to-points} -To transform raster to points, we can use `rasterio.features.shapes`, as in conversion to polygons, only with the addition of the `.centroid` method to go from polygons to their centroids. However, to avoid dissolving nearby pixels, we will actually convert a raster with consecutive IDs, then extract the "true" values by point (it is not strictly necessary in this example, since the values of `elev.tif` are all unique): +To transform raster to points, we can use `rasterio.features.shapes`, same as in conversion to polygons (@sec-raster-to-points), only with the addition of the `.centroid` method to go from polygons to their centroids. +However, to avoid dissolving nearby pixels, we will actually convert a raster with consecutive IDs, then extract the "true" values by point (it is not strictly necessary in this example, since the values of `elev.tif` are all unique). + +First, we create an `ndarray` with consecutive IDs, matching the shape of `elev.tif` raster values: ```{python} -# Prepare IDs array r = src_elev.read(1) ids = r.copy() ids = np.arange(0, r.size).reshape(r.shape).astype(np.int32) ids ``` +Next, we use the `rasterio.features.shapes` function to create a point layer with the raster cell IDs: + ```{python} -# IDs raster to points shapes = rasterio.features.shapes(ids, transform=src_elev.transform) pol = list(shapes) geom = [shapely.geometry.shape(i[0]).centroid for i in pol] geom = gpd.GeoSeries(geom, crs=src_elev.crs) result = gpd.GeoDataFrame(geometry=geom) +result ``` +Finally, we extract (@sec-extraction-to-points) the `elev.tif` raster values to points, technically finalizing the raster-to-points conversion: + ```{python} -# Extract values to points result['value'] = rasterstats.point_query( result, r, @@ -818,7 +832,7 @@ result['value'] = rasterstats.point_query( ) ``` -The result is shown in @fig-raster-to-points: +The input raster and the resulting point layer are shown in @fig-raster-to-points: ```{python} #| label: fig-raster-to-points @@ -841,7 +855,9 @@ rasterio.plot.show(src_elev, cmap='Greys', ax=ax); ### Raster to contours {#sec-raster-to-contours} -Another common type of spatial vectorization is the creation of contour lines representing lines of continuous height or temperatures (isotherms) for example. We will use a real-world digital elevation model (DEM) because the artificial raster elev produces parallel lines (task for the reader: verify this and explain why this happens). Plotting contour lines is straightforward, using the `contour=True` option of `rasterio.plot.show` (@fig-raster-contours1): +Another common type of spatial vectorization is the creation of contour lines representing lines of continuous height or temperatures (isotherms) for example. +We will use a real-world digital elevation model (DEM) because the artificial raster `elev.tif` produces parallel lines (task for the reader: verify this and explain why this happens). +Plotting contour lines is straightforward, using the `contour=True` option of `rasterio.plot.show` (@fig-raster-contours1): ```{python} #| label: fig-raster-contours1 @@ -858,45 +874,86 @@ rasterio.plot.show( ); ``` -Unfortunately, `rasterio` does not provide any way of extracting the contour lines in the form of a vector layer, for uses other than plotting. There are two possible workarounds: +Unfortunately, `rasterio` does not provide any way of extracting the contour lines in the form of a vector layer, for uses other than plotting. + +There are two possible workarounds: 1. Using `gdal_contour` on the [command line](https://gdal.org/programs/gdal_contour.html) (see below), or through its Python interface [`osgeo`](https://gis.stackexchange.com/questions/360431/how-can-i-create-contours-from-geotiff-and-python-gdal-rasterio-etc-into-sh) 1. Writing a custom function to export contour coordinates generated by, e.g., [`matplotlib`](https://www.tutorialspoint.com/how-to-get-coordinates-from-the-contour-in-matplotlib) or [`skimage`](https://gis.stackexchange.com/questions/268331/how-can-i-extract-contours-from-a-raster-with-python) -We hereby demonstrate the first and easiest approach, using `gdal_contour`. Although we deviate from exclusively using the Python language, the benefit of `gdal_contour` is the proven algorithm, customized to spatial data, and with many relevant options. `gdal_contour` (along with other GDAL programs) should already be installed on your system since this is a dependency of `rasterio`. For example, generating 50 $m$ contours of the `dem.tif` file can be done as follows: +We hereby demonstrate the first approach, using `gdal_contour`. +Although we deviate from the Python-focused approach towards more direct interaction with GDAL, the benefit of `gdal_contour` is the proven algorithm, customized to spatial data, and with many relevant options. +Both the `gdal_contour` program (along with other GDAL programs) and its `osgeo` Python wrapper, should already be installed on your system since GDAL is a dependency of `rasterio`. +Using the command line pathway, generating 50 $m$ contours of the `dem.tif` file can be done as follows: ```{python} #| eval: false -os.system('gdal_contour -a elev data/dem.tif output/dem_contour.gpkg -i 50.0') +os.system('gdal_contour -a elev data/dem.tif output/dem_contour1.gpkg -i 50.0') ``` -Note that we ran the `gdal_contour` command through `os.system`, in order to remain in the Python environment. You can also run the standalone command in the command line interface you are using, such as the Anaconda Prompt: +Note that here we ran the `gdal_contour` command through `os.system`, in order to remain in the Python environment. +You can also run the standalone command in the command line interface you are using, such as the Anaconda Prompt: ```sh -gdal_contour -a elev data/dem.tif output/dem_contour.gpkg -i 50.0 +gdal_contour -a elev data/dem.tif output/dem_contour1.gpkg -i 50.0 ``` Like all GDAL programs, `gdal_contour` works with files. Here: * The input is the `data/dem.tif` file -* The result is exported to the `output/dem_contour.gpkg` file +* The result is exported to the `output/dem_contour1.gpkg` file + +To illustrate the result, let's read the resulting `dem_contour2.gpkg` layer back into the Python environment. +Note that the layer contains an arrtibute named `'elev'` (as specified using `-a elev`) with the contour elevation values: + +```{python} +contours1 = gpd.read_file('output/dem_contour1.gpkg') +contours1 +``` -To illustrate the result, let's read the result back into the Python environment. Note that the layer contains an arrtibute named `elev` (as specified using `-a elev`) with the contour elevation values: +Alternatively, we can use the `osgeo` interface to GDAL to do the same, using the [`osgeo.gdal.ContourGenerate`](https://gdal.org/api/python/osgeo.gdal.html#osgeo.gdal.ContourGenerate) wrapper of `gdal_contour`. +Theis makes the code much longer (then just one command, as above), but avoids explicitly running a separate program (`gdal_contour`) outside of the Python environment (since we focus on `rasterio`, explaining the code is beyond the scope of this book): ```{python} -contours = gpd.read_file('output/dem_contour.gpkg') -contours +input_r = osgeo.gdal.Open('data/dem.tif', osgeo.gdalconst.GA_ReadOnly) +input_r1 = input_r.GetRasterBand(1) +sr = osgeo.osr.SpatialReference(input_r.GetProjection()) +ogr_ds = osgeo.ogr.GetDriverByName('GPKG').CreateDataSource('output/dem_contour2.gpkg') +contour_shp = ogr_ds.CreateLayer('contours', sr) +field_defn = osgeo.ogr.FieldDefn('ID', osgeo.ogr.OFTInteger) +contour_shp.CreateField(field_defn) +field_defn = osgeo.ogr.FieldDefn('elev', osgeo.ogr.OFTReal) +contour_shp.CreateField(field_defn) +osgeo.gdal.ContourGenerate(input_r1, 50, 0, [], 0, 0, contour_shp, 0, 1) +ogr_ds = None ``` -Here is a plot of the contour layer in `dem_contour.gpkg` (@fig-raster-contours2): +Here is the identical result of the second approach, named `dem_contour2.gpkg`: + +```{python} +contours2 = gpd.read_file('output/dem_contour2.gpkg') +contours2 +``` + +Here is a plot of the identical contour layers `dem_contour1.gpkg` and `dem_contour2.gpkg` created using the two methods (@fig-raster-contours2): ```{python} #| label: fig-raster-contours2 -#| fig-cap: Raster contours calculated with the `gdal_contour` program +#| fig-cap: Contours of the `dem.tif` raster +#| layout-ncol: 2 +#| fig-subcap: +#| - Using the command line `gdal_contour` program +#| - Using the `osgeo.gdal.ContourGenerate` wrapper + +# Using 'gdal_contour' +fig, ax = plt.subplots() +rasterio.plot.show(src_dem, ax=ax) +contours1.plot(ax=ax, edgecolor='black'); +# Using 'osgeo.gdal.ContourGenerate' fig, ax = plt.subplots() rasterio.plot.show(src_dem, ax=ax) -contours.plot(ax=ax, edgecolor='black'); +contours2.plot(ax=ax, edgecolor='black'); ``` ## Distance to nearest geometry {#sec-distance-to-nearest-geometry} diff --git a/output/dem_contour.gpkg-shm b/output/dem_contour.gpkg-shm deleted file mode 100644 index fe9ac284..00000000 Binary files a/output/dem_contour.gpkg-shm and /dev/null differ diff --git a/output/dem_contour.gpkg-wal b/output/dem_contour.gpkg-wal deleted file mode 100644 index e69de29b..00000000 diff --git a/output/dem_contour.gpkg b/output/dem_contour1.gpkg similarity index 99% rename from output/dem_contour.gpkg rename to output/dem_contour1.gpkg index 3a84adac..8df9f152 100644 Binary files a/output/dem_contour.gpkg and b/output/dem_contour1.gpkg differ diff --git a/output/dem_contour2.gpkg b/output/dem_contour2.gpkg new file mode 100644 index 00000000..c65adc49 Binary files /dev/null and b/output/dem_contour2.gpkg differ