.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples/hydamo_network.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_examples_hydamo_network.py: Processing Elevation Data Along Ugrid1d Networks ================================================ In this example, we'll fetch line data describing ditches in the south of the Netherlands and demonstrate several xugrid techniques for working with 1D networks: * Converting geometries to Ugrid1d networks * Intersecting edges with external line features * Refining networks by inserting vertices * Topology-aware interpolation (nearest neighbor and Laplace) * Converting node data to edge data * Network visualization The dataset contains three components: * Hydro-objects: center lines representing the ditches. * Profile points: elevation measurements sampled perpendicular to each center line, describing the cross-sectional profile. * Profile lines: perpendicular transects connecting the measured points. We'll do some basic processing on these data. Our goal is to get an estimate of the bed elevation along the center lines of the ditches. .. GENERATED FROM PYTHON SOURCE LINES 26-34 .. code-block:: Python import matplotlib.patches as patches import matplotlib.pyplot as plt import numpy as np import shapely import xugrid as xu .. GENERATED FROM PYTHON SOURCE LINES 35-36 The example data consists of three separate GeoDataFrames: .. GENERATED FROM PYTHON SOURCE LINES 36-39 .. code-block:: Python objects, points, profiles = xu.data.hydamo_network() .. GENERATED FROM PYTHON SOURCE LINES 40-42 We will take a look at the center lines and the transects. To get a better idea of the data, we will also zoom in on a 100 by 100 meter window: .. GENERATED FROM PYTHON SOURCE LINES 42-57 .. code-block:: Python xy = (140_270.0, 393_140.0) dx = dy = 100.0 fig, (ax0, ax1) = plt.subplots(ncols=2, figsize=(10, 5)) objects.plot(ax=ax0, column="type", legend=True) profiles.plot(ax=ax0, color="red") ax0.add_patch(patches.Rectangle(xy, dx, dy, fill=False)) objects.plot(ax=ax1, column="type") profiles.plot(ax=ax1, color="red") points.plot(ax=ax1, color="black") ax1.set_xlim(xy[0], xy[0] + dx) ax1.set_ylim(xy[1], xy[1] + dy) .. image-sg:: /examples/images/sphx_glr_hydamo_network_001.png :alt: hydamo network :srcset: /examples/images/sphx_glr_hydamo_network_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none (393140.0, 393240.0) .. GENERATED FROM PYTHON SOURCE LINES 58-62 For spatial analysis, we require sufficient spatial detail -- e.g. to compare with a raster digital elevation map (DEM) or for it to serve as model input. We may use shapely to limit the maximum length of linear elements to 25 m. Next, we will generate a Ugrid1d network from the resulting GeoDataFrame. .. GENERATED FROM PYTHON SOURCE LINES 62-66 .. code-block:: Python discretized = shapely.segmentize(objects.geometry, 25.0) network = xu.Ugrid1d.from_shapely(discretized) .. GENERATED FROM PYTHON SOURCE LINES 67-69 Let's see whether the conversion to a network topology went well, and what the result of segmentizing is. .. GENERATED FROM PYTHON SOURCE LINES 69-74 .. code-block:: Python fig, ax = plt.subplots() network.plot(ax=ax) ax.scatter(*network.node_coordinates.T) .. image-sg:: /examples/images/sphx_glr_hydamo_network_002.png :alt: hydamo network :srcset: /examples/images/sphx_glr_hydamo_network_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 75-86 We can see a fairly consistent distribution of nodes. In the next steps, we would like to associate elevation data with nodes of the network. We take the following steps: * Organize the points per transects and determine its lowest elevation. * Intersect the center lines with the transects. * Insert the locations of these intersections as nodes into our network. * Associate the elevations with the inserted nodes. Finding the lowest elevation can be done with a pandas groupby. We then add this data to the profile transects, then plot it to check it. .. GENERATED FROM PYTHON SOURCE LINES 86-91 .. code-block:: Python profiles = profiles.set_index("line_id") profiles["elevation"] = points.groupby("line_id")["elevation"].mean() profiles.plot(column="elevation", legend=True) .. image-sg:: /examples/images/sphx_glr_hydamo_network_003.png :alt: hydamo network :srcset: /examples/images/sphx_glr_hydamo_network_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 92-100 Not much pattern: a mean elevation centered around 11.5 m above mean sea level and some local variation. Looks like we've empirically verified that this part of the Netherlands is indeed rather flat! For the next step, we will intersect the network with the perpendicular transect lines, and insert the intersections as nodes. The :meth:`xugrid.intersect_edges` method accepts linework as a numpy array of coordinates. .. GENERATED FROM PYTHON SOURCE LINES 100-104 .. code-block:: Python edges = shapely.get_coordinates(profiles.geometry).reshape((-1, 2, 2)) edge_index, core_index, intersections = network.intersect_edges(edges) .. GENERATED FROM PYTHON SOURCE LINES 105-106 Let's verify the intersection locations. .. GENERATED FROM PYTHON SOURCE LINES 106-115 .. code-block:: Python fig, ax = plt.subplots() network.plot(ax=ax, zorder=1) profiles.plot(ax=ax, zorder=1) ax.scatter(*network.node_coordinates.T) ax.scatter(*intersections.T) ax.set_xlim(xy[0], xy[0] + dx) ax.set_ylim(xy[1], xy[1] + dy) .. image-sg:: /examples/images/sphx_glr_hydamo_network_004.png :alt: hydamo network :srcset: /examples/images/sphx_glr_hydamo_network_004.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none (393140.0, 393240.0) .. GENERATED FROM PYTHON SOURCE LINES 116-118 This looks reasonable. We will now insert the intersections in the network, create a UgridDataArray with the known elevations at the intersections. .. GENERATED FROM PYTHON SOURCE LINES 118-130 .. code-block:: Python refined, insert_index = network.refine_by_vertices( intersections, return_index=True, tolerance=1e-6 ) data = np.full(refined.n_node, np.nan) data[insert_index] = profiles["elevation"].to_numpy()[edge_index] uda = xu.UgridDataArray.from_data(data, refined, facet="node") fig, ax = plt.subplots() refined.plot(ax=ax, zorder=1) uda.ugrid.plot(ax=ax) .. image-sg:: /examples/images/sphx_glr_hydamo_network_005.png :alt: hydamo network :srcset: /examples/images/sphx_glr_hydamo_network_005.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 131-148 Our end goal is to have elevation data for all edges of the the network. To do so, we can interpolate along the edges, replacing the NoData (NaN) values. Xugrid supports :meth:`xugrid.UgridDataArrayAccessor.interpolate_na()` and :meth:`xugrid.UgridDataArrayAccessor.laplace_interpolate()` for these goals. Unlike general spatial interpolation methods, these interpolation methods take the topology into account. For example, the nearest distance here is defined along the edges; if there is no connection between two nodes, they does not count towards each other's nearest neighbors. Laplace interpolation here smoothly fills gaps by making each unknown value the average of its connected neighbors along the network, similar to how water levels would equilibrate if flowing through the ditches (with a linear hydraulic resistance). As a final step, we will associate the elevations with the edges, by assigning for each edge the average of its nodes. .. GENERATED FROM PYTHON SOURCE LINES 148-155 .. code-block:: Python nearest = uda.ugrid.interpolate_na() smooth = uda.ugrid.laplace_interpolate() edge_nearest = nearest.ugrid.to_edge().mean("nmax") edge_smooth = smooth.ugrid.to_edge().mean("nmax") .. GENERATED FROM PYTHON SOURCE LINES 156-157 Let's plot both, and compare the differences. .. GENERATED FROM PYTHON SOURCE LINES 157-168 .. code-block:: Python fig, axes = plt.subplots(ncols=3, figsize=(18, 5)) edge_nearest.ugrid.plot(ax=axes[0], linewidth=3) edge_smooth.ugrid.plot(ax=axes[1], linewidth=3) (edge_smooth - edge_nearest).ugrid.plot(ax=axes[2], linewidth=3) for ax in axes: ax.set_aspect(1.0) axes[0].set_title("Nearest neighbor interpolation") axes[1].set_title("Laplace interpolation") axes[2].set_title("Difference") .. image-sg:: /examples/images/sphx_glr_hydamo_network_006.png :alt: Nearest neighbor interpolation, Laplace interpolation, Difference :srcset: /examples/images/sphx_glr_hydamo_network_006.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Text(0.5, 1.0, 'Difference') .. GENERATED FROM PYTHON SOURCE LINES 169-173 The actual data here do little to illustrate the interpolation methods! Let's generate some synthetic data with a clearer pattern. We will assume that the land (and bed) elevation slopes downward in a north-easterly direction with a gradient of 5%. .. GENERATED FROM PYTHON SOURCE LINES 173-188 .. code-block:: Python refine_intersections = refined.node_coordinates[insert_index] xmax, ymax = refine_intersections.max(axis=0) synthetic_elevation = ( xmax - refine_intersections[:, 0] + ymax - refine_intersections[:, 1] ) * 0.05 data = np.full(refined.n_node, np.nan) data[insert_index] = synthetic_elevation uda = xu.UgridDataArray.from_data(data, refined, facet="node") fig, ax = plt.subplots() refined.plot(ax=ax, zorder=1, linewidth=0.5) uda.ugrid.plot(ax=ax) .. image-sg:: /examples/images/sphx_glr_hydamo_network_007.png :alt: hydamo network :srcset: /examples/images/sphx_glr_hydamo_network_007.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 189-194 Now let's compare the two interpolation methods. Nearest neighbor interpolation assigns each gap the value of its closest known neighbor along the network. Laplace interpolation creates a smoother solution where each unknown value becomes the average of its connected neighbors— similar to how water levels would equilibrate if flowing through the ditches. .. GENERATED FROM PYTHON SOURCE LINES 194-213 .. code-block:: Python nearest = uda.ugrid.interpolate_na() smooth = uda.ugrid.laplace_interpolate() edge_nearest = nearest.ugrid.to_edge().mean("nmax") edge_smooth = smooth.ugrid.to_edge().mean("nmax") fig, axes = plt.subplots(ncols=3, figsize=(18, 5)) edge_nearest.ugrid.plot(ax=axes[0], linewidth=3) edge_smooth.ugrid.plot(ax=axes[1], linewidth=3) (edge_smooth - edge_nearest).ugrid.plot(ax=axes[2], linewidth=3) for ax in axes: ax.set_aspect(1.0) axes[0].set_title("Nearest neighbor interpolation") axes[1].set_title("Laplace interpolation") axes[2].set_title("Difference") .. image-sg:: /examples/images/sphx_glr_hydamo_network_008.png :alt: Nearest neighbor interpolation, Laplace interpolation, Difference :srcset: /examples/images/sphx_glr_hydamo_network_008.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Text(0.5, 1.0, 'Difference') .. GENERATED FROM PYTHON SOURCE LINES 214-218 We can now clearly see how reaches without data are filled "along" the network topology. The difference plot shows that Laplace interpolation produces smoother transitions, while nearest neighbor creates more abrupt changes where the nearest known value switches from one transect to another. .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 1.701 seconds) .. _sphx_glr_download_examples_hydamo_network.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: hydamo_network.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: hydamo_network.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: hydamo_network.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_