使用xarray计算NetCDF栅格逐元最大值遇错求助

问题需求 #

计算目录中所有NetCDF文件对应像素的逐元最大值（忽略空值，等效ArcGIS Cell Statistics工具的最大值功能），后续需基于周批量最大值数据集，统计全局最大值、平均值及非NaN值的数量。

现有代码问题 #

原代码使用exec动态创建变量，后续调用exec时返回None，导致ds1和ds2均为NoneType，触发np.fmax的类型错误。错误信息如下：

Traceback (most recent call last):
  File "D:\xarray_stuff\DataProcess.py", line 287, in <module>
    ds_ge_0_final = np.fmax(ds1, ds2)
TypeError: '>=' not supported between instances of 'NoneType' and 'NoneType'

正确实现方案 #

核心改进点 #

• 放弃exec动态创建变量，改用列表存储数据集，更安全且易维护
• 利用xarray原生的concat合并多数据集，通过max(dim=..., skipna=True)直接计算逐元最大值（自动忽略NaN）
• 后续统计需求可直接通过xarray的聚合函数实现

完整代码示例 #

import numpy as np
import xarray as xr
import os

# 数据目录
filepath = 'weekBatch/'

# 用列表存储处理后的数据集，替代exec动态创建变量
ds_ge_0_list = []
ds_ge_15_list = []

for filename in os.listdir(filepath):
    # 拼接完整文件路径
    full_path = os.path.join(filepath, filename)
    ds = xr.open_dataset(full_path)
    
    # 筛选MESH>=0的记录，仅保留MESH变量（减少内存占用）
    ds_ge_0 = ds.where(ds.MESH >= 0.0)['MESH']
    ds_ge_0_list.append(ds_ge_0)
    
    # 筛选MESH>=15的记录
    ds_ge_15 = ds.where(ds.MESH >= 15.0)['MESH']
    ds_ge_15_list.append(ds_ge_15)

# 合并所有数据集，新增一个维度（命名为'file'）
combined_ge_0 = xr.concat(ds_ge_0_list, dim='file')
# 计算逐元最大值，skipna=True自动忽略NaN
ds_ge_0_final = combined_ge_0.max(dim='file', skipna=True)

# 同理处理MESH>=15的场景
combined_ge_15 = xr.concat(ds_ge_15_list, dim='file')
ds_ge_15_final = combined_ge_15.max(dim='file', skipna=True)

周数据集后续统计 #

假设已生成包含周维度（如week）的数据集weekly_ds，可直接用xarray聚合函数完成统计：

# 周数据集全局最大值
weekly_global_max = weekly_ds.max(dim=['week', 'latitude', 'longitude'], skipna=True)
# 周数据集全局平均值
weekly_global_mean = weekly_ds.mean(dim=['week', 'latitude', 'longitude'], skipna=True)
# 统计每个像素的非NaN周数据数量
weekly_pixel_count = weekly_ds.count(dim='week')
# 全局非NaN值总数量
weekly_total_count = weekly_ds.count(dim=['week', 'latitude', 'longitude'])

替代工具推荐 #

如果xarray使用仍有障碍，可考虑以下工具：

• GDAL Command Line：用gdal_merge.py结合-n nan和-max参数实现逐元最大值计算
• CDO (Climate Data Operators)：专门处理气候数据的命令行工具，cdo ensmax命令可直接计算多文件逐元最大值
• PyVista：针对栅格数据的可视化与数值计算工具，支持NetCDF格式

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