保姆级教程：用Python和netCDF4处理气象数据，绘制南京周边温度垂直廓线图

Python气象数据处理实战：从netCDF到温度垂直廓线图

气象数据分析是环境科学、气候研究等领域的重要基础技能。对于刚接触气象数据处理的Python用户来说，面对原始的netCDF格式数据往往会感到无从下手。本文将带你完整走一遍从数据读取、处理到可视化的全流程，重点讲解如何用Python绘制专业的大气温度垂直廓线图。

1. 环境准备与数据理解

在开始之前，我们需要确保工作环境配置正确。气象数据处理通常需要几个核心Python库：

pip install netCDF4 numpy matplotlib

netCDF(Network Common Data Form)是气象领域广泛使用的数据格式，它采用类似于字典的结构存储多维数组数据。一个典型的netCDF气象数据文件可能包含以下变量维度：

提示：使用ncdump -h filename.nc命令可以快速查看netCDF文件的结构信息，无需实际读取数据。

2. 数据读取与空间筛选

让我们从读取netCDF文件开始。以下代码展示了如何安全地打开文件并检查其内容：

import numpy as np from netCDF4 import Dataset # 使用上下文管理器确保文件正确关闭 with Dataset("era5_temperature.nc", "r") as nc_data: # 查看文件中所有变量 print(nc_data.variables.keys()) # 获取基本维度信息 lon = nc_data.variables["longitude"][:] lat = nc_data.variables["latitude"][:] levels = nc_data.variables["level"][:] times = nc_data.variables["time"][:] # 打印维度信息 print(f"经度范围: {lon.min()}~{lon.max()}") print(f"纬度范围: {lat.min()}~{lat.max()}") print(f"可用气压层: {levels}")

在实际分析中，我们通常只需要特定区域和时间的数据。以下是空间筛选的关键步骤：

1. 确定目标区域：假设我们需要南京周边600公里范围内的数据
2. 计算经纬度范围：地球表面1度约111公里，因此600公里≈5.4度
3. 找到最近格点：使用numpy的argmin方法

# 南京大致坐标 (118.8°E, 32.1°N) target_lon, target_lat = 118.8, 32.1 # 找到最近格点索引 lon_idx = np.argmin(np.abs(lon - target_lon)) lat_idx = np.argmin(np.abs(lat - target_lat)) # 计算600公里范围内的索引 lon_range = np.where((lon >= target_lon-5.4) & (lon <= target_lon+5.4))[0] lat_range = np.where((lat >= target_lat-5.4) & (lat <= target_lat+5.4))[0]

3. 时间维度处理与数据提取

气象数据通常包含多个时间步长，我们需要准确选择目标时间点。ERA5再分析数据的time变量通常以"小时数自1900-01-01"的格式存储：

import datetime # 假设我们要分析2023年4月17日12:00的数据 target_time = datetime.datetime(2023,4,17,12) # 计算与基准时间的差值 base_time = datetime.datetime(1900,1,1) delta = target_time - base_time hours_since = delta.total_seconds()/3600 # 在时间变量中找到最接近的索引 time_idx = np.argmin(np.abs(times - hours_since))

提取目标区域和时间的温度数据：

# 提取温度数据 (注意维度顺序通常是[time,level,lat,lon]) temperature = nc_data.variables["t"][time_idx, :, :, :] # 筛选目标区域 region_temp = temperature[:, lat_range, :][:, :, lon_range] # 计算区域平均温度廓线 vertical_profile = np.mean(region_temp, axis=(1,2))

4. 绘制专业温度垂直廓线图

温度垂直廓线图是分析大气热力结构的重要工具。使用matplotlib绘制时需要注意几个专业细节：

import matplotlib.pyplot as plt plt.figure(figsize=(8,10)) # 绘制温度廓线 plt.plot(vertical_profile, levels, color='red', linewidth=2, marker='o', markersize=5) # 专业图表设置 plt.xlabel('Temperature (K)', fontsize=12) plt.ylabel('Pressure Level (hPa)', fontsize=12) plt.title('Vertical Temperature Profile\nNanjing Area, 2023-04-17 12:00', fontsize=14, pad=20) # 反转y轴表示高度增加 plt.gca().invert_yaxis() # 设置对数坐标 (大气压力通常在对数尺度下显示) plt.yscale('log') # 添加网格线 plt.grid(True, which='both', linestyle='--', alpha=0.5) # 自定义yticks标签 plt.yticks([1000, 850, 700, 500, 300, 200, 100], ['1000', '850', '700', '500', '300', '200', '100']) plt.tight_layout() plt.savefig('temperature_profile.png', dpi=300) plt.show()

对于更复杂的分析，我们可能想比较不同纬度位置的温度廓线：

# 选择几个代表性纬度 selected_lats = [lat_range[0], lat_range[len(lat_range)//2], lat_range[-1]] plt.figure(figsize=(10,12)) for i, lat_idx in enumerate(selected_lats): # 提取该纬度所有经度的平均温度 lat_profile = np.mean(region_temp[:, i, :], axis=1) plt.plot(lat_profile, levels, label=f'Lat: {lat[lat_idx]:.1f}°N', linewidth=2) # 其他专业设置 plt.legend(fontsize=12, framealpha=0.9) plt.yscale('log') plt.gca().invert_yaxis() plt.grid(True, which='both', linestyle='--', alpha=0.5)

5. 高级分析与可视化技巧

5.1 添加标准大气参考线

为了更好评估温度异常，可以添加标准大气温度廓线作为参考：

# 国际标准大气温度模型 (简化版) def isa_temperature(altitude_km): if altitude_km <= 11: return 288.15 - 6.5 * altitude_km # 对流层递减率 elif altitude_km <= 20: return 216.65 # 平流层下层等温 else: return 216.65 + (altitude_km-20) # 平流层上层增温 # 将气压层转换为近似高度 (简化计算) levels_km = -7 * np.log(levels/1000) # 计算ISA温度 isa_temp = [isa_temperature(h) for h in levels_km] # 添加到图中 plt.plot(isa_temp, levels, 'k--', label='ISA Reference')

5.2 多变量综合分析

结合风场数据可以分析温度场与环流的关系。假设我们有U/V风分量数据：

u_wind = nc_data.variables["u"][time_idx, :, lat_range, lon_range] v_wind = nc_data.variables["v"][time_idx, :, lat_range, lon_range] # 计算风速和风向 wind_speed = np.sqrt(u_wind**2 + v_wind**2) wind_direction = np.arctan2(v_wind, u_wind) * 180/np.pi # 创建带风矢量的温度剖面图 fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16,10)) # 温度剖面 ax1.plot(vertical_profile, levels, 'r-') ax1.set_yscale('log') ax1.invert_yaxis() # 风速剖面 ax2.plot(wind_speed.mean(axis=(1,2)), levels, 'b-') ax2.set_yscale('log') ax2.invert_yaxis()

5.3 使用Cartopy添加地理背景

对于空间分布分析，Cartopy库可以添加专业地图背景：

import cartopy.crs as ccrs import cartopy.feature as cfeature # 创建地图投影 proj = ccrs.PlateCarree() # 绘制850hPa温度场 level_idx = np.where(levels == 850)[0][0] temp_850 = temperature[level_idx, :, :] fig = plt.figure(figsize=(12,8)) ax = fig.add_subplot(111, projection=proj) # 添加地图要素 ax.add_feature(cfeature.COASTLINE) ax.add_feature(cfeature.BORDERS, linestyle=':') ax.add_feature(cfeature.LAKES, alpha=0.5) # 绘制温度场 contour = ax.contourf(lon, lat, temp_850, levels=20, cmap='RdBu_r', transform=proj) # 添加色标 plt.colorbar(contour, label='Temperature (K)') # 标记南京位置 ax.plot(target_lon, target_lat, 'ro', markersize=8, transform=proj) ax.text(target_lon+1, target_lat, 'Nanjing', fontsize=12, transform=proj)

6. 常见问题与性能优化

处理大规模气象数据时可能会遇到性能瓶颈。以下是几个实用技巧：

• 分块处理：对于TB级数据，使用netCDF4的chunking功能

# 优化后的数据读取方式 with Dataset("large_file.nc") as nc: temp = nc.variables["t"] # 每次只读取一个时间步 for i in range(len(times)): temp_slice = temp[i,:,:,:]

• 使用Dask并行处理：

import dask.array as da # 创建dask数组 temp_dask = da.from_array(nc_data.variables["t"], chunks=(1,10,100,100)) # 并行计算区域平均 region_avg = temp_dask[:, lat_range, :][:, :, lon_range].mean(axis=(2,3))

• 内存映射：避免一次性加载大文件

# 使用内存映射模式 ds = Dataset("big_data.nc", "r", mmap=True)

• 选择性读取：只提取需要的变量

# 只读取温度变量 with Dataset("data.nc") as nc: temp = nc.variables["t"][...]

处理气象数据时，我经常发现数据质量控制至关重要。一个实用的做法是在绘图前添加基本的数据检查：

# 数据有效性检查 if np.ma.is_masked(vertical_profile).any(): print("警告：数据包含缺失值！") # 温度合理性检查 if (vertical_profile < 180).any() or (vertical_profile > 350).any(): print("警告：发现异常温度值！")

对于需要长期运行的分析任务，建议添加日志记录：

import logging logging.basicConfig(filename='weather_analysis.log', level=logging.INFO, format='%(asctime)s - %(levelname)s - %(message)s') try: # 数据处理代码 logging.info("开始处理数据文件") except Exception as e: logging.error(f"处理失败: {str(e)}")