1000hPa气象数据输出全为NaN问题求助及原因排查

问题说明 #

处理气象场数据时，读取1000hPa层级数据输出全为NaN，但切换至850hPa等其他层级可正常获取有效数据。

原始代码 #

import numpy as np
import cartopy.io.shapereader as shpreader
import xarray as xr
import xarray as xr
def load_data(time_array, var_name, ds):
    data = sum(ds[var_name].loc[time,1000, 25:36, 115:128] for time in time_array) / len(time_array)
    return data

time_array = [
'2010-07-16T03',
'2011-06-09T04',
'2011-06-09T20',
'2011-08-18T10',
'2014-07-12T06',
'2016-06-19T10',
'2016-06-20T19',
'2016-06-21T08',
'2017-04-06T08',
'2017-06-29T03',
'2017-06-30T02',
'2017-08-17T08',
'2017-08-18T03',
'2017-08-19T03',
'2017-09-24T18']

ds1 = xr.open_dataset('D:/filter_nc/minor3hbig.nc')
T1 = load_data(time_array, 't', ds1)
print(T1)

输出结果 #

<xarray.DataArray 't' (lat: 45, lon: 53)>
    array([[nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan],
           ...,
           [nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan]], dtype=float32)
    Coordinates:
        level    int32 1000
       * lon      (lon) float32 115.0 115.2 115.5 115.8 ... 127.2 127.5 127.8 128.0
       * lat      (lat) float32 25.0 25.25 25.5 25.75 26.0 ... 35.25 35.5 35.75 36.0
         time     datetime64[ns] 2017-09-24T18:00:00

排查步骤 #

1. 验证数据集1000hPa层级完整性
   先检查整个数据集的1000hPa层级是否本身全为NaN：

   print(ds1['t'].sel(level=1000).isnull().all())
   

   如果返回True，说明数据源中1000hPa层级无有效数据。

2. 检查单个时间点的1000hPa数据
   逐个验证指定时间点是否存在有效数据：

   for t in time_array:
       has_valid = ds1['t'].sel(time=t, level=1000).notnull().any()
       print(f"时间{t}的1000hPa数据是否有效：{has_valid}")
   

   若部分时间点返回False，说明这些时刻的1000hPa数据缺失。

3. 确认层级索引匹配
   检查level维度的实际值和类型，避免因类型不匹配导致索引错误：

   print(ds1['level'])
   

   如果层级存储为浮点数（如1000.0），用整数1000索引会导致匹配失败，返回全NaN。

4. 检查区域筛选合理性
   确认指定的经纬度范围（25-36°N，115-128°E）在1000hPa层级是否有数据：

   print(ds1['t'].sel(level=1000, lat=slice(25,36), lon=slice(115,128)).notnull().any())
   

   部分地区海拔高于1000hPa，地面气压低于1000hPa，该区域的1000hPa层级自然无数据。

解决方法 #

1. 数据源问题
   若数据集本身缺失1000hPa数据，更换包含该层级有效数据的数据源，或确认原始数据的预处理流程是否有误。

2. 时间点过滤
   自动过滤无有效数据的时间点后再计算平均：

   def load_data(time_array, var_name, ds, level):
       data_subset = ds[var_name].sel(time=time_array, level=level, lat=slice(25,36), lon=slice(115,128))
       # 筛选有有效数据的时间点
       valid_times = data_subset.notnull().any(dim=['lat','lon']).where(lambda x: x, drop=True).time
       if len(valid_times) == 0:
           print("无有效数据的时间点")
           return data_subset.mean(dim='time')
       # 计算平均
       return data_subset.sel(time=valid_times).mean(dim='time')
   

3. 修正层级索引
   使用sel方法替代loc，自动匹配层级标签，避免类型或值不匹配问题：

   # 替换原代码中的loc索引为sel
   ds[var_name].sel(time=time, level=1000, lat=slice(25,36), lon=slice(115,128))
   

   若层级是浮点数，直接传入1000.0即可。

4. 调整区域范围
   若目标区域海拔普遍高于1000hPa，缩小经纬度范围至低海拔区域，或确认该区域是否真的存在1000hPa气象数据。

优化后代码示例 #

import numpy as np
import xarray as xr

def load_data(time_array, var_name, ds, level):
    # 使用sel方法安全索引
    data_subset = ds[var_name].sel(time=time_array, level=level, lat=slice(25, 36), lon=slice(115, 128))
    
    # 检查全局数据情况
    if data_subset.isnull().all():
        print(f"{level}hPa层级在指定时间和区域内无有效数据")
        return data_subset.mean(dim='time')
    
    # 过滤无有效数据的时间点
    valid_times = data_subset.notnull().any(dim=['lat','lon']).where(lambda x: x, drop=True).time
    missing_count = len(time_array) - len(valid_times)
    if missing_count > 0:
        print(f"过滤了{missing_count}个无有效数据的时间点")
    
    # 计算平均
    return data_subset.sel(time=valid_times).mean(dim='time')

time_array = [
'2010-07-16T03',
'2011-06-09T04',
'2011-06-09T20',
'2011-08-18T10',
'2014-07-12T06',
'2016-06-19T10',
'2016-06-20T19',
'2016-06-21T08',
'2017-04-06T08',
'2017-06-29T03',
'2017-06-30T02',
'2017-08-17T08',
'2017-08-18T03',
'2017-08-19T03',
'2017-09-24T18']

ds1 = xr.open_dataset('D:/filter_nc/minor3hbig.nc')

# 先做基础检查
print("1000hPa全局数据是否全NaN：", ds1['t'].sel(level=1000).isnull().all())

T1 = load_data(time_array, 't', ds1, 1000)
print(T1)

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