Issue #1768 pandas3 support (#1771)

Fixes #1768 

# Description
Fixes broken pieces of the code when updating to pandas 3, which [has
the following major
changes](https://pandas.pydata.org/community/blog/pandas-3.0.html) which
also affect iMOD Python somewhat:

- _Dedicated string data type by default: string columns are now
inferred as the new str dtype instead of object, providing better
performance and type safety_
- _Consistent copy/view behaviour with Copy-on-Write (CoW) (a.k.a.
getting rid of the SettingWithCopyWarning): more predictable and
consistent behavior for all operations, with improved performance
through avoiding unnecessary copies_
- _New default resolution for datetime-like data: no longer defaulting
to nanoseconds, but generally microseconds (or the resolution of the
input), when constructing datetime or timedelta data (avoiding
out-of-bounds errors for dates with a year before 1678 or after 2262)_

Because of these changes, some parts of the code had to be modified
slightly to get the tests to work. I think with the latter change we
could simplify the code base a lot, as all logic to deal with datetimes
going beyond the year 2262 and prior to 1678 wouldn't be necessary
anymore. This was quite a headache in the past, and I think the choice
for microseconds by default will make our lives significantly easier. [I
created an issue for
this](#1773).

In detail this PR alters the following:

- Pin pandas version to 3.* in pixi.toml
- Update setup and asserts in quite some tests to work with pandas 3.0.
Most of these tests will not work with pandas 2.0 anymore, as we now
assert for pandas 3.0 behavior
- Update examples to work with pandas 3.0
- Fix statements where a view was altered (e.g. altering ``.values``
attribute), this is now forbidden by pandas
- Update to new string handling, sometimes I had to internally enforce
object dtype to get the old behavior. In other cases we could check for
strings by checking ``dtype=="object"`` anymore. Luckily there is
``pd.api.types.is_string_dtype``, which both works for pandas 2 as well
as pandas 3.
- Add unittest for ``to_pandas_datetime_series`` to check whether pandas
behavior is similar in this regard.

# Checklist
<!---
Before requesting review, please go through this checklist:
-->

- [x] Links to correct issue
- [x] Update changelog, if changes affect users
- [x] PR title starts with ``Issue #nr``, e.g. ``Issue #737``
- [x] Unit tests were added
- [ ] **If feature added**: Added/extended example
- [ ] **If feature added**: Added feature to API documentation
- [ ] **If pixi.lock was changed**: Ran `pixi run generate-sbom` and
committed changes

Author: JoerivanEngelen

docs/api/changelog.rst

@@ -44,6 +44,8 @@ Fixed
 - Fixed bug where :func:`imod.evaluate.convert_pointwaterhead_freshwaterhead`
   produced incorrect results when point water heads were below elevation levels
   for unstructured grids.
+- Support pandas 3.0.
+
 
 Changed
 ~~~~~~~

examples/imod-wq/Henry-wq.py

@@ -102,7 +102,7 @@
 )
 m["oc"] = imod.wq.OutputControl(save_head_idf=True, save_concentration_idf=True)
 m.create_time_discretization(
-    additional_times=pd.date_range("2000-01-01", "2001-01-01", freq="M")
+    additional_times=pd.date_range("2000-01-01", "2001-01-01", freq="ME")
 )
 
 # %%

examples/mf6/circle_transport.py

@@ -194,7 +194,7 @@
 # stress period is a year long. We can then use the "last" keyword in the output
 # control to save the output.
 
-simtimes = pd.date_range(start="2000-01-01", end="2030-01-01", freq="As")
+simtimes = pd.date_range(start="2000-01-01", end="2030-01-01", freq="YS")
 simulation.create_time_discretization(additional_times=simtimes)
 
 # %%

examples/mf6/hondsrug.py

@@ -336,7 +336,7 @@ def outer_edge(da):
 # resampled to a yearly step by using the xarray function
 # `Dataset.resample <http://xarray.pydata.org/en/stable/generated/xarray.Dataset.resample.html#xarray.Dataset.resample>`_.
 
-rch_trans_yr = rch_trans.resample(time="A", label="left").mean()
+rch_trans_yr = rch_trans.resample(time="YS", label="left").mean()
 rch_trans_yr
 
 # %%

imod/evaluate/constraints.py

@@ -126,7 +126,8 @@ def stability_constraint_advection(front, lower, right, top_bot, porosity=0.3, R
     dt = 1.0 / (1.0 / dt_x + 1.0 / dt_y + 1.0 / dt_z)
 
     dt_xyz = xr.concat(
-        (dt_x, dt_y, dt_z), dim=pd.Index(["x", "y", "z"], name="direction")
+        (dt_x, dt_y, dt_z),
+        dim=pd.Index(["x", "y", "z"], name="direction", dtype="object"),
     )
     return dt, dt_xyz
 
@@ -254,9 +255,11 @@ def _get_stage_name(sid):
                 if not drop_allnan or not dt.isnull().all():
                     results.append(dt)
                     resultids.append(comb)
-    dt_all = xr.concat(
-        results, pd.Index(resultids, name="combination"), coords="minimal"
-    )
+    # Set index to object dtype to work around xarray concat issue where
+    # StringDtype could not be interpreted as a data type with pandas 3.0 (as
+    # np.dtype is called.)
+    id_index = pd.Index(resultids, name="combination", dtype="object")
+    dt_all = xr.concat(results, id_index, coords="minimal")
 
     # overall dt
     dt_min = dt_all.min(dim="combination")

imod/formats/gen/gen.py

@@ -267,7 +267,7 @@ def read_binary(path: Union[str, Path]) -> "geopandas.GeoDataFrame":  # type: ig
     else:
         df = pd.DataFrame()
     df["feature_type"] = feature_type
-    df["feature_type"] = df["feature_type"].replace(GENTYPE_TO_NAME)
+    df["feature_type"] = df["feature_type"].replace(GENTYPE_TO_NAME).astype(str)
 
     geometry = []
     for ftype, geom in zip(feature_type, xy):

imod/formats/ipf.py

@@ -458,7 +458,8 @@ def write_assoc(path, df, itype=1, nodata=1.0e20, assoc_columns=None):
     # The reason is that datetime columns are converted to string as well
     # and then quoted. This causes trouble with some iMOD(batch) functions.
     for column in df.columns:
-        if df.loc[:, column].dtype == np.dtype("O"):
+        # Test for strings compatible with pandas 2 and 3
+        if pd.api.types.is_string_dtype(df[column].dtype):
             df.loc[:, column] = df.loc[:, column].astype(str)
             df.loc[:, column] = '"' + df.loc[:, column] + '"'
 
@@ -509,7 +510,8 @@ def write(path, df, indexcolumn=0, assoc_ext="txt", nodata=1.0e20):
     # The reason is that datetime columns are converted to string as well
     # and then quoted. This causes trouble with some iMOD(batch) functions.
     for column in df.columns:
-        if df.loc[:, column].dtype == np.dtype("O"):
+        # Test for strings compatible with pandas 2 and 3
+        if pd.api.types.is_string_dtype(df[column].dtype):
             df.loc[:, column] = df.loc[:, column].astype(str)
             df.loc[:, column] = '"' + df.loc[:, column] + '"'
 

imod/mf6/multimodel/exchange_creator.py

@@ -290,10 +290,12 @@ def rearrange_connected_cells(self):
 
         label_decreasing = df["cell_label1"] > df["cell_label2"]
 
-        colnames = ["cell_idx1", "cell_idx2", "cell_label1", "cell_label2"]
-        colnames_reversed = ["cell_idx2", "cell_idx1", "cell_label2", "cell_label1"]
+        if label_decreasing.any():
+            colnames = ["cell_idx1", "cell_idx2", "cell_label1", "cell_label2"]
+            colnames_reversed = ["cell_idx2", "cell_idx1", "cell_label2", "cell_label1"]
 
-        decreasing_connections = df.loc[label_decreasing, colnames].values
-        df.loc[label_decreasing, colnames_reversed] = decreasing_connections
+            df_decreasing = df.loc[label_decreasing, colnames]
+            df_decreasing.columns = colnames_reversed
+            df.loc[label_decreasing, colnames_reversed] = df_decreasing
 
         self._connected_cells = df

imod/testing.py

@@ -16,7 +16,7 @@ def assert_frame_equal(left: pd.DataFrame, right: pd.DataFrame, **kwargs):
     def always_int64(df):
         df = df.copy()
         for column, dtype in df.dtypes.items():
-            if np.issubdtype(dtype, np.integer):
+            if pd.api.types.is_integer_dtype(dtype):
                 df[column] = df[column].astype(np.int64)
         return df
 

imod/tests/test_evaluate/test_constraints.py

@@ -94,13 +94,13 @@ def test_intra_cell_boundary_conditions(test_da1):
     riv1drn = test_da1 * (0.3 * 1.0) / min(100.0, 150) * (1.0 - 0.0)
     dt_min_ref = np.minimum(ghbdrn, riv1drn)
 
-    assert dt_min_ref.equals(dt_min)
-    assert dt_all.equals(
-        xr.concat(
-            (ghbdrn, riv1drn), pd.Index(["ghb-drn", "riv_0-drn"], name="combination")
-        )
+    expected_index = pd.Index(
+        ["ghb-drn", "riv_0-drn"], name="combination", dtype="object"
     )
 
+    assert dt_min_ref.equals(dt_min)
+    assert dt_all.equals(xr.concat((ghbdrn, riv1drn), expected_index))
+
 
 def test_intra_cell_boundary_conditions_thickness_zero(test_da1):
     top_bot = xr.Dataset({"top": test_da1 * -1.0, "bot": test_da1 * -1.0})