Python气象数据处理实战从netCDF到温度垂直廓线图气象数据分析是环境科学、气候研究等领域的重要基础技能。对于刚接触气象数据处理的Python用户来说面对原始的netCDF格式数据往往会感到无从下手。本文将带你完整走一遍从数据读取、处理到可视化的全流程重点讲解如何用Python绘制专业的大气温度垂直廓线图。1. 环境准备与数据理解在开始之前我们需要确保工作环境配置正确。气象数据处理通常需要几个核心Python库pip install netCDF4 numpy matplotlibnetCDF(Network Common Data Form)是气象领域广泛使用的数据格式它采用类似于字典的结构存储多维数组数据。一个典型的netCDF气象数据文件可能包含以下变量维度变量名描述典型维度顺序longitude经度坐标一维数组latitude纬度坐标一维数组time时间坐标一维数组level气压层高度一维数组t温度数据[time,level,lat,lon]提示使用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})在实际分析中我们通常只需要特定区域和时间的数据。以下是空间筛选的关键步骤确定目标区域假设我们需要南京周边600公里范围内的数据计算经纬度范围地球表面1度约111公里因此600公里≈5.4度找到最近格点使用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_lon5.4))[0] lat_range np.where((lat target_lat-5.4) (lat target_lat5.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, colorred, linewidth2, markero, markersize5) # 专业图表设置 plt.xlabel(Temperature (K), fontsize12) plt.ylabel(Pressure Level (hPa), fontsize12) plt.title(Vertical Temperature Profile\nNanjing Area, 2023-04-17 12:00, fontsize14, pad20) # 反转y轴表示高度增加 plt.gca().invert_yaxis() # 设置对数坐标 (大气压力通常在对数尺度下显示) plt.yscale(log) # 添加网格线 plt.grid(True, whichboth, linestyle--, alpha0.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, dpi300) 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, :], axis1) plt.plot(lat_profile, levels, labelfLat: {lat[lat_idx]:.1f}°N, linewidth2) # 其他专业设置 plt.legend(fontsize12, framealpha0.9) plt.yscale(log) plt.gca().invert_yaxis() plt.grid(True, whichboth, linestyle--, alpha0.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--, labelISA 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, projectionproj) # 添加地图要素 ax.add_feature(cfeature.COASTLINE) ax.add_feature(cfeature.BORDERS, linestyle:) ax.add_feature(cfeature.LAKES, alpha0.5) # 绘制温度场 contour ax.contourf(lon, lat, temp_850, levels20, cmapRdBu_r, transformproj) # 添加色标 plt.colorbar(contour, labelTemperature (K)) # 标记南京位置 ax.plot(target_lon, target_lat, ro, markersize8, transformproj) ax.text(target_lon1, target_lat, Nanjing, fontsize12, transformproj)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, mmapTrue)选择性读取只提取需要的变量# 只读取温度变量 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(filenameweather_analysis.log, levellogging.INFO, format%(asctime)s - %(levelname)s - %(message)s) try: # 数据处理代码 logging.info(开始处理数据文件) except Exception as e: logging.error(f处理失败: {str(e)})
保姆级教程：用Python和netCDF4处理气象数据，绘制南京周边温度垂直廓线图
Python气象数据处理实战从netCDF到温度垂直廓线图气象数据分析是环境科学、气候研究等领域的重要基础技能。对于刚接触气象数据处理的Python用户来说面对原始的netCDF格式数据往往会感到无从下手。本文将带你完整走一遍从数据读取、处理到可视化的全流程重点讲解如何用Python绘制专业的大气温度垂直廓线图。1. 环境准备与数据理解在开始之前我们需要确保工作环境配置正确。气象数据处理通常需要几个核心Python库pip install netCDF4 numpy matplotlibnetCDF(Network Common Data Form)是气象领域广泛使用的数据格式它采用类似于字典的结构存储多维数组数据。一个典型的netCDF气象数据文件可能包含以下变量维度变量名描述典型维度顺序longitude经度坐标一维数组latitude纬度坐标一维数组time时间坐标一维数组level气压层高度一维数组t温度数据[time,level,lat,lon]提示使用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})在实际分析中我们通常只需要特定区域和时间的数据。以下是空间筛选的关键步骤确定目标区域假设我们需要南京周边600公里范围内的数据计算经纬度范围地球表面1度约111公里因此600公里≈5.4度找到最近格点使用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_lon5.4))[0] lat_range np.where((lat target_lat-5.4) (lat target_lat5.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, colorred, linewidth2, markero, markersize5) # 专业图表设置 plt.xlabel(Temperature (K), fontsize12) plt.ylabel(Pressure Level (hPa), fontsize12) plt.title(Vertical Temperature Profile\nNanjing Area, 2023-04-17 12:00, fontsize14, pad20) # 反转y轴表示高度增加 plt.gca().invert_yaxis() # 设置对数坐标 (大气压力通常在对数尺度下显示) plt.yscale(log) # 添加网格线 plt.grid(True, whichboth, linestyle--, alpha0.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, dpi300) 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, :], axis1) plt.plot(lat_profile, levels, labelfLat: {lat[lat_idx]:.1f}°N, linewidth2) # 其他专业设置 plt.legend(fontsize12, framealpha0.9) plt.yscale(log) plt.gca().invert_yaxis() plt.grid(True, whichboth, linestyle--, alpha0.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--, labelISA 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, projectionproj) # 添加地图要素 ax.add_feature(cfeature.COASTLINE) ax.add_feature(cfeature.BORDERS, linestyle:) ax.add_feature(cfeature.LAKES, alpha0.5) # 绘制温度场 contour ax.contourf(lon, lat, temp_850, levels20, cmapRdBu_r, transformproj) # 添加色标 plt.colorbar(contour, labelTemperature (K)) # 标记南京位置 ax.plot(target_lon, target_lat, ro, markersize8, transformproj) ax.text(target_lon1, target_lat, Nanjing, fontsize12, transformproj)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, mmapTrue)选择性读取只提取需要的变量# 只读取温度变量 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(filenameweather_analysis.log, levellogging.INFO, format%(asctime)s - %(levelname)s - %(message)s) try: # 数据处理代码 logging.info(开始处理数据文件) except Exception as e: logging.error(f处理失败: {str(e)})
