Support for monthly/yearly data

[TLMueller]

Hi everyone,

I'm sorry for this rather general question, but I could not find an answer to it in the manual.

Is it currently possible to use pastas for monthly or yearly data? Because theoretically the intervalls are not equidistant and in the docs it is pointed out several times, that the stress needs to be equidistant.
I already tried it with monthly testdata, however I ran into the following problem:

ml1 = ps.Model(head, name="model1")

# Define precipitation stress model
stress_prec = ps.StressModel(stress=prec,rfunc=ps.Gamma(),up=True,name='prec')
# Define evapotranspiration stress model
stress_evap = ps.StressModel(stress=evap,rfunc=ps.Gamma(),up=False,name='evap')

ml1.add_stressmodel(stress_prec)
ml1.add_stressmodel(stress_evap)
ml1.solve(noise=False)

When I try to create a model the solve() step throws out an error, telling me that the time series cannot be extended into the past. I guess this is because for some reason the solve method automatically tries to create a warmup period. Setting warmup=None didn't help. If I set warmup=0, then I get an error telling me, that the simulation contains NaN values. However I made sure, that there are no NaNs in my timeseries, so this is a bit confusing to me.

I tried to use the predefined settings for the stresses, like settings='prec' and settings='evap' and then it worked. But I would like to avoid using settings for now to get a better understanding of how everything works and obviously there is something happening in the background, waht I don't quite understand.

Some extra info:
The get_stressmodel_settings function showed, that pastas seems to think, that my frequency is daily (freq='D'), depite it being monthly.
ps.validate_stress(prec) throws no errors, so according to the validate function, my timeseries is equidistant, which in theory it sould not.

If anyone can help, that would be greatly appreciated.

[mbakker7]

Just to make sure that I understand: You have monthly stresses and want to simulate with monthly time steps. That means you have head measurements only at the end of each month. Is that indeed what you want to do? Or do you want to enter monthly stresses (say monthly average values) and then simulate heads on a daily basis?

[TLMueller]

My head measurements are monthly. The real world measurements are at varying days within the month, but for modeling purposes I put them at the end of the month. If I have more than one measurement per month, I averegage them and put the value at the end of the month.

The precipitation data is also monthly (monthly sum) [edited].

The evaporation data is daily. However, to my understanding the stresses need to have the same resolution. So I aggregated my evaporation to monthly sums.

So yes: I want to simulate monthly time steps based on monthly stresses.

Thanks already for your quick reply :)

Edit: Sorry, the prec data is also monthly.

[dbrakenhoff]

Good questions @TLMueller. Is it an option for you to resample your data to a regular (near-monthly) timestep? E.g. 30 days? That way you meet pastas equidistant stresses requirement, but the result also fits the resolution of your data. The function ps.ts.timestep_weighted_resample() might be helpful to resample your monthly series to 30D series. That method is described in this notebook: https://pastas.readthedocs.io/en/latest/examples/prepare_timeseries.html

Another option is maybe to keep the daily calculation time step with daily precipitation and evaporation stresses. And resample any stresses you have on a monthly basis to daily values, for example by dividing the monthly stress over the number of days in that month. That way you won't have to resample your head time series at all, allowing you to calibrate the model on all your observations. It could be worth trying, but this might not work very well if your data is really too sparse.

If I set warmup=0, then I get an error telling me, that the simulation contains NaN values. However I made sure, that there are no NaNs in my timeseries, so this is a bit confusing to me.

This is something we have to look into, as I think we should provide a more clear error message for this case. A quick guess: it seems as if there might be a period in which there are head observations but no precipitation data? In solve() pastas is extending the precipitation time series (but not necessarily to fill the warmup period, because warmup=0), but to match the tmin/tmax? This introduces a NaN that cannot be filled because no filling logic is provided.

[TLMueller]

Thanks again for the super quick reply!

So if I understand you correct, monthly (or yearly) stress time series is indeed not supported at the time and instead it is advised to instead use a 30 (or 365) day intervall, correct?

I will try to resample my data to a 30 day basis as you advised.

A quick guess: it seems as if there might be a period in which there are head observations but no precipitation data? In solve() pastas is extending the precipitation time series (but not necessarily to fill the warmup period, because warmup=0), but to match the tmin/tmax? This introduces a NaN that cannot be filled because no filling logic is provided.

During preprocessing I adjusted all time series to have the same length, so this should not be a problem.

In another trial, I used the dataset Dr. Collenteur provided in the "Workshop on modeling hydraulic head time series with Pastas" Notebook 1 in Aachen at the Fachsektion Hydrogeologie e.V. a few weeks ago. There the the precipitation time series started 1980 and the head started 1985. I aggregated the timeseries to monthly values and the same error occured. The solve() Function used a tmin=1990, so there was no need to extend the time series.

In both cases the error message references a date, which coincides with the starting date minus the default warmup perios (3650 days).

[dbrakenhoff]

During preprocessing I adjusted all time series to have the same length, so this should not be a problem.

In both cases the error message references a date, which coincides with the starting date minus the default warmup perios (3650 days).

Ok, intriguing. Are you ok with sharing your data and code so we can take a look? You can also email it to me if you'd rather not post it here.

[TLMueller]

Just sent you a Mail. Any further discussion can be continued here of course, so that others can benefit from it as well :)

[dbrakenhoff]

I got your email and had a look at your code. The reason you're seeing an error is that the internal computation time step is still daily ("D"). This is the default setting in pastas.

When running solve, you see this warning, which indicates that pastas is trying to up-sample the stresses but it could not, because the settings were not defined (=None):

To summarize, this is what is happening:

1. Your model is initialized with a default time step of 1 day.
2. Your time series are all already nicely resampled to month-end values. When you solve, the model looks at the first and last date for which all stresses are available to compute tmin/tmax.
3. It then creates a simulation index based on tmin, tmax and the computation frequency which is 1 day.
4. Next it tries to simulate on a daily basis, but it cannot resample the stress to this daily time step because no settings were provided for how to do that. This generates the warning messages.
5. It then performs convolution between a monthly stress time series and a block response (which has unit days). This does provide a result, even though it is objectively a bit weird to convolve a daily response function and monthly stress.
6. When it adds the result to the simulation, NaNs remain because the simulation can only be computed on the final day of the month.
7. Finally, we trap the error to warn the user that what they're doing is not correct.

The model timestep can be set with:

ml = ps.Model(head, freq="30D")

Note if you try to put in freq="M" you will get an error message indicating that monthly freq is not a fixed-frequency.

And to answer your question from before, yes it is recommended to use a fixed frequency like "30D" or "365D" for monthly or yearly models.