xarray数据集删除NaN无效，滚动均值平滑后求梯度最小值报错

我有一个包含Range和time坐标的xarray数据集，需要为每个time时刻找到backscatter梯度最小对应的Range。执行过程中触发ValueError: All-NaN slice encountered错误，怀疑是滚动均值平滑数据导致的，尚未确认。

处理步骤：

• 使用da.differentiate(coord = 'Range')创建backscatter梯度的新数据变量；
• 通过min_height = ds['bs_grad'].argmin(dim='Range').values和min_height = ds['Range'].isel(Range=min_height)获取梯度最小值的位置，此步骤触发错误。尝试沿time和Range维度单独调用dropna()均未解决问题。

正常运行示例 #

import pandas as pd
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt

# create a sample DataArray
ds = xr.DataArray([[7, 6, 8, np.nan, 9],
                [9, 6, 2, np.nan, 4],
                [3, 4, 1, np.nan, 1]],
    dims=['Range', 'time'],
    coords={
        'Range': [10, 20, 30],
        'time': ['2022-01-01', '2022-01-02', '2022-01-03', '2022-01-04', '2022-01-05']
    }
)


plt.pcolormesh(ds['time'], ds['Range'], ds, shading = 'auto')
plt.colorbar()
plt.gcf().autofmt_xdate()
plt.title('xarray')
plt.show()

ds['bs_grad'] = ds.differentiate(coord = 'Range')

plt.pcolormesh(ds['time'], ds['Range'], ds['bs_grad'], shading = 'auto')
plt.colorbar()
plt.gcf().autofmt_xdate()
plt.title('gradient')
plt.show()


ds = ds.dropna(dim = 'time', how = 'all')  
ds = ds.dropna(dim = 'Range', how = 'all')  

plt.pcolormesh(ds['time'], ds['Range'], ds['bs_grad'], shading = 'auto')
plt.colorbar()
plt.gcf().autofmt_xdate()
plt.title('grad, nan dropped')
plt.show()

## find the Range of the minimum backscatter for each time step
min_Range = ds['bs_grad'].argmin(dim='Range').values  # get the index of the minimum backscatter
min_Range = ds['Range'].isel(Range=min_Range)  # extract the corresponding Range


plt.scatter(ds['time'], min_Range),
plt.gcf().autofmt_xdate()
plt.show()

复现问题示例 #

import pandas as pd
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt

min_date = "2023-03-01"
max_date = "2023-04-20"
dates = pd.date_range(min_date, max_date)

# Define the dimensions
Range = np.arange(0, 1000, 50)
time = dates

# Create the data array with variables for backscatter, temperature, and humidity
data = xr.DataArray(
    np.random.rand(len(Range), len(time)),  # random data for demonstration purposes
    dims=("Range", "time"),
    coords={"Range": Range, "time": time},
    attrs={"long_name": "example data array"},
)

backscatter = xr.DataArray(
    np.random.rand(len(Range), len(time)),  # random data for demonstration purposes
    dims=("Range", "time"),
    coords={"Range": Range, "time": time},
    attrs={"long_name": "backscatter", "units": "dB"},
)

temperature = xr.DataArray(
    np.random.rand(len(Range), len(time))*20,  # random data for demonstration purposes
    dims=("Range", "time"),
    coords={"Range": Range, "time": time},
    attrs={"long_name": "temperature", "units": "K"},
)

humidity = xr.DataArray(
    np.random.rand(len(Range), len(time))*100,  # random data for demonstration purposes
    dims=("Range", "time"),
    coords={"Range": Range, "time": time},
    attrs={"long_name": "humidity", "units": "%"},
)

# Combine the data arrays into a single xarray dataset
ds = xr.Dataset(
    {"data": data, "backscatter": backscatter, "temperature": temperature, "humidity": humidity}
)

plt.pcolormesh(ds['time'], ds['Range'], ds['backscatter'], shading = 'auto', vmin = 0)
plt.colorbar()
plt.gcf().autofmt_xdate()
plt.xlabel('time')
plt.ylabel('Range')
plt.show()

## Smooth the data with a rolling mean 
ds['backscatter'] = ds['backscatter'].rolling(
    Range = 5, center=True).mean().rolling(
    time = 5, center=True).mean()

plt.pcolormesh(ds['time'], ds['Range'], ds['backscatter'], shading = 'auto', vmin = 0)
plt.colorbar()
plt.gcf().autofmt_xdate()
plt.xlabel('time')
plt.ylabel('Range')
plt.show()

## Artificially add some nan columns, rows, and min gradient
## Set a clear min gradient
ds = ds.where(ds['Range'] != 800, other= -20000)


## Set values equal to nan to siumulate problem I'm running into
ds = ds.where(ds['time'] != pd.to_datetime('2023-04-01'), np.nan)
ds = ds.where(ds['time'] != pd.to_datetime('2023-04-03'), np.nan)
ds = ds.where(ds['time'] != pd.to_datetime('2023-03-12'), np.nan)

ds = ds.where(ds['Range'] != 550, np.nan)
ds = ds.where(ds['Range'] != 300, np.nan)

plt.pcolormesh(ds['time'], ds['Range'], ds['backscatter'], shading = 'auto', vmin = 0)
plt.colorbar()
plt.gcf().autofmt_xdate()
plt.xlabel('time')
plt.ylabel('Range')
plt.show()

## Calculate the gradient
ds['bs_grad'] = ds['backscatter'].differentiate(coord = 'Range')

plt.pcolormesh(ds['time'], ds['Range'], ds['bs_grad'], shading = 'auto')
plt.colorbar()
plt.gcf().autofmt_xdate()
plt.xlabel('time')
plt.ylabel('Range')
plt.show()

## Drop nan values
ds = ds.dropna(dim = 'time', how = 'all') ## Works
ds = ds.dropna(dim = 'Range', how = 'all') ## Does nothing ? 

plt.pcolormesh(ds['time'], ds['Range'], ds['bs_grad'], shading = 'auto')
plt.colorbar()
plt.gcf().autofmt_xdate()
plt.xlabel('time')
plt.ylabel('Range')
plt.show()

## Try to find location of the minimum values as in example 1.
min_height = ds['bs_grad'].argmin(dim='Range').values  # get the index of the minimum backscatter
display(min_height)
min_height = ds['Range'].isel(Range=min_height)  # extract the corresponding height
display(min_height)

plt.scatter(ds['time'], min_height),
plt.gcf().autofmt_xdate()
plt.show()

解决思路与方案 #

错误原因 #

dropna(dim='time', how='all')只会删除所有变量在该time时刻全为NaN的切片，但如果数据集里的其他变量（如temperature、humidity）在该time维度有有效值，就不会被删除，导致对应的bs_grad仍然全NaN，调用argmin时触发错误。

修复步骤 #

1. 针对bs_grad单独过滤无效time切片
   筛选出bs_grad在Range维度上不全为NaN的time值，再过滤数据集：

   # 找出bs_grad在Range维度上存在有效值的time
   valid_time = ds['bs_grad'].notnull().any(dim='Range')
   # 过滤并删除无效的time切片
   ds = ds.where(valid_time, drop=True)
   

2. 优化滚动平滑后的NaN处理
   滚动均值会在边界生成NaN，加上手动设置的NaN，容易出现全NaN切片。可以在平滑后先对backscatter做过滤：

   # 平滑后保留backscatter不全为NaN的time和Range
   ds = ds.where(ds['backscatter'].notnull().any(dim='Range'), drop=True)
   ds = ds.where(ds['backscatter'].notnull().any(dim='time'), drop=True)
   

3. 修复后的关键代码片段
   在复现示例的计算梯度之后添加过滤步骤：

   ## Calculate the gradient
   ds['bs_grad'] = ds['backscatter'].differentiate(coord = 'Range')
   
   # 过滤bs_grad全为NaN的time切片
   valid_time = ds['bs_grad'].notnull().any(dim='Range')
   ds = ds.where(valid_time, drop=True)
   
   ## 查找最小值位置
   min_height = ds['bs_grad'].argmin(dim='Range').values
   min_height = ds['Range'].isel(Range=min_height)
   
   plt.scatter(ds['time'], min_height)
   plt.gcf().autofmt_xdate()
   plt.show()
   

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