用Python处理气象数据：从NetCDF文件到南京周边温度垂直廓线图（附完整代码）

用Python处理气象数据：从NetCDF文件到南京周边温度垂直廓线图（附完整代码）

气象数据分析是理解大气现象的重要手段，而Python凭借其丰富的数据处理库和可视化工具，已成为气象科研和业务工作中的首选语言。本文将手把手教你如何用Python处理NetCDF格式的气象数据，从文件读取到最终绘制出专业的温度垂直廓线图，每个步骤都配有详细解释和可复用的代码模块。

1. 环境准备与数据理解

在开始处理气象数据前，需要确保Python环境中安装了必要的科学计算库。推荐使用Anaconda发行版，它已经集成了我们所需的大部分工具：

conda install numpy matplotlib netCDF4

NetCDF（Network Common Data Form）是气象领域广泛使用的数据格式，它以自描述的方式存储多维数组数据。一个典型的NetCDF气象数据文件可能包含以下变量：

• 经度（longitude）：描述数据点的东西向位置
• 纬度（latitude）：描述数据点的南北向位置
• 时间（time）：记录观测或模拟的时间点
• 高度（level）：表示大气垂直层次
• 温度（t）：存储各时空点的温度值

提示：使用ncdump -h filename.nc命令可以快速查看NetCDF文件的结构信息，无需加载完整数据。

2. 数据加载与初步探索

让我们从加载NetCDF文件开始，逐步理解数据结构：

import numpy as np from netCDF4 import Dataset # 加载NetCDF文件 file_path = "2010_air_12.nc" data = Dataset(file_path, "r") # 查看文件中的变量列表 print("文件包含的变量:", data.variables.keys()) # 获取各维度信息 lon = data.variables['longitude'][:] lat = data.variables['latitude'][:] time = data.variables['time'][:] level = data.variables['level'][:] temp = data.variables['t'][:] print(f"经度范围: {lon.min()}°E ~ {lon.max()}°E") print(f"纬度范围: {lat.min()}°N ~ {lat.max()}°N") print(f"时间点数量: {len(time)}") print(f"垂直层次: {level}")

通过这段代码，我们可以快速了解数据的时空覆盖范围和垂直分辨率。例如，输出可能显示：

3. 时空范围筛选

气象数据通常覆盖大范围区域和长时间序列，但我们的分析可能只需要特定时空范围内的数据。以下代码演示如何筛选南京周边600公里范围内的数据：

# 定义目标位置和时间 target_lon = 119 # 南京经度 target_lat = 32 # 南京纬度 target_date = "2010-04-17 12:00" # 目标日期时间 # 计算600公里对应的经纬度范围（约5.4度） radius_deg = 5.4 # 找到最接近目标点的索引 lon_idx = np.argmin(np.abs(lon - target_lon)) lat_idx = np.argmin(np.abs(lat - target_lat)) # 获取600公里范围内的索引 lon_mask = (lon >= (target_lon - radius_deg)) & (lon <= (target_lon + radius_deg)) lat_mask = (lat >= (target_lat - radius_deg)) & (lat <= (target_lat + radius_deg)) # 筛选数据 subset_lon = lon[lon_mask] subset_lat = lat[lat_mask] subset_temp = temp[:, :, lat_mask, :][:, :, :, lon_mask]

注意：地球表面1度约对应111公里，但这个换算在靠近极地地区会有变化。对于中纬度地区，这个近似足够精确。

4. 温度垂直廓线绘制

温度垂直廓线图是分析大气热力结构的重要工具。下面我们绘制南京周边各纬度位置在同一经度上的温度垂直分布：

import matplotlib.pyplot as plt # 选择固定经度（取范围内的第一个经度点） fixed_lon_idx = 0 # 准备绘图 fig, ax = plt.subplots(figsize=(10, 8)) # 为每个纬度绘制一条廓线 for i in range(len(subset_lat)): # 获取该纬度上所有高度的温度数据 temp_profile = subset_temp[0, :, i, fixed_lon_idx] # 绘制廓线 ax.plot(temp_profile, level, label=f'Lat: {subset_lat[i]:.1f}°N', linewidth=2) # 美化图形 ax.set_xlabel('Temperature (K)', fontsize=12) ax.set_ylabel('Pressure Level (hPa)', fontsize=12) ax.set_title('Temperature Vertical Profiles around Nanjing\n2010-04-17 12:00', fontsize=14, pad=20) ax.invert_yaxis() # 反转y轴，使高度从下往上增加 ax.grid(True, linestyle='--', alpha=0.6) ax.legend(bbox_to_anchor=(1.05, 1), loc='upper left') plt.tight_layout() plt.savefig('temperature_profiles.png', dpi=300, bbox_inches='tight') plt.show()

这段代码会生成一张包含多条温度廓线的专业图表，每条线代表不同纬度位置在同一经度上的温度垂直分布。关键绘图技巧包括：

• 坐标轴反转：大气压力随高度增加而降低，反转y轴符合气象学惯例
• 图例外置：避免遮挡数据曲线
• 网格线：辅助读者准确读取数值
• 高分辨率输出：适合学术报告和论文使用

5. 进阶分析与可视化增强

基础廓线图已经能揭示很多信息，但我们可以通过以下增强手段提取更多洞见：

5.1 温度异常分析

计算各高度层相对于区域平均温度的异常：

# 计算区域平均温度廓线 mean_temp_profile = subset_temp[0].mean(axis=(1, 2)) # 计算各点的温度异常 temp_anomaly = subset_temp[0] - mean_temp_profile[:, np.newaxis, np.newaxis] # 绘制异常图 plt.figure(figsize=(12, 6)) contour = plt.contourf(subset_lon, level, temp_anomaly.mean(axis=2).T, levels=20, cmap='RdBu_r') plt.colorbar(contour, label='Temperature Anomaly (K)') plt.xlabel('Longitude (°E)') plt.ylabel('Pressure Level (hPa)') plt.title('Temperature Anomaly Relative to Regional Mean') plt.gca().invert_yaxis()

5.2 三维可视化

使用mplot3d工具包创建三维温度场可视化：

from mpl_toolkits.mplot3d import Axes3D # 准备网格数据 X, Y = np.meshgrid(subset_lon, level) fig = plt.figure(figsize=(12, 10)) ax = fig.add_subplot(111, projection='3d') # 绘制表面图 surf = ax.plot_surface(X, Y, subset_temp[0, :, 5, :].T, cmap='coolwarm', linewidth=0, antialiased=False) ax.set_xlabel('Longitude (°E)') ax.set_ylabel('Pressure Level (hPa)') ax.set_zlabel('Temperature (K)') ax.set_title('3D Temperature Field') fig.colorbar(surf, shrink=0.5, aspect=5)

5.3 风场叠加分析

如果有风场数据，可以将其叠加到温度场分析中：

# 假设我们已经加载了u和v风分量 u_wind = data.variables['u'][0] # 经向风 v_wind = data.variables['v'][0] # 纬向风 # 绘制温度填色和风矢量 plt.figure(figsize=(12, 8)) plt.contourf(subset_lon, level, subset_temp[0, :, :, 5], 20, cmap='coolwarm') plt.colorbar(label='Temperature (K)') # 每隔5个点绘制一个风矢量 skip = 5 plt.quiver(subset_lon[::skip], level[::skip], u_wind[::skip, ::skip], v_wind[::skip, ::skip], color='black', scale=200) plt.title('Temperature with Wind Vectors') plt.xlabel('Longitude (°E)') plt.ylabel('Pressure Level (hPa)') plt.gca().invert_yaxis()

6. 实用技巧与常见问题

在实际工作中处理气象数据时，有几个实用技巧值得掌握：

6.1 处理大型NetCDF文件

当处理GB级别的大型气象数据集时，内存管理变得尤为重要：

• 分块读取：利用NetCDF4库的分块读取功能

# 只读取需要的变量和范围 temp = data.variables['t'][0, :, 100:200, 50:150]

• 使用Dask：对于超大型数据集，可以使用Dask进行延迟计算

import dask.array as da temp_dask = da.from_array(data.variables['t'][:], chunks=(1, 10, 100, 100))

6.2 时间处理技巧

气象数据中的时间变量通常使用特殊格式（如"hours since 1900-01-01"），需要正确转换：

from netCDF4 import num2date # 转换时间变量为datetime对象 times = num2date(data.variables['time'][:], data.variables['time'].units) # 找到特定日期的索引 target_date = np.datetime64('2010-04-17T12:00') time_idx = np.where(times == target_date)[0][0]

6.3 数据质量控制

气象数据中常包含缺失值或异常值，处理方法是：

# 假设缺失值用-9999表示 temp = data.variables['t'][:] temp = np.ma.masked_equal(temp, -9999) # 计算时自动忽略缺失值 mean_temp = temp.mean(axis=0)

6.4 可视化优化建议

• 色标选择：使用适合气象数据的色标（如'coolwarm'、'viridis'）
• 标注规范：包括单位、数据来源、制图日期等信息
• 多图组合：使用subplots将相关变量放在一起比较

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 6)) # 左图：温度 im1 = ax1.contourf(lon, lat, temp[0, 0], cmap='coolwarm') fig.colorbar(im1, ax=ax1, label='Temperature (K)') ax1.set_title('Surface Temperature') # 右图：相对湿度 im2 = ax2.contourf(lon, lat, rh[0, 0], cmap='Blues') fig.colorbar(im2, ax=ax2, label='Relative Humidity (%)') ax2.set_title('Surface Relative Humidity')