BAVI prediction with multiple sources

typhoon_bavi_gfs_report.py

#!/usr/bin/env python3
from __future__ import annotations

import importlib.util
import math
import os
import re
import sys
import warnings
from datetime import datetime, timedelta, timezone
from pathlib import Path
from typing import Iterable

import cfgrib
import numpy as np
import requests
from PIL import Image, ImageDraw

BASE_DIR = Path("/tmp/hermes_typhoon_bavi")
GFS_DIR = BASE_DIR / "gfs_cron"
GFS_DIR.mkdir(parents=True, exist_ok=True)
# Higher-resolution base map for the wide typhoon viewport. Set before loading
# typhoon_bavi_report.py because that module reads these values at import time.
os.environ.setdefault("TYPHOON_MAP_TILE_ZOOM", "7")
os.environ.setdefault("TYPHOON_MAP_OUTPUT_SCALE", "4")
os.environ.setdefault("TYPHOON_MAP_TILE_URL", "https://basemaps.cartocdn.com/rastertiles/voyager_nolabels/{z}/{x}/{y}.png")
os.environ.setdefault("TYPHOON_MAP_LABEL_LEVEL", "province")
warnings.filterwarnings("ignore", category=FutureWarning, module="cfgrib")
warnings.filterwarnings("ignore", category=FutureWarning, module="cfgrib.xarray_store")
REPORT_SCRIPT = Path(__file__).resolve().with_name("typhoon_bavi_report.py")
NOMADS_FILTER = "https://nomads.ncep.noaa.gov/cgi-bin/filter_gfs_0p25.pl"
LON_MIN, LON_MAX = 105, 135
LAT_MIN, LAT_MAX = 15, 40
FORECAST_STEP_HOURS = (0, 6, 12, 18, 24)
WIND_THRESHOLDS = {
    "radius7": 13.9,
    "radius10": 24.5,
    "radius12": 32.7,
}

spec = importlib.util.spec_from_file_location("typhoon_bavi_report", REPORT_SCRIPT)
if spec is None or spec.loader is None:
    raise RuntimeError(f"cannot load {REPORT_SCRIPT}")
typhoon = importlib.util.module_from_spec(spec)
spec.loader.exec_module(typhoon)


def bj_to_utc_naive(dt: datetime) -> datetime:
    return dt - timedelta(hours=8)


def utc_to_bj_string(dt_utc: datetime) -> str:
    return (dt_utc + timedelta(hours=8)).strftime("%Y-%m-%d %H:%M:%S")


def floor_to_gfs_cycle(now_utc: datetime) -> datetime:
    now_utc = now_utc.replace(minute=0, second=0, microsecond=0)
    return now_utc.replace(hour=(now_utc.hour // 6) * 6)


def candidate_cycles(limit: int = 12) -> Iterable[datetime]:
    cycle = floor_to_gfs_cycle(datetime.now(timezone.utc)).replace(tzinfo=None)
    for i in range(limit):
        yield cycle - timedelta(hours=6 * i)


def gfs_hour_available(hour: int) -> bool:
    if hour < 0 or hour > 384:
        return False
    if hour <= 120:
        return True
    return hour % 3 == 0


def rounded_available_hour(raw_hour: float) -> int:
    hour = int(round(raw_hour))
    if hour <= 120:
        return hour
    return int(round(hour / 3) * 3)


def gfs_url(cycle: datetime, hour: int) -> str:
    date = cycle.strftime("%Y%m%d")
    hh = cycle.strftime("%H")
    return (
        f"{NOMADS_FILTER}?file=gfs.t{hh}z.pgrb2.0p25.f{hour:03d}"
        "&lev_10_m_above_ground=on&lev_mean_sea_level=on"
        "&var_UGRD=on&var_VGRD=on&var_PRMSL=on&subregion="
        f"&leftlon={LON_MIN}&rightlon={LON_MAX}&toplat={LAT_MAX}&bottomlat={LAT_MIN}"
        f"&dir=%2Fgfs.{date}%2F{hh}%2Fatmos"
    )


def fetch_gfs_subset(cycle: datetime, hour: int) -> Path:
    date = cycle.strftime("%Y%m%d")
    hh = cycle.strftime("%H")
    path = GFS_DIR / f"gfs_{date}_{hh}_f{hour:03d}_prmsl_u10v10.grb2"
    response = requests.get(gfs_url(cycle, hour), timeout=180)
    response.raise_for_status()
    content = response.content
    if len(content) < 10_000 or content.lstrip().startswith(b"<"):
        raise RuntimeError(f"GFS subset unavailable for {date} {hh}Z f{hour:03d}")
    path.write_bytes(content)
    return path


def open_gfs_subset(path: Path):
    for idx_path in path.parent.glob(path.name + ".*.idx"):
        try:
            idx_path.unlink()
        except OSError:
            pass
    prmsl = uv10 = None
    for ds in cfgrib.open_datasets(str(path)):
        if "prmsl" in ds:
            prmsl = ds
        if "u10" in ds and "v10" in ds:
            uv10 = ds
    if prmsl is None or uv10 is None:
        raise RuntimeError(f"missing PRMSL/U10/V10 in {path}")
    return prmsl, uv10


def haversine_grid_km(lat_grid: np.ndarray, lon_grid: np.ndarray, lat0: float, lon0: float) -> np.ndarray:
    radius = 6371.0
    phi1 = np.radians(lat0)
    phi2 = np.radians(lat_grid)
    dphi = np.radians(lat_grid - lat0)
    dlambda = np.radians(lon_grid - lon0)
    a = np.sin(dphi / 2) ** 2 + np.cos(phi1) * np.cos(phi2) * np.sin(dlambda / 2) ** 2
    return 2 * radius * np.arcsin(np.sqrt(a))


def bearing_grid_deg(lat_grid: np.ndarray, lon_grid: np.ndarray, lat0: float, lon0: float) -> np.ndarray:
    phi1 = math.radians(lat0)
    phi2 = np.radians(lat_grid)
    dlambda = np.radians(lon_grid - lon0)
    y = np.sin(dlambda) * np.cos(phi2)
    x = math.cos(phi1) * np.sin(phi2) - math.sin(phi1) * np.cos(phi2) * np.cos(dlambda)
    return (np.degrees(np.arctan2(y, x)) + 360) % 360


def find_low_center(prmsl, guess_lat: float, guess_lon: float, radius_km: float) -> dict:
    lats = prmsl.latitude.values.astype(float)
    lons = prmsl.longitude.values.astype(float)
    lon_grid, lat_grid = np.meshgrid(lons, lats)
    distance = haversine_grid_km(lat_grid, lon_grid, guess_lat, guess_lon)
    values = prmsl.prmsl.values.astype(float)
    masked = np.where(distance <= radius_km, values, np.nan)
    if np.all(np.isnan(masked)):
        raise RuntimeError("no GFS low-pressure center found near prior position")
    i, j = np.unravel_index(np.nanargmin(masked), masked.shape)
    return {
        "lat": float(lats[i]),
        "lng": float(lons[j]),
        "prmsl_hpa": float(values[i, j]) / 100.0,
        "distance_from_guess_km": float(distance[i, j]),
    }


def quadrant_radius_string(uv10, center_lat: float, center_lon: float, threshold: float, radius_cap_km: float = 650.0) -> str:
    lats = uv10.latitude.values.astype(float)
    lons = uv10.longitude.values.astype(float)
    lon_grid, lat_grid = np.meshgrid(lons, lats)
    speed = np.hypot(uv10.u10.values.astype(float), uv10.v10.values.astype(float))
    distance = haversine_grid_km(lat_grid, lon_grid, center_lat, center_lon)
    bearing = bearing_grid_deg(lat_grid, lon_grid, center_lat, center_lon)
    valid = (speed >= threshold) & (distance <= radius_cap_km)
    quadrants = [
        ("ne", (bearing >= 0) & (bearing < 90)),
        ("se", (bearing >= 90) & (bearing < 180)),
        ("sw", (bearing >= 180) & (bearing < 270)),
        ("nw", (bearing >= 270) & (bearing < 360)),
    ]
    radii: list[int] = []
    for _, quad_mask in quadrants:
        selected = distance[valid & quad_mask]
        if selected.size == 0:
            radii.append(0)
        else:
            radii.append(int(round(float(np.nanmax(selected)) / 10.0) * 10))
    return "" if max(radii) <= 0 else "|".join(str(v) for v in radii)


def derive_wind_radii(uv10, center: dict) -> dict:
    return {
        field: quadrant_radius_string(uv10, center["lat"], center["lng"], threshold)
        for field, threshold in WIND_THRESHOLDS.items()
    }


def build_gfs_track(start_point: dict) -> tuple[datetime, list[int], list[dict]]:
    start_bj = typhoon.parse_point_time(start_point.get("time"))
    if start_bj is None:
        raise RuntimeError("official forecast point has no time")
    start_utc = bj_to_utc_naive(start_bj)

    last_error: Exception | None = None
    for cycle in candidate_cycles():
        raw_start_hour = (start_utc - cycle).total_seconds() / 3600.0
        start_hour = rounded_available_hour(raw_start_hour)
        hours = [start_hour + step for step in FORECAST_STEP_HOURS]
        if not all(gfs_hour_available(hour) for hour in hours):
            continue
        try:
            paths = [fetch_gfs_subset(cycle, hour) for hour in hours]
            guess_lat = float(typhoon.normalize(start_point.get("lat")))
            guess_lon = float(typhoon.normalize(start_point.get("lng")))
            points: list[dict] = []
            for idx, (hour, path) in enumerate(zip(hours, paths)):
                prmsl, uv10 = open_gfs_subset(path)
                center = find_low_center(prmsl, guess_lat, guess_lon, radius_km=800 if idx == 0 else 650)
                center["time"] = utc_to_bj_string(cycle + timedelta(hours=hour))
                center["forecast_hour"] = hour
                center.update(derive_wind_radii(uv10, center))
                points.append(center)
                guess_lat = center["lat"]
                guess_lon = center["lng"]
            return cycle, hours, points
        except Exception as exc:  # Try an older available cycle.
            last_error = exc
            continue
    raise RuntimeError(f"unable to fetch usable GFS cycle: {last_error}")


def draw_continuous_dashed_polyline(
    draw: ImageDraw.ImageDraw,
    points: list[tuple[float, float]],
    *,
    fill: tuple[int, int, int, int],
    width: int,
    dash: int,
    gap: int,
) -> None:
    """Draw a dashed polyline without resetting dash phase at every vertex.

    The generic typhoon.draw_dashed_polyline resets on every short circle
    segment; with 5-degree wind-circle vertices that visually becomes a solid
    line. This carries dash phase across vertices so closed wind circles stay
    visibly dashed.
    """
    if len(points) < 2:
        return
    cycle = max(1, dash + gap)
    phase = 0.0
    for start, end in zip(points, points[1:]):
        x1, y1 = start
        x2, y2 = end
        length = math.hypot(x2 - x1, y2 - y1)
        if length <= 0:
            continue
        ux = (x2 - x1) / length
        uy = (y2 - y1) / length
        pos = 0.0
        while pos < length:
            cycle_pos = phase % cycle
            in_dash = cycle_pos < dash
            remaining = (dash - cycle_pos) if in_dash else (cycle - cycle_pos)
            seg_end = min(length, pos + max(0.1, remaining))
            if in_dash:
                draw.line(
                    ((x1 + ux * pos, y1 + uy * pos), (x1 + ux * seg_end, y1 + uy * seg_end)),
                    fill=fill,
                    width=width,
                )
            advanced = seg_end - pos
            phase += advanced
            pos = seg_end


def draw_gfs_overlay(base_image_path: str, bounds: dict, latest: dict, official_point: dict, gfs_points: list[dict], provider_tracks: list[tuple[str, list[dict]]], cycle: datetime) -> str:
    image = Image.open(base_image_path).convert("RGBA")

    font = typhoon.scaled_font(image, 14)
    gfs_path_color = (210, 40, 255, 250)

    path_xys: list[tuple[float, float]] = []
    for point in gfs_points:
        xy = typhoon.point_to_image_xy(f"{point['lng']:.4f}", f"{point['lat']:.4f}", bounds, image.size)
        if xy:
            path_xys.append(xy)

    styles = (
        ("radius7", (0, 180, 255, 38), (0, 210, 255, 230)),
        ("radius10", (255, 210, 60, 34), (255, 230, 60, 245)),
        ("radius12", (255, 90, 90, 42), (255, 90, 90, 250)),
    )
    for field, fill, outline in styles:
        geom = typhoon.wind_swath_geometry(gfs_points, field)
        if geom is None or getattr(geom, "is_empty", True):
            continue
        typhoon.draw_geometry_swath(
            image,
            bounds,
            geom,
            fill=fill,
            outline=outline,
            width=typhoon.scaled_px(image, 1),
            dashed=True,
        )

    overlay = Image.new("RGBA", image.size, (0, 0, 0, 0))
    draw = ImageDraw.Draw(overlay)

    if len(path_xys) >= 2:
        for start, end in zip(path_xys, path_xys[1:]):
            typhoon.draw_dashed_line(
                draw,
                start,
                end,
                fill=gfs_path_color,
                width=typhoon.scaled_px(image, 2),
                dash=typhoon.scaled_px(image, 8),
                gap=typhoon.scaled_px(image, 6),
            )

    for idx, (point, xy) in enumerate(zip(gfs_points, path_xys)):
        x, y = xy
        radius = typhoon.scaled_px(image, 5 if idx == len(gfs_points) - 1 else 3)
        draw.ellipse((x - radius, y - radius, x + radius, y + radius), fill=gfs_path_color, outline=(255, 255, 255, 255), width=typhoon.scaled_px(image, 1))
        if idx in (0, len(gfs_points) - 1):
            label = typhoon.short_point_time_label({"time": point["time"]})
            text_xy = (x + typhoon.scaled_px(image, 10), y - typhoon.scaled_px(image, 20))
            draw.text(
                text_xy,
                label,
                fill=(0, 0, 0, 128),
                font=font,
                stroke_width=max(1, typhoon.scaled_px(image, 1)),
                stroke_fill=(255, 255, 255, 128),
            )

    route_palette = [
        (255, 70, 70, 245),
        (50, 180, 255, 245),
        (255, 190, 40, 245),
        (110, 230, 120, 245),
    ]
    merged_entries: list[tuple[str, tuple[int, int, int, int]]] = []
    for idx, (name, track_points) in enumerate(provider_tracks[:4]):
        if track_points:
            merged_entries.append((name, route_palette[idx % len(route_palette)]))
    merged_entries.extend([
        ("GFS数值场路径", gfs_path_color),
        ("7级风圈", (0, 210, 255, 230)),
        ("10级风圈", (255, 230, 60, 245)),
        ("12级风圈", (255, 90, 90, 250)),
    ])
    header = f"预报/GFS：{cycle.strftime('%m-%d %HZ')}"
    label_width = draw.textbbox((0, 0), header, font=font)[2]
    entry_width = max(draw.textbbox((0, 0), name, font=font)[2] for name, _ in merged_entries)
    box_left = typhoon.scaled_px(image, 8)
    box_width = max(label_width + typhoon.scaled_px(image, 32), entry_width + typhoon.scaled_px(image, 68), typhoon.scaled_px(image, 300))
    box_height = typhoon.scaled_px(image, 34) + typhoon.scaled_px(image, 22) * len(merged_entries)
    box_top = image.height - box_height - typhoon.scaled_px(image, 8)
    draw.rounded_rectangle(
        (box_left, box_top, box_left + box_width, box_top + box_height),
        radius=typhoon.scaled_px(image, 10),
        fill=(0, 0, 0, 128),
        outline=(255, 255, 255, 120),
        width=typhoon.scaled_px(image, 1),
    )
    draw.text((box_left + typhoon.scaled_px(image, 10), box_top + typhoon.scaled_px(image, 6)), header, fill=(255, 255, 255, 255), font=font)
    for idx, (name, color) in enumerate(merged_entries):
        yy = box_top + typhoon.scaled_px(image, 31 + 22 * idx)
        draw.line((box_left + typhoon.scaled_px(image, 12), yy + typhoon.scaled_px(image, 8), box_left + typhoon.scaled_px(image, 42), yy + typhoon.scaled_px(image, 8)), fill=color, width=typhoon.scaled_px(image, 4))
        draw.text((box_left + typhoon.scaled_px(image, 52), yy), name, fill=(255, 255, 255, 255), font=font)

    image.alpha_composite(overlay)
    output_path = str(Path(base_image_path).with_name(Path(base_image_path).stem + "_gfs.png"))
    image.save(output_path, dpi=(200, 200))
    return output_path


def media_path_from_report(report: str) -> str | None:
    match = re.search(r"MEDIA:(\S+)", report)
    return match.group(1) if match else None


def report_without_media(report: str) -> str:
    return re.sub(r"\n?MEDIA:\S+", "", report).rstrip()


def gfs_summary(cycle: datetime, points: list[dict]) -> str:
    lines = [f"GFS数值场预测（{cycle.strftime('%Y-%m-%d %HZ')} 起报；PRMSL低压中心 + 10米风阈值）："]
    for point in points:
        r7 = point.get("radius7") or "-"
        r10 = point.get("radius10") or "-"
        r12 = point.get("radius12") or "-"
        lines.append(
            f"- {typhoon.short_point_time_label(point)}: {point['lng']:.2f}°E, {point['lat']:.2f}°N; "
            f"{point['prmsl_hpa']:.1f}hPa; GFS风圈 7级={r7} 10级={r10} 12级={r12}"
        )
    return "\n".join(lines)


def main() -> int:
    os.environ.setdefault("TYPHOON_DRAW_CONSENSUS", "0")
    os.environ.setdefault("TYPHOON_SUPPRESS_TRACK_LEGEND", "1")
    base_report = typhoon.build_report()
    base_media = media_path_from_report(base_report)
    if not base_media or not Path(base_media).exists():
        print(base_report)
        return 0

    found = typhoon.find_typhoon()
    if not found:
        print(base_report)
        return 0
    tfid = typhoon.normalize(found.get("tfid") or found.get("value"))
    info = typhoon.fetch_json(f"/TyphoonInfo/{tfid}")
    points = info.get("points") if isinstance(info, dict) and isinstance(info.get("points"), list) else []
    if not points:
        print(base_report)
        return 0
    latest = points[-1]
    _, china_points = typhoon.select_forecast_provider(latest)
    if not china_points:
        print(base_report)
        return 0
    official_start = china_points[-1]

    cloud = typhoon.fetch_json("/LeastCloud/?type=24")
    bounds = {
        "min_lng": typhoon.normalize(cloud.get("minLng")),
        "max_lng": typhoon.normalize(cloud.get("maxLng")),
        "min_lat": typhoon.normalize(cloud.get("minLat")),
        "max_lat": typhoon.normalize(cloud.get("maxLat")),
    }
    if not all(bounds.values()):
        print(base_report)
        print("\nGFS数值场预测：云图边界缺失，未叠加。")
        return 0

    try:
        cycle, _, gfs_points = build_gfs_track(official_start)
        gfs_media = draw_gfs_overlay(base_media, bounds, latest, official_start, gfs_points, typhoon.list_forecast_providers(latest), cycle)
        print(report_without_media(base_report))
        print("\n" + gfs_summary(cycle, gfs_points))
        print(f"MEDIA:{gfs_media}")
        return 0
    except Exception as exc:
        print(report_without_media(base_report))
        print(f"\nGFS数值场预测：抓取或推断失败：{exc}")
        print(f"MEDIA:{base_media}")
        return 0