python Chapter 6:雨滴速度-直径-数浓度分布图

一、数组数浓度计算：

1.时间戳＋32×32数组数浓度

2.整个降水过程32×32数浓度均值

#!/usr/bin/env python
# -*- coding:utf-8 -*-
"""
@author: Suyue
@file: calculate_nt_only.py
@time: 2026/09/07
@desc: 只计算32×32矩阵中每个直径-速度对应的Nt
       输出：
         1. Nt_particles.txt —— 时间戳 + 32×32 Nt数组（每个时间戳）
         2. Nt_mean.txt      —— 32×32 均值数组（所有时间戳逐元素平均）
       数字格式：0 输出 0，非 0 保留2位小数（不用科学计数法）
"""
import numpy as np
import re
import os


def read_cleaned_data(cleaned_file):
    """
    读取清洗后的数据文件
    格式：时间戳 + 32行数据 + 空行
    返回：时间戳列表 和 数据矩阵列表
    矩阵每一行是直径通道（32行），每一列是速度通道（32列）
    """
    timestamps = []
    data_matrices = []
    current_data = []
    current_timestamp = None

    timestamp_pattern = re.compile(r'^\d{4}-\d{2}-\d{2} \d{2}:\d{2}:\d{2}$')

    print(f"正在读取文件: {cleaned_file}")

    with open(cleaned_file, 'r') as f:
        for line in f:
            line = line.strip()

            if timestamp_pattern.match(line):
                if current_data and len(current_data) == 32 and current_timestamp:
                    data_matrices.append(np.array(current_data, dtype=float))
                    timestamps.append(current_timestamp)

                current_timestamp = line
                current_data = []

            elif line and any(c.isdigit() for c in line):
                try:
                    row_data = list(map(float, line.split()))
                    if len(row_data) == 32:
                        current_data.append(row_data)
                except:
                    pass

            elif not line and current_data:
                if len(current_data) == 32 and current_timestamp:
                    data_matrices.append(np.array(current_data, dtype=float))
                    timestamps.append(current_timestamp)
                    current_timestamp = None
                    current_data = []

    if current_data and len(current_data) == 32 and current_timestamp:
        data_matrices.append(np.array(current_data, dtype=float))
        timestamps.append(current_timestamp)

    print(f"读取完成: 共{len(timestamps)}个时间戳，每个矩阵形状: {data_matrices[0].shape if data_matrices else '无数据'}")

    if len(timestamps) != len(data_matrices):
        print(f"警告: 时间戳数量({len(timestamps)})与数据矩阵数量({len(data_matrices)})不匹配")
        min_len = min(len(timestamps), len(data_matrices))
        timestamps = timestamps[:min_len]
        data_matrices = data_matrices[:min_len]

    return timestamps, data_matrices


def calculate_nt_matrix(data_32x32, velocities_ms):
    """
    计算32×32矩阵中每个直径-速度对应的Nt
    Nt_ij = n_ij / (A * delta_t * V_j)
    """
    A_m2 = 0.0054   # 采样面积 m²
    delta_t = 60    # 采样时间 s

    nt_matrix = np.zeros((32, 32), dtype=float)

    for i in range(32):          # 行 = 直径通道
        for j in range(32):      # 列 = 速度通道
            V_j = velocities_ms[j]
            if V_j <= 0:
                continue
            n_ij = data_32x32[i, j]
            nt_matrix[i, j] = n_ij / (A_m2 * delta_t * V_j)

    return nt_matrix


def format_number(x):
    """
    数字格式化：
    - 0 直接输出 0
    - 非 0 保留2位小数（不用科学计数法）
    """
    if x == 0:
        return "0"
    return f"{x:.2f}"


def main():
    input_file = 'D:/lianxi/50934-20200813083300-20200813153959-0-clean.txt'
    output_particles = 'D:/lianxi/Nt_particles.txt'
    output_mean = 'D:/lianxi/Nt_mean.txt'

    velocities_ms = np.array([
        0.05, 0.15, 0.25, 0.35, 0.45, 0.55, 0.65, 0.75, 0.85, 0.95,
        1.1, 1.3, 1.5, 1.7, 1.9, 2.2, 2.6, 3, 3.4, 3.8,
        4.4, 5.2, 6, 6.8, 7.6, 8.8, 10.4, 12, 13.6, 15.2, 17.6, 20.8
    ])

    if not os.path.exists(input_file):
        print(f"错误: 输入文件不存在: {input_file}")
        return

    print("=" * 60)
    print("32×32 Nt 计算程序")
    print("=" * 60)

    timestamps, data_matrices = read_cleaned_data(input_file)

    if not timestamps:
        print("错误: 没有读取到有效数据")
        return

    print(f"\n开始计算Nt，共 {len(timestamps)} 个时间戳...")

    nt_matrices = []   # 保存每个时间戳的Nt矩阵

    for idx, (timestamp, data_matrix) in enumerate(zip(timestamps, data_matrices)):
        if idx % 200 == 0:
            print(f"  处理到第 {idx + 1}/{len(timestamps)} 个时间戳")

        nt_matrix = calculate_nt_matrix(data_matrix, velocities_ms)
        nt_matrices.append(nt_matrix)

    print("计算完成")

    # ===== 输出文件1：时间戳 + 32×32 Nt数组 =====
    print(f"\n正在保存: {output_particles}")
    with open(output_particles, 'w', encoding='utf-8') as f:
        for timestamp, nt_matrix in zip(timestamps, nt_matrices):
            f.write(timestamp + '\n')
            for i in range(32):
                row_str = ' '.join(format_number(nt_matrix[i, j]) for j in range(32))
                f.write(row_str + '\n')
            f.write('\n')   # 时间戳块之间空一行（跟原始数据格式一致）
    print(f"已保存: {output_particles}")

    # ===== 输出文件2：32×32均值数组 =====
    nt_mean_matrix = np.mean(np.array(nt_matrices), axis=0)   # 逐元素平均

    print(f"\n正在保存: {output_mean}")
    with open(output_mean, 'w', encoding='utf-8') as f:
        for i in range(32):
            row_str = ' '.join(format_number(nt_mean_matrix[i, j]) for j in range(32))
            f.write(row_str + '\n')
    print(f"已保存: {output_mean}")

    # ===== 预览 =====
    print("\n" + "=" * 60)
    print("Nt_mean.txt 预览 (前5行):")
    print("=" * 60)
    with open(output_mean, 'r', encoding='utf-8') as f:
        for i, line in enumerate(f):
            if i < 5:
                print(line.strip())
            else:
                break

    print("\n" + "=" * 60)
    print(f"时间点数量: {len(nt_matrices)}")
    print(f"均值矩阵形状: {nt_mean_matrix.shape}")
    print(f"均值矩阵范围: {nt_mean_matrix.min():.2f} - {nt_mean_matrix.max():.2f}")
    print("=" * 60)


if __name__ == "__main__":
    main()

 

 二、速度-直径-数浓度分布图（雨滴）

1.上述32×32数组数浓度均值

2.32×32数组数浓度均值取对数

3.画图横轴为直径、纵轴速度、值数浓度对数

所以代码需要先取对数，再画图

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import BoundaryNorm, ListedColormap

# =========================
# 1. 读取你的 32x32 数浓度矩阵
# =========================
file_path = "D:/lianxi/Nt_mean.txt"

data = np.loadtxt(file_path, delimiter=None)

print("数据形状:", data.shape)
assert data.shape == (32, 32), "数据必须是 32x32"

# =========================
# 2. 先取 log10：>0 的取对数，=0 的设为 NaN（显示为白色背景）
# =========================
data_log = np.full_like(data, np.nan, dtype=float)
mask = data > 0
data_log[mask] = np.log10(data[mask])

np.set_printoptions(precision=4, suppress=True, linewidth=200)
print("\nlog10 后的数据 data_log =")
print(data_log)

log_txt_path = "D:/lianxi/Nt_mean_log10.txt"
np.savetxt(log_txt_path, data_log, fmt="%.6f")
print(f"\n已保存 log10 结果到: {log_txt_path}")

# =========================
# 3. 定义通道坐标
# =========================
diameters_mm = np.array([
    0.062, 0.187, 0.312, 0.437, 0.562, 0.687, 0.812, 0.937,
    1.064, 1.193, 1.388, 1.651, 1.918, 2.188, 2.463, 2.882,
    3.457, 4.050, 4.665, 5.303, 6.194, 7.321, 8.447, 9.573, 10.699,
    12.389, 14.641, 16.894, 19.146, 21.399, 24.214, 27.593
])

velocities_ms = np.array([
    0.05, 0.15, 0.25, 0.35, 0.45, 0.55, 0.65, 0.75, 0.85, 0.95,
    1.1, 1.3, 1.5, 1.7, 1.9, 2.2, 2.6, 3, 3.4, 3.8,
    4.4, 5.2, 6, 6.8, 7.6, 8.8, 10.4, 12, 13.6, 15.2, 17.6, 20.8
])

# =========================
# 4. 自定义颜色分级
# =========================
upper = max(12, np.nanmax(data_log) + 1)
bounds = [0, 1, 2, 3, 4, 5, 6, upper]

colors = [
    'darkblue',
    'skyblue',
    'green',
    'yellow',
    'orange',
    'red',
    'purple'
]

cmap = ListedColormap(colors)
cmap.set_bad(color='white')   # NaN 显示为白色
norm = BoundaryNorm(bounds, cmap.N)

# =========================
# 5. 画二维分布图
# =========================
plt.rcParams['font.family'] = 'DejaVu Sans'
plt.rcParams['axes.unicode_minus'] = False

fig, ax = plt.subplots(figsize=(9, 7))
fig.patch.set_facecolor('white')
ax.set_facecolor('white')

mesh = ax.pcolormesh(
    diameters_mm,
    velocities_ms,
    data_log.T,
    shading='auto',
    cmap=cmap,
    norm=norm
)

# =========================
# 6. 色标
# =========================
cbar = plt.colorbar(mesh, ax=ax, ticks=[0, 1, 2, 3, 4, 5, 6, 7])
cbar.set_label('log10(Number concentration)')
cbar.ax.set_yticklabels([
    '0-1',
    '1-2',
    '2-3',
    '3-4',
    '4-5',
    '5-6',
    '6+',
    ''
])

# =========================
# 7. 叠加经典落速曲线
# =========================
D_fit = np.linspace(diameters_mm.min(), diameters_mm.max(), 500)
V_fit = 9.65 - 10.3 * np.exp(-0.6 * D_fit)

ax.plot(D_fit, V_fit, color='black', linewidth=2.5,
        label='Vt = 9.65 - 10.3exp(-0.6D)')

# =========================
# 8. 坐标轴与标签
# =========================
ax.set_xlabel('Diameter (mm)')
ax.set_ylabel('Fall speed (m/s)')
ax.set_title('Raindrop size–velocity distribution (log10 concentration)')

ax.set_xlim(diameters_mm.min(), diameters_mm.max())
ax.set_ylim(velocities_ms.min(), velocities_ms.max())

ax.legend(loc='upper left')

plt.tight_layout()

# =========================
# 9. 保存图片
# =========================
png_path = "D:/lianxi/raindrop_distribution_log10.png"
plt.savefig(png_path, dpi=300, facecolor='white')
print(f"已保存图片到: {png_path}")

plt.show()