import os
import re
import numpy as np
import matplotlib.pyplot as plt
from datetime import datetime, timedelta
import matplotlib.dates as mdates
import matplotlib.colors as mcolors
from matplotlib.collections import LineCollection
def plot_cloud_series(folder_path, start_time, end_time, output_dir='.'):
"""
绘制全天云层时间序列图，直接显示所有云底云顶高度，并在下方添加反射率因子图
参数:
folder_path (str): 包含.npy文件的文件夹路径
start_time (str): 起始时间(格式: YYYYMMDDHHMMSS)
end_time (str): 结束时间(格式: YYYYMMDDHHMMSS)
output_dir (str): 输出图片的目录(默认为当前目录)
返回:
None (结果保存为图片文件)
"""
def parse_cloud_value(val_str):
"""解析云高数值"""
if val_str.strip() == '/////':
return np.nan
try:
return float(val_str)
except ValueError:
return np.nan
def decode_dBZ(encoded_str):
# 提取符号位和后4位数字.针对存储的dBZ
sign_bit = encoded_str[0]
abs_value_str = encoded_str[1:]
# 转换为整数并除以100还原浮点数
abs_value = int(abs_value_str) / 100
# 根据符号位决定正负
if sign_bit == '0':
return -abs_value
else:
return abs_value
# 创建输出目录
os.makedirs(output_dir, exist_ok=True)
# 转换时间范围
start_dt = datetime.strptime(start_time, "%Y%m%d%H%M%S")
end_dt = datetime.strptime(end_time, "%Y%m%d%H%M%S")
# 收集所有有效云层数据
cloud_data = [] # 存储所有云层数据
maxdBZ_data = [] # 存储最大dBZ
meandBZ_data = [] # 存储平均DBZ
file_times = [] # 存储文件时间
file_pattern = re.compile(r'.*_(\d{14})_.*\.npy$')
# 统计信息
total_files = 0
valid_layers = 0
for filename in os.listdir(folder_path):
if filename.endswith('.npy'):
match = file_pattern.match(filename)
if match:
file_time_str = match.group(1)
try:
file_time = datetime.strptime(file_time_str, "%Y%m%d%H%M%S")
except ValueError:
continue
if start_dt <= file_time <= end_dt:
total_files += 1
file_path = os.path.join(folder_path, filename)
try:
# 加载数据，确保是一维数组
data = np.load(file_path).flatten()
data = data[0].split(',')
except:
continue
# 记录当前文件的所有云层数据
file_clouds = {'time': file_time, 'layers': []}
# 解析5层云
for layer_idx in range(5):
base_idx = 5 + layer_idx * 2
top_idx = 6 + layer_idx * 2
base = parse_cloud_value(data[base_idx])
top = parse_cloud_value(data[top_idx])
# 跳过无效数据
if base > 0 and top > base and top < 30000:
# 转换为千米
base_km = base / 1000.0
top_km = top / 1000.0
thickness_km = top_km - base_km
file_clouds['layers'].append({
'layer': layer_idx + 1, # 云层编号(1-5)
'base': base_km,
'top': top_km,
'thickness': thickness_km
})
valid_layers += 1
maxdBZ, meandBZ = data[15], data[16]
maxdBZ_data.append(decode_dBZ(maxdBZ))
meandBZ_data.append(decode_dBZ(meandBZ))
# 如果有有效云层数据，则添加到总数据中
if file_clouds['layers']:
cloud_data.append(file_clouds)
file_times.append(file_time)
print(f"处理完成: 共扫描{total_files}个文件, 找到{valid_layers}个有效云层")
if not cloud_data:
print("在指定时间范围内未找到有效云层数据")
return
# 按时间排序
sorted_indices = np.argsort(file_times)
cloud_data = [cloud_data[i] for i in sorted_indices]
file_times = [file_times[i] for i in sorted_indices]
maxdBZ_data = [maxdBZ_data[i] for i in sorted_indices] # 同样排序
meandBZ_data = [meandBZ_data[i] for i in sorted_indices] # 同样排序
# 创建图形和子图
fig = plt.figure(figsize=(20, 12))
# 使用GridSpec创建自定义布局
gs = plt.GridSpec(3, 2, figure=fig, height_ratios=[3, 1, 0.1], width_ratios=[20, 1])
# 云层图 (占据第一行两列)
ax1 = fig.add_subplot(gs[0, 0])
# 反射率图 (占据第二行第一列)
ax2 = fig.add_subplot(gs[1, 0], sharex=ax1)
# 颜色条轴 (占据第二行第二列)
cax = fig.add_subplot(gs[1, 1])
# 为不同云层设置不同颜色
layer_colors = ['blue', 'green', 'red', 'purple', 'orange']
# 存储每层云的数据点，用于后续绘制折线图
all_bases = {1: [], 2: [], 3: [], 4: [], 5: []}
all_tops = {1: [], 2: [], 3: [], 4: [], 5: []}
all_times = {1: [], 2: [], 3: [], 4: [], 5: []}
# 绘制每个数据点的云层
for file_data in cloud_data:
time_val = mdates.date2num(file_data['time'])
for layer in file_data['layers']:
layer_num = layer['layer']
base_val = layer['base']
top_val = layer['top']
# 保存数据点用于折线图
all_times[layer_num].append(time_val)
all_bases[layer_num].append(base_val)
all_tops[layer_num].append(top_val)
# 绘制每层云的云底和云顶折线图
for layer_num in range(1, 2):
if all_times[layer_num]:
# 按时间排序
sorted_idx = np.argsort(all_times[layer_num])
sorted_times = np.array(all_times[layer_num])[sorted_idx]
sorted_bases = np.array(all_bases[layer_num])[sorted_idx]
sorted_tops = np.array(all_tops[layer_num])[sorted_idx]
# 分段绘制折线 - 当时间间隔超过1小时时断开
segments = []
current_segment = []
# 处理第一个点
current_segment.append((sorted_times[0], sorted_bases[0]))
# 遍历所有点，分段处理
for i in range(1, len(sorted_times)):
time_diff = (sorted_times[i] - sorted_times[i - 1]) * 24 # 转换为小时
# 如果时间间隔超过0.5小时，结束当前段并开始新段
if time_diff > 0.25:
if current_segment:
segments.append(current_segment)
current_segment = [(sorted_times[i], sorted_bases[i])]
else:
current_segment.append((sorted_times[i], sorted_bases[i]))
# 添加最后一个段
if current_segment:
segments.append(current_segment)
# 绘制云底线段
for seg_idx, segment in enumerate(segments):
seg_times, seg_bases = zip(*segment)
# 只在第一个段添加标签，避免图例重复
label = f'云层{layer_num} 云底' if seg_idx == 0 else None
ax1.plot(seg_times, seg_bases,
color=layer_colors[layer_num - 1],
linestyle='-', linewidth=1.5,
label=label)
# 设置云层图属性
ax1.set_ylim([0, 10])
ax1.yaxis.set_major_locator(plt.MultipleLocator(1))
ax1.grid(True, linestyle='--', alpha=0.6, which='both')
date_str = start_dt.strftime('%Y年%m月%d日')
ax1.set_title(f'{date_str} 全天云层分布与最大反射率因子', fontsize=16)
ax1.set_ylabel('高度 (km)', fontsize=14)
ax1.legend(loc='upper right', fontsize=10)
# ===================================================
# 绘制最大反射率因子
# ===================================================
# 创建时间点数组（转换为matplotlib数值）
times_num = mdates.date2num(file_times)
# 创建颜色映射：从蓝色（负值）到红色（正值）
cmap = plt.get_cmap('coolwarm')
norm = plt.Normalize(-50, 50)
# 为每个线段创建点
points = np.array([times_num, maxdBZ_data]).T.reshape(-1, 1, 2)
segments = np.concatenate([points[:-1], points[1:]], axis=1)
# 创建线集合，根据反射率值着色
lc = LineCollection(segments, cmap=cmap, norm=norm)
lc.set_array(np.array(maxdBZ_data))
lc.set_linewidth(2)
line = ax2.add_collection(lc)
# 添加颜色条（放在右侧）
cbar = plt.colorbar(line, cax=cax)
cbar.set_label('反射率因子 (dBZ)', fontsize=12)
# 设置反射率图属性
ax2.set_ylim([-50, 50])
ax2.set_ylabel('最大反射率因子 (dBZ)', fontsize=14)
ax2.grid(True, linestyle='--', alpha=0.6)
# 添加阈值线
ax2.axhline(y=0, color='gray', linestyle='-', alpha=0.5)
ax2.axhline(y=20, color='green', linestyle='--', alpha=0.7, label='中雨阈值')
ax2.axhline(y=40, color='red', linestyle='--', alpha=0.7, label='大雨阈值')
ax2.legend(loc='upper right', fontsize=10)
# 设置时间轴（放在上方）
hours = mdates.HourLocator(interval=1)
ax1.xaxis.set_major_locator(hours)
ax1.xaxis.set_major_formatter(mdates.DateFormatter('%H:%M'))
# 将时间刻度显示在中间
ax1.xaxis.tick_bottom()
ax1.xaxis.set_label_position('bottom')
# ax1.xaxis.tick_top()
# ax1.xaxis.set_label_position('top')
# 隐藏ax2的时间刻度（因为共享x轴）
plt.setp(ax2.get_xticklabels(), visible=False)
# 设置时间范围
ax1.set_xlim([mdates.date2num(start_dt), mdates.date2num(end_dt)])
# 自动调整布局
plt.tight_layout()
# 保存图片
start_date = start_dt.strftime("%Y%m%d")
end_date = end_dt.strftime("%Y%m%d")
output_filename = os.path.join(output_dir, f'clouds_{start_date}_{end_date}_with_dBZ.png')
plt.savefig(output_filename, dpi=150, bbox_inches='tight')
plt.close('all')
print(f"云层与反射率因子图已保存至: {output_filename}")
# 使用示例
if __name__ == "__main__":
# 示例参数 - 根据实际情况修改
folder = r"E:\HYData\腾冲云高\20250513\202505\20250513"
start = "20250513000000"
end = "20250513235959" # 改为全天时间
output = "./cloud_plots"
plot_cloud_series(folder, start, end, output)同时，在上述代码基础上，我有一个.csv文件，有3列，time列是日期，类似于2025/5/13 0:00:11，MaxRange/m列是最大距离，单位是米，VisiAvg/m是能见度，单位是米，请你读出里面的数据，并同步画在云层图中，也就是和云底数据一起画在第一张子图中，并标注出来自激光雷达，最终你一共画出两张子图。横轴24:00也需要显示，代码中只显示到23:00,请你修改代码。
扫一扫
3605
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
