Détecteur de parapente — paraglider_lib

Librairie Python légère pour détecter si une trace GPS correspond à un vol de parapente.

Principe

Lors d'un vol en parapente, les vecteurs vitesses (Vx, Vy) forment un anneau circulaire. La librairie :

1. Calcule les vecteurs vitesse entre points GPS consécutifs
2. Filtre les valeurs aberrantes (IQR)
3. Fitte un cercle par la méthode de Taubin (robuste, sans biais)
4. Vérifie que le rayon et la régularité du cercle correspondent à un parapente

Installation

pip install numpy

Seuils (modifiables en haut de paraglider_lib.py)

API

GPSPoint

from paraglider_lib import GPSPoint

pt = GPSPoint(lat=45.632, lon=6.847, timestamp=1725187200.0)

analyser(points)

R_kmh, sigma_pct, is_para = analyser(points)

Entrée
- points : liste de GPSPoint — n'importe quelle taille (minimum ~30 points)

Sorties

Critères de détection :
- R_MIN_KMH ≤ R_kmh ≤ R_MAX_KMH
- sigma_pct < MAX_RESIDUAL_PCT * 100

Exemple rapide

from paraglider_lib import analyser, GPSPoint

points = [
    GPSPoint(45.6320, 6.8470, 1725187200.0),
    GPSPoint(45.6322, 6.8475, 1725187201.0),
    # ... au moins ~30 points
]

R_kmh, sigma_pct, is_para = analyser(points)

print(f"Rayon    : {R_kmh:.1f} km/h")
print(f"sigma/R  : {sigma_pct:.1f} %")
print(f"Parapente: {is_para}")

Voir exemple.py pour un exemple complet sur fichier IGC.

Recommandations d'usage

• Fréquence GPS : fonctionne à 1 Hz (idéal) ou avec un point toutes les 5–6 s (traces 2009)
• Nombre de points : au minimum 30, idéalement 300–600 (5–10 minutes à 1 Hz)
• Ignorer le décollage : supprimer la première minute évite les biais liés au sol
• Pas de lissage temporel : ne pas interpoler ou moyenner les points avant d'appeler analyser() — cela détruirait la structure circulaire dans l'espace des vitesses

exemple.py

"""
Exemple d'utilisation de paraglider_lib sur un fichier IGC.

Usage :
    python exemple.py trace.igc
"""

import re
import sys
from pathlib import Path

from paraglider_lib import GPSPoint, analyser

# ---------------------------------------------------------------------------
# Parser IGC minimal
# ---------------------------------------------------------------------------

_RE_B = re.compile(
    r"^B(\d{2})(\d{2})(\d{2})"   # HHMMSS
    r"(\d{2})(\d{5})([NS])"       # DDMMmmm N/S
    r"(\d{3})(\d{5})([EW])"       # DDDMMmmm E/W
    r"A",                          # Fix valide
)


def _coord(ddmm: str, hemi: str) -> float:
    deg = int(ddmm[:2 if len(ddmm) == 7 else 3])
    minutes = int(ddmm[2 if len(ddmm) == 7 else 3:]) / 1000.0
    v = deg + minutes / 60.0
    return -v if hemi in ("S", "W") else v


def parse_igc(path: Path) -> list[GPSPoint]:
    """Parse un fichier IGC et retourne une liste de GPSPoint."""
    import datetime

    flight_date = None
    points: list[GPSPoint] = []

    with open(path, encoding="latin-1", errors="replace") as f:
        for line in f:
            line = line.strip()
            if line.startswith("HFDTE") or line.startswith("HFDTEDATE"):
                m = re.search(r"(\d{2})(\d{2})(\d{2})", line)
                if m:
                    dd, mm, yy = int(m.group(1)), int(m.group(2)), int(m.group(3))
                    try:
                        flight_date = datetime.date(2000 + yy, mm, dd)
                    except ValueError:
                        flight_date = datetime.date(2000, 1, 1)
                continue

            m = _RE_B.match(line)
            if not m:
                continue

            hh, mi, ss = int(m.group(1)), int(m.group(2)), int(m.group(3))
            lat = _coord(m.group(4) + m.group(5), m.group(6))
            lon = _coord(m.group(7) + m.group(8), m.group(9))

            date = flight_date or datetime.date(2000, 1, 1)
            ts = datetime.datetime(
                date.year, date.month, date.day, hh, mi, ss,
                tzinfo=datetime.timezone.utc,
            ).timestamp()

            points.append(GPSPoint(lat=lat, lon=lon, timestamp=ts))

    # Gestion du passage minuit
    for i in range(1, len(points)):
        while points[i].timestamp < points[i - 1].timestamp:
            points[i] = GPSPoint(points[i].lat, points[i].lon, points[i].timestamp + 86400)

    return points


# ---------------------------------------------------------------------------
# Programme principal
# ---------------------------------------------------------------------------

def main():
    if len(sys.argv) < 2:
        print("Usage : python exemple.py <fichier.igc>")
        sys.exit(1)

    igc_path = Path(sys.argv[1])
    if not igc_path.exists():
        print(f"Fichier introuvable : {igc_path}")
        sys.exit(1)

    print(f"\nFichier : {igc_path.name}")

    # 1. Charger les points GPS
    points = parse_igc(igc_path)
    print(f"Points GPS chargés : {len(points)}")

    if not points:
        print("Aucun point valide.")
        sys.exit(1)

    # 2. Ignorer la première minute (décollage), prendre les 600 suivants
    t0 = points[0].timestamp
    points = [p for p in points if p.timestamp >= t0 + 60][:600]
    print(f"Points analysés    : {len(points)} (après suppression de la 1re minute)")

    # 3. Appel de la librairie
    R_kmh, sigma_pct, is_para = analyser(points)

    # 4. Affichage des résultats
    print(f"\n--- Résultats ---")
    print(f"Rayon du cercle    : {R_kmh:.1f} km/h")
    print(f"Ecart-type résidus : {sigma_pct:.1f} %")
    print(f"Verdict            : {'✓ PARAPENTE DÉTECTÉ' if is_para else '✗ non détecté'}")


if __name__ == "__main__":
    main()

paraglider_lib.py

"""
Librairie de détection de parapente par analyse des vecteurs vitesse GPS.

Utilisation :
    from paraglider_lib import analyser, GPSPoint

    points = [GPSPoint(lat, lon, timestamp), ...]
    R_kmh, sigma_pct, is_para = analyser(points)
"""

from __future__ import annotations

import math
from typing import NamedTuple

import numpy as np


# ---------------------------------------------------------------------------
# Seuils de détection
# ---------------------------------------------------------------------------

R_MIN_KMH        = 25.0   # rayon minimal du cercle de vitesse (km/h)
R_MAX_KMH        = 55.0   # rayon maximal du cercle de vitesse (km/h)
MAX_RESIDUAL_PCT = 0.20   # écart-type des résidus < 20 % du rayon


# ---------------------------------------------------------------------------
# Type d'entrée
# ---------------------------------------------------------------------------

class GPSPoint(NamedTuple):
    lat:       float   # degrés décimaux
    lon:       float   # degrés décimaux
    timestamp: float   # timestamp Unix (secondes)


# ---------------------------------------------------------------------------
# Fonction principale
# ---------------------------------------------------------------------------

def analyser(points: list[GPSPoint]) -> tuple[float, float, bool]:
    """
    Analyse une liste de points GPS et détermine si la trajectoire
    correspond à un vol en de parapente.

    Paramètres
    ----------
    points : liste de GPSPoint (lat, lon, timestamp)
        Peut contenir n'importe quel nombre de points (pas forcément 600).

    Retourne
    --------
    R_kmh      : rayon du cercle de vitesse fitté (km/h)
    sigma_pct  : écart-type des résidus / rayon, en %
    is_para    : True si R_MIN_KMH ≤ R ≤ R_MAX_KMH et σ/R < MAX_RESIDUAL_PCT
    """
    vels = _velocity_vectors(points)
    if len(vels) < 10:
        return 0.0, 0.0, False

    # Filtrage IQR : on garde les vecteurs dont la norme est dans [Q25, Q75]
    speeds = np.array([math.hypot(vx, vy) for vx, vy in vels])
    q25, q75 = float(np.percentile(speeds, 25)), float(np.percentile(speeds, 75))
    vels_iqr = [(vx, vy) for (vx, vy), s in zip(vels, speeds) if q25 <= s <= q75]
    if len(vels_iqr) < 10:
        return 0.0, 0.0, False

    _, _, R, std = _taubin_fit(vels_iqr)
    R_kmh     = R * 3.6
    sigma_pct = (std / R * 100.0) if R > 0 else float("inf")

    speed_ok  = R_MIN_KMH <= R_kmh <= R_MAX_KMH
    circle_ok = std < MAX_RESIDUAL_PCT * R

    return R_kmh, sigma_pct, speed_ok and circle_ok


# ---------------------------------------------------------------------------
# Fonctions internes
# ---------------------------------------------------------------------------

_EARTH_M_DEG = 111_319.9   # mètres par degré de latitude
_MIN_SPEED   = 2.0          # vitesse sol minimale retenue (m/s)


def _velocity_vectors(points: list[GPSPoint]) -> list[tuple[float, float]]:
    """Calcule les vecteurs vitesse (Vx, Vy) en m/s entre points consécutifs."""
    vels = []
    for i in range(len(points) - 1):
        p0, p1 = points[i], points[i + 1]
        dt = p1.timestamp - p0.timestamp
        if dt <= 0.0:
            continue
        lat_mid = math.radians((p0.lat + p1.lat) / 2.0)
        dx = (p1.lon - p0.lon) * _EARTH_M_DEG * math.cos(lat_mid)
        dy = (p1.lat - p0.lat) * _EARTH_M_DEG
        vx, vy = dx / dt, dy / dt
        if math.hypot(vx, vy) >= _MIN_SPEED:
            vels.append((vx, vy))
    return vels


def _taubin_fit(vels: list[tuple[float, float]]) -> tuple[float, float, float, float]:
    """Fit de cercle par la méthode de Taubin. Retourne (cx, cy, R, std_residuals)."""
    pts  = np.asarray(vels, dtype=np.float64)
    mean = pts.mean(axis=0)
    u    = pts - mean
    ux, uy = u[:, 0], u[:, 1]
    z   = ux ** 2 + uy ** 2

    Mxx, Myy = np.mean(ux ** 2), np.mean(uy ** 2)
    Mxy       = np.mean(ux * uy)
    Mxz, Myz  = np.mean(ux * z), np.mean(uy * z)
    Mzz, Mz   = np.mean(z ** 2), np.mean(z)

    Cov_xy = Mxx * Myy - Mxy ** 2
    Var_z  = Mzz - Mz ** 2
    A2  = 4.0 * Cov_xy - 3.0 * Mz ** 2 - Mzz
    A1  = Var_z * Mz + 4.0 * Cov_xy * Mz - Mxz ** 2 - Myz ** 2
    A0  = (Mxz * (Mxz * Myy - Myz * Mxy)
           + Myz * (Myz * Mxx - Mxz * Mxy)
           - Var_z * Cov_xy)

    x = 0.0
    for _ in range(99):
        f  = A0 + x * (A1 + x * (A2 + 4.0 * x ** 2))
        Df = A1 + x * (2.0 * A2 + 16.0 * x ** 2)
        if abs(Df) < 1e-14:
            break
        x_new = x - f / Df
        if abs(x_new - x) < 1e-12:
            x = x_new
            break
        x = x_new

    DET = x ** 2 - x * Mz + Cov_xy
    if abs(DET) < 1e-10:
        return _kasa_fit(vels)

    Xc = (Mxz * (Myy - x) - Myz * Mxy) / (2.0 * DET)
    Yc = (Myz * (Mxx - x) - Mxz * Mxy) / (2.0 * DET)
    cx, cy = Xc + mean[0], Yc + mean[1]
    R = math.sqrt(Xc ** 2 + Yc ** 2 + Mz)

    dist = np.sqrt((pts[:, 0] - cx) ** 2 + (pts[:, 1] - cy) ** 2)
    return cx, cy, R, float(np.std(dist - R))


def _kasa_fit(vels: list[tuple[float, float]]) -> tuple[float, float, float, float]:
    """Fit de cercle par la méthode de Kåsa (repli si Taubin dégénère)."""
    pts = np.asarray(vels, dtype=np.float64)
    x, y = pts[:, 0], pts[:, 1]
    A = np.column_stack([2 * x, 2 * y, np.ones(len(x))])
    b = x ** 2 + y ** 2
    res, *_ = np.linalg.lstsq(A, b, rcond=None)
    cx, cy = res[0], res[1]
    R = math.sqrt(res[2] + cx ** 2 + cy ** 2)
    dist = np.sqrt((x - cx) ** 2 + (y - cy) ** 2)
    return cx, cy, R, float(np.std(dist - R))