Setting up the time axis properly

[davidorme]

@vgro and I had a chat recently about time axis dates and inputs and outputs that raised the issue of calendar aware dates in the VE see here. We've been ducking this forever but I think we're at a stage where it is finally blocking things we want to do.

Please could you have a read and leave any thoughts. I don't think this would be a difficult change - and I think it probably is on the critical path - but I want input before we turn it into an actual issue for action.

Important

Executive summary

It wouldn't be hard to pass actual calendar aware periods and update intervals to models - although the changes needed in models that already subdivide intervals internally might not be trivial. The current sticking point is that the models are not currently written to definitively return outputs that calibrate to the interval length. And fixing that is a longer issue.

---

At the moment, we're using pint to configure the run length and update interval.

• That's fine for giving us a number of intervals and the length of each interval, ...
• but it's useless for tying the simulation clearly to a calendar-aware set of dates. The pint units for years and months are "average" duration, so are floats:

In [3]: pint.Quantity(1, "month").to("days")
Out[3]: <Quantity(30.4375, 'day')>
In [4]: pint.Quantity(1, "year").to("days")
Out[4]: <Quantity(365.25, 'day')>

At the moment, we're simply ducking this and updates iterate over the expected number of intervals just by index [0, 1, ..., N].

That's a problem (with reasons this might not be urgent) because:

1. It's a horrible, clearly low quality state of affairs that we should be embarrassed by.

   That @davidorme is embarrassed by.

2. We can't define the time coordinates of the time update axis properly so can't validate the time axis across inputs.

   Well, OK, but we've survived so far without validation - as long as the number of intervals along the time axis in inputs is fine it works. Models know how long an interval is, and that is the main thing.

3. We have already had to implement hacks to work around this issue:

   virtual_ecosystem/virtual_ecosystem/models/hydrology/hydrology_model.py

   Lines 443 to 452 in 0f4f9ed

   	# Determine number of days
   	days_float: float = self.model_timing.update_interval_seconds / 86400
   	days: int = int(days_float // 1)
   	# Check if the number of days is exact and warn if not
   	if not np.allclose(days_float % 1, 0):
   	LOGGER.warning(
   	f"Update interval is not a whole number of days ({days_float}),"
   	f" partitioning inputs among {days} days."
   	)

   Sure, but those hacks work! Does a recurrent hack warning in every log matter?

4. It makes it impossible to link inputs or outputs at different temporal frequencies. If we have a central update frequency of - say - monthly values but a model wants to use daily inputs, then we have to have calendar-aware update periods to be able to tie update index [0] to '2000-01' in the monthly inputs and '2000-01-01/2000-01-31' in the daily inputs. I think the only way we could duck this would be to insist on 360 day years (12 x 30) - some climate models do this - but it's going to be a hard blocker for most users.

   We already have an current issue about outputting data that needs date aware coordinates: Write out daily hydrology values for calibration validation #1176

   This is hard to duck - it depends on whether daily inputs are on the path to the short-term minimum viable product. The hydrology model at least wants daily inputs and deciduous plants are going to need to know what the date is.

What do we care about

There are many Python packages providing "calendar aware" dates.

• Classic datetime like packages: datetime itself and improvements (Arrow, Pendulum, whenever)
• Array/time series like datetimes (numpy with np.datetime64 and then xarray and pandas which wrap more or less wrap numpy)

These differ in how fancy they get about things like leap seconds, daylight savings, timezones etc. I think we largely don't care about any of those: VE will run in small, single locations on at least yearly scales where we don't want to track saving times. So, something like pandas provides functionality that provides calendar aware arrays of date ranges without us having to reinvent any wheels.

Possible solution

This pandas code allows us to generate a sequence of calendar aware periods, from reasonably user friendly inputs. Honestly, it would be easier to configure end_time rather than run_length, but from where we are this works:

In [2]: import pandas as pd
   ...: 
   ...: start_date = "2000-02-01"
   ...: update_interval = "2M"
   ...: run_length = "2Y"
   ...: 
   ...: # Calculate the requested end time from the run length
   ...: simulation_span = pd.period_range(pd.Timestamp(start_date), periods=1, freq=run_length)
   ...: span_start_time = simulation_span[0].start_time
   ...: span_end_time = simulation_span[0].end_time
   ...: 
   ...: # Calculate the intervals
   ...: simulation_intervals = pd.period_range(
   ...:     start=span_start_time, end=span_end_time, freq=update_interval
   ...: )
   ...: 
   ...: # Get the interval durations
   ...: interval_starts = simulation_intervals.start_time.to_numpy().astype("datetime64[s]")
   ...: interval_ends = simulation_intervals.end_time.to_numpy().astype("datetime64[s]")
   ...: interval_durations = (interval_ends - interval_starts) + 1
   ...: 
In [3]: interval_durations
Out[3]: 
array([5184000, 5270400, 5270400, 5356800, 5270400, 5270400, 5097600,
       5270400, 5270400, 5356800, 5270400, 5270400],
      dtype='timedelta64[s]')
In [4]: simulation_intervals[0]
Out[4]: Period('2000-01', '2M')
In [5]: simulation_intervals[0].start_time
Out[5]: Timestamp('2000-01-01 00:00:00')
In [6]: simulation_intervals[0].end_time
Out[6]: Timestamp('2000-02-29 23:59:59.999999999')

Explanations

• This is using pandas.Period not pandas.Timestamp. That's partly because the data do represent periods not point values, but also because the timestamp implementation is edging back towards the fancier end of datetime handling with all sorts of business days, timezone etc. malarky.

• The pandas.Period implementation is a lot simpler - each item is explicitly a set number of *days, weeks, months or years (I think it would handle faster frequencies fine, but we aren't going there anytime soon). The possible letter codes we'd want to use are therefore nD, nM and nY for days, months and years (see https://pandas.pydata.org/pandas-docs/stable/user_guide/timeseries.html#period-aliases for details) and those seem pretty user friendly.

• We might have wanted weeks but pandas anchors weeks to a calendar week with a particular start day (e.g Monday) rather than just 7 days. We can just use 7D instead.

Limitation

Months and Years automatically align to calendar months and years. So if you use 2020-01-15 as a start date with an update interval of 1M , the period you get back is 2020-01-01 to '2020-01-31 23:59:59.99999. It's the same for years 2019-06-23 plus 1Y could be a couple of different things because of the leap year.

This is reasonable behaviour - it is clear what a month from 2020-01-01 is. If you add a month onto 2020-01-15 then the number of days you should add is ambiguous. Same with years - because of leap years. The restriction for users would be - if you want to use monthly or annual steps then the start date needs to be the start of a calendar month and years, not some weird mid period start.

That seems like a fringe use case anyway. I don't doubt we could work around it, but we'd have to come up with arbitrary rules for resolving those ambiguities. We might circle back to this later, but something like this fixes this issue for now.

So... thoughts?

[arne-exe]

These things are new to me but it's interesting to see the thought process behind the decisions for implementing time in the model.
If it adds anything to the discussion, I read through the above and was able to understand and follow the reasoning.

Months and Years automatically align to calendar months and years. So if you use 2020-01-15 as a start date with an update interval of 1M , the period you get back is 2020-01-01 to '2020-01-31 23:59:59.99999. It's the same for years 2019-06-23 plus 1Y could be a couple of different things because of the leap year.

This is reasonable behaviour - it is clear what a month from 2020-01-01 is. If you add a month onto 2020-01-15 then the number of days you should add is ambiguous. Same with years - because of leap years. The restriction for users would be - if you want to use monthly or annual steps then the start date needs to be the start of a calendar month and years, not some weird mid period start.

Agree this sounds reasonable.

[vgro]

@vgro and I had a chat recently about time axis dates and inputs and outputs that raised the issue of calendar aware dates in the VE see here. We've been ducking this forever but I think we're at a stage where it is finally blocking things we want to do.

Please could you have a read and leave any thoughts. I don't think this would be a difficult change - and I think it probably is on the critical path - but I want input before we turn it into an actual issue for action.

I agree that it is where we want to go in the long-term, and I also agree that it is not mission critical for the minimum viable product; getting all the features in place is I believe more urgent and totally doable with the current time index setup.

At the moment, we're using pint to configure the run length and update interval.

* That's fine for giving us a number of intervals and the length of each interval, ...

* but it's **useless** for tying the simulation clearly to a calendar-aware set of dates. The `pint` units for years and months are "average" duration, so are floats:

In [3]: pint.Quantity(1, "month").to("days")
Out[3]: <Quantity(30.4375, 'day')>
In [4]: pint.Quantity(1, "year").to("days")
Out[4]: <Quantity(365.25, 'day')>

At the moment, we're simply ducking this and updates iterate over the expected number of intervals just by index [0, 1, ..., N].

That's a problem (with reasons this might not be urgent) because:

1. It's a horrible, clearly low quality state of affairs that we should be embarrassed by.
   _That @davidorme is embarrassed by._

Given that everything runs with monthly means at the moment, I don't think this is a reason to be embarrassed ;-)

2. We can't define the time coordinates of the time update axis properly so can't validate the time axis across inputs.
   _Well, OK, but we've survived so far without validation - as long as the number of intervals along the time axis in inputs is fine it works. Models know how long an interval is, and that is the main thing._

I think we could validate the correct number of time steps just like the correct number of grid cells without having exact dates. Users could theoretically make up data, for example to do paleo experiments or just hypothetical climate scenarios -> then they have to make up dates.

3. We have  already had to implement hacks to work around this issue: https://github.com/ImperialCollegeLondon/virtual_ecosystem/blob/0f4f9ed39896693f966b14fdbc56d53d0378d592/virtual_ecosystem/models/hydrology/hydrology_model.py#L443-L452
   
   _Sure, but those hacks work! Does a recurrent hack warning in every log matter?_

This hydrology section is to make it flexible for running the full VE at different timescales which again is not something that we will do in the near future. I would rather focus on improving the rainfall generator than reading in daily rainfall data at this stage.

4. It makes it impossible to link inputs or outputs at different temporal frequencies. If we have a central update frequency of - say - monthly  values but a model wants to use daily inputs, then we _have_ to have calendar-aware update periods to be able to tie update index [0] to '2000-01' in the monthly inputs and '2000-01-01/2000-01-31' in the daily inputs. I think the only way we could duck this would be to insist on 360 day years (12 x 30) - some climate models do this - but it's going to be a hard blocker for most users.

I think we could get around this by introducing a time axis which could contain the daily values, or were you thinking of having lots of nans in the other variables?

   We already have an current issue about outputting data that needs date aware coordinates: [Write out daily hydrology values for calibration validation #1176](https://github.com/ImperialCollegeLondon/virtual_ecosystem/issues/1176)
   _This is hard to duck - it depends on whether daily inputs are on the path to the short-term minimum viable product. The hydrology model at least wants daily inputs and deciduous plants are going to need to know what the date is._

Again, the hydrology daily input is not mission critical and for the validation the daily output would be a nice to have but could be worked around. Not that I don't want the time awareness, but this comes with all sorts of other problems...

I also think it would only make sense to have the daily inputs if the plants run at the same daily/sub-daily timescales - as long as this connection is not made there is no need to be overly detailed with the input, it might even give a false sense of certainty.

What do we care about

There are many Python packages providing "calendar aware" dates.

* Classic `datetime` like packages: `datetime` itself and improvements (Arrow, Pendulum, whenever)

* Array/time series like datetimes (`numpy` with `np.datetime64` and then `xarray` and `pandas` which wrap more or less wrap `numpy`)

These differ in how fancy they get about things like leap seconds, daylight savings, timezones etc. I think we largely don't care about any of those: VE will run in small, single locations on at least yearly scales where we don't want to track saving times. So, something like pandas provides functionality that provides calendar aware arrays of date ranges without us having to reinvent any wheels.

Possible solution

This pandas code allows us to generate a sequence of calendar aware periods, from reasonably user friendly inputs. Honestly, it would be easier to configure end_time rather than run_length, but from where we are this works:

In [2]: import pandas as pd
   ...: 
   ...: start_date = "2000-02-01"
   ...: update_interval = "2M"
   ...: run_length = "2Y"
   ...: 
   ...: # Calculate the requested end time from the run length
   ...: simulation_span = pd.period_range(pd.Timestamp(start_date), periods=1, freq=run_length)
   ...: span_start_time = simulation_span[0].start_time
   ...: span_end_time = simulation_span[0].end_time
   ...: 
   ...: # Calculate the intervals
   ...: simulation_intervals = pd.period_range(
   ...:     start=span_start_time, end=span_end_time, freq=update_interval
   ...: )
   ...: 
   ...: # Get the interval durations
   ...: interval_starts = simulation_intervals.start_time.to_numpy().astype("datetime64[s]")
   ...: interval_ends = simulation_intervals.end_time.to_numpy().astype("datetime64[s]")
   ...: interval_durations = (interval_ends - interval_starts) + 1
   ...: 
In [3]: interval_durations
Out[3]: 
array([5184000, 5270400, 5270400, 5356800, 5270400, 5270400, 5097600,
       5270400, 5270400, 5356800, 5270400, 5270400],
      dtype='timedelta64[s]')
In [4]: simulation_intervals[0]
Out[4]: Period('2000-01', '2M')
In [5]: simulation_intervals[0].start_time
Out[5]: Timestamp('2000-01-01 00:00:00')
In [6]: simulation_intervals[0].end_time
Out[6]: Timestamp('2000-02-29 23:59:59.999999999')

Explanations

* This is using `pandas.Period` not `pandas.Timestamp`. That's partly because the data _do represent  periods_ not point values, but also because the timestamp implementation is edging back towards the fancier end of datetime handling with all sorts of business days, timezone etc. malarky.

* The `pandas.Period` implementation is a lot simpler - each item is explicitly a set number of *_days, weeks, months or years_ (I think it would handle faster frequencies fine, but we aren't going there anytime soon). The possible letter codes we'd want to use are therefore `nD`, `nM` and `nY` for days, months and years (see https://pandas.pydata.org/pandas-docs/stable/user_guide/timeseries.html#period-aliases for details) and those seem pretty user friendly.

* We might have wanted weeks but `pandas` anchors weeks to a calendar week with a particular start day (e.g Monday) rather than just 7 days. We can just use `7D` instead.

I have no preference what to use, I am not so familiar with either at the moment. I think we might want to have a think about the general time steps we allow in the model. Are we going to make this a monthly model, to allow for the animals to make sense, and have some components always on a shorter timescale? Or are we giving the user the option to mix and match timescales based on what they want to focus on?

Limitation

Months and Years automatically align to calendar months and years. So if you use 2020-01-15 as a start date with an update interval of 1M , the period you get back is 2020-01-01 to '2020-01-31 23:59:59.99999. It's the same for years 2019-06-23 plus 1Y could be a couple of different things because of the leap year.

This is reasonable behaviour - it is clear what a month from 2020-01-01 is. If you add a month onto 2020-01-15 then the number of days you should add is ambiguous. Same with years - because of leap years. The restriction for users would be - if you want to use monthly or annual steps then the start date needs to be the start of a calendar month and years, not some weird mid period start.

That seems like a fringe use case anyway. I don't doubt we could work around it, but we'd have to come up with arbitrary rules for resolving those ambiguities. We might circle back to this later, but something like this fixes this issue for now.

So... thoughts?

In summary,

• I agree that this is where we want to go, maybe not a priority for the MVP given that some key functionalities are still missing.
• If we chose to implement time awareness, there will be quite a lot of structural changes required to the abiotic_simple,abiotic, and hydrology model. This is ok, but again, takes time away from other MVP jobs.
• I suggest that we first have a chat about the general VE and individual model time steps flexibility - do we still even want to allow for that? It is sort of built in but makes everything a bit clunky and probably never used; might be better to invest energy in the proper integration.

[TaranRallings]

Following @vgro , I agree this is something we want but also that it is not necessary until we have a feature complete MVP.

The use case for the animal model is that we have some processes that function on seasonality. We may, in the future, get mechanistic seasonality cues but for now having a calendar date to work off would simplify things significantly. It is also the case that a human-facing disturbance system would require calendar aware time for setting things like harvest and planting schedules.

The pros and cons of months and years look acceptable to me.

[davidorme]

Trying to summarise a more wide-ranging IRL conversation about timing. @vgro, @sallymatson and @jacobcook1995 feel free to edit this message!

We started from @vgro's question:

I suggest that we first have a chat about the general VE and individual model time steps flexibility - do we still even want to allow for that? It is sort of built in but makes everything a bit clunky and probably never used; might be better to invest energy in the proper integration.

The consensus was that we do still want to allow for general flexibility, but with the main motivation being that we want to be able to assess the calibration of the overall model with different intervals. We're not really so interested in the currently formal coded limits of the individual models, but how does the model work with steps in the ~14-60 day range. The coded model-specific update interval ranges are still useful information, but we probably want users to focus on a rough month long time step, whilst allowing them to explore how changing that slightly affects predictions.

So - given that we want to be able to vary the update interval - what is the way forward.

• The proposal above works - and comments so far are generally positive - but having calendar-aware intervals leads to (mildly) variable length intervals. That leads to scaling issues for predictions - should values be converted to rates per day rather than totals per interval to iron out these differences (worst case is a ~10% change in total value from Jan (31 days) to Feb (28 days) monthly values with identical rates). So it's desirable, easy to implement the indexing and passing dates in, but making the models work with variable length months is a much deeper problem.

• A simpler approach is to simply enforce equal intervals as a fixed number of days using a 360 day year. We're basically ditching dates here completely and use sampling to draw input data at particular frequencies. As long as we have an update interval that is a divisor of 360, we'd still be tied to fixed points in the annual cycle, so can play with update frequencies a bit.

  Here, we're basically using the existing indexing system for time but just ditching any pretence of using pint to indicate runtimes: you configured 120 time steps of 15 days. It's up to users to sample their inputs such that the input data covers 5 annual cycles (120 * 15 / 360 = 5.0).

[jacobcook1995]

Think this is a good summary.

My inclination would be to go with the second option as it is simpler (we really don't need to introduce more things to potentially trip ourselves up at this point), unless someone has a immediate use case that requires calendar awareness

[davidorme]

I think we'll need to be very clear about how data prep should work - there are lots of potential issues in resampling data that is on a calendar into this framework. Or at least lots of ways it could get over-complicated.

To be clear, I don't like the second option. But.. I think @jacobcook1995's point is telling: unless we are sure the outputs are calendar-aware and won't vary purely with individual update interval lengths then it's the easiest way for now. Honestly, I think we'd be in a better place if we'd fixed this core issue earlier and written the models to use calendar-calibrated outputs. I think we should still do that at some point because I think our target users are typically going to be trying to make it work with calendar data and are going to find that assumption unpleasant and off-putting.

We also have @TaranRallings point:

It is also the case that a human-facing disturbance system would require calendar aware time for setting things like harvest and planting schedules.

Well, again, it's much more intuitive. But we can simply say that planting occurs starting on the 5th index and then every 12 indices etc. I don't like that either but reviewing the outputs to check they are stable with respect to interval time is not a thing we want to do now.

[dalonsoa]

I was going to comment on #1176, but it seems that this discussion is at the root of the problem, so will do so here.

Correct me if I'm wrong, but my impression is that a calendar is a somewhat arbitrary human concept. Neither animals, plants or the weather cares much a about calendars, just about repetitive cycles, and these are split into equal length intervals. So, for the calculations to be "correct", I think they should use a time schema that ensure these cycles. Models will be simpler this way and any option that goes in this direction will be preferable.

Now, humans do use calendars - some do, at least :P - and therefore input data is likely to be based on calendar-aware information, and to be useful, output data should also be provided in a calendar aware format, most of the time.

From my point of view, what is needed here is an adaptor in the input and output layers that make this transformation between calendar-aware data and equal interval data, and vice-versa, in such a way that works for users, making their life a bit easier when compiling the inputs and analysing the outputs, but works well for the models.

This might be going in a different direction that what you had in mind during the discussion, but I feel this is what is lacking at the moment.

[davidorme]

Essentially, yes, the calendar is a construct but they provide the coordinates attached to most time series data. Not all - some climate models use fixed length 360 day years - and if we end up doing paleo-ecology then I don't think those inputs are likely to be generated on a calendar! As @TaranRallings says, we also have to support use-cases that have human-based timings they want to use in the model (and those are a pretty core part of the mission).

I think we could create an adaptor - and there are already things that might help,such as xarray.Dataset.convert_calendar.

But - I think I'd argue that writing that adaptor is a thing that further bakes in an existing structure that I don't think we want to keep. I think ultimately it makes the model better if it is e.g. outputting daily rates rather than per-update totals. Those are ecologically much easier to interpret and can be compared much more easily across models with different frequencies. That does require a root and branch review of the models though, so falls off our current critical path.

[dalonsoa]

Fine by me :)

[vgro]

Would there be an option to have the starting month? I realize now that I will need to base some calculations on the month and latitude for day length. I can make a 12h day assumption for the tropics as a placeholder but eventually need to know at least if the simulation start in January or July...

[davidorme]

We already have that - we have a start date. That's the last point we go anywhere near a calendar 😄

[davidorme]

There's some day length stuff and solar elevation in pyrealm if that is useful.

[davidorme]

@jacobcook1995 @TaranRallings @vgro @sallymatson

I was just thinking about @jacobcook1995's point about output variation due to interval length. Was out running in the rain along a canal and thinking about precipitation and suddenly wondered: don't we already have a problem here because monthly input datasets using current data (like CRU etc - not climate forecasts or paleo models) are going to have variations in totals because of monthly variation. And if the data preparation team are summing values from daily inputs to monthly totals (without having to do mad calendar corrections -which was @dalonsoa's point about an interface) then we have the same problem.

How many of our input variables are "totals over the update interval"? I'm not saying that we suddenly switch to implementing this, but I wonder if our outputs are actually any less calendar biased now than they would be if we went calendar aware.

[jacobcook1995]

The external data that litter and soil take in are all estimates of initial state rather than changes over time. All the variables that relate to changes over time (i.e. addition of carbon from plants to soil) are calculated by other models. The internals of the models use rate parameters estimated from data, but these would always be expressed relative to time units of fixed length (i.e. seconds or days).

I guess this is just because soil and litter don't use monthly input datasets though, do many of the other models use monthly input datasets?

[davidorme]

Plants just uses the average microclimate and the average solar flux in µmol m-2 s-1 and scales productivity up to the update interval length. We have rates that are applied at the update interval timescale but they are configured per day and rescaled. We might need to make sure they are scaled per interval, rather than once at the start.

[TaranRallings]

Animal input data is just functional group distributions so we largely dodge this issue. Even the rates are per-day such that handling variable length intervals shouldn't be a big problem on my end.

[davidorme]

Linking in an issue from the data science team that illustrates why what we currently have is a problem:

https://github.com/ImperialCollegeLondon/ve_data_science/issues/161

[nickwctan]

Closely following on this discussion.

This maybe a bit off topic though. I would like to add:

As a user, I found it a bit counter intuitive that in the exported output file, the "time" column (obtained from timestamp) is in POSIX date format when it is a 30 days fixed interval. Personally, I think having "Day since 0" or both "Day since 0" and "POSIX" makes more sense.

If I were to use VE as a researcher, I will be looking out for daily rate outputs to design my studies or interpret my results. Was there a reason why we do not use day as timestamp interval for input and output data?

[davidorme]

Days since start is something we could do, but it does make it more awkward to tie the model to annual cycles. I think that most use cases initially are going to be taking recent data which is typically arranged on calendar dates. So its trying to find the best way to handle time.