Refactor time handling utilities now that we have pandas 3

[JoerivanEngelen]

Pandas 2 takes nanosecond by default as its base for datetimes. This means that the years that can stored in this format are limited between 1678 and 2261. We ran into these limits more often than one would expect beforehand, and thus have spent quite some time in the past adding functionality to work around this, to check for this, or to deal with this. For example, including support for cftime.

Pandas 3 uses microsecond by default as its base for datetimes, which has a year range roughly inbetween ~290,000 BC to ~290,000 AD. iMOD5 supports years ranging from 1900 AD up to 9999 AD, so well within limits of what is supported by a microseconds base. For paleo-reconstructions where million years are simulated, this still isn't enough, though I don't foresee that usecase any time soon. In short: Our lives are made a lot easier now.

If we drop support of pandas 2, we can simplify the codebase as follows:

• the calls to _check_year need to be removed. These now block users from setting up simulations beyond the year 2261 with pandas 3.
• support for use_cftime can be dropped, reducing quite some if...else clauses

Before we drop this stuff, we will first enter a transition phase where the following is required:

• add some logic to _check_year that accounts for the base of the datetime that is used, so that users using pandas 3.0 don't run into errors and throws a warning in case pandas 2 is used that behavior will change.
• throw a DeprecationWarning whenever use_cftime is used and some instructions how to migrate to pandas 3.0

[Huite]

That's excellent news!
I think this mostly (exclusively?) affects the iMOD formats? Because that's where we do the time parsing.

Ideally, we'd be able to use cftime in a MODFLOW 6 setting. But that should be just a case of the user providing CFTimeIndex objects in their xarray objects (not sure about wells).

[JoerivanEngelen]

Ideally, we'd be able to use cftime in a MODFLOW 6 setting. But that should be just a case of the user providing CFTimeIndex objects in their xarray objects (not sure about wells).

What would be the usecases for providing an CFTimeIndex compared to a pandas datetime with microseconds as base? Only things I can think of:

• You want to conduct a simulation beyond the years ~290,000 BC to ~290,000 AD. I think this is a very rare research question. Only for this research have I ever gone before 290,000 BC, for my other paleo-reconstructions the 290,000 BC would have been sufficient. A glacial cycle lasts roughly 125,000 years, so you can capture 2 cycles with this range, which is sufficient for most fresh-salt problems.
• You are forced to use some obscure calendar and don't want headaches converting datetimes yourself.

I don't think these 2 usecases would warrant the extra effort of supporting CFTimeIndex as well, but perhaps I'm overestimating the amount of work required to support CFTimeIndex. Or I've missed an important usecase.

[Huite]

Yes, I think paleo reconstructions would be the (only) application. Climate models routinely use 360_day, noleap, or all_leap calendars -- hence the existence of cftime package in the first place! I do foresee cases where we use climate forcing and having to convert between calendars is quite error prone. I agree the use case is rare. But most of the pain was caused by parsing and converting the iMOD formats date formats, even though it was a single (truncated) calendar, because of the nanosecond resolution. Supporting MF6 & cftime seems very little effort and might save someone significant hassle or confusion in the future.

The effort should be minimal since xarray is the data structure and we're looking either at a pandas TimeIndex or a CFTimeIndex. As far as I can tell, we do very little processing (fortunately):

def timestep_duration(times: np.ndarray, use_cftime: bool):
...
    for start, end in zip(times[:-1], times[1:]):
        timedelta = end - start
        duration = timedelta.days + timedelta.seconds / 86400.0
        timestep_duration.append(duration)
...

This works because a cftime difference also returns a datetime.timedelta, so actually the use_cftime check can be removed.

It's called here:

    @standard_log_decorator()
    def create_time_discretization(self, additional_times, validate: bool = True):
        """
        ...
        """
        self.use_cftime = any(
            model._use_cftime()
            for model in self.values()
            if isinstance(model, Modflow6Model)
        )

        times = [
            imod.util.time.to_datetime_internal(time, self.use_cftime)
            for time in additional_times
        ]
        for model in self.values():
            if isinstance(model, Modflow6Model):
                times.extend(model._yield_times())

        # np.unique also sorts
        times = np.unique(np.hstack(times))

        duration = imod.util.time.timestep_duration(times, self.use_cftime)
        # Generate time discretization, just rely on default arguments
        # Probably won't be used that much anyway?
        timestep_duration = xr.DataArray(
            duration, coords={"time": np.array(times)[:-1]}, dims=("time",)
        )
        self["time_discretization"] = imod.mf6.TimeDiscretization(
            timestep_duration=timestep_duration, validate=validate
        )

This can likely be simplified. I think there needs only be a check whether the date time formats are consistent across the models and datasets. The model._use_cftime does the right thing:

    def _use_cftime(self):
        """
        Also checks if datetime types are homogeneous across packages.
        """
        types = [
            type(pkg.dataset["time"].values[0])
            for pkg in self.values()
            if "time" in pkg.dataset.coords
        ]
        set_of_types = set(types)
        # Types will be empty if there's no time dependent input
        if len(set_of_types) == 0:
            return False
        else:  # there is time dependent input
            if not len(set_of_types) == 1:
                raise ValueError(
                    f"Multiple datetime types detected: {set_of_types}"
                    "Use either cftime or numpy.datetime64[ns]."
                )
            # Since we compare types and not instances, we use issubclass
            if issubclass(types[0], cftime.datetime):
                return True
            elif issubclass(types[0], np.datetime64):
                return False
            else:
                raise ValueError("Use either cftime or numpy.datetime64[ns].")

Though perhaps this could be encoded in some Schema, I don't know.

The isodate should work too, as is:

def iso8601(datetime):
    datetype = type(datetime)
    if issubclass(datetype, np.datetime64):
        return np.datetime_as_string(datetime)
    elif issubclass(datetype, cftime.datetime):
        return datetime.isoformat()
    else:
        raise TypeError(f"Expected np.datetime64 or cftime.datetime, got {datetype}")

Likely the only real change required is here, in assigning a date time index to the output:

def assign_datetime_coords(
    da: GridDataArray,
    simulation_start_time: np.datetime64,
    time_unit: Optional[str] = "d",
) -> GridDataArray:
    if "time" not in da.coords:
        raise ValueError(
            "cannot convert time column, because a time column could not be found"
        )

    time = pd.Timestamp(simulation_start_time) + pd.to_timedelta(
        da["time"], unit=time_unit
    )
    return da.assign_coords(time=time)

As this explicitly calls pd.Timestamps while we might dispatch on the start time instead.

Basically my reasoning is:

• The iMOD format logic can be simplified (no need to check years, etc.)
• The MODFLOW 6 timedis logic can be simplified too, the date time APIs are sufficiently alike.
• A typecheck via Schema for the time axis may be worthwhile: should be a datetype (pandas, numpy, cftime are the options, I guess? maybe pandas and numpy are the same?)
• Model already checks for consistent types; simulation should do the same. Maybe pull it into a separate utility.
• Extra dispatch in output handling: None -> nothing happens, cftime -> CFTimeCalendar or appropriate type, everything else defaults to pandas (e.g. strings).

So yes, if you're touching it anyway, it doesn't seem like a lot of work to allow cftime as well.