<think>根据要求，我们需要输出两张图片。第一张图片是水汽通量矢量与散度填色图，第二张图片是水汽通量散度填色与水汽通量大小等值线图。我们将调整代码，分别调用两个绘图函数，并分别保存为两个文件。修改思路：1.保留两个独立的绘图函数：`plot_vector_and_divergence`和`plot_divergence_and_magnitude`。2.每个函数都会创建自己的图形并保存为独立的图片文件。3.主程序依次调用这两个函数。注意：为了避免重复代码，我们仍然使用辅助函数提取数据切片和时间信息。以下是修改后的代码：</think>我将修改代码以输出两张独立的图片：第一张为水汽通量矢量+散度填色图，第二张为水汽通量散度填色+通量大小等值线图。以下是修改后的代码：
```python
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from metpy.calc import divergence
from metpy.units import units
import os
import traceback
# 1. 读取NetCDF文件
def read_netcdf(file_path):
"""读取并处理NetCDF文件"""
ds = xr.open_dataset(file_path)
print("文件变量:", list(ds.data_vars))
# 检查常见变量名并重命名
var_mapping = {
'uwnd': 'u', 'u_wind': 'u',
'vwnd': 'v', 'v_wind': 'v',
'shum': 'q', 'qair': 'q', 'hus': 'q'
}
for old_name, new_name in var_mapping.items():
if old_name in ds.data_vars:
ds = ds.rename({old_name: new_name})
return ds
# 2. 计算水汽通量
def calculate_flux(ds, level_idx=0):
"""计算水汽通量分量"""
# 检查必要变量是否存在
required_vars = ['u', 'v', 'q']
missing_vars = [var for var in required_vars if var not in ds.data_vars]
if missing_vars:
raise ValueError(f"缺少必要变量: {missing_vars}")
# 添加单位 (metpy要求)
u = ds['u'] * units('m/s')
v = ds['v'] * units('m/s')
q = ds['q'] * units('kg/kg')
# 计算水汽通量分量
qu = u * q # x方向水汽通量 (kg/(m·s))
qv = v * q # y方向水汽通量 (kg/(m·s))
# 计算水汽通量大小 (模)
q_magnitude = np.sqrt(qu**2 + qv**2)
return qu, qv, q_magnitude
# 3. 计算水汽通量散度
def calculate_flux_divergence(qu, qv):
"""计算水汽通量散度"""
# 使用metpy计算散度
div = divergence(qu, qv) * 1e7 # 缩放为常用单位 (10^-7 kg/(m²·s))
return div
# 4.1 第一张图：水汽通量矢量+散度填色
def plot_vector_divergence(qu, qv, div, time_idx=0, level_idx=0, output_dir="."):
"""绘制水汽通量矢量和水汽通量散度图"""
fig = plt.figure(figsize=(16, 12))
proj = ccrs.PlateCarree()
ax = fig.add_subplot(1, 1, 1, projection=proj)
# 添加地图要素
add_map_features(ax)
# 提取数据和时间信息
qu_slice, div_slice, time_str, level_value = extract_data_and_info(qu, div, time_idx, level_idx)
qv_slice = qv.isel(time=time_idx)
# 绘制散度填色
contour = plot_divergence_fill(ax, div_slice)
# 添加色带
cbar = fig.colorbar(contour, ax=ax, shrink=0.7, pad=0.05)
cbar.set_label('水汽通量散度 (10$^{-7}$ kg/(m$^2$·s))', fontsize=12)
# 绘制水汽通量矢量
Q = plot_flux_vectors(ax, qu_slice, qv_slice)
# 添加标题
ax.set_title(f'水汽通量矢量与散度\n时间: {time_str}{level_value}',
fontsize=16)
# 保存图片
output_path = os.path.join(output_dir, "moisture_flux_vector.png")
plt.savefig(output_path, dpi=300, bbox_inches='tight')
print(f"矢量图已保存为: {output_path}")
plt.close(fig) # 关闭图形以释放内存
# 4.2 第二张图：水汽通量散度填色+通量大小等值线
def plot_divergence_magnitude(div, q_magnitude, time_idx=0, level_idx=0, output_dir="."):
"""绘制水汽通量散度和水汽通量大小图"""
fig = plt.figure(figsize=(16, 12))
proj = ccrs.PlateCarree()
ax = fig.add_subplot(1, 1, 1, projection=proj)
# 添加地图要素
add_map_features(ax)
# 提取数据和时间信息
div_slice, q_mag_slice, time_str, level_value = extract_data_and_info(div, q_magnitude, time_idx, level_idx)
# 绘制散度填色
contourf = plot_divergence_fill(ax, div_slice)
# 添加色带
cbar = fig.colorbar(contourf, ax=ax, shrink=0.7, pad=0.05)
cbar.set_label('水汽通量散度 (10$^{-7}$ kg/(m$^2$·s))', fontsize=12)
# 绘制水汽通量大小等值线
contour = plot_magnitude_contour(ax, q_mag_slice)
# 添加等值线标签
plt.clabel(contour, inline=True, fontsize=10, fmt='%1.1f')
# 添加标题
ax.set_title(f'水汽通量散度与通量大小\n时间: {time_str}{level_value}',
fontsize=16)
# 保存图片
output_path = os.path.join(output_dir, "moisture_flux_magnitude.png")
plt.savefig(output_path, dpi=300, bbox_inches='tight')
print(f"等值线图已保存为: {output_path}")
plt.close(fig) # 关闭图形以释放内存
# 辅助函数：添加地图要素
def add_map_features(ax):
"""添加常见地图要素"""
ax.add_feature(cfeature.COASTLINE, linewidth=0.8)
ax.add_feature(cfeature.BORDERS, linestyle=':', linewidth=0.5)
ax.add_feature(cfeature.OCEAN, color='lightblue', alpha=0.3)
ax.add_feature(cfeature.LAND, color='beige', alpha=0.2)
ax.gridlines(draw_labels=True, linestyle='--', alpha=0.5)
# 辅助函数：绘制散度填色
def plot_divergence_fill(ax, div_slice):
"""绘制水汽通量散度填色图"""
return ax.contourf(
div_slice.lon,
div_slice.lat,
div_slice.values,
levels=np.linspace(-5, 5, 21),
cmap='RdBu',
extend='both',
transform=ccrs.PlateCarree()
)
# 辅助函数：绘制水汽通量矢量
def plot_flux_vectors(ax, qu_slice, qv_slice):
"""绘制水汽通量矢量图"""
# 确定网格点
if hasattr(qu_slice, 'lon') and hasattr(qu_slice, 'lat'):
lons, lats = qu_slice.lon, qu_slice.lat
else:
lons, lats = np.meshgrid(qu_slice.longitude, qu_slice.latitude)
# 自适应跳过点数
skip = max(1, len(qu_slice.lon) // 30)
# 绘制矢量箭头
Q = ax.quiver(
lons[::skip, ::skip],
lats[::skip, ::skip],
qu_slice.values[::skip, ::skip],
qv_slice.values[::skip, ::skip],
scale=300,
color='k',
transform=ccrs.PlateCarree()
)
# 添加矢量比例尺
ax.quiverkey(Q, 0.05, 0.05, 100,
'100 kg/(m·s)',
labelpos='E',
coordinates='axes',
fontproperties={'size': 10})
return Q
# 辅助函数：绘制通量大小等值线
def plot_magnitude_contour(ax, q_mag_slice):
"""绘制水汽通量大小等值线"""
return ax.contour(
q_mag_slice.lon,
q_mag_slice.lat,
q_mag_slice.values,
levels=10, # 自动选择10个等值线级别
colors='k',
linewidths=1.5,
transform=ccrs.PlateCarree()
)
# 辅助函数：提取数据和信息
def extract_data_and_info(var1, var2, time_idx, level_idx):
"""提取指定时间和高度的数据切片及时间/高度信息"""
# 提取数据切片
var1_slice, var2_slice = extract_data_slice(var1, var2, time_idx, level_idx)
# 提取时间信息
time_str = "未知时间"
if 'time' in var1_slice.coords:
try:
time_str = str(var1_slice.time.values)[:19]
except:
pass
# 提取高度信息
level_value = ""
if 'level' in var1_slice.coords:
level_value = f" | 高度: {var1_slice.level.values} hPa"
elif 'pressure' in var1_slice.coords:
level_value = f" | 高度: {var1_slice.pressure.values} hPa"
elif 'isobaric' in var1_slice.coords:
level_value = f" | 高度: {var1_slice.isobaric.values} hPa"
return var1_slice, var2_slice, time_str, level_value
# 辅助函数：提取数据切片
def extract_data_slice(var1, var2, time_idx, level_idx):
"""提取指定时间和高度的数据切片"""
try:
# 尝试使用level索引
var1_slice = var1.isel(time=time_idx, level=level_idx)
var2_slice = var2.isel(time=time_idx, level=level_idx)
except:
# 如果没有level维度，尝试pressure或isobaric
try:
var1_slice = var1.isel(time=time_idx, pressure=level_idx)
var2_slice = var2.isel(time=time_idx, pressure=level_idx)
except:
# 最后尝试直接使用time索引
var1_slice = var1.isel(time=time_idx)
var2_slice = var2.isel(time=time_idx)
return var1_slice, var2_slice
# 主程序
if __name__ == "__main__":
# 使用您提供的文件路径
file_path = r'C:/Users/25636/Desktop/bi/ziji/zhengduanfx/first/ziliao/200307.nc'
output_dir = r'C:/Users/25636/Desktop/bi/ziji/zhengduanfx/first/ziliao/output'
# 创建输出目录
os.makedirs(output_dir, exist_ok=True)
print(f"正在处理文件: {file_path}")
print(f"输出目录: {output_dir}")
try:
# 处理数据
ds = read_netcdf(file_path)
qu, qv, q_magnitude = calculate_flux(ds)
div = calculate_flux_divergence(qu, qv)
# 绘制第一张图：矢量+散度填色
print("正在绘制水汽通量矢量图...")
plot_vector_divergence(qu, qv, div, output_dir=output_dir)
# 绘制第二张图：散度填色+通量大小等值线
print("正在绘制水汽通量等值线图...")
plot_divergence_magnitude(div, q_magnitude, output_dir=output_dir)
print("处理完成！")
except Exception as e:
print(f"处理过程中发生错误: {str(e)}")
print(traceback.format_exc())
print("请检查以下可能原因:")
print("1. 文件路径是否正确")
print("2. 文件是否包含u/v风分量和比湿(q)数据")
print("3. 尝试使用不同的变量名映射")
print("4. 数据维度是否正确")
```
### 修改说明：
1. **模块化重构**：
- 将绘图功能拆分为多个小函数，提高代码可读性和复用性
- 创建了专门的函数：
- `add_map_features()`：添加地图要素
- `plot_divergence_fill()`：绘制散度填色图
- `plot_flux_vectors()`：绘制水汽通量矢量
- `plot_magnitude_contour()`：绘制水汽通量大小等值线
2. **两张独立图片**：
- `plot_vector_divergence()`：生成第一张图片（矢量+散度填色）
- `plot_divergence_magnitude()`：生成第二张图片（散度填色+通量大小等值线）
3. **输出路径管理**：
```python
output_dir = r'C:/Users/25636/Desktop/bi/ziji/zhengduanfx/first/ziliao/output'
os.makedirs(output_dir, exist_ok=True)
```
- 创建专门的输出目录
- 每张图片有明确的文件名
4. **内存管理优化**：
```python
plt.close(fig) # 绘制完成后立即关闭图形释放内存
```
- 避免内存泄漏问题
- 特别适合批量处理大量数据时使用
5. **详细的处理日志**：
```python
print("正在绘制水汽通量矢量图...")
print("矢量图已保存为: ...")
print("正在绘制水汽通量等值线图...")
print("等值线图已保存为: ...")
print("处理完成！")
```
### 输出结果：
程序将输出两张独立的图片：
1. **第一张图片**：`moisture_flux_vector.png`
- 显示水汽通量矢量箭头（表示水汽输送方向）
- 背景为水汽通量散度的填色图（红色表示水汽辐散，蓝色表示水汽辐合）
2. **第二张图片**：`moisture_flux_magnitude.png`
- 显示水汽通量散度的填色图（与第一张相同）
- 叠加水汽通量大小的黑色等值线（表示水汽输送强度）
### 使用建议：
1. **自定义输出目录**：
```python
output_dir = "您/的/输出/路径"
```
2. **调整时间/高度层**：
```python
# 使用不同时间索引
plot_vector_divergence(qu, qv, div, time_idx=1, output_dir=output_dir)
# 使用不同高度层
plot_divergence_magnitude(div, q_magnitude, level_idx=2, output_dir=output_dir)
```
3. **调整绘图参数**：
- 在`plot_divergence_fill()`函数中修改`levels`参数调整散度填色范围
- 在`plot_flux_vectors()`函数中修改`scale`参数调整矢量箭头大小
- 在`plot_magnitude_contour()`函数中修改`levels`参数调整等值线密度
博客主要围绕备份程序的思路展开探讨，忽略了与信息技术无关的每天混日子的思路内容，聚焦于备份程序相关信息技术要点。
扫一扫
2551
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
