Исходный код soniks_client.waterfall.plot

"""Разбор сырого ``.dat`` водопада и отрисовка PNG для портала.

Здесь живёт весь matplotlib. Фигура создаётся явно (``Figure`` +
``FigureCanvasAgg``), без ``pyplot``: глобальный менеджер фигур не переживает
параллельные постобработки в воркерах планировщика.
"""

import json
import os
from datetime import UTC, datetime, timedelta

import matplotlib.colors as mcolors
import matplotlib.dates as mdates
import numpy as np
from matplotlib.axes import Axes
from matplotlib.backends.backend_agg import FigureCanvasAgg
from matplotlib.dates import date2num
from matplotlib.figure import Figure
from matplotlib.gridspec import GridSpec

from core.configs import logger, settings
from core.exceptions import EmptyWaterfallError, TimestampError, WaterfallFileError
from soniks_client.waterfall.deviation import _estimate_deviation_v7

# ═══════════════════════════════════════════════════════════════════════════
# Константы отрисовки
# ═══════════════════════════════════════════════════════════════════════════
WATERFALL_POWER_GAMMA = 1.4




class Waterfall:
    """Водопад одного наблюдения: разбор сырого ``.dat``, отрисовка и анализ.
    
    Raises:
        WaterfallFileError: Файл не найден или не читается.
        EmptyWaterfallError: Файл прочитан, но отсчётов нет.
    """

    def __init__(self, datafile_path: str, output_path: str):
        self.datafile_path = datafile_path
        self.output_path = output_path

        self.wf_data = None
        self.metadata = None

        # Заполняется в plot(); до вызова анализа нет, и get_signal_metadata()
        # обязан отдать None, а не пустую метадату.
        self.analysis: dict | None = None

        self._read_data()

        if not self.wf_data['data'].size:
            logger.error("Пустой файл для построения водопада: %s", self.datafile_path)
            raise EmptyWaterfallError()



    def plot(
        self,
        vmin: float | None = None,
        vmax: float | None = None,
    ) -> None:
        """Отрисовка водопада с анализом сигнала, маркерами и перекрестием"""

        est = None
        try:
            est = _estimate_deviation_v7(self.wf_data, self.metadata)
            self.analysis = est or {}
            if est:
                logger.info(
                    f"Анализ сигнала: {est.get('modulation')} | "
                    f"Δf={(est.get('delta_f_hz') or 0):.2f} Hz | "
                    f"SNR={(est.get('snr_db') or 0):.1f} dB | "
                    f"Метод: {est.get('delta_f_method', 'N/A')}"
                )
        except Exception as e:
            logger.warning(f"Ошибка оценки девиации: {e}")
            self.analysis = {}

        tmin = np.min(self.wf_data['data']['tabs'] / 10**6)
        tmax = np.max(self.wf_data['data']['tabs'] / 10**6)

        try:
            # Используем микросекунды (%S.%f), как в Mod 07d75e7
            timefmt = settings.waterfall.TIMESTAMP_FORMAT
            if '%f' not in timefmt:
                timefmt = '%Y-%m-%dT%H:%M:%S.%fZ'  # Fallback на точный формат

            # Метка в файле водопада всегда в UTC; делаем datetime aware, чтобы
            # ось и маркеры кадров (тоже aware) не смешивали naive и aware.
            t_ref = datetime.strptime(
                self.metadata['timestamp'], timefmt
            ).replace(tzinfo=UTC)
        except ValueError as e:
            logger.error("Некорректный формат временной метки в файле водопада")
            raise TimestampError() from e

        dt_min = t_ref + timedelta(seconds=tmin)
        dt_max = t_ref + timedelta(seconds=tmax)

        # Вычисляем середину наблюдения по времени (для горизонтальной линии)
        dt_mid = dt_min + (dt_max - dt_min) / 2.0

        spec = self.wf_data['data']['spec']

        # ── Обрезка краев SDR (rolloff) ──
        _freq_crop = self.wf_data['freq']
        fmin = float(_freq_crop[0] / 1000.0)
        fmax = float(_freq_crop[-1] / 1000.0)

        # ── Удаление DC spike ──
        spec_plot = spec.copy()
        center_bin = spec_plot.shape[1] // 2
        # Соседи берутся на ±2 бина: при более узком спектре center_bin + 2
        # выходит за массив, а center_bin - 2 молча заворачивается в хвост.
        if spec_plot.shape[1] >= 5:
            for dr in [-1, 0, 1]:
                spec_plot[:, center_bin + dr] = (spec_plot[:, center_bin - 2] + spec_plot[:, center_bin + 2]) / 2.0
        else:
            logger.debug(
                "Спектр из %s бинов: удаление DC-пика пропущено",
                spec_plot.shape[1],
            )

        # ── Адаптивные vmin/vmax с гамма-коррекцией ──
        if vmin is None or vmax is None:
            if settings.waterfall.AUTORANGE:
                # np.isfinite() убирает -inf пустых бинов, WATERFALL__THRESHOLD
                # — заведомо нефизические уровни. Порог именно в дополнение к
                # isfinite: сам по себе он отсекал реальный шум при низком
                # усилении SDR, поэтому дефолт (-200 дБ) намеренно щадящий.
                _flat = spec_plot[
                    np.isfinite(spec_plot) & (spec_plot > settings.waterfall.THRESHOLD)
                ]
                if _flat.size > settings.waterfall.MIN_VALID_SAMPLES:
                    # Робастная оценка уровня шума через сигма-клиппинг
                    _valid = _flat.copy()
                    for _ in range(3):
                        _med = float(np.median(_valid))
                        _mad = float(np.median(np.abs(_valid - _med))) * 1.4826
                        if _mad <= 0:
                            break
                        _clip = np.abs(_valid - _med) < 3.0 * _mad
                        if np.sum(_clip) < 5:
                            break
                        _valid = _valid[_clip]

                    if _valid.size > 5:
                        _noise = float(np.mean(_valid))
                        _noise_std = float(np.std(_valid))
                    else:
                        _noise = float(np.median(_flat))
                        _noise_std = max(float(np.std(_flat)), 10.0)

                    vmin = _noise - 0.5 * _noise_std
                    vmax = max(float(np.percentile(_flat, 99.5)), _noise + 1.5 * _noise_std)

                    # Гарантия: vmin строго меньше vmax
                    if vmax <= vmin:
                        vmin = float(np.percentile(_flat, 1.0))
                        vmax = float(np.percentile(_flat, 99.0))

                    logger.debug(
                        f"Waterfall AUTORANGE vmin={vmin:.1f} vmax={vmax:.1f} "
                        f"(noise={_noise:.1f} dB, std={_noise_std:.1f} dB, "
                        f"valid={_valid.size}/{_flat.size})"
                    )
                else:
                    vmin = settings.waterfall.DEFAULT_MIN_VALUE
                    vmax = settings.waterfall.DEFAULT_MAX_VALUE
                    logger.warning(
                        f"Waterfall: недостаточно данных для AUTORANGE "
                        f"(finite values: {_flat.size}). Использованы значения по умолчанию."
                    )
            else:
                # Если AUTORANGE=False, используем фиксированные значения из конфига
                vmin = settings.waterfall.MIN_VALUE
                vmax = settings.waterfall.MAX_VALUE
                logger.debug(
                    f"Waterfall FIXED vmin={vmin} vmax={vmax}"
                )

        norm = mcolors.PowerNorm(gamma=WATERFALL_POWER_GAMMA, vmin=vmin, vmax=vmax, clip=True)

        # Отрисовка графика. Фигура строится явно, без pyplot: постпроцессинг
        # идёт в воркер-потоках APScheduler, а глобальный менеджер фигур один
        # на процесс — два наложившихся прохода портили бы графики друг другу.
        fig = Figure(figsize=(settings.waterfall.WIDTH, settings.waterfall.HEIGHT))
        FigureCanvasAgg(fig)
        gs = GridSpec(**settings.waterfall.GRIDSPEC)
        ax = fig.add_subplot(gs[0])
        cax = fig.add_subplot(gs[1])

        im = ax.imshow(
            spec_plot, origin='lower', aspect='auto', interpolation='bilinear',
            extent=[fmin, fmax, date2num(dt_min), date2num(dt_max)],
            norm=norm,
            cmap='jet'
        )
        # Настройка осей
        # tz=UTC задан явно: иначе matplotlib берёт rcParams['timezone'] и
        # подписи оси начинают зависеть от настроек хоста.
        ax.yaxis_date(tz=UTC)
        ax.yaxis.set_major_locator(
            mdates.MinuteLocator(interval=settings.waterfall.INTERVAL_IN_MINUTES, tz=UTC)
        )
        ax.yaxis.set_major_formatter(
            mdates.DateFormatter(settings.waterfall.TIME_FORMAT, tz=UTC)
        )
        ax.set_xlabel("Частота (кГц)")
        ax.set_ylabel("Время (UTC)")

        cbar = fig.colorbar(im, aspect=50, cax=cax)
        cbar.set_label("Мощность (дБ)", labelpad=4)

        # ── Перекрестие: Центр по частоте (0 кГц) и Центр по времени ──
        self._plot_crosshair(ax, dt_mid)

        # ── Маркеры декодированных кадров (Красные/Оранжевые) ──
        self._plot_decoded_frames_markers(ax, dt_min, dt_max)

        # ── Аннотация девиации Δf ──
        if self.analysis.get("f_pos_hz") is not None:
            self._plot_deviation_annotation(ax, self.analysis)

        ax.set_ylim(date2num(dt_min), date2num(dt_max))

        # Сохранение с метаданными
        metadata = {
            'satnogs:wf-dat': json.dumps({k: str(v) for k, v in self.metadata.items()}),
            'satnogs:wf-signal': json.dumps(self._format_signal_metadata())
        }
        fig.savefig(self.output_path, metadata=metadata, bbox_inches="tight")


    def _plot_crosshair(self, ax: Axes, dt_mid: datetime):
        """Отрисовка тонкой красной линии по центру частоты (0 кГц) и времени"""
        if not settings.waterfall.HELP_LINES:
            return

        ax.axvline(x=0.0, color='red', linewidth=0.4, alpha=0.5, linestyle='-', zorder=10)
        ax.axhline(y=date2num(dt_mid), color='red', linewidth=0.4, alpha=0.5, linestyle='-', zorder=10)

    def _plot_decoded_frames_markers(self, ax: Axes, dt_min: datetime, dt_max: datetime):
        """Чтение decoded_frames.jsonl и отрисовка маркеров (Mod 215bf7e)"""
        try:
            data_dir = os.path.dirname(os.path.abspath(self.datafile_path))
            log_file = os.path.join(data_dir, 'decoded_frames.jsonl')
            if not os.path.isfile(log_file):
                return

            ts_flowgraph, ts_grsatellites = [], []
            with open(log_file, encoding='utf-8') as f:
                for line in f:
                    entry = json.loads(line.strip())
                    ts_str = entry.get('timestamp')
                    if not ts_str:
                        continue

                    try:
                        ts = datetime.fromisoformat(ts_str).astimezone(UTC)
                        # Защита от 1970 года
                        if not (dt_min - timedelta(minutes=5) < ts < dt_max + timedelta(minutes=5)):
                            continue

                        if entry.get('source') == 'gr-satellites':
                            ts_grsatellites.append(ts)
                        else:
                            ts_flowgraph.append(ts)
                    except ValueError:
                        continue

            trans = ax.get_yaxis_transform()
            if ts_flowgraph:
                ax.plot([0.015]*len(ts_flowgraph), [date2num(ts) for ts in ts_flowgraph], 
                        marker='o', color='red', markersize=6, linestyle='None', transform=trans, clip_on=False, zorder=1000)
            if ts_grsatellites:
                ax.plot([0.025]*len(ts_grsatellites), [date2num(ts) for ts in ts_grsatellites], 
                        marker='o', color='#FF8C00', markersize=5, linestyle='None', transform=trans, clip_on=False, zorder=1001)
        except Exception as e:
            logger.debug(f"Не удалось отрисовать маркеры кадров: {e}")

    def _plot_deviation_annotation(self, ax: Axes, est: dict):
        """Отрисовка аннотации Δf под графиком"""
        f_pos_khz = est["f_pos_hz"] / 1000.0
        f_neg_khz = est["f_neg_hz"] / 1000.0
        delta_f_khz = est["delta_f_hz"] / 1000.0
        center_khz = (f_pos_khz + f_neg_khz) / 2.0
        ann_color = '#C0392B' if est.get("confidence", 0) > 0.5 else '#D4760A'
        trans = ax.get_xaxis_transform()
        y0, th = -0.035, 0.014

        for fx in (f_neg_khz, f_pos_khz):
            ax.plot([fx, fx], [y0, y0 - th], color=ann_color, linewidth=1.4, transform=trans, clip_on=False, zorder=100)

        ax.annotate('', xy=(f_neg_khz, y0 - th*0.5), xytext=(f_pos_khz, y0 - th*0.5),
                    xycoords=trans, textcoords=trans,
                    arrowprops={'arrowstyle': '<->', 'color': ann_color, 'lw': 1.0}, clip_on=False, zorder=100)

        ax.text(center_khz, y0 - th - 0.012, f'\u0394f \u2248 {delta_f_khz:.2f} kHz',
                color=ann_color, fontsize=6.5, fontweight='bold', ha='center', va='top', transform=trans, clip_on=False,
                bbox={'boxstyle': 'round,pad=0.2', 'facecolor': 'white', 'edgecolor': ann_color, 'linewidth': 0.6, 'alpha': 0.92})

    def _format_signal_metadata(self) -> dict:
        """Форматирование результатов анализа для встраивания в PNG"""
        e = self.analysis
        def _fmt(v, d=2): return round(float(v), d) if v is not None else None
        return {
            'modulation': e.get('modulation', 'unknown'),
            'deviation_hz': _fmt(e.get('delta_f_hz'), 1),
            'deviation_confidence': _fmt(e.get('confidence')),
            'snr_db': _fmt(e.get('snr_db')),
            'bw_99_hz': _fmt(e.get('bw_99_hz'), 1),
        }

    # ══════════════════════════════════════════════════════════════════════
    # Signal analysis interface for metadata export
    # ══════════════════════════════════════════════════════════════════════



    def get_signal_metadata(self) -> dict | None:
        """Return formatted signal metadata for API submission.

        Returns dict with string-formatted values matching the server schema,
        or None if analysis was not performed.
        """
        if self.analysis is None:
            return None

        r = self.analysis

        def _fmt_hz(val) -> str | None:
            if val is None:
                return None
            return f"{val:.1f} Hz"

        def _fmt_db(val) -> str | None:
            if val is None:
                return None
            return f"{val:.1f} dB"

        def _fmt_mhz(val) -> str | None:
            if val is None:
                return None
            return f"{val / 1e6:.6f} MHz"

        def _fmt_ppm(val) -> str | None:
            if val is None:
                return None
            return f"{val:.3f} ppm"

        signal = {
            "snr_db": _fmt_db(r.get("snr_db")),
            "snr_quality": f"{r.get('snr_quality', 0):.2f}" if r.get("snr_quality") is not None else None,
            "noise_db": _fmt_db(r.get("noise_db")),
            "noise_std_db": _fmt_db(r.get("noise_std_db")),
            "peak_db": _fmt_db(r.get("peak_db")),
            "peak_freq_hz": _fmt_hz(r.get("peak_f_hz")),
            "modulation": r.get("modulation"),
            "n_spec_peaks": r.get("n_spec_peaks"),
            "valley_depth_db": _fmt_db(r.get("valley_depth_db")),
            "deviation_hz": _fmt_hz(r.get("delta_f_hz")),
            "deviation_uncertainty_hz": _fmt_hz(r.get("delta_f_uncertainty_hz")),
            "deviation_confidence": f"{r.get('confidence', 0):.2f}" if r.get("confidence") is not None else None,
            "deviation_method": r.get("delta_f_method"),
            "f_low_hz": _fmt_hz(r.get("f_neg_hz")),
            "f_high_hz": _fmt_hz(r.get("f_pos_hz")),
            "bw_3db_hz": _fmt_hz(r.get("bw_3db_hz")),
            "bw_10db_hz": _fmt_hz(r.get("bw_10db_hz")),
            "bw_20db_hz": _fmt_hz(r.get("bw_20db_hz")),
            "bw_26db_hz": _fmt_hz(r.get("bw_26db_hz")),
            # Ключ на портале — bw99_hz, а оценщики пишут его как bw_99_hz;
            # bw99_hz в их словаре — legacy-ключ, которому никто не присваивает
            # значения.
            "bw99_hz": _fmt_hz(r.get("bw_99_hz")),
            "sigma_rms_hz": _fmt_hz(r.get("sigma_rms_hz")),
            "n_bursts": r.get("n_bursts"),
            "n_valid_rows": r.get("n_valid_rows"),
            # Опора цели упёрлась в край полосы: сигнал шире водопада либо
            # приёмник перегружен. Bool, не строка — порталу считать долю.
            "saturated": r.get("saturated"),
            "carrier_offset_hz": _fmt_hz(r.get("center_freq_offset_hz")),
            "ppm_error": _fmt_ppm(r.get("ppm_error")),
            "drift_ppb": f"{r.get('ppb_error', 0):.1f}" if r.get("ppb_error") is not None else None,
            "measured_freq_mhz": _fmt_mhz(r.get("measured_freq_hz")),
            "method": r.get("note"),
        }

        # Геометрия водопада — из заголовка .dat, а не из analysis: заголовок
        # есть всегда, оценщик при nchan < 16 выходит рано. Полоса и ширина
        # бина зависят от режима (Фаза 3.3 в docs/roadmap-network.md), а по
        # правилу 21 старые станции шлют разнородные водопады вечно — значит
        # сравнивать уровни шума между проходами можно только зная геометрию.
        # noise_dbhz = noise_db − 10·lg(bin_hz): для белого шума мощность в
        # бине ∝ N0·bin_hz при прямоугольном окне и масштабе 1/fft_size
        # waterfall_sink. Смещение Max Hold (10·lg(H_n) по nfft_per_row) не
        # вычитается: режим sink в заголовке не записан.
        samp_rate = float(self.metadata["samp_rate"])
        nchan = int(self.metadata["nchan"])
        bin_hz = samp_rate / nchan if nchan else None
        signal["samp_rate_hz"] = _fmt_hz(samp_rate)
        signal["bin_hz"] = _fmt_hz(bin_hz)
        signal["nchan"] = nchan
        signal["nfft_per_row"] = int(self.metadata["nfft_per_row"])
        if bin_hz and r.get("noise_db") is not None:
            signal["noise_dbhz"] = f"{r['noise_db'] - 10 * np.log10(bin_hz):.1f} dB/Hz"

        return {k: v for k, v in signal.items() if v is not None}


    @property
    def _frequencies(self) -> np.ndarray:
        return self.wf_data['freq'] / 1000.0

    def _read_data(self) -> None:
        """Чтение данных с адаптацией под внутренний формат DSP-алгоритмов"""
        try:
            with open(self.datafile_path, "rb") as datafile:
                metadata = {
                    'timestamp': np.fromfile(datafile, dtype="|S32", count=1)[0].decode("utf-8"),
                    'nchan': np.fromfile(datafile, dtype=">i4", count=1)[0],
                    'samp_rate': np.fromfile(datafile, dtype=">i4", count=1)[0],
                    'nfft_per_row': np.fromfile(datafile, dtype=">i4", count=1)[0],
                    'center_freq': np.fromfile(datafile, dtype=">f4", count=1)[0],
                    'endianness': np.fromfile(datafile, dtype="<i4", count=1)[0],
                }
                dtype_prefix = '<' if metadata['endianness'] else '>'
                data_dtypes = np.dtype([
                    ('tabs', dtype_prefix + 'i8'),
                    ('spec', dtype_prefix + 'f4', (metadata['nchan'],))
                ])
                spectral_data = np.fromfile(datafile, dtype=data_dtypes)

                nint = spectral_data.shape[0]
                wf_data = {
                    'data': spectral_data,
                    'trel': np.arange(nint) * metadata['nfft_per_row'] * metadata['nchan'] / float(metadata['samp_rate']),
                    'freq': np.linspace(
                        -0.5 * metadata['samp_rate'], 0.5 * metadata['samp_rate'],
                        metadata['nchan'], endpoint=False
                    ),
                    'compressed': None
                }

                self.metadata = metadata
                self.wf_data = wf_data

        except FileNotFoundError as e:
            logger.error("Файл водопада не найден: %s", self.datafile_path)
            raise WaterfallFileError() from e
        except OSError as e:
            logger.error("Ошибка чтения файла водопада: %s", self.datafile_path)
            raise WaterfallFileError() from e
        except (IndexError, ValueError, KeyError, UnicodeDecodeError) as e:
            # Обрезанный или битый заголовок: np.fromfile(...)[0] даёт IndexError,
            # неUTF-8 метка времени — UnicodeDecodeError. Без этой ветки они
            # уходили наверх мимо контракта исключений водопада.
            logger.error(
                "Повреждён файл водопада: %s (%s)", self.datafile_path, e
            )
            raise WaterfallFileError() from e