从下载到出图：手把手教你用Python处理ERA5再分析数据（以地表温度为例）
从下载到出图：Python处理ERA5地表温度数据的完整指南
气象数据科学正成为跨学科研究的重要工具，而ERA5作为欧洲中期天气预报中心（ECMWF）的第五代再分析数据集，为研究者提供了覆盖全球的高质量气候数据。本文将带你从零开始，使用Python处理ERA5地表温度数据，最终生成专业级可视化图表。
1. 环境准备与数据获取
处理ERA5数据前，需要配置合适的Python环境。推荐使用Anaconda创建独立环境：
conda create -n era5 python=3.9 conda activate era5 conda install -c conda-forge xarray netCDF4 cartopy matplotlib pandas从Copernicus Climate Data Store（CDS）获取数据时，需要先注册账号并安装CDS API客户端：
# 安装CDS API pip install cdsapi # 在~/.cdsapirc文件中配置API密钥 url: https://cds.climate.copernicus.eu/api/v2 key: UID:API-key获取地表温度数据的API请求示例：
import cdsapi c = cdsapi.Client() c.retrieve( 'reanalysis-era5-single-levels', { 'product_type': 'reanalysis', 'variable': '2m_temperature', 'year': '2022', 'month': ['01', '02'], 'day': ['01', '02', '03'], 'time': ['00:00', '06:00', '12:00', '18:00'], 'format': 'netcdf', }, 'era5_temp.nc')2. 数据加载与初步探索
使用xarray加载下载的NetCDF文件：
import xarray as xr ds = xr.open_dataset('era5_temp.nc') print(ds)典型输出结构：
<xarray.Dataset> Dimensions: (longitude: 1440, latitude: 721, time: 24) Coordinates: * longitude (longitude) float32 0.0 0.25 0.5 ... 359.25 359.5 359.75 * latitude (latitude) float32 90.0 89.75 89.5 ... -89.75 -90.0 * time (time) datetime64[ns] 2022-01-01T00:00:00 ... 2022-01-03T18:00:00 Data variables: t2m (time, latitude, longitude) float32 ... Attributes: Conventions: CF-1.6 history: 2023-04-15 12:34:56 GMT by grib_to_netcdf-2.25.1: /opt/ecmw...关键数据操作：
- 查看变量属性：
ds['t2m'].attrs - 选择特定时间点：
ds.sel(time='2022-01-01T12:00:00') - 提取区域数据：
ds.sel(latitude=slice(50, 20), longitude=slice(70, 140))
3. 数据预处理与转换
ERA5数据通常需要以下预处理步骤：
单位转换（开尔文转摄氏度）
ds['t2m_celsius'] = ds['t2m'] - 273.15 ds['t2m_celsius'].attrs = {'units': '°C', 'long_name': '2m temperature in Celsius'}处理缺失值
# 检查缺失值 print(ds['t2m'].isnull().sum()) # 填充或插值 ds_filled = ds['t2m'].fillna(-9999) # 或用插值方法时间维度处理
# 转换为Pandas时间索引 time_series = ds['time'].to_pandas() # 重采样为日平均 daily_mean = ds['t2m'].resample(time='1D').mean()空间统计计算
计算区域平均温度：
# 简单区域平均 region_mean = ds['t2m_celsius'].sel( latitude=slice(40, 20), longitude=slice(110, 120) ).mean(dim=['latitude', 'longitude']) # 加权平均（考虑纬度影响） weights = np.cos(np.deg2rad(ds.latitude)) weighted_mean = (ds['t2m_celsius'] * weights).mean(dim=['latitude', 'longitude']) / weights.mean()4. 高级分析与可视化
时间序列分析
import matplotlib.pyplot as plt # 简单时间序列图 region_mean.plot(figsize=(12, 6)) plt.title('Daily Mean Temperature (40-20°N, 110-120°E)') plt.ylabel('Temperature (°C)') plt.grid() plt.show()空间分布图
使用Cartopy创建专业地图：
import cartopy.crs as ccrs import cartopy.feature as cfeature # 创建地图 fig = plt.figure(figsize=(15, 10)) ax = plt.axes(projection=ccrs.PlateCarree()) # 添加地理特征 ax.add_feature(cfeature.LAND) ax.add_feature(cfeature.OCEAN) ax.add_feature(cfeature.COASTLINE) ax.add_feature(cfeature.BORDERS, linestyle=':') # 绘制温度数据 temp = ds['t2m_celsius'].sel(time='2022-01-01T12:00:00') contour = ax.contourf(temp.longitude, temp.latitude, temp, levels=20, cmap='coolwarm') plt.colorbar(contour, label='Temperature (°C)') # 设置标题和网格 ax.set_title('Surface Temperature at 2m on 2022-01-01 12:00 UTC') ax.gridlines(draw_labels=True) plt.show()温度异常分析
# 计算月平均气候态 climatology = ds['t2m_celsius'].groupby('time.month').mean(dim='time') # 计算异常 anomalies = ds['t2m_celsius'].groupby('time.month') - climatology # 可视化异常 anomalies.sel(time='2022-01-01T12:00:00').plot( figsize=(15, 8), cmap='RdBu_r', vmin=-10, vmax=10, cbar_kwargs={'label': 'Temperature Anomaly (°C)'} ) plt.title('Temperature Anomaly on 2022-01-01 12:00 UTC') plt.show()5. 实战技巧与常见问题解决
内存优化技巧
处理全球高分辨率数据时，内存管理至关重要：
# 分块处理大型数据集 ds_chunked = xr.open_dataset('era5_temp.nc', chunks={'time': 10}) # 使用Dask进行延迟计算 mean_temp = ds_chunked['t2m'].mean(dim='time').compute()性能优化对比
| 方法 | 执行时间 | 内存使用 | 适用场景 |
|---|---|---|---|
| 直接加载 | 快 | 高 | 小数据集 |
| 分块加载 | 中 | 中 | 大数据集 |
| 按需加载 | 慢 | 低 | 特定区域分析 |
常见错误处理
CDS API请求失败：
- 检查API密钥配置
- 确认数据请求参数符合CDS要求
- 添加适当的等待时间避免请求限制
内存不足错误：
# 使用更高效的数据类型 ds['t2m'] = ds['t2m'].astype('float32') # 选择特定时间段或区域 small_ds = ds.sel(time=slice('2022-01-01', '2022-01-10'))投影转换问题：
# 正确设置Cartopy投影 proj = ccrs.LambertConformal(central_longitude=100, central_latitude=30) fig, ax = plt.subplots(subplot_kw={'projection': proj})
6. 扩展应用与自动化
批量处理多年份数据
years = range(2010, 2023) for year in years: filename = f'era5_temp_{year}.nc' if not os.path.exists(filename): c.retrieve( 'reanalysis-era5-single-levels', { 'product_type': 'reanalysis', 'variable': '2m_temperature', 'year': str(year), 'month': ['01', '02'], 'day': ['01', '02', '03'], 'time': ['00:00', '06:00', '12:00', '18:00'], 'format': 'netcdf', }, filename)创建温度变化动画
import matplotlib.animation as animation fig, ax = plt.subplots(figsize=(12, 6), subplot_kw={'projection': ccrs.PlateCarree()}) def update(frame): ax.clear() temp_frame = ds['t2m_celsius'].isel(time=frame) contour = ax.contourf(temp_frame.longitude, temp_frame.latitude, temp_frame, levels=20, cmap='coolwarm') ax.set_title(f'Surface Temperature at {str(ds.time[frame].values)[:19]}') ax.coastlines() return contour ani = animation.FuncAnimation(fig, update, frames=len(ds.time), interval=200) ani.save('temperature_animation.mp4', writer='ffmpeg', dpi=300)集成机器学习分析
from sklearn.ensemble import RandomForestRegressor from sklearn.model_selection import train_test_split # 准备特征和目标变量 X = ds[['t2m', 'time']].to_dataframe().dropna() y = X.pop('t2m') # 时间特征工程 X['hour'] = X.index.hour X['day_of_year'] = X.index.dayofyear # 训练测试分割 X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2) # 训练模型 model = RandomForestRegressor(n_estimators=100) model.fit(X_train, y_train) # 评估 print(f'Test R^2: {model.score(X_test, y_test):.2f}')