create time bounds after time averaging

[siligam]

Missing time bounds (time_bnds) variable after data is resampled to new frequency (i.e., time averaging step).

After introducing this feature, the netcdf file should roughly look like this

netcdf file_with_time_bounds {
dimensions:
	time = UNLIMITED ; // (2 currently)
	lat = 3 ;
	lon = 4 ;
	bnds = 2 ;   # <-- this entry
variables:
	double time(time) ;
		time:bounds = "time_bnds" ;      # <-- this entry
		time:units = "days since 2018" ;
		time:calendar = "360_day" ;
		time:axis = "T" ;
		time:long_name = "time" ;
		time:standard_name = "time" ;
	double time_bnds(time, bnds) ;      # <-- this entry

[pgierz]

There is already a step for "variable level attributes", where the top part could go, but I would rather group the related code together with actually inserting the time_bnds variable and create a separate step for adding time bounds specifically, or more generally, adding bounds (also lon/lat/other).

What isn't immediately intuitive is which DataArray in the Dataset needs to be modified for extra variable-level attributes. You could have:

def my_bounds_step(data: DataArray, rule: Rule):
    # [FIXME] How to know that we are doing time specifically?
    # Maybe:
    # 1. Read from the rule's drv/table info?
    # 2. Read from the data's attrs? 
    # 3. Read from the data's coords?
    data.attrs['bounds'] = 'time_bnds'
    return data

However, if you have a Dataset it isn't obvious:

def my_bounds_step(data: Dataset, rule: Rule):
    da = data.variables.get(...)  # [FIXME] Which variable?
    da.attrs['bounds'] = 'time_bnds' # [FIXME] Same questions as above
    ds = da.to_dataset()  # [FIXME] Is this needed?
    return ds

And eventually you will need to have a Dataset structure, since you are not only passing around the main variable, but also now the bound variables.