Python气象数据处理实战：从NetCDF文件解析到南京温度垂直廓线可视化
气象数据分析是环境科学、地理信息系统等领域的重要基础技能。对于刚接触气象数据处理的Python用户来说，NetCDF格式文件往往让人望而生畏。本文将手把手带你完成从数据读取到可视化的全流程，重点解决实际工作中的痛点问题。
1. 理解NetCDF气象数据格式
NetCDF（Network Common Data Form）是气象领域广泛使用的科学数据存储格式。与CSV或Excel不同，NetCDF采用多维数组结构存储数据，特别适合处理具有多个维度（如时间、高度、经纬度）的气象要素。
典型的ERA5再分析数据包含以下维度变量：
| 变量名 | 描述 | 单位 | 典型维度顺序 |
|---|---|---|---|
| longitude | 经度坐标 | 度 | 独立维度 |
| latitude | 纬度坐标 | 度 | 独立维度 |
| time | 时间坐标 | 小时 | 独立维度 |
| level | 气压层高度 | hPa | 独立维度 |
| t | 温度变量 | K | time×level×lat×lon |
常见痛点：初次接触时容易混淆维度顺序，特别是在不同数据源中维度排列可能不同。建议先用以下代码检查数据结构：
from netCDF4 import Dataset nc_file = "era5_temperature.nc" data = Dataset(nc_file, "r") print(data.variables.keys()) # 查看所有变量 print(data.variables['t'].dimensions) # 查看温度变量的维度顺序2. 精准定位南京地区的温度数据
提取特定位置的气象数据需要解决三个关键问题：
- 时间索引定位：ERA5数据通常采用UTC时间，需要注意时区转换
- 空间索引匹配：经纬度网格可能不完全匹配目标位置
- 空间范围确定：是否需考虑周边区域平均值
2.1 时间维度处理
假设我们需要获取2023年4月17日12:00（北京时间，即UTC+8的04:00）的数据：
import numpy as np from datetime import datetime # 转换北京时间到UTC target_time = datetime(2023, 4, 17, 12) # 北京时间12:00 utc_time = datetime(2023, 4, 17, 4) # 对应UTC时间04:00 # 在NetCDF时间变量中查找最接近的时间点 nc_times = data.variables['time'][:] time_units = data.variables['time'].units # 如"hours since 1900-01-01 00:00:00"注意：不同数据源的时间基准可能不同，务必检查time变量的units属性
2.2 空间维度处理
南京市中心约位于118.8°E，32.1°N。由于数据网格可能不完全匹配，需要找到最接近的网格点：
def find_nearest_index(array, value): """找到数组中与给定值最接近的索引""" return np.argmin(np.abs(array - value)) lon = data.variables['longitude'][:] # 通常范围是0-360或-180到180 lat = data.variables['latitude'][:] # 通常从北到南排列 # 处理经度范围差异 if lon.max() > 180: # 0-360范围 nanjing_lon = 118.8 else: # -180-180范围 nanjing_lon = 118.8 if 118.8 <= 180 else 118.8 - 360 nanjing_lat = 32.1 lon_idx = find_nearest_index(lon, nanjing_lon) lat_idx = find_nearest_index(lat, nanjing_lat)3. 提取垂直温度廓线数据
获取温度垂直廓线需要考虑气压层维度。ERA5通常提供从地表到高空的多层数据：
levels = data.variables['level'][:] # 气压层，单位hPa temperature = data.variables['t'] # 温度，单位K # 获取特定时间、位置的温度垂直廓线 time_idx = 42 # 假设已经确定的时间索引 vertical_temp = temperature[time_idx, :, lat_idx, lon_idx] # 转换为摄氏度（如原始数据为开尔文） vertical_temp_c = vertical_temp - 273.15高级技巧：如需考虑周边区域平均值，可以扩展空间范围：
# 定义搜索半径（约50km） radius_deg = 0.45 # 约50km lon_mask = (lon >= nanjing_lon - radius_deg) & (lon <= nanjing_lon + radius_deg) lat_mask = (lat >= nanjing_lat - radius_deg) & (lat <= nanjing_lat + radius_deg) # 提取区域平均温度廓线 regional_temp = np.mean(temperature[time_idx, :, lat_mask, :][:, :, lon_mask], axis=(1,2))4. 专业级温度廓线可视化
使用matplotlib绘制温度垂直廓线图时，需要注意气象绘图的专业规范：
import matplotlib.pyplot as plt from matplotlib.ticker import MultipleLocator plt.style.use('seaborn') # 使用美观的样式 fig, ax = plt.subplots(figsize=(8, 10)) # 绘制温度曲线 ax.plot(vertical_temp_c, levels, color='red', linewidth=2, marker='o', markersize=6) # 设置坐标轴 ax.set_xlabel('Temperature (°C)', fontsize=12) ax.set_ylabel('Pressure Level (hPa)', fontsize=12) ax.set_title('Nanjing Temperature Profile\n2023-04-17 12:00 CST', fontsize=14, pad=20) # 专业气象图的y轴设置 ax.set_yscale('log') # 对数坐标更符合大气结构 ax.invert_yaxis() # 高压在下，低压在上 ax.yaxis.set_major_locator(MultipleLocator(100)) # 每100hPa一个主刻度 ax.grid(True, which='both', linestyle='--', alpha=0.6) # 添加高度参考线 for level in [850, 700, 500, 300, 200]: ax.axhline(y=level, color='gray', linestyle=':', alpha=0.5) plt.tight_layout() plt.savefig('nanjing_temp_profile.png', dpi=300, bbox_inches='tight') plt.show()关键改进点：
- 使用对数坐标表示气压高度，更符合实际大气结构
- 添加了标准等压线参考（850hPa、700hPa等）
- 优化了图形尺寸和字体大小，适合学术使用
- 提高了输出图像的分辨率（300dpi）
5. 处理常见问题与异常情况
实际工作中经常会遇到各种数据问题，以下是几种典型场景的解决方案：
5.1 数据缺失处理
# 检查缺失值 missing_value = temperature.missing_value if hasattr(temperature, 'missing_value') else None if missing_value is not None: vertical_temp = np.ma.masked_equal(vertical_temp, missing_value) print(f"发现{np.sum(vertical_temp.mask)}个缺失值") # 简单插值处理 from scipy.interpolate import interp1d valid_mask = ~vertical_temp.mask if hasattr(vertical_temp, 'mask') else slice(None) valid_levels = levels[valid_mask] valid_temp = vertical_temp[valid_mask] if len(valid_levels) < len(levels): interp_func = interp1d(valid_levels, valid_temp, bounds_error=False, fill_value="extrapolate") vertical_temp = interp_func(levels)5.2 多时间点对比分析
比较不同时间的温度廓线有助于分析日变化：
# 获取同一天00:00和12:00的数据 time_indices = [36, 42] # 假设对应00:00和12:00的索引 times_of_day = ['00:00', '12:00'] fig, ax = plt.subplots(figsize=(8, 10)) for idx, time_label in zip(time_indices, times_of_day): temp_profile = temperature[idx, :, lat_idx, lon_idx] - 273.15 ax.plot(temp_profile, levels, label=f'{time_label} CST', linewidth=2) ax.set(xlabel='Temperature (°C)', ylabel='Pressure Level (hPa)', title='Nanjing Temperature Profile Comparison\n2023-04-17') ax.invert_yaxis() ax.set_yscale('log') ax.legend() ax.grid()5.3 批量处理多个位置
对于区域研究，可能需要处理多个站点的数据：
stations = { 'Nanjing': (118.8, 32.1), 'Shanghai': (121.5, 31.2), 'Hefei': (117.3, 31.8) } fig, ax = plt.subplots(figsize=(10, 10)) for name, (st_lon, st_lat) in stations.items(): lon_idx = find_nearest_index(lon, st_lon) lat_idx = find_nearest_index(lat, st_lat) temp_profile = temperature[time_idx, :, lat_idx, lon_idx] - 273.15 ax.plot(temp_profile, levels, label=name) ax.set(xlabel='Temperature (°C)', ylabel='Pressure Level (hPa)', title='Temperature Profile Comparison\n2023-04-17 12:00 CST') ax.invert_yaxis() ax.set_yscale('log') ax.legend() ax.grid()6. 进阶技巧与性能优化
当处理长时间序列或高分辨率数据时，需要考虑内存和计算效率问题。
6.1 分块处理大数据文件
# 使用xarray处理大型NetCDF文件 import xarray as xr # 延迟加载数据，不立即读入内存 ds = xr.open_dataset('large_era5_data.nc', chunks={'time': 10}) # 只加载南京附近区域 nanjing_slice = ds.sel( longitude=slice(118, 120), latitude=slice(33, 31) # 注意纬度顺序 ) # 计算季节平均 seasonal_mean = nanjing_slice['t'].groupby('time.season').mean(dim='time')6.2 使用Dask加速计算
import dask.array as da # 将NetCDF数据转为Dask数组 dask_temp = da.from_array(temperature, chunks=(24, 10, 100, 100)) # 并行计算区域平均 def compute_regional_mean(temp_array, lon_range, lat_range): lon_mask = (lon >= lon_range[0]) & (lon <= lon_range[1]) lat_mask = (lat >= lat_range[0]) & (lat <= lat_range[1]) return temp_array[:, :, lat_mask, :][:, :, :, lon_mask].mean(axis=(2,3)) # 南京周边区域 nanjing_region = ((118, 120), (31, 33)) result = compute_regional_mean(dask_temp, *nanjing_region) # 触发实际计算 regional_mean = result.compute()6.3 结果保存与共享
将处理后的数据保存为新NetCDF文件：
from netCDF4 import Dataset import os output_file = 'nanjing_temp_profiles.nc' if os.path.exists(output_file): os.remove(output_file) with Dataset(output_file, 'w') as dst: # 创建维度 dst.createDimension('level', len(levels)) dst.createDimension('time', None) # 无限制长度 # 创建变量 level_var = dst.createVariable('level', 'f4', ('level',)) level_var[:] = levels level_var.units = 'hPa' time_var = dst.createVariable('time', 'f8', ('time',)) time_var[:] = nc_times[time_indices] time_var.units = time_units temp_var = dst.createVariable('temp', 'f4', ('time', 'level')) temp_var[:, :] = temperature[time_indices, :, lat_idx, lon_idx] temp_var.units = 'K' temp_var.long_name = 'Air temperature' # 添加全局属性 dst.title = "Nanjing temperature profiles extracted from ERA5" dst.source = "ERA5 reanalysis data" dst.location = f"lon={lon[lon_idx]:.2f}, lat={lat[lat_idx]:.2f}"
