Как получить максимальное время разрыва в наборе данных xarray

Question

У меня есть набор данных xarray с 3 измерениями (широта, долгота, время) для каждой переменной. У меня есть 720 значений для лата, 1440 для лона и 13140 для времени. Для каждого пикселя у меня есть некоторые промежутки во временном ряду, и я пытаюсь узнать, каково максимальное и среднее время промежутков. Поскольку это большой набор данных, я пытался обрабатывать его по годам.

Выходные данные для набора данных за 1 год (имя указано в коде):

#the dataset for 1 year:
<xarray.Dataset>
Dimensions:         (lat: 720, lon: 1440, time: 365)
Coordinates:
* lat             (lat) float32 89.875 89.625 89.375 89.125 88.875 88.625    ...
* lon             (lon) float32 -179.875 -179.625 -179.375 -179.125 ...
* time            (time) datetime64[ns] 1981-04-06 1981-01-18 1981-09-29 ...
Data variables:
t0              (time, lat, lon) datetime64[ns] dask.array<shape=(365, 720, 1440), chunksize=(1, 720, 1440)>
sm              (time, lat, lon) float32 dask.array<shape=(365, 720, 1440), chunksize=(1, 720, 1440)> 

Я пробовал этот код с циклом для каждого года:

# create dataset of nan to then fill it with the values
var=np.zeros((36,720,1440))*np.NaN
lat = combined.lat.values
lon = combined.lon.values
time_na = time # each year
diff_day = xr.Dataset(
    data_vars={'max':    (('time','lat', 'lon'), var),'mean':    (('time','lat', 'lon'), var)},
    coords={'time': time_na, 'lat': lat, 'lon':lon})

for t,name in tqdm(enumerate(filenames)): #loop for each year
  filename_year = glob(name+'/*.nc') # read all the files for the year
  combined = xr.open_mfdataset(filename_year,concat_dim='time',autoclose =True, decode_times=True)
  combined = combined.sortby(combined['time'],ascending=True) # otherwise the time is not montonic

  # calculation pixel by pixel
  for i in range(len(combined.lat)):
    for j in range(len(combined.lon)):
        if len(combined.time.values[np.isfinite(combined.sm.values[:,i,j])])>1 : # avoid cases where it's a list of nan 
            # the idea is to make the diff of time between finite (not finite values correspond to the gap) values.
            diff_day['max'].loc[t,i,j] = np.diff(combined.time.values[np.isfinite(combined.sm.values[:,i,j])]).astype('timedelta64[D]').max()/ np.timedelta64(1, 'D')
            diff_day['mean'].loc[t,i,j] = np.diff(combined.time.values[np.isfinite(combined.sm.values[:,i,j])]).astype('timedelta64[D]').mean()/ np.timedelta64(1, 'D')

Этот код работает, но время процесса слишком велико. Мне интересно узнать, есть ли более простой способ сделать это. Спасибо

Answer 1

Если вы хотите получить среднее число значений NaN, вам поможет что-то такое простое, как da.isnull().mean(dim='time'). Но получить среднюю и максимальную длину непрерывных блоков NaNs - более сложный алгоритмический вопрос, чем простой процедурный вопрос xarray.

Я уверен, что есть много способов сделать это, но тот, который я придумал, был таким:

Сначала создайте массив такой же формы, что и ваши данные, который просто увеличивается по временному измерению:

In [10]: arange = xr.ones_like(da) * np.arange(len(da.time))

В данных игрушки, которые я сделал для этого, временные ряды каждой ячейки выглядят так:

In [11]: arange.isel(lat=0, lon=0).plot()

Затем создайте аналогичный массив, но с периодами, которые содержат константу для каждого блока NaN:

In [12]: cumulative_nans = (arange.where(da.notnull()).ffill(dim='time')).fillna(0)

В каждой ячейке этого массива есть ступени лестницы для каждого блока NaN:

In [13]: cumulative_nans.isel(lat=0, lon=0).plot()

Теперь вы можете вычесть эти два, чтобы получить массив, в котором значение в каждой ячейке является счетчиком с совокупным числом NaN в этом блоке:

In [14]: time_series_of_cumulative_nan_blocks = (arange - cumulative_nans)

В каждой ячейке:

In [15]: time_series_of_cumulative_nan_blocks.isel(lat=0, lon=0).plot()

Вы можете легко рассчитать максимальное значение для этого:

In [16]: max_nan_duration = time_series_of_cumulative_nan_blocks.max(dim='time')

Среднее значение сложнее. Мы можем использовать переход от одного временного шага к следующему, чтобы отфильтровать данные, чтобы включить только точки, где отклоняется следующая ячейка, например, где у нас есть пик:

In [17]: nan_block_length_peaks_only = (
    time_series_of_cumulative_nan_blocks
    .where(
        time_series_of_cumulative_nan_blocks
        .diff(dim='time', label='lower')
        < 0)

В каждой ячейке третье число ограничено набором точек:

In [18]: nan_block_length_peaks_only.isel(lat=0, lon=0).plot(marker='.')

Это значение может быть усреднено, чтобы найти среднюю продолжительность:

In [19]: mean_nan_duration = nan_block_length_peaks_only.mean(dim='time')

Это должно иметь значительные преимущества в производительности по сравнению с циклическим просмотром каждой ячейки и периода времени в наборе данных, главным образом потому, что оно основано на скомпилированных векторизованных функциях xarray, а не на циклах python. Это вычисление также может быть выполнено для всего набора данных с использованием dask , что может привести к дальнейшему увеличению, в зависимости от ваших настроек.