filter on the time dimension with a large dataset

[miniufo]

I am interested in the lntime accessor that could filter along the time dimension. So a long-standing problem (at least for me) is how to do filtering along time dimension in a large dataset while the time dimension is chunked?

For example, I have the daily AVISO sea-surface height dataset over several years, which is chunked along time (2D lat-lon data in a single file per day):

<xarray.Dataset> Size: 364GB
Dimensions:         (time: 1096, lat: 1440, nv: 2, lon: 2880)
Coordinates:
  * time            (time) datetime64[ns] 9kB 2015-01-01 ... 2017-12-31
  * lat             (lat) float32 6kB -89.94 -89.81 -89.69 ... 89.69 89.81 89.94
  * lon             (lon) float32 12kB -179.9 -179.8 -179.7 ... 179.8 179.9
  * nv              (nv) int32 8B 0 1
Data variables: (12/14)
    crs             (time) int32 4kB -2147483647 -2147483647 ... -2147483647
    lat_bnds        (time, lat, nv) float32 13MB dask.array<chunksize=(1, 1440, 2), meta=np.ndarray>
    lon_bnds        (time, lon, nv) float32 25MB dask.array<chunksize=(1, 2880, 2), meta=np.ndarray>
    sla             (time, lat, lon) float64 36GB dask.array<chunksize=(1, 1440, 2880), meta=np.ndarray>
    err_sla         (time, lat, lon) float64 36GB dask.array<chunksize=(1, 1440, 2880), meta=np.ndarray>
    ugosa           (time, lat, lon) float64 36GB dask.array<chunksize=(1, 1440, 2880), meta=np.ndarray>
    ...              ...
    err_vgosa       (time, lat, lon) float64 36GB dask.array<chunksize=(1, 1440, 2880), meta=np.ndarray>
    adt             (time, lat, lon) float64 36GB dask.array<chunksize=(1, 1440, 2880), meta=np.ndarray>
    ugos            (time, lat, lon) float64 36GB dask.array<chunksize=(1, 1440, 2880), meta=np.ndarray>
    vgos            (time, lat, lon) float64 36GB dask.array<chunksize=(1, 1440, 2880), meta=np.ndarray>
    flag_ice        (time, lat, lon) float64 36GB dask.array<chunksize=(1, 1440, 2880), meta=np.ndarray>
    tpa_correction  (time) float64 9kB dask.array<chunksize=(1,), meta=np.ndarray>
Attributes: (12/42)
    Conventions:                     CF-1.6
    Metadata_Conventions:            Unidata Dataset Discovery v1.0
    cdm_data_type:                   Grid
    comment:                         Sea Surface Height measured by Altimetry...
    contact:                         servicedesk.cmems@mercator-ocean.eu
    creator_email:                   servicedesk.cmems@mercator-ocean.eu
    ...                              ...
    geospatial_vertical_units:       m
    time_coverage_duration:          P1D
    time_coverage_resolution:        P1D
    time_coverage_end:               2015-01-01T12:00:00Z
    time_coverage_start:             2014-12-31T12:00:00Z
    platform:                        Cryosat-2, OSTM/Jason-2, Haiyang-2A, Altika

then how to do filtering along time (for example to extract the signal of a certain bandwidth) on a machine that cannot load the whole dataset into memory?

Just want to know if lenapy has some better way to do this.

[hulecom]

I would say from personal experience that the best way to time filter such datasets is, if possible, to rechunk along lon and lat dimensions. The chunks along time dimension might caused some issues with the filtering operators.

I did some tests with lntime.filter and some large synthetic data. It seems to works fine even if its quite long on my side (but chunk optimization depending on the computing environment might improve that).

import lenapy
import numpy as np
import pandas as pd
import dask.array as da
import xarray as xr
import matplotlib.pyplot as plt

# coordinates
time = pd.date_range("2015-01-01", periods=1096, freq="D")  # ~3 years
lat = np.linspace(-90, 89.99, 144000)

# Build a synthetic (time, lat) signal with a seasonal cycle + latitude structure + noise
t = np.arange(time.size)
seasonal = np.sin(2*np.pi*t/365.0)[:, None] * np.cos(np.deg2rad(lat))[None, :]
seasonal = da.from_array(seasonal, chunks=(-1, 1000))

noise = da.random.normal(0, 0.1, size=(time.size, lat.size), chunks=(-1, 1000))
data = seasonal + noise

# assemble DataArray (time, lat)
da_tlat = xr.DataArray(
    data,
    coords={"time": time, "lat": lat},
    dims=("time", "lat"),
    name="synthetic_signal",
).chunk({"time": -1, "lat": 1000})


lat0 = 45.0

# plot before filter
plt.figure()
da_tlat.sel(lat=lat0, method="nearest").plot()


da_tlat_filt = da_tlat.lntime.filter(filter_name='lanczos',cutoff=36,order=2)

#plot time series after filter
da_tlat_filt.sel(lat=lat0, method="nearest").plot()

plt.title(f"Synthetic signal at lat=45° (pre and post-filter)")
plt.xlabel("Time")
plt.show()

As you can see with issue #8, we are also currently working on the refactor and time improvement of the climatology estimation for time-series.

There are also some optimization works that can be done on the .lntime.filter part with better arrangement of the chunks during the operation I think, it is an area of improvement for the future.

[miniufo]

Great you are working on improving this. The multidimensional dataset is usually stored in (t, y, x)-sequence (x-dim is continuous on the disk), which makes it more difficult to process along time dimension even if you re-chunked as above. So for an efficient processing along time (like filtering), one would expect the dataset is stored in (y, x, t)-sequence. Based on my experience, such a conversion of the data file, from (t, y, x) to (y, x, t) is also time-consuming and inefficient.

I would close this now, as this is not part of the review. But it would be great if you have better solution.