noaa-goes-visualization/suvi/measures.py

"""Per-frame measurements the judgment is made from, taken once per frame.

Reading and decoding the archive is the expensive part of judging it (~2 days for
2023-01..2025-03); comparing numbers is not.  So ``filter_FITS.py measure`` reads
every frame once and stores what these functions return, and the rules that
compare a frame with its neighbours run later, from the stored numbers, in
minutes -- which is what lets their limits be tuned against review sheets.

Everything here is a plain measurement with no verdict.  Images are worked on at
MEASURE_SIZE (a quarter of the native 1280), radii come from the frame's own
DIAM_SUN, and positions are in units of the solar radius.
"""

import math

import cv2 as cv
import numpy as np

#: Working resolution: a 386px radius becomes ~97px, ample for medians, a limb
#: position to ~1% of the radius, and a double exposure's echo.
MEASURE_SIZE = 320
#: Angular sectors of the disc for regional brightness.  Eight 45-degree sectors
#: resolve the partial dimming of a grazing eclipse (one side of the disc behind
#: Earth's atmosphere) that a whole-disc number dilutes.
SECTORS = 8
#: Radii, in solar radii, the disc brightness is taken over: inside the limb, so
#: the limb's own brightening and a displaced disc's edge stay out, and clear of
#: the centre, where every sector would share the same pixels.
SECTOR_BAND = (0.2, 0.9)
#: Where limb edges are looked for, and how far a centre may sit from the frame's.
#: Beyond 1.45R the instrument's field-of-view edge would win the vote.
EDGE_LIMIT = 1.45
CENTRE_LIMIT = 0.8
#: Lags searched for a double exposure's echo, in solar radii: beyond a feature's
#: own width, within the disc's span.
ECHO_LAGS = (0.08, 1.6)

_EARTH_RADIUS_KM = 6378.137
_SUN_RADIUS_KM = 695700.0
_AU_KM = 149597870.7


def _working(image, size=MEASURE_SIZE):
    """Area-averaged radiance at the working size, NaN as zero, negatives clipped
    after averaging -- the order the position limits were calibrated with: noise
    around zero averages out instead of being rectified into a faint glow."""
    finite = np.nan_to_num(np.asarray(image, dtype=np.float32), nan=0.0, posinf=0.0,
                           neginf=0.0)
    return np.clip(cv.resize(finite, (size, size), interpolation=cv.INTER_AREA), 0.0, None)


def _grid(size=MEASURE_SIZE):
    yy, xx = np.mgrid[:size, :size].astype(np.float32)
    centre = size / 2.0
    return yy - centre, xx - centre


def working_radius(image, diam_sun, size=MEASURE_SIZE):
    """The solar radius in working pixels, from the header's full-resolution
    DIAM_SUN, or None when there is none."""
    if not diam_sun or not np.isfinite(diam_sun) or diam_sun <= 0:
        return None
    return float(diam_sun) / 2.0 * size / np.asarray(image).shape[0]


def disc_brightness(linear, radius):
    """Median radiance of the disc and of each of SECTORS angular sectors.

    Medians, not means: a flare kernel lifts a frame's mean (and IMG_MEAN) by tens
    of percent while barely moving its median, so a frame beside a flare does not
    look dim against it.  `linear` is the working-size radiance image.
    """
    dy, dx = _grid(linear.shape[0])
    r = np.hypot(dy, dx) / radius
    angle = np.mod(np.arctan2(dy, dx), 2 * np.pi)
    low, high = SECTOR_BAND
    disc = r < high
    out = {"disc_median": float(np.median(linear[disc])) if disc.any() else None}
    ring = (r >= low) & (r < high)
    width = 2 * np.pi / SECTORS
    for k in range(SECTORS):
        mask = ring & (angle >= k * width) & (angle < (k + 1) * width)
        out[f"sector_{k}"] = float(np.median(linear[mask])) if mask.any() else None
    return out


def _edges(log_image, radius):
    """Edge strength of the log image, capped at its 99th percentile and zero
    beyond EDGE_LIMIT, so no handful of edges (a flare kernel's) can dominate."""
    dy, dx = _grid(log_image.shape[0])
    gx = cv.Sobel(log_image, cv.CV_32F, 1, 0, ksize=3)
    gy = cv.Sobel(log_image, cv.CV_32F, 0, 1, ksize=3)
    strength = np.hypot(gx, gy)
    inside = np.hypot(dy, dx) < EDGE_LIMIT * radius
    strength[~inside] = 0.0
    cap = float(np.percentile(strength[inside], 99)) if inside.any() else 0.0
    return np.minimum(strength, cap) / cap if cap > 0 else strength


def ring_position(log_image, radius):
    """Where the disc is: the centre of the ring of radius `radius` that the frame's
    edges fit best, as an offset from the frame centre in solar radii.

    An intensity centroid is dragged by active regions, and the limb fitted angle
    by angle fails where too few angles show an edge.  Correlating the whole edge
    image with a ring uses every limb pixel at once and lets no bright region
    vote alone.  On the labelled set clean frames, X-flares included, all stay
    within 0.08R; 98% of displaced discs lie beyond it.
    """
    size = log_image.shape[0]
    edges = _edges(log_image, radius)
    dy, dx = _grid(size)
    r = np.hypot(dy, dx)
    ring = (np.abs(r - radius) < 1.0).astype(np.float32)
    match = np.real(np.fft.ifft2(np.fft.fft2(edges) * np.conj(np.fft.fft2(ring))))
    match = np.fft.fftshift(match)                 # zero shift at the centre pixel
    match[r > CENTRE_LIMIT * radius] = -np.inf
    y, x = np.unravel_index(int(np.argmax(match)), match.shape)
    oy, ox = (y - size / 2.0) / radius, (x - size / 2.0) / radius
    return {"ring_dx": float(ox), "ring_dy": float(oy), "ring_offset": float(math.hypot(ox, oy))}


def echo(log_image, radius):
    """How strongly the frame repeats itself at some offset: a double exposure is
    the image plus a shifted copy, so its structure's autocorrelation has a second
    peak at the shift.  Returns the strongest off-centre autocorrelation (zero lag
    = 1) and its lag in solar radii.

    The structure is the log image minus its smooth part, clipped so a flare
    kernel cannot carry the measure; a flare's diffraction spikes still repeat,
    which is why 94A and 131A flares reach ~0.4 and the limit is set per band.
    """
    size = log_image.shape[0]
    dy, dx = _grid(size)
    r = np.hypot(dy, dx)
    structure = log_image - cv.GaussianBlur(log_image, (0, 0), 4.0)
    inside = r < EDGE_LIMIT * radius
    structure[~inside] = 0.0
    if not inside.any():
        return {"echo": None, "echo_lag": None}
    limit = float(np.percentile(np.abs(structure[inside]), 99.5))
    structure = np.clip(structure, -limit, limit)
    window = np.outer(np.hanning(size), np.hanning(size)).astype(np.float32)
    spectrum = np.fft.fft2(structure * window)
    auto = np.fft.fftshift(np.real(np.fft.ifft2(spectrum * np.conj(spectrum))))
    zero = auto[size // 2, size // 2]
    if zero <= 0:
        return {"echo": None, "echo_lag": None}
    low, high = ECHO_LAGS
    band = (r > low * radius) & (r < high * radius)
    candidates = np.where(band, auto / zero, -np.inf)
    y, x = np.unravel_index(int(np.argmax(candidates)), candidates.shape)
    return {"echo": float(candidates[y, x]), "echo_lag": float(r[y, x] / radius)}


def flare_index(linear, radius, median):
    """The disc's 99.99th-percentile radiance over its median: ~10 for a quiet Sun,
    hundreds to thousands under a flare kernel."""
    dy, dx = _grid(linear.shape[0])
    disc = np.hypot(dy, dx) < SECTOR_BAND[1] * radius
    if not disc.any() or not median or median <= 0:
        return {"flare_index": None}
    return {"flare_index": float(np.percentile(linear[disc], 99.99) / median)}


def pixel_measures(image, diam_sun):
    """Every pixel measure for one decoded frame; {} when it has no usable radius."""
    if image is None or np.asarray(image).ndim != 2:
        return {}
    radius = working_radius(image, diam_sun)
    if radius is None:
        return {}
    linear = _working(image)
    log_image = cv.GaussianBlur(np.log1p(linear), (0, 0), 1.0)
    out = disc_brightness(linear, radius)
    out.update(ring_position(log_image, radius))
    out.update(echo(log_image, radius))
    out.update(flare_index(linear, radius, out["disc_median"]))
    return out


# ------------------------------------------------------------------ Earth geometry


def _julian_date(unix_seconds):
    return unix_seconds / 86400.0 + 2440587.5


def sun_direction(unix_seconds):
    """Unit vector to the Sun, equatorial coordinates of date, and the distance in
    km: the Astronomical Almanac's low-precision formula, good to ~0.01 degrees."""
    n = _julian_date(unix_seconds) - 2451545.0
    mean_longitude = math.radians((280.460 + 0.9856474 * n) % 360.0)
    anomaly = math.radians((357.528 + 0.9856003 * n) % 360.0)
    ecliptic = mean_longitude + math.radians(1.915 * math.sin(anomaly)
                                             + 0.020 * math.sin(2 * anomaly))
    obliquity = math.radians(23.439 - 0.0000004 * n)
    distance = (1.00014 - 0.01671 * math.cos(anomaly)
                - 0.00014 * math.cos(2 * anomaly)) * _AU_KM
    direction = np.array([math.cos(ecliptic),
                          math.cos(obliquity) * math.sin(ecliptic),
                          math.sin(obliquity) * math.sin(ecliptic)])
    return direction, distance


def _sidereal_angle(unix_seconds):
    """Greenwich mean sidereal time, radians."""
    days = _julian_date(unix_seconds) - 2451545.0
    return math.radians((280.46061837 + 360.98564736629 * days) % 360.0)


def earth_tangent_km(obsgeo_m, unix_seconds):
    """How high above Earth's surface the line of sight to the Sun's nearest edge
    passes, in km; None when the Sun is on the far side of the satellite from Earth.

    Around local midnight in the eclipse seasons the Sun drops behind Earth as
    seen from geostationary orbit.  NOAA's ECLIPSE flag marks only the shadow
    itself; before it the disc grazes the atmosphere, which absorbs EUV, and one
    side of it dims over several frames.  `obsgeo_m` is the satellite position
    (OBSGEO-X/Y/Z, metres, Earth-fixed).
    """
    theta = _sidereal_angle(unix_seconds)
    x, y, z = (float(v) / 1000.0 for v in obsgeo_m)
    satellite = np.array([x * math.cos(theta) - y * math.sin(theta),
                          x * math.sin(theta) + y * math.cos(theta), z])
    direction, distance = sun_direction(unix_seconds)
    towards = direction * distance - satellite
    span = float(np.linalg.norm(towards))
    unit = towards / span
    along = -float(satellite @ unit)               # to the point nearest Earth
    if along <= 0:
        return None
    closest = float(np.linalg.norm(satellite + along * unit)) - _EARTH_RADIUS_KM
    return closest - along * (_SUN_RADIUS_KM / span)


def header_measures(cards):
    """Measures from the header alone: the Earth tangent height."""
    cards = cards or {}
    position = [cards.get(k) for k in ("OBSGEO-X", "OBSGEO-Y", "OBSGEO-Z")]
    when = cards.get("DATE-OBS") or cards.get("DATE-BEG")
    if not when or not all(isinstance(v, (int, float)) for v in position):
        return {}
    from datetime import datetime, timezone

    try:
        moment = datetime.fromisoformat(str(when).replace("Z", "")).replace(
            tzinfo=timezone.utc).timestamp()
    except ValueError:
        return {}
    height = earth_tangent_km(position, moment)
    return {} if height is None else {"earth_tangent_km": float(height)}