如何在Xarray数据集中获取像素首次满足阈值的时间并解决类型冲突？

我有一个记录火山灰累积量的Xarray数据集，结构如下：

xarray.DataArray 'mass' (time: 86, y: 115, x: 132)>
array([[[0.        , 0.        , 0.        , ..., 0.        ,
         0.        , 0.        ],
        [0.        , 0.        , 0.        , ..., 0.        ,
         0.        , 0.        ],
        [0.        , 0.        , 0.        , ..., 0.        ,
         0.        , 0.        ],
        ...,
        [7.86871581, 8.03121152, 8.1963989 , ..., 8.52833785,
         8.34455418, 8.16493449],
        [7.71566969, 7.87296236, 8.03278228, ..., 8.34824074,
         8.1713868 , 7.99838556],
        [7.56522864, 7.71749057, 7.87212488, ..., 8.17195143,
         8.00172886, 7.83507221]]])
Coordinates:
  * time     (time) datetime64[ns] 2021-09-19 2021-09-20 ... 2021-12-13
    band     int64 1
  * y        (y) int64 3161827 3161927 3162027 ... 3173027 3173127 3173227
  * x        (x) int64 213152 213252 213352 213452 ... 226052 226152 226252

这个变量代表火山喷发期间地面火山灰的累积量，数值从0开始单调递增。我需要为每个像素计算首次达到累积量≥1的时间，最终返回一个维度为(time: 1, y: 115, x: 132)的DataArray，每个像素对应对应的日期。

我知道可以沿time维度循环实现，但想找更优雅高效的方法。尝试用da.where和da.idxmin实现时：

da.where(da>=1).idxmin(dim='time')

触发了数据类型冲突错误：

TypeError: The DTypes <class 'numpy.dtype[int64]'> and <class 'numpy.dtype[datetime64]'> do not have a common DType. For example they cannot be stored in a single array unless the dtype is `object`.

目前临时解决方法是把time坐标转换成距起始时间的天数（int类型），但希望有更优雅的方案。

错误重现步骤 #

使用本地文件da.nc，执行以下代码：

import xarray as xr

da = xr.open_dataset('~/Desktop/tmp/da.nc')
da = da.rename({'__xarray_dataarray_variable__':'thickness'})

da['thickness'].where(da['thickness']>=1).idxmin(dim='time').plot()

相关预处理代码 #

# 计算累积和
da = da.cumsum(dim='time')

# 设置日期
date = ['2021-09-19','2021-10-11','2021-11-26','2021-12-13']
date = pd.DataFrame(date,columns=['time'])
date['time'] = pd.to_datetime(date['time'])

# 插值到日分辨率
dateD = date.set_index('time').resample('d').sum().reset_index()
da = da.interp(time=dateD['time'])

# 转换为数据集并移除band维度
da = da.drop_vars("band")
DA = da.to_dataset(name='thickness').squeeze()

出现这个错误的核心原因是：da.where(da>=1)会将不满足条件的值设为NaN，而idxmin默认返回的是满足条件的最小值对应的索引位置（int类型），但我们需要直接映射到datetime类型的time坐标，两者类型不匹配导致冲突。

下面是两种优雅的解决思路：

方法1：通过索引映射获取日期 #

先找到每个像素首次满足条件的时间索引，再用索引从time坐标中提取对应日期：

import xarray as xr
import pandas as pd

# 生成布尔掩码：累积量≥1的位置标记为True
mask = da['thickness'] >= 1

# 沿time维度找到第一个True的位置（argmax返回第一个最大值的索引，True=1、False=0，刚好对应首次达标时间）
first_idx = mask.argmax(dim='time')

# 用索引提取对应的time值，得到每个像素的首次达标日期
first_date = da['time'].isel(time=first_idx)

# 按需求调整维度为(time:1, y, x)
first_date = first_date.expand_dims('time', axis=0)

这种方法利用了数据单调递增的特性，argmax能精准定位到第一个满足条件的时间点，完全避免类型转换问题。

方法2：填充无效值后再取索引 #

通过将不满足条件的位置填充为极大值，再用idxmin找到有效索引，最后映射回datetime：

import numpy as np

# 生成有效掩码
valid_mask = da['thickness'] >= 1

# 将无效位置填充为int64的最大值，确保idxmin返回第一个有效位置的索引
filled = xr.where(valid_mask, da['time'].astype(np.int64), np.iinfo(np.int64).max).idxmin(dim='time')

# 将索引转换回datetime类型
first_date = pd.to_datetime(filled)

# 重构为带坐标的DataArray，并调整维度
first_date = xr.DataArray(first_date, dims=['y','x'], coords={'y':da['y'], 'x':da['x']})
first_date = first_date.expand_dims('time', axis=0)

结果验证 #

因为数据是单调递增的，两种方法都能准确找到每个像素首次达到阈值的时间。最终得到的first_date就是维度为(time:1, y:115, x:132)的DataArray，每个像素对应首次达标日期。

内容的提问来源于stack exchange，提问作者e5k