从ERA5小时数据到日均数据：一个高效批量处理的Python实践

1. 从小时数据到日均数据的必要性

处理气象数据时，我们常常会遇到数据量过大的问题。以ERA5-Land小时数据为例，一个月的全球数据就可能达到几十GB。这种高频数据虽然信息丰富，但直接使用往往效率低下，尤其当我们需要分析长期趋势或进行空间分析时。

我曾在处理一个气候分析项目时，电脑因为内存不足崩溃了三次。后来发现，将小时数据聚合为日均数据后，文件体积缩小了24倍，处理速度提升了近30倍。这就是为什么我们需要掌握高效的数据聚合方法——它不仅能节省存储空间，还能大幅提升后续分析效率。

日均数据的另一个优势是能平滑短期波动。比如研究农作物生长与温度的关系时，日均温度比小时温度更具参考价值。当然，聚合方式要根据具体需求选择，常见的有日均值、日最大值、日最小值等。

2. 环境准备与数据理解

2.1 必备工具安装

工欲善其事，必先利其器。我们需要准备以下Python库：

• xarray：处理多维数组数据的利器
• rasterio：地理空间数据读写专家
• numpy：科学计算基础包

安装命令很简单：

pip install xarray rasterio numpy

如果你使用Anaconda，也可以用conda安装：

conda install -c conda-forge xarray rasterio

2.2 认识ERA5数据格式

ERA5数据通常以NetCDF格式存储，这是一种自描述的科学数据格式。用xarray打开一个示例文件看看结构：

import xarray as xr data = xr.open_dataset('example.nc') print(data)

你会看到类似这样的输出：

Dimensions: (time: 24, latitude: 721, longitude: 1440) Coordinates: * time (time) datetime64[ns] 2023-01-01 ... 2023-01-01T23:00:00 * latitude (latitude) float64 90.0 89.75 89.5 ... -89.75 -90.0 * longitude (longitude) float64 0.0 0.25 0.5 ... 359.5 359.75 Data variables: t2m (time, latitude, longitude) float32 ...

这里t2m就是我们要处理的2米高度气温数据，它有24个时间点（每小时一个），721个纬度点和1440个经度点。

3. 核心处理流程详解

3.1 批量读取与日均计算

处理大量文件时，我们需要考虑内存管理。直接一次性读取所有文件可能会导致内存溢出。我的经验是逐个文件处理：

import os import xarray as xr input_folder = 'ERA5_hourly' output_folder = 'ERA5_daily' os.makedirs(output_folder, exist_ok=True) for filename in os.listdir(input_folder): if filename.endswith('.nc'): filepath = os.path.join(input_folder, filename) # 使用xarray的延迟加载特性 with xr.open_dataset(filepath) as ds: # 计算日均值 daily_mean = ds['t2m'].mean(dim='time') # 从文件名提取日期 date_part = filename.split('_')[-1].split('.')[0] output_path = os.path.join(output_folder, f'daily_{date_part}.nc') # 保存结果 daily_mean.to_netcdf(output_path)

这个脚本的关键点：

1. 使用with语句确保文件正确关闭
2. 利用xarray的延迟加载，避免过早加载全部数据
3. 保持原始数据的坐标信息

3.2 转换为GeoTIFF格式

地理信息系统(GIS)通常使用GeoTIFF格式。转换时需要特别注意坐标系统：

import rasterio from rasterio.transform import from_origin def save_geotiff(data, lons, lats, output_path): # 计算分辨率和边界 res = lons[1] - lons[0] lon_min, lon_max = lons.min(), lons.max() lat_min, lat_max = lats.min(), lats.max() # 创建仿射变换 transform = from_origin(lon_min, lat_max, res, res) # 写入GeoTIFF with rasterio.open( output_path, 'w', driver='GTiff', height=data.shape[0], width=data.shape[1], count=1, dtype=str(data.dtype), crs='EPSG:4326', # WGS84坐标系统 transform=transform ) as dst: dst.write(data, 1)

使用时只需传入numpy数组和对应的经纬度坐标：

save_geotiff(daily_mean.values, lons, lats, 'output.tif')

4. 性能优化技巧

4.1 并行处理加速

当文件数量很多时，可以使用多进程加速。Python的multiprocessing模块是个不错的选择：

from multiprocessing import Pool def process_file(filename): # 处理单个文件的代码 pass if __name__ == '__main__': nc_files = [f for f in os.listdir('input') if f.endswith('.nc')] with Pool(processes=4) as pool: # 使用4个进程 pool.map(process_file, nc_files)

在我的测试中，4进程处理100个文件的时间从约15分钟缩短到了4分钟。注意不要设置过多进程，否则可能因I/O瓶颈反而降低效率。

4.2 内存优化策略

处理全球高分辨率数据时，内存可能成为瓶颈。可以采用分块处理策略：

# 使用xarray的分块加载 ds = xr.open_dataset('large.nc', chunks={'time': 24, 'latitude': 100, 'longitude': 100}) # 分块计算日均值 daily_mean = ds['t2m'].chunk({'time': 24}).mean(dim='time')

这种方法特别适合处理多年数据。我曾用这个方法处理了10年的ERA5数据，内存使用始终保持在2GB以下。

5. 常见问题与解决方案

5.1 时间坐标处理

ERA5数据的时间坐标有时会带来困扰。比如遇到闰秒或时区问题时，可以这样标准化：

import pandas as pd # 确保时间坐标是规范的datetime64 ds['time'] = pd.to_datetime(ds['time'].values)

如果数据跨越夏令时变更，建议统一使用UTC时间以避免混淆。

5.2 缺失值处理

气象数据常有缺失值（如海洋区域的陆地变量）。处理时需要注意：

# 计算日均值时跳过缺失值 daily_mean = ds['t2m'].mean(dim='time', skipna=True) # 或者用插值填充 filled = ds['t2m'].interpolate_na(dim='time', method='linear')

在保存为GeoTIFF时，最好明确指定缺失值的编码：

with rasterio.open(...) as dst: dst.write(data.filled(np.nan), 1) # 用NaN表示缺失 dst.nodata = np.nan

6. 完整案例演示

让我们看一个从原始数据到最终成果的完整流程。假设我们有2023年1月的前三天的每小时数据：

1. 首先创建处理脚本process_era5.py：

import os import xarray as xr import numpy as np import rasterio from rasterio.transform import from_origin from datetime import datetime def process_day(input_path, output_dir): """处理单日数据""" with xr.open_dataset(input_path) as ds: # 计算日均 daily = ds['t2m'].mean(dim='time') # 准备输出路径 date_str = datetime.strptime(ds.time.dt.strftime('%Y%m%d').values[0], '%Y%m%d') output_path = os.path.join(output_dir, f't2m_daily_{date_str}.tif') # 保存为GeoTIFF lons, lats = ds.longitude.values, ds.latitude.values res = lons[1] - lons[0] transform = from_origin(lons.min(), lats.max(), res, res) with rasterio.open( output_path, 'w', driver='GTiff', height=len(lats), width=len(lons), count=1, dtype=str(daily.dtype), crs='EPSG:4326', transform=transform ) as dst: dst.write(daily.values, 1) if __name__ == '__main__': input_dir = 'era5_hourly' output_dir = 'era5_daily' os.makedirs(output_dir, exist_ok=True) for filename in sorted(os.listdir(input_dir)): if filename.endswith('.nc'): process_day(os.path.join(input_dir, filename), output_dir)

2. 运行脚本：

python process_era5.py

3. 检查输出：

import matplotlib.pyplot as plt import rasterio with rasterio.open('era5_daily/t2m_daily_20230101.tif') as src: data = src.read(1) plt.imshow(data, cmap='coolwarm') plt.colorbar(label='Temperature (K)') plt.title('Daily Mean Temperature 2023-01-01') plt.show()

这个流程在我的ThinkPad P53上处理一天数据只需约5秒，内存占用不超过500MB。对于批量处理多年数据，可以考虑加入进度条显示：

from tqdm import tqdm files = [f for f in os.listdir(input_dir) if f.endswith('.nc')] for filename in tqdm(sorted(files), desc='Processing'): process_day(os.path.join(input_dir, filename), output_dir)

7. 进阶应用建议

掌握了基础处理方法后，可以考虑以下进阶方向：

1. 异常值检测：在聚合前加入质量控制步骤

# 排除超出合理范围的值 valid_t2m = ds['t2m'].where((ds['t2m'] > 200) & (ds['t2m'] < 350))

2. 时空聚合：同时聚合空间和时间维度

# 计算区域日均值 region_mean = ds['t2m'].mean(dim=['time', 'longitude', 'latitude'])

3. 自定义聚合函数：不只是平均值

# 计算日温度范围 daily_range = ds['t2m'].max(dim='time') - ds['t2m'].min(dim='time')

4. 元数据保留：确保结果数据包含完整的属性信息

daily_mean.attrs.update({ 'description': 'Daily mean 2m temperature', 'units': 'K', 'processing_history': 'Generated from ERA5 hourly data' })

在实际项目中，我建议将处理流程封装成函数或类，方便复用。比如可以创建一个ERA5Processor类，封装所有处理逻辑，支持不同的聚合方法和输出格式。