Merge branch '329-imod-visualize-cross_section-does-not-support-quasi-3d-layering' of https://gitlab.com/deltares/imod/imod-python into 329-imod-visualize-cross_section-does-not-support-quasi-3d-layering

Author: Joost Delsman

docs/api/changelog.rst

@@ -6,15 +6,93 @@ All notable changes to this project will be documented in this file.
 The format is based on `Keep a Changelog`_, and this project adheres to
 `Semantic Versioning`_.
 
+
 [Unreleased]
 ------------
 
 Fixed
 ~~~~~
 
+- :meth:`imod.mf6.Simulation._validate` now print all validation errors for all
+  models and packages in one message.
+
+
+[0.13.1] - 2023-05-05
+---------------------
+
+Added
+~~~~~
+
+- :class:`imod.mf6.SpecificStorage` and :class:`imod.mf6.StorageCoefficient` 
+  now have a ``save_flow`` argument.
+
+Fixed
+~~~~~
+
+- :func:`imod.mf6.open_cbc` can now read storage fluxes without error.
+
+
+[0.13.0] - 2023-05-02
+---------------------
+
+Added
+~~~~~
+
+- :class:`imod.mf6.OutputControl` now takes parameters ``head_file``,
+  ``concentration_file``, and ``budget_file`` to specify where to store
+  MODFLOW6 output files.
+- :func:`imod.util.from_mdal_compliant_ugrid2d` to "restack" the variables that
+  have have been "unstacked" in :func:`imod.util.mdal_compliant_ugrid2d`.
+- Added support for the Modflow6 Lake package
+- :func:`imod.select.points_in_bounds`, :func:`imod.select.points_indices`,
+  :func:`imod.select.points_values` now support unstructured grids.
+- Added support for the Modflow 6 Lake package: :class:`imod.mf6.Lake`,
+  :class:`imod.mf6.LakeData`, :class:`imod.mf6.OutletManning`, :class:`OutletSpecified`,
+  :class:`OutletWeir`. See the examples for an application of the Lake package.
+
+Fixed
+~~~~~
+
+- :meth:`imod.wq.bas.BasicFlow.thickness` returns a DataArray with the correct
+  dimension order again. This confusingly resulted in an error when writing the
+  :class:`imod.wq.btn.BasicTransport` package.
+- Fixed bug in :class:`imod.mf6.dis.StructuredDiscretization` and
+  :class:`imod.mf6.dis.VerticesDiscretization` where 
+  ``inactive bottom above active cell`` was incorrectly raised.
+
+[0.12.0] - 2023-03-17
+---------------------
+
+Added
+~~~~~
+
+- :func:`imod.prj.read_projectfile` to read the contents of a project file into
+  a Python dictionary.
+- :func:`imod.prj.open_projectfile_data` to read/open the data that is pointed
+  to in a project file.
+- :func:`imod.gen.read_ascii` to read the geometry stored in ASCII text .gen files.
+- :class:`imod.mf6.hfb.HorizontalFlowBarrier` to support Modflow6's HFB
+  package, works well with `xugrid.snap_to_grid` function.
+
+Fixed
+~~~~~
+
 - :func:`imod.evaluate.budget.flow_velocity` now properly computes velocity by
   dividing by the porosity. Before, this function computed the Darcian velocity.
 
+Changed
+~~~~~~~
+
+- :func:`imod.ipf.save` will error on duplicate IDs for associated files if a
+  ``"layer"`` column is present. As a dataframe is automatically broken down
+  into a single IPF per layer, associated files for the first layer would be
+  overwritten by the second, and so forth.
+- :meth:`imod.wq.Well.save` will now write time varying data to associated
+  files for extration rate and concentration.
+- Choosing ``method="geometric_mean"`` in the Regridder will now result in NaN
+  values in the regridded result if a geometric mean is computed over negative
+  values; in general, a geometric mean should only be computed over physical
+  quantities with a "true zero" (e.g. conductivity, but not elevation).
 
 [0.11.6] - 2023-02-01
 ---------------------
@@ -27,7 +105,6 @@ Added
   which can be used to provide keyword arguments to the writing of the Modflow 6
   Simulation.
 
-
 Fixed
 ~~~~~
 
@@ -42,8 +119,8 @@ Fixed
 
 - :meth:`imod.mf6.Modflow6Simulation.write` with ``binary=False`` no longer
   results in invalid MODFLOW6 input for 2D grid data, such as DIS top.
-- :meth:`imod.flow.ImodflowModel.write` no longer writes incorrect projectfiles
-  for non-grid values with a time and layer dimension.
+- :meth:`imod.flow.ImodflowModel.write` no longer writes incorrect project
+  files for non-grid values with a time and layer dimension.
 - :func:`imod.evaluate.interpolate_value_boundaries`: Fix edge case when
   successive values in z direction are exactly equal to the boundary value.
 

docs/api/io.rst

@@ -1,4 +1,4 @@
-.. currentmodule:: imod
+.. currentmodule:: imod.formats
 
 Input/output
 ------------
@@ -24,3 +24,9 @@ Input/output
 
     gen.read
     gen.write
+    gen.read_ascii
+    
+    prj.read_projectfile
+    prj.open_projectfile_data
+    prj.read_timfile
+    prj.convert_to_disv

docs/api/mf6.rst

@@ -57,6 +57,9 @@ Flow Packages
     Drainage
     Evapotranspiration
     GeneralHeadBoundary
+    HorizontalFlowBarrierHydraulicCharacteristic
+    HorizontalFlowBarrierMultiplier
+    HorizontalFlowBarrierResistance
     InitialConditions
     NodePropertyFlow
     Recharge

docs/faq/how-do-i/modification.rst

@@ -148,7 +148,29 @@ Generally, raster data is y-descending, so ``ymax`` comes before ``ymin``:
 .. code-block:: python
 
     da_selection = da.sel(x=slice(xmin, xmax), y=slice(ymax, ymin))
-  
+    
+Get a single value (e.g. a summary statistic)
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+When computing an aggregation such as a mean, xarray will return a single
+valued DataArray. To get an ordinary Python scalar (for example a single
+integer or float), use ``.item()``:
+
+.. code-block:: python
+   
+    single_value = da.mean().item()
+ 
+Increase the extent of a raster
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Use another raster with appropriate extent, and use ``align``:
+
+.. code-block:: python
+
+    small_aligned, big_aligned = xr.align(small, big, join="outer")
+    
+Make sure the cell size is the same, or the result will be non-equidistant.
+
 Create an empty raster
 ~~~~~~~~~~~~~~~~~~~~~~
 

docs/installation.rst

@@ -63,7 +63,7 @@ Spatial operations:
 Geospatial libaries (optional):
 
 * `geopandas <https://geopandas.org/en/stable/>`__
-* `pygeos <https://pygeos.readthedocs.io/en/stable/>`__
+* `shapely <https://shapely.readthedocs.io/en/stable/manual.html>`__
 * `pyproj <https://pyproj4.github.io/pyproj/stable/>`__
 * `rasterio <https://rasterio.readthedocs.io/en/latest/>`__
 

examples/metaswap/metaswap_example.py

@@ -102,7 +102,7 @@
 # Initial conditions
 # ``````````````````
 
-gwf_model["ic"] = mf6.InitialConditions(head=0.5)
+gwf_model["ic"] = mf6.InitialConditions(start=0.5)
 
 # %%
 # Output Control

examples/mf6/Henry.py

@@ -105,7 +105,7 @@
     transient=False,
     convertible=0,
 )
-gwf_model["ic"] = imod.mf6.InitialConditions(head=0.0)
+gwf_model["ic"] = imod.mf6.InitialConditions(start=0.0)
 gwf_model["oc"] = imod.mf6.OutputControl(save_head="last", save_budget="last")
 
 # %%

examples/mf6/circle.py

@@ -114,7 +114,7 @@
 gwf_model["chd"] = imod.mf6.ConstantHead(
     constant_head, print_input=True, print_flows=True, save_flows=True
 )
-gwf_model["ic"] = imod.mf6.InitialConditions(head=0.0)
+gwf_model["ic"] = imod.mf6.InitialConditions(start=0.0)
 gwf_model["npf"] = imod.mf6.NodePropertyFlow(
     icelltype=icelltype,
     k=k,
@@ -218,7 +218,7 @@
 
 ds = ds.ugrid.assign_edge_coords()
 fig, ax = plt.subplots()
-head.isel(time=0, layer=0).ugrid.plot(ax=ax)
+head.isel(time=0, layer=0).compute().ugrid.plot(ax=ax)
 ds.isel(time=0, layer=0).plot.quiver(
     x="mesh2d_edge_x", y="mesh2d_edge_y", u="u", v="v", color="white"
 )
@@ -227,3 +227,28 @@
 # %%
 # As would be expected from our model input, we observe circular groundwater
 # mounding and increasing flows as we move from the center to the exterior.
+
+# %%
+# Slice the model domain
+# ----------------------
+#
+# We may also quickly setup a smaller model. We'll select half of the original
+# domain.
+
+half_simulation = simulation.clip_box(x_max=0.0)
+
+# %%
+# Let's run the model, read the results, and visualize.
+
+modeldir = imod.util.temporary_directory()
+half_simulation.write(modeldir)
+half_simulation.run()
+head = imod.mf6.out.open_hds(
+    modeldir / "GWF_1/GWF_1.hds",
+    modeldir / "GWF_1/disv.disv.grb",
+)
+
+fig, ax = plt.subplots()
+head.isel(time=0, layer=0).compute().ugrid.plot(ax=ax)
+ax.set_aspect(1)
+# %%

examples/mf6/circle_transport.py

@@ -152,7 +152,7 @@
     print_flows=True,
     save_flows=True,
 )
-gwf_model["ic"] = imod.mf6.InitialConditions(head=0.0)
+gwf_model["ic"] = imod.mf6.InitialConditions(start=0.0)
 gwf_model["npf"] = imod.mf6.NodePropertyFlow(
     icelltype=icelltype,
     k=k,

examples/mf6/ex01_twri.py

@@ -117,7 +117,7 @@
     print_flows=True,
     save_flows=True,
 )
-gwf_model["ic"] = imod.mf6.InitialConditions(head=0.0)
+gwf_model["ic"] = imod.mf6.InitialConditions(start=0.0)
 gwf_model["npf"] = imod.mf6.NodePropertyFlow(
     icelltype=icelltype,
     k=k,

examples/mf6/example_1d_transport.py

@@ -80,7 +80,7 @@ def create_transport_model(flowmodel, speciesname, dispersivity, retardation, de
         decay_sorbed=decay,
         bulk_density=rhobulk,
         distcoef=kd,
-        decay_order="first",
+        first_order_decay=True,
         sorption=sorption,
     )
 

examples/mf6/hondsrug.py

@@ -24,7 +24,7 @@
 import matplotlib.pyplot as plt
 import numpy as np
 import pandas as pd
-import scipy.ndimage.morphology
+import scipy.ndimage
 import xarray as xr
 
 import imod
@@ -208,7 +208,7 @@
 
 def outer_edge(da):
     data = da.copy()
-    from_edge = scipy.ndimage.morphology.binary_erosion(data)
+    from_edge = scipy.ndimage.binary_erosion(data)
     is_edge = (data == 1) & (from_edge == 0)
     return is_edge.astype(bool)
 
@@ -226,7 +226,7 @@ def outer_edge(da):
 # Using the previously created function and the 2d template,
 # the outer edge is defined for this example.
 
-edge = outer_edge(xr.full_like(like_2d.drop("layer"), 1))
+edge = outer_edge(xr.full_like(like_2d.drop_vars("layer"), 1))
 
 # %%
 # Adding information to the CHD package
@@ -489,15 +489,18 @@ def outer_edge(da):
     {"layer": layer},
     ("layer",),
 )
-times_sto = [
-    np.datetime64("2009-12-30T23:59:59.000000000"),
-    np.datetime64("2009-12-31T00:00:00.000000000"),
-    np.datetime64("2010-12-31T00:00:00.000000000"),
-    np.datetime64("2011-12-31T00:00:00.000000000"),
-    np.datetime64("2012-12-31T00:00:00.000000000"),
-    np.datetime64("2013-12-31T00:00:00.000000000"),
-    np.datetime64("2014-12-31T00:00:00.000000000"),
-]
+times_sto = np.array(
+    [
+        "2009-12-30T23:59:59.00",
+        "2009-12-31T00:00:00.00",
+        "2010-12-31T00:00:00.00",
+        "2011-12-31T00:00:00.00",
+        "2012-12-31T00:00:00.00",
+        "2013-12-31T00:00:00.00",
+        "2014-12-31T00:00:00.00",
+    ],
+    dtype="datetime64[ns]",
+)
 
 transient = xr.DataArray(
     [False, True, True, True, True, True, True], {"time": times_sto}, ("time",)
@@ -513,7 +516,11 @@ def outer_edge(da):
 # specific yield and if the layers are convertible.
 
 gwf_model["sto"] = imod.mf6.SpecificStorage(
-    specific_storage=ss, specific_yield=sy, transient=transient, convertible=0
+    specific_storage=ss,
+    specific_yield=sy,
+    transient=transient,
+    convertible=0,
+    save_flows=True,
 )
 
 # %%