OneScience-Group/FourCastNet

1import numpy as np2import matplotlib.pyplot as plt3import os4import sys5import glob6import h5py7from datetime import datetime8from tqdm import tqdm9from onescience.utils.fcn.YParams import YParams10from matplotlib import rcParams11 12# rcParams['font.family'] = 'serif'13# rcParams['font.serif'] = ['DejaVu Serif']14rcParams['mathtext.fontset'] = 'stix'15rcParams['axes.linewidth'] = 0.916rcParams['xtick.major.width'] = 0.917rcParams['ytick.major.width'] = 0.918 19 20def get_metadata(data_dir, channels):21    """从新版 h5 attrs 中读取变量列表和 time_step"""22    h5_files = sorted(glob.glob(os.path.join(data_dir, "data", "*.h5")))23    with h5py.File(h5_files[0], "r") as f:24        ds = f["fields"]25        all_variables = [v.decode() if isinstance(v, bytes) else v for v in ds.attrs["variables"]]26        time_step = int(ds.attrs["time_step"])27 28    channel_indices = [all_variables.index(v) for v in channels]29 30    total_files = [f for f in os.listdir('./result/output/') if f.endswith('.npy')]31    total_files.sort()32    return total_files, channel_indices, time_step33 34 35def filename_to_index(filename, time_step):36    """将 YYYYMMDDHH 格式的文件名转换为年度 h5 文件中的时间步索引"""37    dt = datetime.strptime(filename, "%Y%m%d%H")38    year_start = datetime(dt.year, 1, 1)39    hours = (dt - year_start).total_seconds() / 360040    return int(hours / time_step)41 42 43def get_result(total_files, channel_indices, time_step, data_dir, clim_mean):44    channel_rmse = np.zeros(len(channel_indices))45    channel_acc = np.zeros(len(channel_indices))46    clim_mean = clim_mean[0, :, :, :]47    if not os.path.exists('./result/rmse.npy') or not os.path.exists('result/acc.npy'):48        numerator = np.zeros(len(channel_indices))49        pred_sq_sum = np.zeros(len(channel_indices))50        label_sq_sum = np.zeros(len(channel_indices))51        for file in tqdm(total_files, unit="files"):52            fname = file[:-4]  # 去掉 .npy53            year = fname[:4]54            t_idx = filename_to_index(fname, time_step)55            with h5py.File(os.path.join(data_dir, 'data', f'{year}.h5'), "r") as f:56                label = f["fields"][t_idx]  # [C, H, W]57                label = label[channel_indices]58                label = label[:, :-1, :]59            pred = np.load(f'result/output/{file}').squeeze()60 61            label_anom = label - clim_mean62            pred_anom = pred - clim_mean63            # 累加64            numerator += np.sum(pred_anom * label_anom, axis=(1, 2))65            pred_sq_sum += np.sum(pred_anom ** 2, axis=(1, 2))66            label_sq_sum += np.sum(label_anom ** 2, axis=(1, 2))67 68            channel_rmse += np.sqrt(np.mean((label - pred) ** 2, axis=(1, 2)))69        channel_rmse /= len(total_files)70        channel_acc = numerator / (np.sqrt(pred_sq_sum * label_sq_sum) + 1e-8)71        np.save('./result/acc.npy', channel_acc)72        np.save('./result/rmse.npy', channel_rmse)73 74 75def show_result():76    channel_rmse = np.load('./result/rmse.npy')77    channel_acc = np.load('./result/acc.npy')78 79    channels = [cfg_data.dataset.channels[i] for i in range(len(channel_indices))]80    w = 24  # 最长 channel 名宽度81 82    # 表头83    print(f"┌{'─' * (w + 2)}┬{'─' * 14}┬{'─' * 14}┐")84    print(f"│ {'Channel':<{w}} │ {'RMSE':>12} │ {'ACC':>12} │")85    print(f"├{'─' * (w + 2)}┼{'─' * 14}┼{'─' * 14}┤")86    # 数据行87    for i, ch in enumerate(channels):88        print(f"│ {ch:<{w}} │ {channel_rmse[i]:>12.4f} | {channel_acc[i]:>12.4f} |")89    print(f"├{'─' * (w + 2)}┼{'─' * 14}┼{'─' * 14}┤")90    print(f"│ {'Average':<{w}} │ {np.mean(channel_rmse):>12.4f} │ {np.mean(channel_acc):>12.4f} │")91    print(f"└{'─' * (w + 2)}┴{'─' * 14}┴{'─' * 14}┘")92 93 94def plot(label, pred, var, filename):95    fig, axes = plt.subplots(1, 3, figsize=(15, 4))96 97    xtick_labels = ['180°W', '90°W', '0°', '90°E', '180°E']98    ytick_labels = ['90°S', '45°S', '0°', '45°N', '90°N']99    xticks = np.linspace(0, label.shape[-1] - 1, 5)100    yticks = np.linspace(0, label.shape[-2] - 1, 5)101 102    vmin = min(label.min(), pred.min())103    vmax = max(label.max(), pred.max())104 105    diff = label - pred106    rmse = np.sqrt(np.mean(diff ** 2))107    diff_abs_max = np.abs(diff).max()108 109    plot_configs = [110        {'data': label, 'title': 'Truth', 'cmap': 'viridis', 'vmin': vmin, 'vmax': vmax},111        {'data': pred,  'title': 'Prediction', 'cmap': 'viridis', 'vmin': vmin, 'vmax': vmax},112        {'data': diff,  'title': f'Difference (RMSE={rmse:.2f})', 'cmap': 'RdBu_r', 'vmin': -diff_abs_max, 'vmax': diff_abs_max},113    ]114 115    for ax, cfg in zip(axes, plot_configs):116        im = ax.imshow(cfg['data'], cmap=cfg['cmap'], vmin=cfg['vmin'], vmax=cfg['vmax'])117        ax.set_title(cfg['title'], fontsize=12, pad=4)118        ax.set_xlabel('Longitude')119        ax.set_ylabel('Latitude')120        ax.set_xticks(xticks)121        ax.set_xticklabels(xtick_labels)122        ax.set_yticks(yticks)123        ax.set_yticklabels(ytick_labels)124        plt.colorbar(im, ax=ax, orientation='horizontal')125 126    fig.suptitle(var, fontsize=14, fontweight='bold', y=0.98)127    plt.savefig(filename, dpi=300, bbox_inches='tight')128    plt.close()129 130 131def plot_loss(train_loss, valid_loss):132    mask = ~(np.isnan(train_loss) | np.isnan(valid_loss))133    train_loss = train_loss[mask]134    valid_loss = valid_loss[mask]135 136    fig, ax = plt.subplots(figsize=(5, 3.5))137    colors = {'train': '#2563EB', 'valid': '#EA580C'}138    epochs = np.arange(1, len(train_loss) + 1)139 140    ax.plot(epochs, train_loss, color=colors['train'], linewidth=1.5, label='Train')141    ax.plot(epochs, valid_loss, color=colors['valid'], linewidth=1.5, label='Valid', linestyle='--')142    min_idx = np.argmin(valid_loss)143    ax.scatter(epochs[min_idx], valid_loss[min_idx],144               color=colors['valid'], s=40, zorder=5, edgecolors='white')145    ax.annotate(f'Best: {valid_loss[min_idx]:.3f}',146                xy=(epochs[min_idx], valid_loss[min_idx]),147                xytext=(10, 10), textcoords='offset points', fontsize=8, color=colors['valid'],148                arrowprops=dict(arrowstyle='-', color=colors['valid'], lw=0.5))149 150    ax.set(xlabel='Epoch', ylabel='Loss', xlim=(0, len(train_loss) + 1))151    ax.legend(frameon=False, loc='upper right')152    ax.grid(True, linestyle='--', alpha=0.3)153    ax.spines[['top', 'right']].set_visible(False)154 155    plt.tight_layout()156    plt.savefig('./result/loss.png', dpi=300, bbox_inches='tight')157    plt.close()158 159 160if __name__ == "__main__":161    current_path = os.getcwd()162    sys.path.append(current_path)163    config_file_path = os.path.join(current_path, 'conf/config.yaml')164    cfg = YParams(config_file_path, 'model')165    cfg_data = YParams(config_file_path, "datapipe")166 167    train_loss = np.load('./data/checkpoints/trloss.npy')168    valid_loss = np.load('./data/checkpoints/valoss.npy')169    plot_loss(train_loss, valid_loss)170 171    data_dir = cfg_data.dataset.data_dir172    total_files, channel_indices, time_step = get_metadata(data_dir, cfg_data.dataset.channels)173 174    # Load data & Compute RMSE/ACC per channel175    h5_files = sorted(glob.glob(os.path.join(data_dir, "data", "*.h5")))176    with h5py.File(h5_files[0], "r") as f:177        mu = f["global_means"][:]178    clim_mean = mu[:, channel_indices, :, :]179    get_result(total_files, channel_indices, time_step, data_dir, clim_mean)180    show_result()181 182    ##### 默认绘制 test_time 第一年的第一个时间步，用户可自行指定日期和变量 #####183    test_year = cfg_data.dataset.test_time[0]184    eg_files = [f'{test_year}010206']185    channel_index = [cfg_data.dataset.channels.index(v) for v in ['2m_temperature', 'geopotential_500', 'temperature_500']]186 187    selected_var = [cfg_data.dataset.channels[int(i)] for i in channel_index]188    print(f"seleted date: {eg_files}")189    print(f"selected channels: {selected_var}")190    for file in eg_files:191        year = file[:4]192        t_idx = filename_to_index(file, time_step)193        with h5py.File(os.path.join(data_dir, 'data', f'{year}.h5'), "r") as f:194            label = f["fields"][t_idx]  # [C, H, W]195            label = label[channel_indices]196            label = label[:, :-1, :]197        pred = np.load(f'result/output/{file}.npy').squeeze()198        for i in range(len(selected_var)):199            filename = f'./result/{file}_{selected_var[i]}.png'200            plot(label[channel_index[i]], pred[channel_index[i]], selected_var[i], filename)201            print(f'✅plot {filename}')202