.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "user-guide\07-regridding.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_07-regridding.py: Regridding ========== Introduction ------------ Most MODFLOW 6 packages have spatial data arrays as input. These arrays are discrete: they are defined over the simulation grid and contain values associated to each cell. Regridding these package means: create a new package with spatial data arrays defined over a different grid. Computing what the values in these new arrays should be, is done by xugrid. xugrid will compute a regridded array based on: - the original array - the original discretization (this is described in the coordinates of the original arry) - the new discretization - a regridding method More information on the available regridding methods can be found in the xugrid documentation https://deltares.github.io/xugrid/user_guide.html The regridding method that should be used depends on the property being regridded. For example a point-based property (whose value do not depend intrinsically on the grid block size) such as temperature or density can be regridded by an averaging approach (for upscaling) or sampling (for downscaling). Volume-based properties (whose values do depend on the grid block size) include (water) mass and pore volume of a gridblock, and the regridding method should be chosen to reflect that. Finally regridding methods for conductivity-like properties follow the rules for parallel or serial resistors- at least when the tensor rotation angles are constant or comparable in the involved gridblocks. Note that the different regridding methods may have a different output domain when regridding: if the original array has no-data values in some cells, then the output array may have no-data values as well, and where these end up depends on the chosen regridding method. Also note that regridding is only possible in the xy-plane, and not across the layer dimension. The output array will have the same number of layers as the input array. .. GENERATED FROM PYTHON SOURCE LINES 45-47 .. code-block:: Python :dedent: 1 .. GENERATED FROM PYTHON SOURCE LINES 49-91 Obtaining the final (i)domain ----------------------------- In many real-world models, some cells will be inactive or marked as "vertical passthrough" (VPT) in the idomain array of the simulation. Some packages require that all cells that are inactictive or VPT in idomain are excluded from the package as well. An example is the :class:`imod.mf6.NodePropertyFlow` package: cells that are inactive or VPT in idomain, should not have conductivity data in the npf package. Therefore at the end of the regridding process, a final step consists in enforcing consistency between those of idomain and all the packages. This is a 2-step process: 1) for cells that do not have inputs in crucial packages like npf or storage, idomain will be set to inactive. 2) for cells that are marked as inactive or VPT in idomain, all package inputs will be removed from all the packages This synchronization step between idomain and the packages is automated, and it is carried out when calling regrid_like on the simulation object or the model object. There are 2 caveats: 1) the synchronization between idomain and the package domains is done on the model-level. If we have a simulation containing both a flow model and a transport model then idomain for flow is determined independent from that for transport. These models may therefore end up using different domains (this may lead to undesired results, so a manual synchronization may be necessary between the flow regridding and the transport regridding) This manual synchronization can be done using the "mask_all_packages" function- this function removes input from all packages that are marked as inactive or VPT in the idomain passed to this method. 2) The idomain/packages synchronization step is carried out when regridding a model, and when regridding a model it will use default methods for all the packages. So if you have regridded some packages yourself with non-default methods, then these are not taken into account during this synchonization step. Regridding using default methods -------------------------------- The regrid_like function is available on packages, models and simulations. When the default methods are acceptable, regridding the whole simulation is the most convenient from a user-perspective. .. GENERATED FROM PYTHON SOURCE LINES 91-98 .. code-block:: Python import imod tmpdir = imod.util.temporary_directory() original_simulation = imod.data.hondsrug_simulation(tmpdir / "hondsrug_saved") .. GENERATED FROM PYTHON SOURCE LINES 99-101 To reduce computational overhead for this example, we are going to clip off most of the timesteps. .. GENERATED FROM PYTHON SOURCE LINES 101-104 .. code-block:: Python original_simulation = original_simulation.clip_box(time_max="2010-01-01") .. GENERATED FROM PYTHON SOURCE LINES 105-106 Let's take a look at the original discretization: .. GENERATED FROM PYTHON SOURCE LINES 106-109 .. code-block:: Python original_simulation["GWF"]["dis"] .. raw:: html
StructuredDiscretization
<xarray.Dataset> Size: 11MB
    Dimensions:  (layer: 13, y: 200, x: 500)
    Coordinates:
      * layer    (layer) int32 52B 1 2 3 4 5 6 7 8 9 10 11 12 13
      * y        (y) float64 2kB 5.64e+05 5.64e+05 5.639e+05 ... 5.59e+05 5.59e+05
      * x        (x) float64 4kB 2.375e+05 2.375e+05 2.376e+05 ... 2.5e+05 2.5e+05
        dx       float64 8B ...
        dy       float64 8B ...
    Data variables:
        idomain  (layer, y, x) int32 5MB dask.array<chunksize=(13, 200, 500), meta=np.ndarray>
        top      (y, x) float32 400kB dask.array<chunksize=(200, 500), meta=np.ndarray>
        bottom   (layer, y, x) float32 5MB dask.array<chunksize=(13, 200, 500), meta=np.ndarray>


.. GENERATED FROM PYTHON SOURCE LINES 110-111 We want to regrid this to the following target grid: .. GENERATED FROM PYTHON SOURCE LINES 111-127 .. code-block:: Python import xarray as xr import imod dx = 100 dy = -100 xmin = 237500.0 xmax = 250000.0 ymin = 559000.0 ymax = 564000.0 target_grid = imod.util.empty_2d( dx=dx, dy=dy, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax ) target_grid .. raw:: html 
<xarray.DataArray (y: 50, x: 125)> Size: 50kB
    array([[nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan],
           ...,
           [nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan]], shape=(50, 125))
    Coordinates:
      * y        (y) float64 400B 5.64e+05 5.638e+05 ... 5.592e+05 5.59e+05
      * x        (x) float64 1kB 2.376e+05 2.376e+05 2.378e+05 ... 2.498e+05 2.5e+05
        dx       float64 8B 100.0
        dy       float64 8B -100.0


.. GENERATED FROM PYTHON SOURCE LINES 128-130 This is a grid of nans, and we require a grid of ones, which can create with xarray: .. GENERATED FROM PYTHON SOURCE LINES 130-134 .. code-block:: Python target_grid = xr.ones_like(target_grid) target_grid .. raw:: html 
<xarray.DataArray (y: 50, x: 125)> Size: 50kB
    array([[1., 1., 1., ..., 1., 1., 1.],
           [1., 1., 1., ..., 1., 1., 1.],
           [1., 1., 1., ..., 1., 1., 1.],
           ...,
           [1., 1., 1., ..., 1., 1., 1.],
           [1., 1., 1., ..., 1., 1., 1.],
           [1., 1., 1., ..., 1., 1., 1.]], shape=(50, 125))
    Coordinates:
      * y        (y) float64 400B 5.64e+05 5.638e+05 ... 5.592e+05 5.59e+05
      * x        (x) float64 1kB 2.376e+05 2.376e+05 2.378e+05 ... 2.498e+05 2.5e+05
        dx       float64 8B 100.0
        dy       float64 8B -100.0


.. GENERATED FROM PYTHON SOURCE LINES 135-136 Now regrid the simulation (without recharge): .. GENERATED FROM PYTHON SOURCE LINES 136-138 .. code-block:: Python regridded_simulation = original_simulation.regrid_like("regridded", target_grid) .. GENERATED FROM PYTHON SOURCE LINES 139-140 Let's look at the discretization again: .. GENERATED FROM PYTHON SOURCE LINES 140-142 .. code-block:: Python regridded_simulation["GWF"]["dis"] .. raw:: html
StructuredDiscretization
<xarray.Dataset> Size: 1MB
    Dimensions:  (layer: 13, y: 50, x: 125)
    Coordinates:
      * layer    (layer) int32 52B 1 2 3 4 5 6 7 8 9 10 11 12 13
      * y        (y) float64 400B 5.64e+05 5.638e+05 ... 5.592e+05 5.59e+05
      * x        (x) float64 1kB 2.376e+05 2.376e+05 2.378e+05 ... 2.498e+05 2.5e+05
        dx       float64 8B 100.0
        dy       float64 8B -100.0
    Data variables:
        idomain  (layer, y, x) int64 650kB dask.array<chunksize=(13, 50, 125), meta=np.ndarray>
        top      (y, x) float32 25kB dask.array<chunksize=(50, 125), meta=np.ndarray>
        bottom   (layer, y, x) float32 325kB dask.array<chunksize=(13, 50, 125), meta=np.ndarray>


.. GENERATED FROM PYTHON SOURCE LINES 143-144 All packages have been regridded, for example the NPF package: .. GENERATED FROM PYTHON SOURCE LINES 144-147 .. code-block:: Python regridded_simulation["GWF"]["npf"] .. raw:: html
NodePropertyFlow
<xarray.Dataset> Size: 652kB
    Dimensions:                              (layer: 13, y: 50, x: 125)
    Coordinates:
      * layer                                (layer) int32 52B 1 2 3 4 ... 11 12 13
      * y                                    (y) float64 400B 5.64e+05 ... 5.59e+05
      * x                                    (x) float64 1kB 2.376e+05 ... 2.5e+05
        dx                                   float64 8B 100.0
        dy                                   float64 8B -100.0
    Data variables: (12/22)
        icelltype                            int32 4B ...
        k                                    (layer, y, x) float32 325kB dask.array<chunksize=(13, 50, 125), meta=np.ndarray>
        rewet                                bool 1B False
        rewet_layer                          object 8B None
        rewet_factor                         object 8B None
        rewet_iterations                     object 8B None
        ...                                   ...
        dewatered                            bool 1B True
        perched                              bool 1B True
        save_specific_discharge              bool 1B False
        save_saturation                      bool 1B False
        xt3d_option                          bool 1B False
        rhs_option                           bool 1B False


.. GENERATED FROM PYTHON SOURCE LINES 148-149 Let's make a comparison plot of the hydraulic conductivities: .. GENERATED FROM PYTHON SOURCE LINES 149-165 .. code-block:: Python import matplotlib.pyplot as plt import numpy as np fig, axes = plt.subplots(nrows=2, sharex=True) plot_kwargs = {"colors": "viridis", "levels": np.linspace(0.0, 100.0, 21), "fig": fig} imod.visualize.spatial.plot_map( original_simulation["GWF"]["npf"]["k"].sel(layer=3), ax=axes[0], **plot_kwargs ) imod.visualize.spatial.plot_map( regridded_simulation["GWF"]["npf"]["k"].sel(layer=3), ax=axes[1], **plot_kwargs ) axes[0].set_ylabel("original") axes[1].set_ylabel("regridded") .. image-sg:: /user-guide/images/sphx_glr_07-regridding_001.png :alt: 07 regridding :srcset: /user-guide/images/sphx_glr_07-regridding_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Text(28.472222222222243, 0.5, 'regridded') .. GENERATED FROM PYTHON SOURCE LINES 166-188 Regridding using non-default methods ------------------------------------ When non-default methods are used for one or more packages, these should be regridded separately. In that case, the most convenient approach is likely: - pop the packages that should use non-default methods from the source simulation (the popping is optional, and is only recommended for packages whose presence is not mandatory for validation.) - regrid the source simulation: this takes care of all the packages that should use default methods. - regrid the package(s) where you want to use non-standard rergridding methods indivudually starting from the packages in the source simulation - insert the custom-regridded packages to the regridded simulation (or replace the package regridded with default methods with the one you just regridded with non-default methods if it was not popped) In code, consider an example where we want to regrid the recharge package using non default methods then we would do the following. First we'll load some example simulation. There is a separate example contained in :doc:`/examples/mf6/hondsrug` that you should look at if you are interested in the model building .. GENERATED FROM PYTHON SOURCE LINES 191-194 Set up the input needed for custom regridding. Create a regridder weight-cache. This object can (and should) be reused for all the packages that undergo custom regridding at this stage. .. GENERATED FROM PYTHON SOURCE LINES 195-201 .. code-block:: Python from imod.util.regrid import RegridderWeightsCache regrid_cache = RegridderWeightsCache() regrid_cache .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 202-204 Next, we'll remove the recharge package and obtain it as a variable. We'll do this for later use. .. GENERATED FROM PYTHON SOURCE LINES 205-210 .. code-block:: Python original_rch_package = original_simulation["GWF"].pop("rch") original_rch_package .. raw:: html
Recharge
<xarray.Dataset> Size: 21MB
    Dimensions:        (time: 2, layer: 13, y: 200, x: 500)
    Coordinates:
      * time           (time) datetime64[ns] 16B 2009-12-30T23:59:59 2009-12-31
      * layer          (layer) int32 52B 1 2 3 4 5 6 7 8 9 10 11 12 13
      * y              (y) float64 2kB 5.64e+05 5.64e+05 ... 5.59e+05 5.59e+05
      * x              (x) float64 4kB 2.375e+05 2.375e+05 ... 2.5e+05 2.5e+05
        dx             float64 8B ...
        dy             float64 8B ...
    Data variables:
        rate           (time, layer, y, x) float64 21MB dask.array<chunksize=(2, 13, 200, 500), meta=np.ndarray>
        print_input    bool 1B ...
        print_flows    bool 1B ...
        save_flows     bool 1B ...
        observations   object 8B None
        repeat_stress  object 8B None


.. GENERATED FROM PYTHON SOURCE LINES 211-214 Regrid the recharge package with a custom regridder. In this case we opt for the centroid locator regridder. This regridder is similar to using a "nearest neighbour" lookup. .. GENERATED FROM PYTHON SOURCE LINES 215-228 .. code-block:: Python from imod.common.utilities.regrid import RegridderType from imod.mf6.regrid import RechargeRegridMethod regridder_types = RechargeRegridMethod(rate=(RegridderType.CENTROIDLOCATOR,)) regridded_recharge = original_rch_package.regrid_like( target_grid, regrid_cache=regrid_cache, regridder_types=regridder_types, ) regridded_recharge .. raw:: html
Recharge
<xarray.Dataset> Size: 1MB
    Dimensions:        (layer: 13, time: 2, y: 50, x: 125)
    Coordinates:
      * layer          (layer) int32 52B 1 2 3 4 5 6 7 8 9 10 11 12 13
      * time           (time) datetime64[ns] 16B 2009-12-30T23:59:59 2009-12-31
      * y              (y) float64 400B 5.64e+05 5.638e+05 ... 5.592e+05 5.59e+05
      * x              (x) float64 1kB 2.376e+05 2.376e+05 ... 2.498e+05 2.5e+05
        dx             float64 8B 100.0
        dy             float64 8B -100.0
    Data variables:
        rate           (time, layer, y, x) float64 1MB dask.array<chunksize=(2, 13, 50, 125), meta=np.ndarray>
        print_input    bool 1B False
        print_flows    bool 1B False
        save_flows     bool 1B False
        observations   object 8B None
        repeat_stress  object 8B None
        fixed_cell     bool 1B False


.. GENERATED FROM PYTHON SOURCE LINES 229-230 Next, add the recharge package to the regridded simulation .. GENERATED FROM PYTHON SOURCE LINES 230-235 .. code-block:: Python regridded_simulation["GWF"]["rch"] = regridded_recharge # We can also reattach the original again original_simulation["GWF"]["rch"] = original_rch_package .. GENERATED FROM PYTHON SOURCE LINES 236-244 Comparison with histograms -------------------------- In the next segment we will compare the input of the models on different grids. We advice to always check how your input is regridded. In this example we upscaled grid, many input parameters are regridded with a ``mean`` method. This means that their input range is reduced, which can be seen in tailings in the histograms becoming shorter .. GENERATED FROM PYTHON SOURCE LINES 244-256 .. code-block:: Python def plot_histograms_side_by_side(array_original, array_regridded, title): """This function creates a plot of normalized histograms of the 2 input DataArray. It plots a title above each histogram.""" _, (ax0, ax1) = plt.subplots(1, 2, sharex=True, sharey=True, tight_layout=True) array_original.plot.hist(ax=ax0, bins=25, density=True) array_regridded.plot.hist(ax=ax1, bins=25, density=True) ax0.title.set_text(f"{title} (original)") ax1.title.set_text(f"{title} (regridded)") .. GENERATED FROM PYTHON SOURCE LINES 257-258 Compare constant head arrays. .. GENERATED FROM PYTHON SOURCE LINES 258-264 .. code-block:: Python plot_histograms_side_by_side( original_simulation["GWF"]["chd"]["head"], regridded_simulation["GWF"]["chd"]["head"], "chd head", ) .. image-sg:: /user-guide/images/sphx_glr_07-regridding_002.png :alt: chd head (original), chd head (regridded) :srcset: /user-guide/images/sphx_glr_07-regridding_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 265-266 Compare horizontal hydraulic conductivities. .. GENERATED FROM PYTHON SOURCE LINES 266-271 .. code-block:: Python plot_histograms_side_by_side( original_simulation["GWF"]["npf"]["k"], regridded_simulation["GWF"]["npf"]["k"], "npf k", ) .. image-sg:: /user-guide/images/sphx_glr_07-regridding_003.png :alt: npf k (original), npf k (regridded) :srcset: /user-guide/images/sphx_glr_07-regridding_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 272-273 Compare vertical hydraulic conductivities. .. GENERATED FROM PYTHON SOURCE LINES 273-278 .. code-block:: Python plot_histograms_side_by_side( original_simulation["GWF"]["npf"]["k33"], regridded_simulation["GWF"]["npf"]["k33"], "npf k33", ) .. image-sg:: /user-guide/images/sphx_glr_07-regridding_004.png :alt: npf k33 (original), npf k33 (regridded) :srcset: /user-guide/images/sphx_glr_07-regridding_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 279-280 Compare starting heads. .. GENERATED FROM PYTHON SOURCE LINES 280-286 .. code-block:: Python plot_histograms_side_by_side( original_simulation["GWF"]["ic"]["start"], regridded_simulation["GWF"]["ic"]["start"], "ic start", ) .. image-sg:: /user-guide/images/sphx_glr_07-regridding_005.png :alt: ic start (original), ic start (regridded) :srcset: /user-guide/images/sphx_glr_07-regridding_005.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 287-288 Compare river stages. .. GENERATED FROM PYTHON SOURCE LINES 288-294 .. code-block:: Python plot_histograms_side_by_side( original_simulation["GWF"]["riv"]["stage"], regridded_simulation["GWF"]["riv"]["stage"], "riv stage", ) .. image-sg:: /user-guide/images/sphx_glr_07-regridding_006.png :alt: riv stage (original), riv stage (regridded) :srcset: /user-guide/images/sphx_glr_07-regridding_006.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 295-296 Compare river bottom elevations. .. GENERATED FROM PYTHON SOURCE LINES 296-302 .. code-block:: Python plot_histograms_side_by_side( original_simulation["GWF"]["riv"]["bottom_elevation"], regridded_simulation["GWF"]["riv"]["bottom_elevation"], "riv bottom elevation", ) .. image-sg:: /user-guide/images/sphx_glr_07-regridding_007.png :alt: riv bottom elevation (original), riv bottom elevation (regridded) :srcset: /user-guide/images/sphx_glr_07-regridding_007.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 303-304 Compare riverbed conductance. .. GENERATED FROM PYTHON SOURCE LINES 304-310 .. code-block:: Python plot_histograms_side_by_side( original_simulation["GWF"]["riv"]["conductance"], regridded_simulation["GWF"]["riv"]["conductance"], "riv conductance", ) .. image-sg:: /user-guide/images/sphx_glr_07-regridding_008.png :alt: riv conductance (original), riv conductance (regridded) :srcset: /user-guide/images/sphx_glr_07-regridding_008.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 311-312 Compare recharge rates. .. GENERATED FROM PYTHON SOURCE LINES 312-318 .. code-block:: Python plot_histograms_side_by_side( original_simulation["GWF"]["rch"]["rate"], regridded_simulation["GWF"]["rch"]["rate"], "rch rate", ) .. image-sg:: /user-guide/images/sphx_glr_07-regridding_009.png :alt: rch rate (original), rch rate (regridded) :srcset: /user-guide/images/sphx_glr_07-regridding_009.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 319-320 Compare drainage elevations. .. GENERATED FROM PYTHON SOURCE LINES 320-326 .. code-block:: Python plot_histograms_side_by_side( original_simulation["GWF"]["drn-pipe"]["elevation"], regridded_simulation["GWF"]["drn-pipe"]["elevation"], "drn-pipe elevation", ) .. image-sg:: /user-guide/images/sphx_glr_07-regridding_010.png :alt: drn-pipe elevation (original), drn-pipe elevation (regridded) :srcset: /user-guide/images/sphx_glr_07-regridding_010.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 327-328 Compare drain conductances. .. GENERATED FROM PYTHON SOURCE LINES 328-334 .. code-block:: Python plot_histograms_side_by_side( original_simulation["GWF"]["drn-pipe"]["conductance"], regridded_simulation["GWF"]["drn-pipe"]["conductance"], "drn-pipe conductance", ) .. image-sg:: /user-guide/images/sphx_glr_07-regridding_011.png :alt: drn-pipe conductance (original), drn-pipe conductance (regridded) :srcset: /user-guide/images/sphx_glr_07-regridding_011.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 335-339 Compare simulation outputs -------------------------- Let's compare outputs now. .. GENERATED FROM PYTHON SOURCE LINES 340-349 .. code-block:: Python # Write simulation original_simulation.write(tmpdir / "original") # Run and open heads original_simulation.run() hds_original = original_simulation.open_head() hds_original .. raw:: html 
<xarray.DataArray 'head' (time: 2, layer: 13, y: 200, x: 500)> Size: 21MB
    dask.array<stack, shape=(2, 13, 200, 500), dtype=float64, chunksize=(1, 13, 200, 500), chunktype=numpy.ndarray>
    Coordinates:
      * time     (time) float64 16B 1.157e-05 365.0
      * layer    (layer) int64 104B 1 2 3 4 5 6 7 8 9 10 11 12 13
      * y        (y) float64 2kB 5.64e+05 5.64e+05 5.639e+05 ... 5.59e+05 5.59e+05
      * x        (x) float64 4kB 2.375e+05 2.375e+05 2.376e+05 ... 2.5e+05 2.5e+05
        dx       float64 8B 25.0
        dy       float64 8B -25.0


.. GENERATED FROM PYTHON SOURCE LINES 350-356 .. code-block:: Python regridded_simulation.write(tmpdir / "regridded", validate=False) regridded_simulation.run() hds_regridded = regridded_simulation.open_head() hds_regridded .. raw:: html 
<xarray.DataArray 'head' (time: 2, layer: 13, y: 50, x: 125)> Size: 1MB
    dask.array<stack, shape=(2, 13, 50, 125), dtype=float64, chunksize=(1, 13, 50, 125), chunktype=numpy.ndarray>
    Coordinates:
      * time     (time) float64 16B 1.157e-05 365.0
      * layer    (layer) int64 104B 1 2 3 4 5 6 7 8 9 10 11 12 13
      * y        (y) float64 400B 5.64e+05 5.638e+05 ... 5.592e+05 5.59e+05
      * x        (x) float64 1kB 2.376e+05 2.376e+05 2.378e+05 ... 2.498e+05 2.5e+05
        dx       float64 8B 100.0
        dy       float64 8B -100.0


.. GENERATED FROM PYTHON SOURCE LINES 357-358 Let's make a comparison plot of the regridded heads. .. GENERATED FROM PYTHON SOURCE LINES 358-371 .. code-block:: Python fig, axes = plt.subplots(nrows=2, sharex=True) plot_kwargs = {"colors": "viridis", "levels": np.linspace(0.0, 11.0, 12), "fig": fig} imod.visualize.spatial.plot_map( hds_original.isel(layer=6, time=-1), ax=axes[0], **plot_kwargs ) imod.visualize.spatial.plot_map( hds_regridded.isel(layer=6, time=-1), ax=axes[1], **plot_kwargs ) axes[0].set_ylabel("original") axes[1].set_ylabel("regridded") .. image-sg:: /user-guide/images/sphx_glr_07-regridding_012.png :alt: 07 regridding :srcset: /user-guide/images/sphx_glr_07-regridding_012.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Text(28.472222222222243, 0.5, 'regridded') .. GENERATED FROM PYTHON SOURCE LINES 372-447 A note on regridding conductivity --------------------------------- By default, K and K22 are regrid with a different method than K33. Namely K and K22 are regridded by default with a geometric mean (Bierkens, 1998), whereas K33 is regrid with an arithmetic mean. Note that default regridding methods were chosen assuming that K and K22 are roughly horizontal and K33 roughly vertical. But this may not be the case if the input arrays angle2 and/or angle3 have large values. Furthermore it is possible to only provide one array for the hydraulic conductivity, K. In that case, the same k values are used in all directions (K, K22, and K33), and the model is called isotropic. It is recommended to introduce K33 as a separate array in the source model even if it is isotropic when regridding simulations, as K33 is regridded with a different method. Reference: Bierkens, M., van der Gaast, J. Upscaling hydraulic conductivity: theory and examples from geohydrological studies. *Nutrient Cycling in Agroecosystems* **50**, 193–207 (1998). https://doi.org/10.1023/A:1009740328153 Regridding boundary conditions ------------------------------ Special care must be taken when regridding boundary conditions, and it is recommended that users verify the balance output of a regridded simulation and compare it to the original model. If the regridded simulation is a good representation of the original simulation, the mass contributions on the balance by the different boundary conditions should be comparable in both simulations. To achieve this, it may be necessary to tweak the input or the regridding methods. An example of this is upscaling recharge (so the target grid has coarser cells than the source grid). Its default method is averaging, with the following rules: - if a cell in the source grid is inactive in the source recharge package (meaning no recharge), it will not count when averaging. So if a target cell has partial overlap with one source recharge cell, and the rest of the target cell has no overlap with any source active recharge cell, it will get the recharge of the one cell it has overlap with. But since the target cell is larger, this effectively means the regridded recharge will be more in the regridded simulation than it was in the source simulation - but we do the same regridding this time assigning a zero recharge to cells without recharge then the averaging will take the zero-recharge cells into account and the regridded recharge will be the same as the source recharge. A note on regridding transport ------------------------------ Transport simulations can be unstable if constraints related to the grid Peclet number and the courant number are exceeded. This can easily happen when regridding. It may be necessary to reduce the simulation's time step size especially when downscaling, to prevent numerical issues. Increasing dispersivities or the molecular diffusion constant can also help to stabilize the simulation. Inversely, when upscaling, a larger time step size can be acceptable. Unsupported packages -------------------- Some packages cannot be regridded. This includes the Lake package and the UZF package. Such packages should be removed from the simulation before regridding, and then new packages should be created by the user and then added to the regridded simulation. Listing all default regridding methods -------------------------------------- MODFLOW 6 ^^^^^^^^^ This code snippet prints all default methods: .. GENERATED FROM PYTHON SOURCE LINES 447-502 .. code-block:: Python import inspect import sys from dataclasses import asdict import pandas as pd def collect_regrid_methods(classes: list): """Collect all regrid methods from all list of package classes.""" regrid_method_setup = { "package name": [], "array name": [], "method name": [], "function name": [], } regrid_method_table = pd.DataFrame(regrid_method_setup) counter = 0 for obj in classes: if hasattr(obj, "_regrid_method"): package_name = obj.__name__ regrid_methods = asdict(obj.get_regrid_methods()) for array_name in regrid_methods.keys(): method_name = regrid_methods[array_name][0].name function_name = "" if len(regrid_methods[array_name]) > 0: function_name = regrid_methods[array_name][1] regrid_method_table.loc[counter] = ( package_name, array_name, method_name, function_name, ) counter = counter + 1 # Set multi index to group with packages regrid_method_table = regrid_method_table.set_index(["package name", "array name"]) return regrid_method_table # Get all classes in the imod.mf6 module (e.g. # :class:`imod.mf6.NodePropertyFlow`, :class:`imod.mf6.GroundwaterFlowModel`, # :class:`imod.mf6.River`) mf6_classes = [ obj for _, obj in inspect.getmembers(sys.modules["imod.mf6"], inspect.isclass) ] mf6_regrid_methods = collect_regrid_methods(mf6_classes) # Pandas by default displays at max 60 rows, which this table exceeds. nrows = mf6_regrid_methods.shape[0] # Configure pandas to increase the max amount of displayable rows. pd.set_option("display.max_rows", nrows + 1) # Display rows: mf6_regrid_methods .. raw:: html
method name function name
package name array name
ConstantHead head OVERLAP mean
concentration OVERLAP mean
ibound OVERLAP mode
Dispersion diffusion_coefficient OVERLAP mean
longitudinal_horizontal OVERLAP mean
transversal_horizontal1 OVERLAP mean
longitudinal_vertical OVERLAP mean
transversal_horizontal2 OVERLAP mean
transversal_vertical OVERLAP mean
Drainage elevation OVERLAP mean
conductance RELATIVEOVERLAP conductance
concentration OVERLAP mean
Evapotranspiration surface OVERLAP mean
rate OVERLAP mean
depth OVERLAP mean
proportion_rate OVERLAP mean
proportion_depth OVERLAP mean
GeneralHeadBoundary head OVERLAP mean
conductance RELATIVEOVERLAP conductance
concentration OVERLAP mean
InitialConditions start OVERLAP mean
MobileStorageTransfer porosity OVERLAP mean
decay OVERLAP mean
decay_sorbed OVERLAP mean
bulk_density OVERLAP mean
distcoef OVERLAP mean
sp2 OVERLAP mean
NodePropertyFlow icelltype OVERLAP mode
k OVERLAP geometric_mean
k22 OVERLAP geometric_mean
k33 OVERLAP mean
angle1 OVERLAP mean
angle2 OVERLAP mean
angle3 OVERLAP mean
rewet_layer OVERLAP mean
Recharge rate OVERLAP mean
concentration OVERLAP mean
River stage OVERLAP mean
conductance RELATIVEOVERLAP conductance
bottom_elevation OVERLAP mean
concentration OVERLAP mean
infiltration_factor OVERLAP mean
SpecificStorage convertible OVERLAP mode
specific_storage OVERLAP mean
specific_yield OVERLAP mean
StorageCoefficient convertible OVERLAP mode
storage_coefficient OVERLAP mean
specific_yield OVERLAP mean
StructuredDiscretization top OVERLAP mean
bottom OVERLAP mean
idomain OVERLAP mode
VerticesDiscretization top OVERLAP mean
bottom OVERLAP mean
idomain OVERLAP mode


.. GENERATED FROM PYTHON SOURCE LINES 503-507 MetaSWAP ^^^^^^^^^ Let's list all default regridding methods for all MetaSWAP packages: .. GENERATED FROM PYTHON SOURCE LINES 508-517 .. code-block:: Python msw_classes = [ obj for _, obj in inspect.getmembers(sys.modules["imod.msw"], inspect.isclass) ] msw_regrid_methods = collect_regrid_methods(msw_classes) msw_regrid_methods .. raw:: html
method name function name
package name array name
GridData area RELATIVEOVERLAP conductance
landuse OVERLAP mode
rootzone_depth OVERLAP mean
surface_elevation OVERLAP mean
soil_physical_unit OVERLAP mode
active OVERLAP mode
IdfMapping area RELATIVEOVERLAP conductance
Infiltration infiltration_capacity OVERLAP mean
downward_resistance OVERLAP mean
upward_resistance OVERLAP mean
longitudinal_vertical OVERLAP mean
bottom_resistance OVERLAP mean
extra_storage_coefficient OVERLAP mean
MeteoGrid precipitation OVERLAP mean
evapotranspiration OVERLAP mean
Ponding ponding_depth OVERLAP mean
runon_resistance OVERLAP mean
runoff_resistance OVERLAP mean
ScalingFactors scale_soil_moisture OVERLAP mean
scale_hydraulic_conductivity OVERLAP mean
scale_pressure_head OVERLAP mean
depth_perched_water_table OVERLAP mean
Sprinkling max_abstraction_groundwater OVERLAP mean
max_abstraction_surfacewater OVERLAP mean


.. rst-class:: sphx-glr-timing **Total running time of the script:** (1 minutes 10.381 seconds) .. _sphx_glr_download_user-guide_07-regridding.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 07-regridding.ipynb <07-regridding.ipynb>` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 07-regridding.py <07-regridding.py>` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: 07-regridding.zip <07-regridding.zip>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_