代码实战 | JMA 台风最佳路径可视化教程

66666 2501 022 0001 2501 0 6 WUTIP 20250912

25061300 002 4 183 1084 980 055 90030 0030 90180 0180

import re
from pathlib import Path

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.lines import Line2D
import cartopy.crs as ccrs
import cartopy.feature as cfeature

DATA_PATH = Path('bst_all.txt')
OUTPUT_FIG = Path('jma_2025_typhoon_tracks.png')

HEADER_PATTERN = re.compile(
    r'^66666\s+'
    r'(?P<international_number>\d{4})\s+'
    r'(?P<num_lines>\d{3})\s+'
    r'(?P<tc_id>\d{4})\s+'
    r'(?P<jma_id>\d{4})\s+'
    r'(?P<header_flag>\d)\s+'
    r'(?P<time_interval>\d)\s*'
    r'(?P<name>.*?)\s+'
    r'(?P<revision_date>\d{8})\s*$'
)

GRADE_MAP = {
    2: 'Tropical Depression',
    3: 'Tropical Storm',
    4: 'Severe Tropical Storm',
    5: 'Typhoon',
    6: 'Extratropical Cyclone',
    7: 'Just entering the responsible area of RSMC Tokyo-Typhoon Center',
    9: 'Tropical Cyclone of TS intensity or higher',
}

WIND_DIR_MAP = {
    0: 'None / Radius=0',
    1: 'NE quadrant',
    2: 'E semicircle',
    3: 'SE quadrant',
    4: 'S semicircle',
    5: 'SW quadrant',
    6: 'W semicircle',
    7: 'NW quadrant',
    8: 'N semicircle',
    9: 'Full circle (concentric)',
}


def parse_jma_best_track(file_path: Path, year_prefix: str = '25') -> pd.DataFrame:
    rows = []
    current_storm = None

    for raw_line in file_path.read_text(encoding='utf-8', errors='ignore').splitlines():
        if raw_line.startswith('66666'):
            match = HEADER_PATTERN.match(raw_line)
            current_storm = None
            if match:
                meta = match.groupdict()
                if meta['international_number'].startswith(year_prefix):
                    current_storm = meta
            continue

        if current_storm isNoneornot raw_line.strip():
            continue

        parts = raw_line.split()
        if len(parts) < 6:
            continue

        obs_time = pd.to_datetime(parts[0], format='%y%m%d%H')
        grade_code = int(parts[2])
        latitude = int(parts[3]) / 10.0
        longitude = int(parts[4]) / 10.0
        pressure_hpa = int(parts[5])
        max_wind_kt = int(parts[6]) if len(parts) > 6else np.nan

        # Wind radii parsing (tokens 7–10 present when len(parts) >= 11)
        if len(parts) >= 11:
            r50_dir = int(parts[7][0])
            r50_long = int(parts[7][1:])
            r50_short = int(parts[8])
            r30_dir = int(parts[9][0])
            r30_long = int(parts[9][1:])
            r30_short = int(parts[10])
        else:
            r50_dir = r50_long = r50_short = np.nan
            r30_dir = r30_long = r30_short = np.nan

        rows.append(
            {
                'international_number': current_storm['international_number'],
                'name': current_storm['name'].strip(),
                'time': obs_time,
                'grade_code': grade_code,
                'grade_name': GRADE_MAP.get(grade_code, f'Unknown({grade_code})'),
                'latitude': latitude,
                'longitude': longitude,
                'pressure_hpa': pressure_hpa,
                'max_wind_kt': np.nan if max_wind_kt == 0else max_wind_kt,
                'wind_radius_50kt_dir': r50_dir if r50_dir != 0else np.nan,
                'wind_radius_50kt_long_nm': r50_long if r50_long != 0else np.nan,
                'wind_radius_50kt_short_nm': r50_short if r50_short != 0else np.nan,
                'wind_radius_30kt_dir': r30_dir if r30_dir != 0else np.nan,
                'wind_radius_30kt_long_nm': r30_long if r30_long != 0else np.nan,
                'wind_radius_30kt_short_nm': r30_short if r30_short != 0else np.nan,
            }
        )

    df = pd.DataFrame(rows)
    if df.empty:
        raise ValueError('没有解析到 2025 年台风记录，请检查数据文件格式。')

    return df.sort_values(['international_number', 'time']).reset_index(drop=True)


def summarize_typhoons(df: pd.DataFrame) -> pd.DataFrame:
    summary = (
        df.groupby(['international_number', 'name'], as_index=False)
        .agg(
            start_time=('time', 'min'),
            end_time=('time', 'max'),
            track_points=('time', 'size'),
            min_pressure_hpa=('pressure_hpa', 'min'),
            max_wind_kt=('max_wind_kt', 'max'),
            max_latitude=('latitude', 'max'),
            min_latitude=('latitude', 'min'),
            max_longitude=('longitude', 'max'),
            min_longitude=('longitude', 'min'),
        )
        .sort_values('international_number')
        .reset_index(drop=True)
    )
    return summary

df_2025 = parse_jma_best_track(DATA_PATH, year_prefix='25')
summary_2025 = summarize_typhoons(df_2025)

print(f'共解析出台风 {df_2025["international_number"].nunique()} 个')
print(f'共解析出路径点 {len(df_2025)} 条')

共解析出台风 27 个
共解析出路径点 882 条

summary_2025

def plot_all_typhoon_tracks(df: pd.DataFrame, output_path: Path):
    storm_ids = df['international_number'].drop_duplicates().tolist()
    colors = plt.cm.nipy_spectral(np.linspace(0.05, 0.95, len(storm_ids)))

    fig = plt.figure(figsize=(15, 11))
    ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=180))
    ax.set_extent([100, 260, 0, 60], crs=ccrs.PlateCarree())

    ax.add_feature(cfeature.OCEAN.with_scale('110m'), facecolor='#dceefb', zorder=0)
    ax.add_feature(cfeature.LAND.with_scale('110m'), facecolor='#f7f3e8', edgecolor='none', zorder=0)
    ax.add_feature(cfeature.COASTLINE.with_scale('110m'), linewidth=0.7)
    ax.add_feature(cfeature.BORDERS.with_scale('110m'), linewidth=0.4, linestyle=':')

    gl = ax.gridlines(
        crs=ccrs.PlateCarree(),
        draw_labels=True,
        linewidth=0.5,
        color='gray',
        alpha=0.5,
        linestyle='--'
    )
    gl.top_labels = False
    gl.right_labels = False

    for color, (storm_id, track) in zip(colors, df.groupby('international_number')):
        track = track.sort_values('time')
        storm_name = track['name'].iloc[0]
        label = f'{storm_id} {storm_name}'

        ax.plot(
            track['longitude'],
            track['latitude'],
            color=color,
            linewidth=1.8,
            transform=ccrs.PlateCarree(),
            label=label,
        )
        ax.scatter(
            track['longitude'].iloc[0],
            track['latitude'].iloc[0],
            color=color,
            s=20,
            marker='o',
            transform=ccrs.PlateCarree(),
            zorder=3,
        )
        ax.scatter(
            track['longitude'].iloc[-1],
            track['latitude'].iloc[-1],
            color=color,
            s=24,
            marker='s',
            transform=ccrs.PlateCarree(),
            zorder=3,
        )

    ax.set_title('JMA Best Track of All 2025 Typhoons', fontsize=18, pad=16)
    legend = ax.legend(
        loc='upper center',
        bbox_to_anchor=(0.5, -0.08),
        ncol=3,
        fontsize=8,
        frameon=True,
        title='Storm ID / Name'
    )
    legend.get_title().set_fontsize(9)

    fig.tight_layout()

    if output_path isnotNone:
        fig.savefig(output_path, dpi=180, bbox_inches='tight')
        print(f'图片已保存到：{output_path.resolve()}')

    return fig, ax


fig, ax = plot_all_typhoon_tracks(df_2025, OUTPUT_FIG)
plt.show()

图片已保存到：/home/mw/project/jma_2025_typhoon_tracks.png

GRADE_STYLE = {
    2: {'label': 'TD', 'color': '#9ecae1'},
    3: {'label': 'TS', 'color': '#3182bd'},
    4: {'label': 'STS', 'color': '#31a354'},
    5: {'label': 'TY', 'color': '#de2d26'},
    6: {'label': 'Extratropical', 'color': '#756bb1'},
    7: {'label': 'Entering area', 'color': '#fd8d3c'},
    9: {'label': 'TS or higher', 'color': '#fdae6b'},
}


def plot_typhoon_tracks_by_intensity(df: pd.DataFrame, output_path: Path):
    fig = plt.figure(figsize=(15, 11))
    ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=180))
    ax.set_extent([100, 260, 0, 60], crs=ccrs.PlateCarree())

    ax.add_feature(cfeature.OCEAN.with_scale('110m'), facecolor='#dceefb', zorder=0)
    ax.add_feature(cfeature.LAND.with_scale('110m'), facecolor='#f7f3e8', edgecolor='none', zorder=0)
    ax.add_feature(cfeature.COASTLINE.with_scale('110m'), linewidth=0.7)
    ax.add_feature(cfeature.BORDERS.with_scale('110m'), linewidth=0.4, linestyle=':')

    gl = ax.gridlines(
        crs=ccrs.PlateCarree(),
        draw_labels=True,
        linewidth=0.5,
        color='gray',
        alpha=0.5,
        linestyle='--'
    )
    gl.top_labels = False
    gl.right_labels = False

    for storm_id, track in df.groupby('international_number'):
        track = track.sort_values('time').reset_index(drop=True)

        for i in range(len(track) - 1):
            row = track.iloc[i]
            next_row = track.iloc[i + 1]
            style = GRADE_STYLE.get(int(row['grade_code']), {'label': 'Unknown', 'color': '#636363'})
            ax.plot(
                [row['longitude'], next_row['longitude']],
                [row['latitude'], next_row['latitude']],
                color=style['color'],
                linewidth=1.8,
                alpha=0.95,
                transform=ccrs.PlateCarree(),
                zorder=2,
            )

        point_colors = [GRADE_STYLE.get(int(code), {'label': 'Unknown', 'color': '#636363'})['color'] for code in track['grade_code']]
        ax.scatter(
            track['longitude'],
            track['latitude'],
            c=point_colors,
            s=12,
            edgecolors='none',
            transform=ccrs.PlateCarree(),
            zorder=3,
        )

    observed_codes = sorted(df['grade_code'].dropna().astype(int).unique())
    legend_handles = [
        Line2D([0], [0], color=GRADE_STYLE[code]['color'], lw=2.5, label=f"{code} {GRADE_STYLE[code]['label']}")
        for code in observed_codes if code in GRADE_STYLE
    ]

    ax.set_title('JMA Best Track of All 2025 Typhoons (Colored by Intensity)', fontsize=18, pad=16)
    legend = ax.legend(
        handles=legend_handles,
        loc='upper center',
        bbox_to_anchor=(0.5, -0.08),
        ncol=3,
        fontsize=9,
        frameon=True,
        title='JMA grade code / intensity'
    )
    legend.get_title().set_fontsize(10)

    fig.tight_layout()

    if output_path isnotNone:
        fig.savefig(output_path, dpi=180, bbox_inches='tight')
        print(f'图片已保存到：{output_path.resolve()}')

    return fig, ax


OUTPUT_FIG_INTENSITY = Path('jma_2025_typhoon_tracks_by_intensity.png')
fig, ax = plot_typhoon_tracks_by_intensity(df_2025, OUTPUT_FIG_INTENSITY)
plt.show()

图片已保存到：/home/mw/project/jma_2025_typhoon_tracks_by_intensity.png

from matplotlib.patches import Ellipse
import math


def nm_to_deg(radius_nm: float, lat: float, axis: str = 'ns') -> float:
    """将海里半径转换为地理度数（PlateCarree 数据坐标）。"""
    if axis == 'ns':  # 南北方向
        return radius_nm / 60.0
    else:  # 东西方向，需除以 cos(lat)
        return radius_nm / (60.0 * math.cos(math.radians(lat)))


def plot_wind_radii_for_storm(df: pd.DataFrame, storm_name: str, output_path: Path):
    storm = df[df['name'] == storm_name].sort_values('time').reset_index(drop=True)
    if storm.empty:
        raise ValueError(f'未找到台风 {storm_name}')

    fig = plt.figure(figsize=(14, 12))
    ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=180))

    # 自动设置范围：路径四周留 3° 边距
    lon_min, lon_max = storm['longitude'].min() - 3, storm['longitude'].max() + 3
    lat_min, lat_max = storm['latitude'].min() - 3, storm['latitude'].max() + 3
    ax.set_extent([lon_min, lon_max, lat_min, lat_max], crs=ccrs.PlateCarree())

    ax.add_feature(cfeature.OCEAN.with_scale('110m'), facecolor='#dceefb', zorder=0)
    ax.add_feature(cfeature.LAND.with_scale('110m'), facecolor='#f7f3e8', edgecolor='none', zorder=0)
    ax.add_feature(cfeature.COASTLINE.with_scale('110m'), linewidth=0.7)
    ax.add_feature(cfeature.BORDERS.with_scale('110m'), linewidth=0.4, linestyle=':')

    gl = ax.gridlines(
        crs=ccrs.PlateCarree(),
        draw_labels=True,
        linewidth=0.5,
        color='gray',
        alpha=0.5,
        linestyle='--'
    )
    gl.top_labels = False
    gl.right_labels = False

    # 绘制台风路径
    ax.plot(
        storm['longitude'],
        storm['latitude'],
        color='#2c2c2c',
        linewidth=1.5,
        marker='o',
        markersize=3,
        transform=ccrs.PlateCarree(),
        zorder=4,
        label='Track'
    )

    # 在每个时次绘制风圈
    for _, row in storm.iterrows():
        lat = row['latitude']
        lon = row['longitude']

        # 30 kt 风圈（蓝色）
        r30 = row['wind_radius_30kt_long_nm']
        if pd.notna(r30) and r30 > 0:
            # 用 Ellipse 近似圆，修正 cos(lat) 压缩
            width = 2 * nm_to_deg(r30, lat, axis='ew')
            height = 2 * nm_to_deg(r30, lat, axis='ns')
            ellipse = Ellipse(
                xy=(lon, lat),
                width=width,
                height=height,
                angle=0,
                facecolor='#3182bd',
                edgecolor='#2166ac',
                alpha=0.18,
                linewidth=0.8,
                transform=ccrs.PlateCarree(),
                zorder=2,
            )
            ax.add_patch(ellipse)

        # 50 kt 风圈（红色）
        r50 = row['wind_radius_50kt_long_nm']
        if pd.notna(r50) and r50 > 0:
            width = 2 * nm_to_deg(r50, lat, axis='ew')
            height = 2 * nm_to_deg(r50, lat, axis='ns')
            ellipse = Ellipse(
                xy=(lon, lat),
                width=width,
                height=height,
                angle=0,
                facecolor='#de2d26',
                edgecolor='#a50f15',
                alpha=0.22,
                linewidth=0.8,
                transform=ccrs.PlateCarree(),
                zorder=3,
            )
            ax.add_patch(ellipse)

    # 标注台风名称（路径中点附近）
    mid_idx = len(storm) // 2
    ax.text(
        storm.iloc[mid_idx]['longitude'] + 0.8,
        storm.iloc[mid_idx]['latitude'] + 0.8,
        storm_name,
        fontsize=14,
        fontweight='bold',
        color='#2c2c2c',
        transform=ccrs.PlateCarree(),
        zorder=5,
    )

    ax.set_title(
        f'JMA Best Track of {storm_name} (2025) with 30 kt / 50 kt Wind Radii',
        fontsize=16,
        pad=14
    )

    # 自定义图例
    from matplotlib.patches import Patch
    legend_handles = [
        Line2D([0], [0], color='#2c2c2c', marker='o', markersize=5, linewidth=1.5, label='Best track'),
        Patch(facecolor='#3182bd', edgecolor='#2166ac', alpha=0.4, label='30 kt wind radius'),
        Patch(facecolor='#de2d26', edgecolor='#a50f15', alpha=0.5, label='50 kt wind radius'),
    ]
    ax.legend(
        handles=legend_handles,
        loc='upper center',
        bbox_to_anchor=(0.5, -0.06),
        ncol=3,
        fontsize=10,
        frameon=True,
    )

    fig.tight_layout()

    return fig, ax


OUTPUT_FIG_RADII = Path('jma_2025_krosa_wind_radii.png')
fig, ax = plot_wind_radii_for_storm(df_2025, 'KROSA', OUTPUT_FIG_RADII)
plt.show()