.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "user-guide\00-index-examples.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_00-index-examples.py: Quick start =========== These are the examples in the index page of the iMOD Python user guide. The examples are used to demonstrate the features of iMOD Python and how to use them. This file is meant to test if the index page is working correctly and to provide a quick overview of the examples available in the user guide. .. GENERATED FROM PYTHON SOURCE LINES 12-13 Let's first create some raster data to work with. .. GENERATED FROM PYTHON SOURCE LINES 13-26 .. code-block:: Python import xarray as xr import xugrid as xu import imod elevation_uda = xu.data.elevation_nl() # Drop unnecessary coords. These coords are also stored in elevation_uda.ugrid.grid elevation_uda = elevation_uda.drop_vars(["mesh2d_face_x", "mesh2d_face_y"]) part = elevation_uda.ugrid.sel( x=slice(100_000.0, 200_000.0), y=slice(400_000.0, 500_000.0) ) .. GENERATED FROM PYTHON SOURCE LINES 27-28 Make structured grid .. GENERATED FROM PYTHON SOURCE LINES 28-39 .. code-block:: Python resolution = 500.0 xmin, ymin, xmax, ymax = part.ugrid.grid.bounds structured_grid = imod.util.empty_2d(resolution, xmin, xmax, resolution, ymin, ymax) # Interpolate regridder = xu.BarycentricInterpolator(part, structured_grid) interpolated_elevation = regridder.regrid(part) # Plot interpolated_elevation.plot.imshow() .. image-sg:: /user-guide/images/sphx_glr_00-index-examples_001.png :alt: dy = -500.0, dx = 500.0 :srcset: /user-guide/images/sphx_glr_00-index-examples_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 40-41 Save to GeoTIFF .. GENERATED FROM PYTHON SOURCE LINES 41-47 .. code-block:: Python # Create a temporary directory to store the data tmp_dir = imod.util.path.temporary_directory() tmp_dir.mkdir(parents=True, exist_ok=True) interpolated_elevation.rio.to_raster(tmp_dir / "elevation.tif") .. GENERATED FROM PYTHON SOURCE LINES 48-51 Seamlessly integrate your GIS rasters or meshes with MODFLOW 6, by using `xarray`_ and `xugrid`_ arrays, for structured and unstructured grids respectively, to create grid-based model packages. .. GENERATED FROM PYTHON SOURCE LINES 51-73 .. code-block:: Python # Open Geotiff with elevation data as xarray DataArray elevation = imod.rasterio.open(tmp_dir / "elevation.tif") elevation.load() # Create idomain grid layer_template = xr.DataArray([1, 1, 1], dims=("layer",), coords={"layer": [1, 2, 3]}) idomain = layer_template * xr.ones_like(elevation).astype(int) # Deactivate cells with NoData idomain = idomain.where(elevation.notnull(), 0) # Compute bottom elevations of model layers layer_thickness = xr.DataArray( [10.0, 20.0, 10.0], dims=("layer",), coords={"layer": [1, 2, 3]} ) bottom = elevation - layer_thickness.cumsum(dim="layer") # Create MODFLOW 6 DIS package dis_pkg = imod.mf6.StructuredDiscretization( idomain=idomain, top=elevation, bottom=bottom.transpose("layer", "y", "x") ) .. GENERATED FROM PYTHON SOURCE LINES 74-76 Assign wells based on x, y coordinates and filter screen depths, instead of layer, row and column: .. GENERATED FROM PYTHON SOURCE LINES 76-113 .. code-block:: Python import pandas as pd # Specify well locations x = [150_200.0, 160_800.0] y = [450_300.0, 460_200.0] # Specify well screen depths screen_top = [0.0, 0.0] screen_bottom = [-4.0, -10.0] # Specify flow rate, which changes over time. weltimes = pd.date_range("2000-01-01", "2000-01-03", freq="2D") well_rates_period1 = [0.5, 1.0] well_rates_period2 = [2.5, 3.0] rate = xr.DataArray( [well_rates_period1, well_rates_period2], coords={"time": weltimes}, dims=("time", "index"), ) # Now construct the Well package wel_pkg = imod.mf6.Well( x=x, y=y, rate=rate, screen_top=screen_top, screen_bottom=screen_bottom ) # iMOD Python will take care of the rest and assign the wells to the correct model # layers upon writing the model. It will furthermore distribute well rates based # on transmissivities. To verify how wells will be assigned to model layers before # writing the entire simulation, you can use the following command: wel_mf6_pkg = wel_pkg.to_mf6_pkg(idomain, elevation, bottom, k=1.0) # Wells have been distributed across two model layers based on screen depths. print(wel_mf6_pkg["cellid"]) # Well rates have been distributed based on screen overlap print(wel_mf6_pkg["rate"]) .. rst-class:: sphx-glr-script-out .. code-block:: none Size: 96B array([[ 1, 104, 106], [ 1, 84, 127], [ 2, 104, 106], [ 2, 84, 127]]) Coordinates: * ncellid (ncellid) int64 32B 0 1 2 3 * dim_cellid (dim_cellid) Size: 64B array([[0.28073287, 1.40366435], [0.58867207, 1.7660162 ], [0.21926713, 1.09633565], [0.41132793, 1.2339838 ]]) Coordinates: * ncellid (ncellid) int64 32B 0 1 2 3 x (ncellid) float64 32B 1.502e+05 1.608e+05 1.502e+05 1.608e+05 y (ncellid) float64 32B 4.503e+05 4.602e+05 4.503e+05 4.602e+05 * time (time) datetime64[ns] 16B 2000-01-01 2000-01-03 .. GENERATED FROM PYTHON SOURCE LINES 114-118 MODFLOW 6 requires that all stress periods are defined in the time discretization package. However, usually boundary conditions are defined at insconsistent times. iMOD Python can help you to create a time discretization package that is consistent, based on all the unique times assigned to the boundary conditions. .. GENERATED FROM PYTHON SOURCE LINES 118-132 .. code-block:: Python simulation = imod.mf6.Modflow6Simulation("example") simulation["gwf"] = imod.mf6.GroundwaterFlowModel() simulation["gwf"]["dis"] = dis_pkg simulation["gwf"]["wel"] = wel_pkg # Create a time discretization based on the times assigned to the packages. # Specify the end time of the simulation as one of the additional_times simulation.create_time_discretization(additional_times=["2000-01-07"]) # Note that timesteps in well package are also inserted in the time # discretization print(simulation["time_discretization"].dataset) .. rst-class:: sphx-glr-script-out .. code-block:: none Size: 48B Dimensions: (time: 2) Coordinates: * time (time) datetime64[ns] 16B 2000-01-01 2000-01-03 Data variables: timestep_duration (time) float64 16B 2.0 4.0 n_timesteps int64 8B 1 timestep_multiplier float64 8B 1.0 .. GENERATED FROM PYTHON SOURCE LINES 133-138 Regrid MODFLOW 6 models to different grids, even from structured to unstructured grids. iMOD Python takes care of properly scaling the input parameters. You can also configure scaling methods yourself for each input parameter, for example when you want to upscale drainage elevations with the minimum instead of the average. .. GENERATED FROM PYTHON SOURCE LINES 138-145 .. code-block:: Python new_unstructured_grid = part sim_regridded = simulation.regrid_like("regridded", new_unstructured_grid) # Notice that discretization has converted to VerticesDiscretization (DISV) print(sim_regridded["gwf"]["dis"]) .. rst-class:: sphx-glr-script-out .. code-block:: none VerticesDiscretization Size: 79kB Dimensions: (layer: 3, mesh2d_nFaces: 1238) Coordinates: * layer (layer) int64 24B 1 2 3 * mesh2d_nFaces (mesh2d_nFaces) int64 10kB 0 1 2 3 4 ... 1234 1235 1236 1237 Data variables: idomain (layer, mesh2d_nFaces) int64 30kB 1 1 1 1 1 1 ... 1 1 1 1 1 1 top (mesh2d_nFaces) float64 10kB 7.708 18.72 ... 12.31 11.52 bottom (layer, mesh2d_nFaces) float64 30kB -2.292 8.724 ... -28.48 .. GENERATED FROM PYTHON SOURCE LINES 146-148 To reduce the size of your model, you can clip it to a bounding box. This is useful for example when you want to create a smaller model for testing purposes. .. GENERATED FROM PYTHON SOURCE LINES 148-151 .. code-block:: Python sim_clipped = simulation.clip_box( x_min=125_000, x_max=175_000, y_min=425_000, y_max=475_000 ) .. GENERATED FROM PYTHON SOURCE LINES 152-155 You can even provide states for the model, which will be set on the model boundaries of the clipped model. Create a grid of zeros, which will be used to set as heads at the boundaries of clipped parts. .. GENERATED FROM PYTHON SOURCE LINES 155-170 .. code-block:: Python head_for_boundary = xr.zeros_like(idomain, dtype=float).where(idomain > 0) states_for_boundary = {"gwf": head_for_boundary} sim_clipped = simulation.clip_box( x_min=125_000, x_max=175_000, y_min=425_000, y_max=475_000, states_for_boundary=states_for_boundary, ) # Notice that a Constant Head (CHD) package has been created for the clipped # model. print(sim_clipped["gwf"]) .. rst-class:: sphx-glr-script-out .. code-block:: none GroundwaterFlowModel( listing_file=None, print_input=False, print_flows=False, save_flows=False, newton=False, under_relaxation=False, ){ 'dis': StructuredDiscretization, 'wel': Well, 'chd_clipped': ConstantHead, } .. GENERATED FROM PYTHON SOURCE LINES 171-173 .. _xarray: http://xarray.pydata.org/ .. _xugrid: https://deltares.github.io/xugrid/ .. rst-class:: sphx-glr-timing **Total running time of the script:** (1 minutes 0.353 seconds) .. _sphx_glr_download_user-guide_00-index-examples.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 00-index-examples.ipynb <00-index-examples.ipynb>` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 00-index-examples.py <00-index-examples.py>` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: 00-index-examples.zip <00-index-examples.zip>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_