.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "user-guide\04-unstructured-grids.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_user-guide_04-unstructured-grids.py: Unstructured Grids ================== MODFLOW 6 supports unstructured grids. Unlike raster data, the connectivity of unstructured grids is irregular. This means that the number of neighboring cells is not constant; for structured grids, every cell has 4 neighbors (in 2D), or 6 neighbors (in 3D). The package we use to handle unstructured grids is called `xugrid `_. It is a Python package that provides a data structure for unstructured grids, and allows for plotting and analysis of these grids. It is built on top of `xarray `_`, so handling data is similar to how you would handle raster data in xarray, with some differences. .. GENERATED FROM PYTHON SOURCE LINES 19-21 Let's first load some sample data from xugrid.data module. The example data is a triangular grid of the elevation of the Netherlands. .. GENERATED FROM PYTHON SOURCE LINES 22-29 .. code-block:: Python import xugrid as xu elevation = xu.data.elevation_nl() elevation .. raw:: html

<xarray.DataArray 'elevation' (mesh2d_nFaces: 5248)> Size: 21kB
    [5248 values with dtype=float32]
    Coordinates:
      * mesh2d_nFaces  (mesh2d_nFaces) int64 42kB 0 1 2 3 4 ... 5244 5245 5246 5247
        mesh2d_face_x  (mesh2d_nFaces) float64 42kB ...
        mesh2d_face_y  (mesh2d_nFaces) float64 42kB ...
    Attributes:
        unit:     m NAP

.. GENERATED FROM PYTHON SOURCE LINES 30-35 Note that this data is stored differently than a raster dataset. We don't see an x- and a y-dimension here. So where is this spatial information? It's in this ``grid`` attribute, which is accesses via the ``ugrid`` accessor. This is a special accessor that provides access to the unstructured grid properties of the dataset. .. GENERATED FROM PYTHON SOURCE LINES 36-40 .. code-block:: Python grid = elevation.ugrid.grid grid .. rst-class:: sphx-glr-script-out .. code-block:: none Size: 255kB Dimensions: (mesh2d_nFaces: 5248, nmax_face: 3, mesh2d_nNodes: 2790) Coordinates: mesh2d_face_x (mesh2d_nFaces) float64 42kB ... mesh2d_face_y (mesh2d_nFaces) float64 42kB ... mesh2d_node_x (mesh2d_nNodes) float64 22kB ... mesh2d_node_y (mesh2d_nNodes) float64 22kB ... Dimensions without coordinates: mesh2d_nFaces, nmax_face, mesh2d_nNodes Data variables: mesh2d int32 4B ... mesh2d_face_nodes (mesh2d_nFaces, nmax_face) float64 126kB 138.0 ... 1.9... Attributes: Conventions: CF-1.8 UGRID-1.0 .. GENERATED FROM PYTHON SOURCE LINES 41-48 The grid has a number of properties, such as the number of cells, the coordinates of the vertices, and the connectivity of the cells. The latter is very important when dealing with unstructured grids, as it defines how the cells are connected to each other. This is different from structured grids, where you don't need to specify this as there connectivity follows from the rows and columns. Let's take a look at the connectivity of this unstructured grid: .. GENERATED FROM PYTHON SOURCE LINES 49-52 .. code-block:: Python grid.format_connectivity_as_dense(grid.face_face_connectivity) .. rst-class:: sphx-glr-script-out .. code-block:: none array([[ 89, -1, -1], [ 244, 362, 367], [ 15, 279, 335], ..., [2041, 3617, 5246], [ 528, 3618, 5245], [ 528, 3617, 3640]], shape=(5248, 3)) .. GENERATED FROM PYTHON SOURCE LINES 53-63 This shows the connectivity of the cells in the grid. Each row represents a cell, and the columns represent the neighboring cells. Because this is a grid with triangles, we can see three columns in the connectivity matrix. -1 indicates that a cell is not connected to a neighbour. For example, the first cell is only connected to one cell. You can imagine that plotting unstructured grids is quite a hassle, as you need to provide quite some information to the plotting function. Fortunately, the xugrid package provides a convenient way to plot unstructured grids. Let's plot the grid: .. GENERATED FROM PYTHON SOURCE LINES 64-67 .. code-block:: Python grid.plot() .. image-sg:: /user-guide/images/sphx_glr_04-unstructured-grids_001.png :alt: 04 unstructured grids :srcset: /user-guide/images/sphx_glr_04-unstructured-grids_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 68-69 We can plot the grid with the data values as well: .. GENERATED FROM PYTHON SOURCE LINES 70-73 .. code-block:: Python elevation.ugrid.plot() .. image-sg:: /user-guide/images/sphx_glr_04-unstructured-grids_002.png :alt: 04 unstructured grids :srcset: /user-guide/images/sphx_glr_04-unstructured-grids_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 74-81 Note that we need to call the ``ugrid`` accessor to plot the data. This is because the data is stored in a different format than raster data, and we need to tell python that xugrid needs to be used for plotting, and not the regular xarray plotting methods. We can make calculations just like we would do with raster data. For example, we can calculate the mean elevation of the Netherlands: .. GENERATED FROM PYTHON SOURCE LINES 82-86 .. code-block:: Python mean_elevation = elevation.mean() mean_elevation .. raw:: html

<xarray.DataArray 'elevation' ()> Size: 8B
    array(7.61234045)
    Attributes:
        unit:     m NAP

.. GENERATED FROM PYTHON SOURCE LINES 87-88 Or if we want to compute a grid of the elevation 100 meters below the surface: .. GENERATED FROM PYTHON SOURCE LINES 89-93 .. code-block:: Python elevation_100m_below = elevation - 100 elevation_100m_below.ugrid.plot() .. image-sg:: /user-guide/images/sphx_glr_04-unstructured-grids_003.png :alt: 04 unstructured grids :srcset: /user-guide/images/sphx_glr_04-unstructured-grids_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 94-98 We can use this data to be the bottom elevation of a MODFLOW 6 layered unstructured discretization (DISV). For that we need to assign a layer coordinate as iMOD Python requires a layer dimension to be present on the bottom and idomain variables. .. GENERATED FROM PYTHON SOURCE LINES 99-102 .. code-block:: Python bottom = elevation_100m_below.expand_dims("layer").assign_coords(layer=[1]) .. GENERATED FROM PYTHON SOURCE LINES 103-107 Finally we need to create an idomain variable that indicates which cells are active. For this example we will assume all cells are active, so we can just create an idomain variable with all ones. We can do this with the `xu.ones_like` function. .. GENERATED FROM PYTHON SOURCE LINES 108-111 .. code-block:: Python idomain = xu.ones_like(bottom).astype(int) .. GENERATED FROM PYTHON SOURCE LINES 112-115 Now we can create a MODFLOW 6 unstructured discretization (DISV) object using the :class:`imod.mf6.VerticesDiscretization` class. We will create a discretization with one model layer. .. GENERATED FROM PYTHON SOURCE LINES 116-121 .. code-block:: Python from imod.mf6 import VerticesDiscretization disv = VerticesDiscretization(top=elevation, bottom=bottom, idomain=idomain) .. GENERATED FROM PYTHON SOURCE LINES 122-128 Note: we just created a model discretization on a triangular grid. This means that a line connecting two cell centers is rarely perfectly orthogonal to the cell edge. This causes mass balance errors. A voronoi grid would be a better choice for a model discretization, as it ensures that the cell centers are perfectly orthogonal to the cell edges. Luckily for us, xugrid has functionality to tesselate a grid into voronoi cells. .. GENERATED FROM PYTHON SOURCE LINES 129-134 .. code-block:: Python voronoi_grid = grid.tesselate_centroidal_voronoi() voronoi_grid.plot() .. image-sg:: /user-guide/images/sphx_glr_04-unstructured-grids_004.png :alt: 04 unstructured grids :srcset: /user-guide/images/sphx_glr_04-unstructured-grids_004.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 135-141 We can use this voronoi grid to create a new MODFLOW 6 unstructured discretization (DISV) object. iMOD Python has regridding functionality to regrid the top, bottom and idomain variables to the new grid. For that, iMOD Python regridding functionality requires a UgridDataArray instead of a Ugrid2d, so we create a UgridDataArray with the voronoi grid. .. GENERATED FROM PYTHON SOURCE LINES 142-148 .. code-block:: Python from imod.util import RegridderWeightsCache, ones_like_ugrid voronoi_uda = ones_like_ugrid(voronoi_grid) disv_voronoi = disv.regrid_like(voronoi_uda, RegridderWeightsCache()) disv_voronoi .. raw:: html

.. GENERATED FROM PYTHON SOURCE LINES 149-154 Other packages like the :class:`imod.mf6.Drainage` package, accept both structured and unstructured grids. The drain elevation requires a layer dimension. We can use the `idomain` to broadcast the data to include the layer dimension. .. GENERATED FROM PYTHON SOURCE LINES 155-159 .. code-block:: Python drain_elevation = disv_voronoi["idomain"] * disv_voronoi["top"] drain_elevation .. raw:: html

<xarray.DataArray (layer: 1, mesh2d_nFaces: 2790)> Size: 22kB
    array([[21.14808083, 64.5098114 , 43.33101654, ...,  4.84251642,
            -0.14340343, 37.67595291]], shape=(1, 2790))
    Coordinates:
      * layer          (layer) int64 8B 1
      * mesh2d_nFaces  (mesh2d_nFaces) int64 22kB 0 1 2 3 4 ... 2786 2787 2788 2789
    Attributes:
        unit:     m NAP

.. GENERATED FROM PYTHON SOURCE LINES 160-164 We also need to create a conductance variable for the :class:`imod.mf6.Drainage` package. The conductance variable is a measure of how much water can flow through a cell. The areas of the cells are used to calculate the conductance, and can be be obtained from its grid. .. GENERATED FROM PYTHON SOURCE LINES 165-169 .. code-block:: Python resistance = 1 # Example resistance value in days conductance = disv_voronoi["idomain"] * voronoi_grid.area / resistance .. GENERATED FROM PYTHON SOURCE LINES 170-171 We can now use these variables to create a :class:`imod.mf6.Drainage` package object. .. GENERATED FROM PYTHON SOURCE LINES 172-181 .. code-block:: Python from imod.mf6 import Drainage drainage = Drainage( elevation=drain_elevation, conductance=conductance, ) drainage .. raw:: html

.. GENERATED FROM PYTHON SOURCE LINES 182-186 That wraps up the example of how to use unstructured grids in iMOD Python. For more information on how to manipulate unstructured grids and how to make selections, please refer to the `xugrid documentation `_. .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 4.590 seconds) .. _sphx_glr_download_user-guide_04-unstructured-grids.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 04-unstructured-grids.ipynb <04-unstructured-grids.ipynb>` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 04-unstructured-grids.py <04-unstructured-grids.py>` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: 04-unstructured-grids.zip <04-unstructured-grids.zip>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_