Python气象数据处理实战从NetCDF文件解析到南京温度垂直廓线可视化气象数据分析是环境科学、地理信息系统等领域的重要基础技能。对于刚接触气象数据处理的Python用户来说NetCDF格式文件往往让人望而生畏。本文将手把手带你完成从数据读取到可视化的全流程重点解决实际工作中的痛点问题。1. 理解NetCDF气象数据格式NetCDFNetwork Common Data Form是气象领域广泛使用的科学数据存储格式。与CSV或Excel不同NetCDF采用多维数组结构存储数据特别适合处理具有多个维度如时间、高度、经纬度的气象要素。典型的ERA5再分析数据包含以下维度变量变量名描述单位典型维度顺序longitude经度坐标度独立维度latitude纬度坐标度独立维度time时间坐标小时独立维度level气压层高度hPa独立维度t温度变量Ktime×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北京时间即UTC8的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°E32.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, colorred, linewidth2, markero, markersize6) # 设置坐标轴 ax.set_xlabel(Temperature (°C), fontsize12) ax.set_ylabel(Pressure Level (hPa), fontsize12) ax.set_title(Nanjing Temperature Profile\n2023-04-17 12:00 CST, fontsize14, pad20) # 专业气象图的y轴设置 ax.set_yscale(log) # 对数坐标更符合大气结构 ax.invert_yaxis() # 高压在下低压在上 ax.yaxis.set_major_locator(MultipleLocator(100)) # 每100hPa一个主刻度 ax.grid(True, whichboth, linestyle--, alpha0.6) # 添加高度参考线 for level in [850, 700, 500, 300, 200]: ax.axhline(ylevel, colorgray, linestyle:, alpha0.5) plt.tight_layout() plt.savefig(nanjing_temp_profile.png, dpi300, bbox_inchestight) plt.show()关键改进点使用对数坐标表示气压高度更符合实际大气结构添加了标准等压线参考850hPa、700hPa等优化了图形尺寸和字体大小适合学术使用提高了输出图像的分辨率300dpi5. 处理常见问题与异常情况实际工作中经常会遇到各种数据问题以下是几种典型场景的解决方案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_errorFalse, fill_valueextrapolate) 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, labelf{time_label} CST, linewidth2) ax.set(xlabelTemperature (°C), ylabelPressure Level (hPa), titleNanjing 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, labelname) ax.set(xlabelTemperature (°C), ylabelPressure Level (hPa), titleTemperature 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( longitudeslice(118, 120), latitudeslice(33, 31) # 注意纬度顺序 ) # 计算季节平均 seasonal_mean nanjing_slice[t].groupby(time.season).mean(dimtime)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 flon{lon[lon_idx]:.2f}, lat{lat[lat_idx]:.2f}