用Python处理气象数据从NetCDF文件到南京周边温度垂直廓线图附完整代码气象数据分析是理解大气现象的重要手段而Python凭借其丰富的数据处理库和可视化工具已成为气象科研和业务工作中的首选语言。本文将手把手教你如何用Python处理NetCDF格式的气象数据从文件读取到最终绘制出专业的温度垂直廓线图每个步骤都配有详细解释和可复用的代码模块。1. 环境准备与数据理解在开始处理气象数据前需要确保Python环境中安装了必要的科学计算库。推荐使用Anaconda发行版它已经集成了我们所需的大部分工具conda install numpy matplotlib netCDF4NetCDFNetwork 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})通过这段代码我们可以快速了解数据的时空覆盖范围和垂直分辨率。例如输出可能显示维度最小值最大值分辨率经度100°E130°E0.5°纬度20°N40°N0.5°时间0小时8759小时每小时高度1000hPa100hPa不等距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, labelfLat: {subset_lat[i]:.1f}°N, linewidth2) # 美化图形 ax.set_xlabel(Temperature (K), fontsize12) ax.set_ylabel(Pressure Level (hPa), fontsize12) ax.set_title(Temperature Vertical Profiles around Nanjing\n2010-04-17 12:00, fontsize14, pad20) ax.invert_yaxis() # 反转y轴使高度从下往上增加 ax.grid(True, linestyle--, alpha0.6) ax.legend(bbox_to_anchor(1.05, 1), locupper left) plt.tight_layout() plt.savefig(temperature_profiles.png, dpi300, bbox_inchestight) 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(axis2).T, levels20, cmapRdBu_r) plt.colorbar(contour, labelTemperature 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, projection3d) # 绘制表面图 surf ax.plot_surface(X, Y, subset_temp[0, :, 5, :].T, cmapcoolwarm, linewidth0, antialiasedFalse) 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, shrink0.5, aspect5)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, cmapcoolwarm) plt.colorbar(labelTemperature (K)) # 每隔5个点绘制一个风矢量 skip 5 plt.quiver(subset_lon[::skip], level[::skip], u_wind[::skip, ::skip], v_wind[::skip, ::skip], colorblack, scale200) 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(axis0)6.4 可视化优化建议色标选择使用适合气象数据的色标如coolwarm、viridis标注规范包括单位、数据来源、制图日期等信息多图组合使用subplots将相关变量放在一起比较fig, (ax1, ax2) plt.subplots(1, 2, figsize(16, 6)) # 左图温度 im1 ax1.contourf(lon, lat, temp[0, 0], cmapcoolwarm) fig.colorbar(im1, axax1, labelTemperature (K)) ax1.set_title(Surface Temperature) # 右图相对湿度 im2 ax2.contourf(lon, lat, rh[0, 0], cmapBlues) fig.colorbar(im2, axax2, labelRelative Humidity (%)) ax2.set_title(Surface Relative Humidity)
用Python处理气象数据：从NetCDF文件到南京周边温度垂直廓线图（附完整代码）
用Python处理气象数据从NetCDF文件到南京周边温度垂直廓线图附完整代码气象数据分析是理解大气现象的重要手段而Python凭借其丰富的数据处理库和可视化工具已成为气象科研和业务工作中的首选语言。本文将手把手教你如何用Python处理NetCDF格式的气象数据从文件读取到最终绘制出专业的温度垂直廓线图每个步骤都配有详细解释和可复用的代码模块。1. 环境准备与数据理解在开始处理气象数据前需要确保Python环境中安装了必要的科学计算库。推荐使用Anaconda发行版它已经集成了我们所需的大部分工具conda install numpy matplotlib netCDF4NetCDFNetwork 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})通过这段代码我们可以快速了解数据的时空覆盖范围和垂直分辨率。例如输出可能显示维度最小值最大值分辨率经度100°E130°E0.5°纬度20°N40°N0.5°时间0小时8759小时每小时高度1000hPa100hPa不等距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, labelfLat: {subset_lat[i]:.1f}°N, linewidth2) # 美化图形 ax.set_xlabel(Temperature (K), fontsize12) ax.set_ylabel(Pressure Level (hPa), fontsize12) ax.set_title(Temperature Vertical Profiles around Nanjing\n2010-04-17 12:00, fontsize14, pad20) ax.invert_yaxis() # 反转y轴使高度从下往上增加 ax.grid(True, linestyle--, alpha0.6) ax.legend(bbox_to_anchor(1.05, 1), locupper left) plt.tight_layout() plt.savefig(temperature_profiles.png, dpi300, bbox_inchestight) 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(axis2).T, levels20, cmapRdBu_r) plt.colorbar(contour, labelTemperature 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, projection3d) # 绘制表面图 surf ax.plot_surface(X, Y, subset_temp[0, :, 5, :].T, cmapcoolwarm, linewidth0, antialiasedFalse) 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, shrink0.5, aspect5)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, cmapcoolwarm) plt.colorbar(labelTemperature (K)) # 每隔5个点绘制一个风矢量 skip 5 plt.quiver(subset_lon[::skip], level[::skip], u_wind[::skip, ::skip], v_wind[::skip, ::skip], colorblack, scale200) 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(axis0)6.4 可视化优化建议色标选择使用适合气象数据的色标如coolwarm、viridis标注规范包括单位、数据来源、制图日期等信息多图组合使用subplots将相关变量放在一起比较fig, (ax1, ax2) plt.subplots(1, 2, figsize(16, 6)) # 左图温度 im1 ax1.contourf(lon, lat, temp[0, 0], cmapcoolwarm) fig.colorbar(im1, axax1, labelTemperature (K)) ax1.set_title(Surface Temperature) # 右图相对湿度 im2 ax2.contourf(lon, lat, rh[0, 0], cmapBlues) fig.colorbar(im2, axax2, labelRelative Humidity (%)) ax2.set_title(Surface Relative Humidity)
