MeteoThom
Tempête Austen, 22-24 octobre 2026

Le code

Les 18 fichiers C et PHP qui fabriquent la page de la tempête, soit 4 879 lignes, du fichier GRIB brut à l'image de la carte. Le calcul est en C (décodeur GRIB2 écrit à la main, suivi de la dépression, cartes), l'orchestration en PHP. Chaque bloc explique ce que fait le fichier, puis montre son code tel qu'il tourne sur le serveur, sous un compte système dédié sans aucun droit particulier.

Les chemins et comptes internes du serveur sont remplacés par des noms génériques (/chemin/vers/…). C11 compilé par GCC 12, PHP 8.4 ; coloration du code faite par l'analyseur de PHP et un découpage maison pour le C. Cette chaîne produit la page de la tempête depuis le jeudi 22 octobre à 17 h 53 (heure de Paris) ; elle remplace une première version en Python dont elle reproduit les résultats à l'identique.

[code 1]

Les réglages et les outils communs

lib.php, 380 lignes

La seule source des réglages : fenêtre de l'épisode (jeudi 20 h → samedi 2 h, heure de Paris), modèles, 38 stations, palette des rafales de Meteociel, et la trajectoire de référence, six positions d'Austen posées à la main qui servent uniquement à la reconnaître parmi les autres dépressions. Le binaire C reçoit ces réglages par un petit fichier texte écrit ici.

On y trouve aussi les règles de sécurité tenues partout : aucune commande passée par un shell (tableau d'arguments), téléchargements en HTTPS seulement avec certificat vérifié et taille plafonnée, jeton Météo-France gardé en mémoire et posé dans l'en-tête, écritures par fichier temporaire puis renommage.

lib.php
Fichier brut
<?php
/**
 * Tempête Austen : réglages et outils communs de la partie PHP (orchestration, téléchargements,
 * données de la page). Le calcul (GRIB2, analyse, cartes) est fait par le binaire C « austen ».
 *
 * Règles de sécurité tenues partout :
 *  - aucune commande passée par un shell : proc_open() reçoit un tableau d'arguments ;
 *  - téléchargements en HTTPS seulement, certificat vérifié, sans redirection, taille plafonnée ;
 *  - le jeton Météo-France est lu en mémoire et posé dans l'en-tête curl, jamais en argument
 *    ni sur disque ; aucun POST /token (il révoquerait le jeton de toutes les collectes) ;
 *  - un fichier d'un dossier inscriptible par le site est lu par « austen catsafe » (ni lien, ni FIFO) ;
 *  - toute écriture publiée passe par un fichier temporaire puis un renommage.
 */
declare(strict_types=1);
 
ini_set('memory_limit', '768M');
 
const HOME_DIR = __DIR__ . '/..';
const BIN = HOME_DIR . '/c/austen';
const UA = 'MeteoThom-Austen/2.0 (+https://meteothom.com)';
const MF_TOKEN = '/chemin/vers/cache-jeton-meteo-france.json';
const OBS_DIR = '/chemin/vers/observations-6min/2026/10';
const VIGILANCE = '/chemin/vers/vigilance.json';
const STATIONS_CSV = '/chemin/vers/stations.csv';
 
/* Dossiers : données de travail et dossier publié, réglables pour une passe à blanc. */
function data_dir(): string { return getenv('AUSTEN_DATA') ?: HOME_DIR . '/shadow/data'; }
function site_dir(): string { return getenv('AUSTEN_SITE') ?: HOME_DIR . '/shadow/site'; }
function pub_dir(): string { return getenv('AUSTEN_PUB') ?: HOME_DIR . '/shadow/pub'; }
 
function utc(string $iso): int
{
    $t = DateTimeImmutable::createFromFormat('!Y-m-d\TH:i\Z', $iso, new DateTimeZone('UTC'));
    if ($t === false) {
        throw new RuntimeException("date invalide : $iso");
    }
    return $t->getTimestamp();
}
 
function iso(int $t): string { return gmdate('Y-m-d\TH:i\Z', $t); }
 
/* ---------------------------------------------------------------- réglages de l'épisode */
function ep0(): int { return utc('2026-10-22T18:00Z'); }   /* jeudi 20 h, heure de Paris */
function ep1(): int { return utc('2026-10-24T00:00Z'); }   /* samedi 2 h */
function tr0(): int { return utc('2026-10-22T00:00Z'); }
function tr1(): int { return utc('2026-10-24T00:00Z'); }
function fetch_end(): int { return utc('2026-10-24T06:00Z'); }
 
const MODELS = [
    'arome'  => ['name' => 'AROME', 'res' => '1,3 km', 'color' => '#d6336c', 'gust' => true,
                 'gnote' => 'rafale maxi horaire', 'org' => 'Météo-France', 'cycle' => 3, 'delay' => 2.5],
    'arpege' => ['name' => 'ARPEGE', 'res' => '10 km', 'color' => '#1c7ed6', 'gust' => true,
                 'gnote' => 'rafale maxi horaire', 'org' => 'Météo-France', 'cycle' => 6, 'delay' => 3.5],
    'iconeu' => ['name' => 'ICON-EU', 'res' => '7 km', 'color' => '#2b8a3e', 'gust' => true,
                 'gnote' => 'rafale maxi horaire', 'org' => 'DWD', 'cycle' => 6, 'delay' => 2.5],
    'ifs'    => ['name' => 'CEP (IFS)', 'res' => '25 km', 'color' => '#e8590c', 'gust' => true,
                 'gnote' => 'rafale maxi de la dernière heure, toutes les 3 h', 'org' => 'ECMWF', 'cycle' => 6, 'delay' => 6.5],
    'gfs'    => ['name' => 'GFS', 'res' => '25 km', 'color' => '#7048e8', 'gust' => true,
                 'gnote' => 'rafale instantanée horaire', 'org' => 'NOAA', 'cycle' => 6, 'delay' => 3.5],
    'aifs'   => ['name' => 'AIFS', 'res' => '25 km', 'color' => '#868e96', 'gust' => false,
                 'gnote' => 'pas de rafale (pression seule)', 'org' => 'ECMWF, IA', 'cycle' => 6, 'delay' => 5.5],
];
 
/* Positions d'Austen posées à la main (consensus des runs du 22/10 00-06 UTC) : elles ne servent
 * qu'à reconnaître Austen parmi les autres dépressions au moment d'accrocher le suivi. */
const REF_TRACK = [
    ['2026-10-22T12:00Z', 48.6, -14.0], ['2026-10-22T18:00Z', 49.6, -9.8], ['2026-10-23T00:00Z', 49.9, -5.9],
    ['2026-10-23T06:00Z', 50.3, -1.7], ['2026-10-23T12:00Z', 51.2, 2.0], ['2026-10-23T18:00Z', 52.2, 5.0],
];
 
/* Stations Météo-France : identifiant, nom court, latitude, longitude. */
const STATIONS = [
    ['29155005', 'Ouessant', 48.473333, -5.057000], ['29168001', 'Pointe du Raz', 48.038833, -4.732000],
    ['29075001', 'Brest', 48.453833, -4.391167], ['29158001', 'Penmarch', 47.797500, -4.374833],
    ['29216001', 'Quimper', 47.973000, -4.160667], ['22113006', 'Lannion', 48.755333, -3.468667],
    ['22168001', "Ploumanac'h", 48.825833, -3.473167], ['56069001', 'Île de Groix', 47.652333, -3.502333],
    ['56185001', 'Lorient', 47.762833, -3.435667], ['56009001', 'Belle-Île', 47.294333, -3.218333],
    ['22016001', 'Île de Bréhat', 48.855167, -3.004500], ['22372001', 'Saint-Brieuc', 48.536500, -2.853167],
    ['56243001', 'Vannes', 47.604500, -2.714167], ['35228001', 'Dinard', 48.584833, -2.076333],
    ['44103001', 'Saint-Nazaire', 47.313833, -2.154500], ['35281001', 'Rennes', 48.068833, -1.734000],
    ['44020001', 'Nantes', 47.150000, -1.608833], ['50020001', 'La Hague', 49.725167, -1.939833],
    ['50218001', 'Granville', 48.834500, -1.613667], ['50196001', 'Gatteville', 49.697000, -1.265333],
    ['14137001', 'Caen', 49.180000, -0.456167], ['61001001', 'Alençon', 48.445500, 0.110167],
    ['72181001', 'Le Mans', 47.945833, 0.194167], ['76552001', 'Cap de la Hève', 49.509000, 0.071167],
    ['76217001', 'Dieppe', 49.932833, 1.090000], ['76116001', 'Rouen', 49.389500, 1.178333],
    ['62826001', 'Le Touquet', 50.514500, 1.622833], ['62160001', 'Boulogne', 50.732500, 1.599667],
    ['80001001', 'Abbeville', 50.136000, 1.834000], ['59183001', 'Dunkerque', 51.055833, 2.339500],
    ['59343001', 'Lille', 50.570000, 3.097500], ['60639001', 'Beauvais', 49.446500, 2.127167],
    ['95527001', 'Roissy', 49.015167, 2.534333], ['75114001', 'Paris', 48.821667, 2.337833],
    ['45055001', 'Orléans', 47.990667, 1.778167], ['37179001', 'Tours', 47.444500, 0.727333],
    ['17323001', 'Chassiron', 46.046833, -1.411667], ['51449002', 'Reims', 49.209833, 4.155333],
];
 
/* Échelle Meteociel des rafales (bornes basses des 30 cases) : légende de la page. */
const GUST_BOUNDS = [0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90,
                     100, 110, 120, 130, 140, 150, 160, 180, 200, 220, 240];
const GUST_COLORS = ['#ffffff', '#97e6ff', '#33ccff', '#0099ff', '#00ff99', '#33cc66', '#66cc33', '#66ff00',
                     '#a4f22f', '#b7cf0e', '#d6f017', '#cc9900', '#ff9900', '#ff9966', '#cc9999', '#cc6633',
                     '#cc3333', '#ff0d0d', '#c60000', '#800000', '#800050', '#a00077', '#cc00cc', '#ff00ff',
                     '#ff40ff', '#ff80ff', '#ffa0ff', '#ffffff', '#e3e3e3', '#c6c6c6'];
 
/* conf.tsv pour le binaire C : la seule source des réglages reste ce fichier PHP. */
function write_conf(): string
{
    $lines = ['ep0 ' . ep0(), 'ep1 ' . ep1(), 'tr0 ' . tr0(), 'tr1 ' . tr1()];
    foreach (REF_TRACK as [$t, $la, $lo]) {
        $lines[] = sprintf('ref %d %.4f %.4f', utc($t), $la, $lo);
    }
    foreach (STATIONS as [$id, , $la, $lo]) {
        $lines[] = sprintf('station %s %.6f %.6f', $id, $la, $lo);
    }
    $p = data_dir() . '/conf.tsv';
    atomic_write($p, implode("\n", $lines) . "\n");
    return $p;
}
 
/* ---------------------------------------------------------------- outils */
function logm(string ...$parts): void
{
    fwrite(STDOUT, gmdate('H:i:s') . ' ' . implode(' ', $parts) . "\n");
}
 
function ensure_dir(string $d): void
{
    if (!is_dir($d) && !mkdir($d, 0755, true) && !is_dir($d)) {
        throw new RuntimeException("création impossible : $d");
    }
}
 
function atomic_write(string $path, string $data): void
{
    ensure_dir(dirname($path));
    $tmp = dirname($path) . '/.' . basename($path) . '.new';
    if (is_link($tmp) || file_put_contents($tmp, $data, LOCK_EX) !== strlen($data)) {
        throw new RuntimeException("écriture impossible : $tmp");
    }
    chmod($tmp, 0644);
    if (!rename($tmp, $path)) {
        throw new RuntimeException("renommage impossible : $path");
    }
}
 
function json_out(string $path, mixed $data): void
{
    atomic_write($path, json_encode($data, JSON_UNESCAPED_SLASHES | JSON_UNESCAPED_UNICODE | JSON_PRESERVE_ZERO_FRACTION | JSON_THROW_ON_ERROR));
}
 
function json_in(string $path): mixed
{
    $s = file_get_contents($path);
    if ($s === false) {
        throw new RuntimeException("lecture impossible : $path");
    }
    return json_decode($s, true, 64, JSON_THROW_ON_ERROR);
}
 
/**
 * Lance une commande sans shell, avec délai maxi. Retourne [code, sortie standard, erreurs].
 */
function run(array $cmd, int $timeout = 900): array
{
    $p = proc_open($cmd, [0 => ['file', '/dev/null', 'r'], 1 => ['pipe', 'w'], 2 => ['pipe', 'w']], $pipes);
    if (!is_resource($p)) {
        throw new RuntimeException('lancement impossible : ' . $cmd[0]);
    }
    stream_set_blocking($pipes[1], false);
    stream_set_blocking($pipes[2], false);
    $out = $err = '';
    $deadline = microtime(true) + $timeout;
    while (true) {
        $r = [$pipes[1], $pipes[2]];
        $w = $e = null;
        if (stream_select($r, $w, $e, 1) === false) {
            break;
        }
        foreach ($r as $s) {
            $chunk = fread($s, 1 << 16);
            if ($chunk !== false) {
                if ($s === $pipes[1]) { $out .= $chunk; } else { $err .= $chunk; }
            }
        }
        if (feof($pipes[1]) && feof($pipes[2])) {
            break;
        }
        if (microtime(true) > $deadline) {
            proc_terminate($p, 9);
            $err .= "\ndélai dépassé";
            break;
        }
    }
    fclose($pipes[1]);
    fclose($pipes[2]);
    $code = proc_close($p);
    return [$code, $out, $err];
}
 
/* Lecture d'un fichier d'un dossier inscriptible par le site : par le binaire C (ni lien, ni FIFO). */
function read_untrusted(string $path, int $max): ?string
{
    [$code, $out] = run([BIN, 'catsafe', $path, (string) $max], 60);
    return $code === 0 ? $out : null;
}
 
/* ---------------------------------------------------------------- HTTP */
function curl_handle(string $url, array $headers = [], ?string $range = null, int $max = 1 << 30): CurlHandle
{
    if (!str_starts_with($url, 'https://')) {
        throw new RuntimeException("URL refusée : $url");
    }
    $h = curl_init($url);
    $buf = '';
    curl_setopt_array($h, [
        CURLOPT_PROTOCOLS_STR => 'https',
        CURLOPT_REDIR_PROTOCOLS_STR => 'https',
        CURLOPT_FOLLOWLOCATION => false,
        CURLOPT_SSL_VERIFYPEER => true,
        CURLOPT_SSL_VERIFYHOST => 2,
        CURLOPT_CONNECTTIMEOUT => 30,
        CURLOPT_TIMEOUT => 900,
        CURLOPT_USERAGENT => UA,
        CURLOPT_HTTPHEADER => $headers,
        CURLOPT_RETURNTRANSFER => false,
        CURLOPT_WRITEFUNCTION => function ($ch, string $data) use (&$buf, $max): int {
            if (strlen($buf) + strlen($data) > $max) {
                return 0;   /* taille maxi dépassée : curl abandonne */
            }
            $buf .= $data;
            return strlen($data);
        },
    ]);
    if ($range !== null) {
        curl_setopt($h, CURLOPT_RANGE, $range);
    }
    $GLOBALS['curl_buffers'][spl_object_id($h)] = &$buf;
    return $h;
}
 
function curl_body(CurlHandle $h): string
{
    $id = spl_object_id($h);
    $b = $GLOBALS['curl_buffers'][$id] ?? '';
    unset($GLOBALS['curl_buffers'][$id]);
    return $b;
}
 
/**
 * Téléchargement direct sur disque (paquets Météo-France de 50 à 80 Mo) : fichier temporaire caché,
 * taille plafonnée, renommé seulement si le code HTTP vaut 200. Retourne le code HTTP (0 = réseau).
 */
function http_download(string $url, array $headers, string $path, int $max): int
{
    if (!str_starts_with($url, 'https://')) {
        throw new RuntimeException("URL refusée : $url");
    }
    $tmp = dirname($path) . '/.' . basename($path) . '.part';
    $code = 0;
    for ($k = 0; $k < 3; $k++) {
        $fh = fopen($tmp, 'wb');
        if ($fh === false) {
            throw new RuntimeException("écriture impossible : $tmp");
        }
        $h = curl_init($url);
        curl_setopt_array($h, [
            CURLOPT_PROTOCOLS_STR => 'https',
            CURLOPT_FOLLOWLOCATION => false,
            CURLOPT_SSL_VERIFYPEER => true,
            CURLOPT_SSL_VERIFYHOST => 2,
            CURLOPT_CONNECTTIMEOUT => 30,
            CURLOPT_TIMEOUT => 900,
            CURLOPT_USERAGENT => UA,
            CURLOPT_HTTPHEADER => $headers,
            CURLOPT_FILE => $fh,
            CURLOPT_MAXFILESIZE_LARGE => $max,
            CURLOPT_NOPROGRESS => false,
            CURLOPT_XFERINFOFUNCTION => fn($c, $dlt, $dl) => $dl > $max ? 1 : 0,
        ]);
        curl_exec($h);
        $code = (int) curl_getinfo($h, CURLINFO_RESPONSE_CODE);
        fclose($fh);
        if ($code === 200) {
            chmod($tmp, 0644);
            rename($tmp, $path);
            return 200;
        }
        @unlink($tmp);
        if ($code !== 0 && $code < 500) {
            return $code;
        }
        sleep(20);
    }
    return $code;
}
 
/** GET simple : [code HTTP, corps]. 0 = échec réseau. Quatre essais sur erreur réseau ou 5xx. */
function http_get(string $url, array $headers = [], ?string $range = null, int $max = 1 << 30): array
{
    for ($k = 0; $k < 4; $k++) {
        $h = curl_handle($url, $headers, $range, $max);
        curl_exec($h);
        $code = (int) curl_getinfo($h, CURLINFO_RESPONSE_CODE);
        $body = curl_body($h);
        if ($code === 200 || $code === 206 || ($code >= 400 && $code < 500)) {
            return [$code, $body];
        }
        sleep(5 * ($k + 1));
    }
    return [$code ?? 0, ''];
}
 
/**
 * Téléchargements en parallèle. $jobs : clé => [url, range|null]. Retourne clé => [code, corps].
 */
function http_multi(array $jobs, int $parallel = 8, int $max = 256 << 20): array
{
    $res = [];
    $queue = $jobs;
    $mh = curl_multi_init();
    $active = [];
    $tries = [];
    $start = function () use (&$queue, &$active, $mh, $max, &$tries): void {
        $key = array_key_first($queue);
        [$url, $range] = $queue[$key];
        unset($queue[$key]);
        $h = curl_handle($url, [], $range, $max);
        curl_multi_add_handle($mh, $h);
        $active[spl_object_id($h)] = [$key, $h, $url, $range];
        $tries[$key] = ($tries[$key] ?? 0) + 1;
    };
    while ($queue && count($active) < $parallel) {
        $start();
    }
    do {
        curl_multi_exec($mh, $running);
        curl_multi_select($mh, 1.0);
        while ($info = curl_multi_info_read($mh)) {
            $h = $info['handle'];
            [$key, , $url, $range] = $active[spl_object_id($h)];
            unset($active[spl_object_id($h)]);
            $code = (int) curl_getinfo($h, CURLINFO_RESPONSE_CODE);
            $body = curl_body($h);
            curl_multi_remove_handle($mh, $h);
            if (($code === 0 || $code >= 500) && $tries[$key] < 4) {
                usleep(2_000_000 * $tries[$key]);   /* serveur en difficulté : on ne le martèle pas */
                $queue[$key] = [$url, $range];
            } else {
                $res[$key] = [$code, $body];
            }
            if ($queue) {
                $start();
            }
        }
    } while ($active || $queue);
    curl_multi_close($mh);
    return $res;
}
 
/* Jeton Météo-France partagé : lu sous le compte des collectes, gardé en mémoire, validé. */
function mf_token(): string
{
    [$code, $out] = run(['<lecture-déléguée>', '/usr/bin/cat', MF_TOKEN], 30);
    $j = $code === 0 ? json_decode($out, true) : null;
    $t = is_array($j) ? ($j['access_token'] ?? '') : '';
    if (!is_string($t) || !preg_match('/^[A-Za-z0-9._-]{20,4096}$/', $t)) {
        throw new RuntimeException('jeton Météo-France illisible');
    }
    return $t;
}
 
function valid_run(string $run): bool
{
    return (bool) preg_match('/^20\d{2}(0[1-9]|1[0-2])(0[1-9]|[12]\d|3[01])([01]\d|2[0-3])$/', $run);
}
 
function run_time(string $run): int
{
    if (!valid_run($run)) {
        throw new RuntimeException("run invalide : $run");
    }
    return gmmktime((int) substr($run, 8, 2), 0, 0, (int) substr($run, 4, 2), (int) substr($run, 6, 2), (int) substr($run, 0, 4));
}
[code 2]

Le téléchargement des modèles

fetch.php, 271 lignes

Pour chaque modèle, le run le plus récent, ou le précédent s'il n'est pas complet. Le Centre européen et la NOAA publient un index de chaque fichier : seuls les octets de la pression et de la rafale sont demandés, en parallèle (curl multi). ICON-EU arrive en petits fichiers bzip2, décompressés avec une taille maximale ; AROME et ARPEGE en paquets du portail Météo-France, écrits directement sur disque.

fetch.php
Fichier brut
<?php
/**
 * Téléchargement des sorties modèles : fetch.php <modèle> [AAAAMMJJHH]
 * Sans run : le plus récent disponible et complet. Écrit les GRIB bruts dans <données>/raw/<modèle>/<run>/
 * et un fichier .complete quand le run est entier. Affiche « RUN <run> » en cas de succès.
 */
declare(strict_types=1);
require __DIR__ . '/lib.php';
 
const AROME_SLICES = ['00H06H', '07H12H', '13H18H', '19H24H', '25H30H', '31H36H', '37H42H', '43H48H', '49H51H'];
const ARPEGE_SLICES = ['000H012H', '013H024H', '025H036H', '037H048H', '049H060H', '061H072H'];
 
function raw_dir(string $model, string $run): string
{
    $d = data_dir() . "/raw/$model/$run";
    ensure_dir($d);
    return $d;
}
 
/** Heures d'échéance utiles : de 0 à la fin de la fenêtre, par pas de $step, plafonnées. */
function hours_needed(int $run, int $step, ?int $cap = null): array
{
    $h = intdiv(fetch_end() - $run, 3600);
    if ($cap !== null) {
        $h = min($h, $cap);
    }
    return $h >= 0 ? range(0, $h, $step) : [];
}
 
/* ---------------------------------------------------------------- Météo-France */
function fetch_mf(string $model, int $run): bool
{
    $base = 'https://public-api.meteofrance.fr/previnum';
    [$slices, $url] = $model === 'arome'
        ? [AROME_SLICES, "$base/DPPaquetAROME/v1/models/AROME/grids/0.025/packages/SP1/productARO"]
        : [ARPEGE_SLICES, "$base/DPPaquetARPEGE/v1/models/ARPEGE/grids/0.1/packages/SP1/productARP"];
    $maxh = intdiv(fetch_end() - $run, 3600);
    $dir = raw_dir($model, gmdate('YmdH', $run));
    $token = mf_token();
    foreach ($slices as $s) {
        [$a] = explode('H', $s);
        if ((int) $a > $maxh) {
            break;
        }
        $out = "$dir/SP1_$s.grib2";
        if (is_file($out)) {
            continue;
        }
        $q = '?referencetime=' . gmdate('Y-m-d\TH:00:00\Z', $run) . "&time=$s&format=grib2";
        $code = http_download($url . $q, ["Authorization: Bearer $token", 'Accept: */*'], $out, 1 << 30);
        if ($code === 401) {   /* jeton renouvelé entre-temps par le collecteur : on relit le cache */
            sleep(30);
            $token = mf_token();
            $code = http_download($url . $q, ["Authorization: Bearer $token", 'Accept: */*'], $out, 1 << 30);
        }
        logm($model, gmdate('d/H', $run), $s, (string) $code);
        if ($code !== 200) {
            /* la dernière tranche (49H51H) n'existe pas pour tous les runs : run accepté sans elle */
            return $code === 404 && $s === end($slices);
        }
        sleep(1);   /* quota du portail : 100 requêtes par minute, partagées avec les autres collectes */
    }
    return true;
}
 
/* ---------------------------------------------------------------- ECMWF (IFS, AIFS) */
function fetch_ecmwf(string $model, int $run): bool
{
    $stream = $model === 'ifs' ? 'ifs' : 'aifs-single';
    $params = $model === 'ifs' ? ['msl', '10fg'] : ['msl'];
    $step = $model === 'ifs' ? 3 : 6;
    $hh = gmdate('H', $run);
    $ymd = gmdate('Ymd', $run);
    $dir = raw_dir($model, gmdate('YmdH', $run));
    $todo = [];
    foreach (hours_needed($run, $step, in_array($hh, ['06', '18'], true) ? 90 : null) as $h) {
        if (!is_file(sprintf('%s/f%03d.grib2', $dir, $h))) {
            $todo[$h] = sprintf('https://data.ecmwf.int/forecasts/%s/%sz/%s/0p25/oper/%s%s0000-%dh-oper-fc', $ymd, $hh, $stream, $ymd, $hh, $h);
        }
    }
    /* 1. index de chaque échéance (une ligne JSON par champ, avec position et longueur en octets) */
    $idx = http_multi(array_map(fn($b) => [$b . '.index', null], $todo), 6, 4 << 20);
    $ranges = [];
    foreach ($todo as $h => $b) {
        [$code, $body] = $idx[$h];
        if ($code !== 200) {
            return false;
        }
        foreach (explode("\n", trim($body)) as $line) {
            $j = json_decode($line, true);
            if (is_array($j) && in_array($j['param'] ?? '', $params, true) && ($j['levtype'] ?? '') === 'sfc') {
                $o = (int) $j['_offset'];
                $ranges["$h:$o"] = [$b . '.grib2', $o . '-' . ($o + (int) $j['_length'] - 1)];
            }
        }
    }
    /* 2. seulement les octets des champs utiles */
    $parts = http_multi($ranges, 8, 64 << 20);
    $byHour = [];
    ksort($parts, SORT_NATURAL);
    foreach ($parts as $k => [$code, $body]) {
        if ($code !== 206 && $code !== 200) {
            return false;
        }
        $byHour[(int) explode(':', $k)[0]][] = $body;
    }
    foreach ($byHour as $h => $bodies) {
        atomic_write(sprintf('%s/f%03d.grib2', $dir, $h), implode('', $bodies));
    }
    logm($model, gmdate('d/H', $run), 'ok', (string) count($todo), 'échéances');
    return true;
}
 
/* ---------------------------------------------------------------- GFS (NOAA, S3) */
function fetch_gfs(string $model, int $run): bool
{
    $hh = gmdate('H', $run);
    $ymd = gmdate('Ymd', $run);
    $dir = raw_dir($model, gmdate('YmdH', $run));
    $todo = [];
    foreach (hours_needed($run, 1, 120) as $h) {
        if (!is_file(sprintf('%s/f%03d.grib2', $dir, $h))) {
            $todo[$h] = sprintf('https://noaa-gfs-bdp-pds.s3.amazonaws.com/gfs.%s/%s/atmos/gfs.t%sz.pgrb2.0p25.f%03d', $ymd, $hh, $hh, $h);
        }
    }
    $idx = http_multi(array_map(fn($u) => [$u . '.idx', null], $todo), 8, 4 << 20);
    $ranges = [];
    foreach ($todo as $h => $u) {
        [$code, $body] = $idx[$h];
        if ($code !== 200) {
            return false;
        }
        /* .idx : « n:octet:date:VAR:niveau:... » ; un champ s'arrête où commence le suivant */
        $lines = explode("\n", trim($body));
        foreach ($lines as $i => $line) {
            if ((str_contains($line, ':PRMSL:mean sea level:') || str_contains($line, ':GUST:surface:')) && isset($lines[$i + 1])) {
                $a = (int) explode(':', $line)[1];
                $b = (int) explode(':', $lines[$i + 1])[1] - 1;
                $ranges["$h:$a"] = [$u, "$a-$b"];
            }
        }
    }
    $parts = http_multi($ranges, 8, 64 << 20);
    ksort($parts, SORT_NATURAL);
    $byHour = [];
    foreach ($parts as $k => [$code, $body]) {
        if ($code !== 206 && $code !== 200) {
            return false;
        }
        $byHour[(int) explode(':', $k)[0]][] = $body;
    }
    foreach ($byHour as $h => $bodies) {
        atomic_write(sprintf('%s/f%03d.grib2', $dir, $h), implode('', $bodies));
    }
    logm($model, gmdate('d/H', $run), 'ok', (string) count($todo), 'échéances');
    return true;
}
 
/* ---------------------------------------------------------------- ICON-EU (DWD) */
function fetch_iconeu(string $model, int $run): bool
{
    $hh = gmdate('H', $run);
    $ymdh = gmdate('YmdH', $run);
    $dir = raw_dir($model, $ymdh);
    $jobs = [];
    foreach (hours_needed($run, 1, 78) as $h) {
        foreach (['pmsl', 'vmax_10m'] as $v) {
            $out = sprintf('%s/%s_%03d.grib2', $dir, $v, $h);
            if (!is_file($out)) {
                $jobs[$out] = [sprintf('https://opendata.dwd.de/weather/nwp/icon-eu/grib/%s/%s/icon-eu_europe_regular-lat-lon_single-level_%s_%03d_%s.grib2.bz2',
                    $hh, $v, $ymdh, $h, strtoupper($v)), null];
            }
        }
    }
    foreach (http_multi($jobs, 8, 64 << 20) as $out => [$code, $body]) {
        if ($code !== 200) {
            return false;
        }
        $raw = bzdecompress_limited($body, 256 << 20);
        if ($raw === null) {
            return false;
        }
        atomic_write($out, $raw);
    }
    logm($model, gmdate('d/H', $run), 'ok', (string) count($jobs), 'fichiers');
    return true;
}
 
/* décompression bzip2 bornée (l'extension bz2 n'étant pas compilée, on passe par bzip2 -dc) */
function bzdecompress_limited(string $data, int $max): ?string
{
    $p = proc_open(['bzip2', '-dc'], [0 => ['pipe', 'r'], 1 => ['pipe', 'w'], 2 => ['file', '/dev/null', 'w']], $pipes);
    if (!is_resource($p)) {
        return null;
    }
    stream_set_blocking($pipes[0], false);   /* jamais bloqué en écriture : on lit la sortie entre deux envois */
    stream_set_blocking($pipes[1], false);
    $out = '';
    $in = $data;
    while (true) {
        if ($in !== '') {
            $w = fwrite($pipes[0], substr($in, 0, 1 << 16));
            $in = $w ? substr($in, $w) : $in;
            if ($in === '') {
                fclose($pipes[0]);
            }
        }
        $chunk = fread($pipes[1], 1 << 16);
        if ($chunk !== false && $chunk !== '') {
            $out .= $chunk;
            if (strlen($out) > $max) {
                proc_terminate($p, 9);
                return null;
            }
        }
        if ($in === '' && feof($pipes[1])) {
            break;
        }
        if ($chunk === '' || $chunk === false) {
            usleep(1000);
        }
    }
    fclose($pipes[1]);
    return proc_close($p) === 0 ? $out : null;
}
 
/* ---------------------------------------------------------------- choix du run */
function candidates(string $model): array
{
    $c = MODELS[$model]['cycle'];
    $t = time() - (int) (MODELS[$model]['delay'] * 3600);
    $t -= $t % ($c * 3600);
    return array_map(fn($k) => $t - $k * $c * 3600, range(0, 3));
}
 
function fetch_run(string $model, int $run): bool
{
    return match ($model) {
        'arome', 'arpege' => fetch_mf($model, $run),
        'ifs', 'aifs' => fetch_ecmwf($model, $run),
        'gfs' => fetch_gfs($model, $run),
        'iconeu' => fetch_iconeu($model, $run),
    };
}
 
function main(array $argv): int
{
    $model = $argv[1] ?? '';
    if (!isset(MODELS[$model])) {
        fwrite(STDERR, "usage : fetch.php <modèle> [AAAAMMJJHH]\n");
        return 2;
    }
    $runs = isset($argv[2]) ? [run_time($argv[2])] : candidates($model);
    foreach ($runs as $run) {
        $d = data_dir() . "/raw/$model/" . gmdate('YmdH', $run);
        $done = data_dir() . "/derived/$model/" . gmdate('YmdH', $run) . '.json';
        if (is_file("$d/.complete") || is_file("$d/.decoded") || is_file($done)) {
            echo 'RUN ', gmdate('YmdH', $run), "\n";
            return 0;
        }
        if (fetch_run($model, $run)) {
            touch("$d/.complete");
            echo 'RUN ', gmdate('YmdH', $run), "\n";
            return 0;
        }
        logm($model, gmdate('YmdH', $run), 'indisponible, run précédent');
    }
    return 1;
}
 
exit(main($argv));
[code 3]

Le décodeur GRIB2 : ce qu'il sait lire

grib2.h, 35 lignes

Les modèles livrent leurs champs au format GRIB2 de l'OMM. Le décodeur est écrit à la main et ne gère que ce dont la chaîne a besoin : grille latitude-longitude régulière, produits instantanés ou statistiques (rafale maxi), bitmap des points sans donnée, emballages simple, complexe avec différenciation spatiale (GFS) et CCSDS (tous les autres). Tout le reste est refusé proprement.

grib2.h
Fichier brut
/* grib2.h : décodeur GRIB édition 2, limité à ce dont la chaîne a besoin.
 *
 * Gère : grille 3.0 (latitude-longitude régulière), produits 4.0 et 4.8, section 6 (bitmap),
 * emballages 5.0 (simple), 5.3 (complexe avec différenciation spatiale, GFS) et 5.42 (CCSDS,
 * décompressé par libaec). Tout le reste est refusé proprement.
 * Chaque longueur lue dans le fichier est bornée par ce qui reste du tampon.
 */
#ifndef AUSTEN_GRIB2_H
#define AUSTEN_GRIB2_H
 
#include <stddef.h>
#include <stdint.h>
 
typedef struct {
    int discipline, category, number;   /* identité du paramètre (tables OMM) */
    int pdt;                            /* modèle de produit : 0 ou 8 */
    int surface, stat_proc;             /* type de surface fixe ; traitement statistique (-1 si aucun) */
    double surface_value;
    int64_t ref_time, valid_time;       /* epoch UTC ; validité = fin d'intervalle pour 4.8 */
    int ni, nj, scan;                   /* points en i (longitude) et j (latitude), mode de balayage */
    double la1, lo1, la2, lo2, di, dj;  /* degrés */
    int drt, bits;                      /* modèle de représentation des données, bits par valeur */
    double *values;                     /* ni*nj dans l'ordre du fichier, NaN = manquant */
} G2Field;
 
/* Décode le message qui commence à *pos (ou le suivant). Retourne 1 si un champ a été lu,
 * 0 en fin de tampon, -1 en cas d'erreur (message dans err). *pos avance au message suivant.
 * want(discipline, category, number, ctx) permet de ne dépaqueter que les champs utiles
 * (NULL = tous) ; un champ ignoré est rendu avec values = NULL. */
typedef int (*g2_filter)(int discipline, int category, int number, void *ctx);
int g2_next(const uint8_t *buf, size_t len, size_t *pos, G2Field *f, g2_filter want, void *ctx,
            char *err, size_t errlen);
void g2_free(G2Field *f);
 
#endif
[code 4]

Le décodeur GRIB2 : le code

grib2.c, 478 lignes

Lecture des sections 0 à 8, entiers gros-boutistes en signe et magnitude (pas en complément à deux), date de validité lue en clair pour une rafale maxi, dépaquetage (R + X × 2^E) / 10^D. L'emballage complexe du GFS est déroulé à la main : valeurs par groupes (référence, largeur, longueur de chaque groupe), puis différenciation spatiale d'ordre 1 ou 2 à défaire. Seule la décompression CCSDS passe par libaec, l'implémentation de référence de ce codec.

Vérification : sur 409 messages des 7 modèles, valeurs et en-têtes identiques à ceux d'eccodes, la bibliothèque du Centre européen (écart nul), compilé avec AddressSanitizer et UndefinedBehaviorSanitizer. Puis environ 7 800 fichiers corrompus au hasard : 5 900 refus propres, aucun plantage. Chaque longueur lue dans le fichier est bornée par ce qui reste du tampon.

grib2.c
Fichier brut
/* grib2.c : décodeur GRIB2 (voir grib2.h). Références : manuel des codes OMM n° 306,
 * FM 92 GRIB édition 2, sections 0 à 8 et modèles 3.0, 4.0, 4.8, 5.0, 5.3, 5.42, 6.0, 7.x. */
#include "grib2.h"
 
#include <libaec.h>
#include <math.h>
#include <stdarg.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
 
#include "util.h"
 
/* --------------------------------------------------------------- lecture d'octets */
/* Les entiers GRIB2 sont gros-boutistes ; les entiers signés sont en signe + magnitude
 * (bit de poids fort = signe), pas en complément à deux. */
static uint32_t u8(const uint8_t *p) { return p[0]; }
static uint32_t u16(const uint8_t *p) { return (uint32_t)p[0] << 8 | p[1]; }
static uint32_t u32(const uint8_t *p) { return (uint32_t)p[0] << 24 | (uint32_t)p[1] << 16 | (uint32_t)p[2] << 8 | p[3]; }
static uint64_t u64(const uint8_t *p) { return (uint64_t)u32(p) << 32 | u32(p + 4); }
static int32_t s16(const uint8_t *p) { int32_t v = (int32_t)(u16(p) & 0x7fff); return (p[0] & 0x80) ? -v : v; }
static int32_t s32(const uint8_t *p) { int32_t v = (int32_t)(u32(p) & 0x7fffffffu); return (p[0] & 0x80) ? -v : v; }
static float f32(const uint8_t *p) { uint32_t b = u32(p); float f; memcpy(&f, &b, 4); return f; }
 
/* entier signé (signe + magnitude) sur n octets, pour les descripteurs du 5.3 */
static int64_t sn(const uint8_t *p, int n)
{
    uint64_t v = 0;
    for (int k = 0; k < n; k++)
        v = v << 8 | p[k];
    uint64_t sign = (uint64_t)1 << (8 * n - 1);
    return (v & sign) ? -(int64_t)(v & (sign - 1)) : (int64_t)v;
}
 
/* --------------------------------------------------------------- lecteur de bits */
typedef struct {
    const uint8_t *p;
    size_t nbits, pos;
} Bits;
 
/* lit n bits (0..32) ; retourne -1 si le tampon est dépassé */
static int bits_get(Bits *b, unsigned n, uint32_t *out)
{
    if (n > 32 || b->pos + n > b->nbits)
        return -1;
    uint64_t v = 0;
    size_t pos = b->pos;
    unsigned left = n;
    while (left) {
        unsigned off = (unsigned)(pos & 7), take = 8 - off;
        if (take > left)
            take = left;
        unsigned c = (b->p[pos >> 3] >> (8 - off - take)) & ((1u << take) - 1);
        v = v << take | c;
        pos += take;
        left -= take;
    }
    b->pos = pos;
    *out = (uint32_t)v;
    return 0;
}
 
static void bits_align(Bits *b) { b->pos = (b->pos + 7) & ~(size_t)7; }
 
static int fail(char *err, size_t errlen, const char *fmt, ...) __attribute__((format(printf, 3, 4)));
static int fail(char *err, size_t errlen, const char *fmt, ...)
{
    va_list ap;
    va_start(ap, fmt);
    vsnprintf(err, errlen, fmt, ap);
    va_end(ap);
    return -1;
}
 
/* durée selon l'unité de la table 4.4, en secondes */
static int64_t unit_seconds(int unit)
{
    switch (unit) {
    case 0: return 60;
    case 1: return 3600;
    case 2: return 86400;
    case 10: return 3 * 3600;
    case 11: return 6 * 3600;
    case 12: return 12 * 3600;
    case 13: return 1;
    default: return -1;
    }
}
 
/* --------------------------------------------------------------- dépaquetage */
typedef struct {
    float R;
    int E, D, bits, drt;
    size_t ndata;               /* valeurs codées (sans les points masqués par le bitmap) */
    /* 5.42 */
    int ccsds_flags, block_size, rsi;
    /* 5.3 */
    int group_split, missing_mgmt, ng, ref_width, bits_width, len_incr, bits_len, order, extra;
    uint32_t ref_len, last_len;
} Pack;
 
static int unpack_simple(const Pack *k, const uint8_t *d, size_t dlen, uint32_t *x, char *err, size_t el)
{
    if (k->bits == 0) {
        memset(x, 0, k->ndata * sizeof *x);
        return 0;
    }
    if (k->bits > 32)
        return fail(err, el, "5.0 : %d bits par valeur", k->bits);
    Bits b = {d, dlen * 8, 0};
    for (size_t n = 0; n < k->ndata; n++)
        if (bits_get(&b, (unsigned)k->bits, &x[n]))
            return fail(err, el, "5.0 : données tronquées");
    return 0;
}
 
static int unpack_ccsds(const Pack *k, const uint8_t *d, size_t dlen, uint32_t *x, char *err, size_t el)
{
    if (k->bits == 0) {
        memset(x, 0, k->ndata * sizeof *x);
        return 0;
    }
    if (k->bits > 32 || k->block_size <= 0 || k->rsi <= 0)
        return fail(err, el, "5.42 : paramètres invalides");
    /* échantillons en sortie : 1, 2 ou 4 octets, ordre natif (on retire les drapeaux MSB et 3BYTE) */
    size_t nb = (size_t)(k->bits + 7) / 8;
    if (nb == 3)
        nb = 4;
    if (k->ndata > SIZE_MAX / nb)
        return fail(err, el, "5.42 : taille");
    uint8_t *out = xmalloc(k->ndata, nb);
    struct aec_stream s;
    memset(&s, 0, sizeof s);
    s.next_in = d;
    s.avail_in = dlen;
    s.next_out = out;
    s.avail_out = k->ndata * nb;
    s.bits_per_sample = (unsigned)k->bits;
    s.block_size = (unsigned)k->block_size;
    s.rsi = (unsigned)k->rsi;
    s.flags = (unsigned)k->ccsds_flags & ~(unsigned)(AEC_DATA_3BYTE | AEC_DATA_MSB);
    int rc = aec_buffer_decode(&s);
    if (rc != AEC_OK || s.total_out != k->ndata * nb) {
        free(out);
        return fail(err, el, "5.42 : décompression CCSDS en échec (%d)", rc);
    }
    for (size_t n = 0; n < k->ndata; n++) {
        if (nb == 1)
            x[n] = out[n];
        else if (nb == 2) {
            uint16_t v;
            memcpy(&v, out + 2 * n, 2);
            x[n] = v;
        } else
            memcpy(&x[n], out + 4 * n, 4);
    }
    free(out);
    return 0;
}
 
/* 5.3 : valeurs en groupes (référence, largeur, longueur par groupe), puis différenciation
 * spatiale d'ordre 1 ou 2 à défaire sur la suite des valeurs non manquantes.
 * x reçoit les entiers (après différenciation), miss[n] = 1 si la valeur est manquante. */
static int unpack_complex(const Pack *k, const uint8_t *d, size_t dlen, int64_t *x, uint8_t *miss, char *err, size_t el)
{
    if (k->order < 1 || k->order > 2 || k->extra < 1 || k->extra > 4)
        return fail(err, el, "5.3 : différenciation non gérée");
    if (k->ng <= 0 || (size_t)k->ng > k->ndata || k->bits > 32 || k->bits_width > 32 || k->bits_len > 32)
        return fail(err, el, "5.3 : groupes invalides");
    size_t need = (size_t)k->extra * (size_t)(k->order + 1);
    if (dlen < need)
        return fail(err, el, "5.3 : descripteurs tronqués");
    int64_t ival1 = sn(d, k->extra), ival2 = k->order == 2 ? sn(d + k->extra, k->extra) : 0;
    int64_t minsd = sn(d + (size_t)k->extra * (size_t)k->order, k->extra);
    Bits b = {d + need, (dlen - need) * 8, 0};
    size_t ng = (size_t)k->ng;
    uint32_t *gref = xmalloc(ng, 4), *gw = xmalloc(ng, 4), *gl = xmalloc(ng, 4);
    int rc = -1;
    for (size_t g = 0; g < ng; g++)
        if (bits_get(&b, (unsigned)k->bits, &gref[g]))
            goto out;
    bits_align(&b);
    for (size_t g = 0; g < ng; g++) {
        if (bits_get(&b, (unsigned)k->bits_width, &gw[g]))
            goto out;
        gw[g] += (uint32_t)k->ref_width;
        if (gw[g] > 32)
            goto out;
    }
    bits_align(&b);
    size_t total = 0;
    for (size_t g = 0; g < ng; g++) {
        uint32_t v;
        if (bits_get(&b, (unsigned)k->bits_len, &v))
            goto out;
        gl[g] = g + 1 == ng ? k->last_len : k->ref_len + v * (uint32_t)k->len_incr;
        total += gl[g];
        if (total > k->ndata)
            goto out;
    }
    if (total != k->ndata)
        goto out;
    bits_align(&b);
    uint32_t allbits = k->bits ? (uint32_t)((1ull << k->bits) - 1) : 0;
    size_t n = 0;
    for (size_t g = 0; g < ng; g++) {
        uint32_t w = gw[g], wall = w ? (uint32_t)((1ull << w) - 1) : 0;
        for (uint32_t m = 0; m < gl[g]; m++, n++) {
            uint32_t v = 0;
            miss[n] = 0;
            if (w) {
                if (bits_get(&b, w, &v))
                    goto out;
                if (k->missing_mgmt >= 1 && v == wall)
                    miss[n] = 1;
                else if (k->missing_mgmt == 2 && v == wall - 1)
                    miss[n] = 1;
                x[n] = (int64_t)v + gref[g];
            } else {
                if (k->missing_mgmt >= 1 && gref[g] == allbits)
                    miss[n] = 1;
                else if (k->missing_mgmt == 2 && gref[g] == allbits - 1)
                    miss[n] = 1;
                x[n] = gref[g];
            }
        }
    }
    /* différenciation spatiale, sur les seules valeurs présentes */
    int64_t p1 = 0, p2 = 0;
    size_t seen = 0;
    for (n = 0; n < k->ndata; n++) {
        if (miss[n])
            continue;
        int64_t v;
        if (seen == 0)
            v = ival1;
        else if (seen == 1 && k->order == 2)
            v = ival2;
        else if (k->order == 1)
            v = x[n] + minsd + p1;
        else
            v = x[n] + minsd + 2 * p1 - p2;
        p2 = p1;
        p1 = v;
        x[n] = v;
        seen++;
    }
    rc = 0;
out:
    free(gref);
    free(gw);
    free(gl);
    if (rc)
        fail(err, el, "5.3 : données tronquées ou incohérentes");
    return rc;
}
 
/* --------------------------------------------------------------- message */
void g2_free(G2Field *f)
{
    free(f->values);
    f->values = NULL;
}
 
int g2_next(const uint8_t *buf, size_t len, size_t *pos, G2Field *f, g2_filter want, void *ctx,
            char *err, size_t el)
{
    memset(f, 0, sizeof *f);
    f->stat_proc = -1;
    /* recherche de "GRIB" (des octets parasites entre messages sont tolérés) */
    size_t p = *pos;
    while (p + 16 <= len && memcmp(buf + p, "GRIB", 4) != 0)
        p++;
    if (p + 16 > len) {
        *pos = len;
        return 0;
    }
    const uint8_t *m = buf + p;
    if (m[7] != 2)
        return *pos = p + 4, fail(err, el, "édition %u non gérée", m[7]);
    uint64_t mlen = u64(m + 8);
    if (mlen < 16 + 4 || mlen > len - p)
        return *pos = len, fail(err, el, "longueur de message invalide");
    *pos = p + (size_t)mlen;
    if (memcmp(m + mlen - 4, "7777", 4) != 0)
        return fail(err, el, "fin de message absente");
    f->discipline = (int)m[6];
 
    const uint8_t *s5 = NULL, *s6 = NULL, *s7 = NULL;
    size_t l5 = 0, l6 = 0, l7 = 0;
    int have3 = 0, have4 = 0;
    size_t q = 16, end = (size_t)mlen - 4;
    while (q < end) {
        if (end - q < 5)
            return fail(err, el, "section tronquée");
        size_t sl = u32(m + q);
        int sn_ = (int)m[q + 4];
        if (sl < 5 || sl > end - q)
            return fail(err, el, "section %d de longueur invalide", sn_);
        const uint8_t *s = m + q;
        switch (sn_) {
        case 1:
            if (sl < 21)
                return fail(err, el, "section 1 courte");
            /* octets 13-19 : année, mois, jour, heure, minute, seconde (numérotation à partir de 1) */
            f->ref_time = epoch_utc((int)u16(s + 12), (int)s[14], (int)s[15], (int)s[16], (int)s[17], (int)s[18]);
            break;
        case 2:
            break;
        case 3: {
            if (sl < 72 || u16(s + 12) != 0)
                return fail(err, el, "grille non gérée (modèle %u)", sl >= 14 ? u16(s + 12) : 0);
            f->ni = (int)u32(s + 30);
            f->nj = (int)u32(s + 34);
            uint32_t basic = u32(s + 38), sub = u32(s + 42);
            double unit = (basic == 0 || basic == 0xffffffffu || sub == 0xffffffffu) ? 1e-6 : (double)basic / sub;
            f->la1 = s32(s + 46) * unit;
            f->lo1 = s32(s + 50) * unit;
            f->la2 = s32(s + 55) * unit;
            f->lo2 = s32(s + 59) * unit;
            f->di = u32(s + 63) * unit;
            f->dj = u32(s + 67) * unit;
            f->scan = (int)s[71];
            if (f->ni <= 0 || f->nj <= 0 || (uint64_t)f->ni * (uint64_t)f->nj > MAX_POINTS)
                return fail(err, el, "grille trop grande");
            if (f->scan & 0x20)
                return fail(err, el, "balayage par colonnes non géré");
            have3 = 1;
            break;
        }
        case 4: {
            if (sl < 34)
                return fail(err, el, "section 4 courte");
            f->pdt = (int)u16(s + 7);
            if (f->pdt != 0 && f->pdt != 8)
                return fail(err, el, "produit 4.%d non géré", f->pdt);
            f->category = (int)s[9];
            f->number = (int)s[10];
            int unit = (int)s[17];
            int64_t ft = (int64_t)u32(s + 18), us = unit_seconds(unit);
            f->surface = (int)s[22];
            int sc = (int8_t)s[23];
            uint32_t sv = u32(s + 24);
            f->surface_value = sv == 0xffffffffu ? NAN : sv * pow(10.0, -((s[23] & 0x80) ? -(sc & 0x7f) : sc));
            if (f->pdt == 0) {
                if (us < 0)
                    return fail(err, el, "unité de temps %d inconnue", unit);
                f->valid_time = f->ref_time + ft * us;
            } else {
                if (sl < 58)
                    return fail(err, el, "section 4.8 courte");
                /* octets 35-41 : fin de l'intervalle statistique, écrite en clair */
                f->valid_time = epoch_utc((int)u16(s + 34), (int)s[36], (int)s[37], (int)s[38], (int)s[39], (int)s[40]);
                f->stat_proc = (int)s[46];
            }
            have4 = 1;
            break;
        }
        case 5:
            s5 = s;
            l5 = sl;
            break;
        case 6:
            s6 = s;
            l6 = sl;
            break;
        case 7:
            s7 = s;
            l7 = sl;
            break;
        default:
            return fail(err, el, "section %d inattendue", sn_);
        }
        q += sl;
    }
    if (!have3 || !have4 || !s5 || !s6 || !s7)
        return fail(err, el, "message incomplet");
    if (want && !want(f->discipline, f->category, f->number, ctx))
        return 1;
 
    /* ---- section 5 */
    Pack k;
    memset(&k, 0, sizeof k);
    if (l5 < 21)
        return fail(err, el, "section 5 courte");
    k.ndata = u32(s5 + 5);
    k.drt = (int)u16(s5 + 9);
    k.R = f32(s5 + 11);
    k.E = s16(s5 + 15);
    k.D = s16(s5 + 17);
    k.bits = (int)u8(s5 + 19);
    f->drt = k.drt;
    f->bits = k.bits;
    size_t npts = (size_t)f->ni * (size_t)f->nj;
    if (k.ndata > npts)
        return fail(err, el, "plus de valeurs que de points");
    if (k.drt == 42) {
        if (l5 < 25)
            return fail(err, el, "section 5.42 courte");
        k.ccsds_flags = (int)s5[21];
        k.block_size = (int)s5[22];
        k.rsi = (int)u16(s5 + 23);
    } else if (k.drt == 3) {
        if (l5 < 49)
            return fail(err, el, "section 5.3 courte");
        k.group_split = (int)s5[21];
        k.missing_mgmt = (int)s5[22];
        k.ng = (int)u32(s5 + 31);
        k.ref_width = (int)s5[35];
        k.bits_width = (int)s5[36];
        k.ref_len = u32(s5 + 37);
        k.len_incr = (int)s5[41];
        k.last_len = u32(s5 + 42);
        k.bits_len = (int)s5[46];
        k.order = (int)s5[47];
        k.extra = (int)s5[48];
        if (k.group_split != 1)
            return fail(err, el, "5.3 : découpage de groupes %d non géré", k.group_split);
    } else if (k.drt != 0) {
        return fail(err, el, "emballage 5.%d non géré", k.drt);
    }
 
    /* ---- section 6 : bitmap */
    const uint8_t *bitmap = NULL;
    if (l6 < 6)
        return fail(err, el, "section 6 courte");
    if (s6[5] == 0) {
        if (l6 - 6 < (npts + 7) / 8)
            return fail(err, el, "bitmap tronqué");
        bitmap = s6 + 6;
    } else if (s6[5] != 255) {
        return fail(err, el, "bitmap %u non géré", s6[5]);
    }
    size_t present = npts;
    if (bitmap) {
        present = 0;
        for (size_t n = 0; n < npts; n++)
            present += (bitmap[n >> 3] >> (7 - (n & 7))) & 1;
    }
    if (present != k.ndata)
        return fail(err, el, "bitmap et nombre de valeurs incohérents");
 
    /* ---- section 7 : données */
    const uint8_t *d = s7 + 5;
    size_t dlen = l7 - 5;
    double scale_e = ldexp(1.0, k.E), scale_d = pow(10.0, -k.D);
    double *vals = xmalloc(k.ndata ? k.ndata : 1, sizeof(double));
    int rc;
    if (k.drt == 3) {
        int64_t *x = xmalloc(k.ndata ? k.ndata : 1, sizeof(int64_t));
        uint8_t *miss = xmalloc(k.ndata ? k.ndata : 1, 1);
        rc = unpack_complex(&k, d, dlen, x, miss, err, el);
        for (size_t n = 0; rc == 0 && n < k.ndata; n++)
            vals[n] = miss[n] ? NAN : ((double)k.R + (double)x[n] * scale_e) * scale_d;
        free(x);
        free(miss);
    } else {
        uint32_t *x = xmalloc(k.ndata ? k.ndata : 1, sizeof(uint32_t));
        rc = k.drt == 42 ? unpack_ccsds(&k, d, dlen, x, err, el) : unpack_simple(&k, d, dlen, x, err, el);
        for (size_t n = 0; rc == 0 && n < k.ndata; n++)
            vals[n] = ((double)k.R + (double)x[n] * scale_e) * scale_d;
        free(x);
    }
    if (rc) {
        free(vals);
        return -1;
    }
    f->values = xmalloc(npts, sizeof(double));
    if (bitmap) {
        size_t c = 0;
        for (size_t n = 0; n < npts; n++)
            f->values[n] = ((bitmap[n >> 3] >> (7 - (n & 7))) & 1) ? vals[c++] : NAN;
    } else {
        memcpy(f->values, vals, npts * sizeof(double));
    }
    free(vals);
    return 1;
}
[code 5]

Du fichier GRIB à la grille de travail

decode.c, 224 lignes

Garde la pression au niveau moyen de la mer (paramètre OMM 0/3/1, ou 0/3/0 au Centre européen, toujours sur la surface 101) et la rafale à 10 m (0/2/22), découpe un cadre de l'Atlantique à l'Europe centrale, remet latitudes et longitudes dans l'ordre croissant et écrit une grille binaire compacte. Un seul message illisible fait refuser le run entier.

decode.c
Fichier brut
/* decode.c : « austen decode <rep_grib> <sortie.grd> »
 * Lit tous les *.grib2 du répertoire, garde la pression au niveau de la mer (hPa) et la rafale à
 * 10 m (km/h), découpe le cadre Atlantique-Europe et range latitudes et longitudes en ordre croissant.
 * Paramètres OMM : pression réduite (0,3,1) ou pression au niveau moyen de la mer (0,3,0, ECMWF),
 * rafale (0,2,22) ; une rafale maxi (traitement statistique 2) l'emporte sur une valeur instantanée. */
#define _DEFAULT_SOURCE
#include <dirent.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
 
#include "grib2.h"
#include "util.h"
 
static const double BOX[4] = {-40.0, 25.0, 33.0, 66.0};   /* lon mini, lon maxi, lat mini, lat maxi */
 
enum { V_MSL, V_GUST };
 
typedef struct {
    int64_t t;
    int stat;
    float *v;
} Slot;
 
typedef struct {
    Slot s[2][MAX_TIMES];
    int n[2];
    int set;                     /* grille fixée par le premier champ */
    int ni, nj, scan;
    double la1, lo1, di, dj;
    int *ix, *jy;                /* indices source retenus, déjà triés */
    Grid g;
} Acc;
 
static int want(int disc, int cat, int num, void *ctx)
{
    (void)ctx;
    return disc == 0 && ((cat == 3 && (num == 0 || num == 1)) || (cat == 2 && num == 22));
}
 
static int cmp_str(const void *a, const void *b) { return strcmp(*(char *const *)a, *(char *const *)b); }
 
typedef struct { double lon; int i; } LonIdx;
static int cmp_lon(const void *a, const void *b)
{
    double x = ((const LonIdx *)a)->lon, y = ((const LonIdx *)b)->lon;
    return x < y ? -1 : x > y;
}
 
static void set_grid(Acc *a, const G2Field *f)
{
    a->ni = f->ni;
    a->nj = f->nj;
    a->scan = f->scan;
    a->la1 = f->la1;
    a->lo1 = f->lo1;
    a->di = f->di;
    a->dj = f->dj;
    if (f->scan & 0x80)
        die("balayage est → ouest non géré");
    /* longitudes ramenées dans [-180, 180), triées */
    LonIdx *li = xmalloc((size_t)f->ni, sizeof *li);
    int nl = 0;
    for (int i = 0; i < f->ni; i++) {
        double lo = fmod(f->lo1 + i * f->di, 360.0);
        if (lo >= 180.0)
            lo -= 360.0;
        if (lo < -180.0)
            lo += 360.0;
        if (lo >= BOX[0] - 1e-9 && lo <= BOX[1] + 1e-9)
            li[nl++] = (LonIdx){lo, i};
    }
    qsort(li, (size_t)nl, sizeof *li, cmp_lon);
    a->ix = xmalloc((size_t)nl, sizeof(int));
    for (int k = 0; k < nl; k++)
        a->ix[k] = li[k].i;
    /* latitudes : sens selon le bit « j positif » du mode de balayage */
    double sgn = (f->scan & 0x40) ? 1.0 : -1.0;
    a->jy = xmalloc((size_t)f->nj, sizeof(int));
    int nlat = 0;
    for (int j = 0; j < f->nj; j++) {
        int jj = sgn > 0 ? j : f->nj - 1 - j;            /* parcours du sud au nord */
        double la = f->la1 + sgn * jj * f->dj;
        if (la >= BOX[2] - 1e-9 && la <= BOX[3] + 1e-9)
            a->jy[nlat++] = jj;
    }
    if (nl < 2 || nlat < 2)
        die("le cadre ne recoupe pas la grille du modèle");
    a->g.nx = nl;
    a->g.ny = nlat;
    a->g.lon0 = li[0].lon;
    a->g.dlon = f->di;
    a->g.lat0 = f->la1 + sgn * a->jy[0] * f->dj;
    a->g.dlat = f->dj;
    free(li);
    a->set = 1;
}
 
static void add(Acc *a, const G2Field *f)
{
    if (!a->set)
        set_grid(a, f);
    if (f->ni != a->ni || f->nj != a->nj || f->scan != a->scan || fabs(f->la1 - a->la1) > 1e-6 ||
        fabs(f->lo1 - a->lo1) > 1e-6 || fabs(f->di - a->di) > 1e-9 || fabs(f->dj - a->dj) > 1e-9)
        die("grille différente d'un message à l'autre");
    int var = f->category == 3 ? V_MSL : V_GUST;
    double scale = var == V_MSL ? 0.01 : 3.6;
    /* même heure de validité : la rafale maxi l'emporte, sinon le dernier message */
    int slot = -1;
    for (int k = 0; k < a->n[var]; k++)
        if (a->s[var][k].t == f->valid_time)
            slot = k;
    if (slot >= 0 && var == V_GUST && a->s[var][slot].stat == 2 && f->stat_proc != 2)
        return;
    if (slot < 0) {
        if (a->n[var] >= MAX_TIMES)
            die("trop d'échéances");
        slot = a->n[var]++;
    } else {
        free(a->s[var][slot].v);
    }
    size_t np = (size_t)a->g.nx * (size_t)a->g.ny;
    float *v = xmalloc(np, sizeof(float));
    for (int y = 0; y < a->g.ny; y++)
        for (int x = 0; x < a->g.nx; x++) {
            double val = f->values[(size_t)a->jy[y] * (size_t)f->ni + (size_t)a->ix[x]];
            v[(size_t)y * (size_t)a->g.nx + (size_t)x] = isnan(val) ? NAN : (float)(val * scale);
        }
    a->s[var][slot] = (Slot){f->valid_time, f->stat_proc, v};
}
 
static int cmp_slot(const void *a, const void *b)
{
    int64_t x = ((const Slot *)a)->t, y = ((const Slot *)b)->t;
    return x < y ? -1 : x > y;
}
 
int cmd_decode(int argc, char **argv)
{
    if (argc != 2)
        die("usage : austen decode <rep_grib> <sortie.grd>");
    DIR *d = opendir(argv[0]);
    if (!d)
        die("répertoire illisible : %s", argv[0]);
    char **names = NULL;
    size_t nn = 0;
    struct dirent *e;
    while ((e = readdir(d))) {
        size_t l = strlen(e->d_name);
        if (e->d_name[0] != '.' && l > 6 && !strcmp(e->d_name + l - 6, ".grib2")) {
            names = xrealloc(names, nn + 1, sizeof *names);
            names[nn++] = xstrdup(e->d_name);
        }
    }
    closedir(d);
    qsort(names, nn, sizeof *names, cmp_str);
    Acc *a = xcalloc(1, sizeof *a);
    int nerr = 0;
    for (size_t k = 0; k < nn; k++) {
        char path[4096];
        if (snprintf(path, sizeof path, "%s/%s", argv[0], names[k]) >= (int)sizeof path)
            die("chemin trop long");
        size_t len;
        uint8_t *buf = read_file(path, (size_t)2 << 30, &len);
        size_t pos = 0;
        char err[256];
        for (;;) {
            G2Field f;
            int rc = g2_next(buf, len, &pos, &f, want, NULL, err, sizeof err);
            if (rc == 0)
                break;
            if (rc < 0) {
                fprintf(stderr, "austen decode : %s : %s\n", names[k], err);
                nerr++;
                continue;
            }
            /* pression : seulement au niveau moyen de la mer (surface 101), jamais la pression au sol */
            if (f.values && !(f.category == 3 && f.surface != 101))
                add(a, &f);
            g2_free(&f);
        }
        free(buf);
        free(names[k]);
    }
    free(names);
    if (nerr)
        die("%d message(s) illisible(s) : run refusé", nerr);
    if (!a->set || a->n[V_MSL] == 0)
        die("aucune pression trouvée");
    Grd out = {0};
    out.g = a->g;
    size_t np = grid_points(&out.g);
    for (int var = 0; var < 2; var++) {
        int n = a->n[var];
        qsort(a->s[var], (size_t)n, sizeof(Slot), cmp_slot);
        int64_t *t = xmalloc((size_t)n, sizeof(int64_t));
        float *v = xmalloc((size_t)n * np, sizeof(float));
        for (int k = 0; k < n; k++) {
            t[k] = a->s[var][k].t;
            memcpy(v + (size_t)k * np, a->s[var][k].v, np * sizeof(float));
            free(a->s[var][k].v);
        }
        if (var == V_MSL) {
            out.nt_msl = n;
            out.t_msl = t;
            out.msl = v;
        } else {
            out.nt_gust = n;
            out.t_gust = t;
            out.gust = v;
        }
    }
    grd_write(argv[1], &out);
    printf("{\"nx\":%d,\"ny\":%d,\"msl\":%d,\"gust\":%d}\n", out.g.nx, out.g.ny, out.nt_msl, out.nt_gust);
    free(a->ix);
    free(a->jy);
    free(a);
    free(out.t_msl);
    free(out.t_gust);
    free(out.msl);
    free(out.gust);
    return 0;
}
[code 6]

Les opérations sur les champs

field.c, 191 lignes

Le cœur numérique, sans bibliothèque : filtre gaussien séparable (bords prolongés), moyenne gaussienne normalisée qui ignore les trous au lieu de les remplir par une valeur fixe, comblement au plus proche, distance de Manhattan, interpolation bilinéaire. Chaque fonction tient en une vingtaine de lignes.

field.c
Fichier brut
#include "field.h"
 
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
 
void conf_read(const char *path, Conf *c)
{
    size_t len;
    char *buf = (char *)read_file(path, 1 << 20, &len);
    memset(c, 0, sizeof *c);
    int lineno = 0;
    char *cur_ = buf;
    for (char *line = next_line(&cur_); line; line = next_line(&cur_)) {
        lineno++;
        char key[16], id[16];
        long long t;
        double a, b;
        if (line[0] == '#' || line[0] == 0)
            continue;
        if (sscanf(line, "%15s", key) != 1)
            die("conf.tsv ligne %d illisible", lineno);
        if (!strcmp(key, "ep0") || !strcmp(key, "ep1") || !strcmp(key, "tr0") || !strcmp(key, "tr1")) {
            if (sscanf(line, "%*s %lld", &t) != 1)
                die("conf.tsv ligne %d : date attendue", lineno);
            int64_t *dst = key[0] == 'e' ? (key[2] == '0' ? &c->ep0 : &c->ep1) : (key[2] == '0' ? &c->tr0 : &c->tr1);
            *dst = t;
        } else if (!strcmp(key, "ref")) {
            if (c->nref >= MAX_REF || sscanf(line, "%*s %lld %lf %lf", &t, &a, &b) != 3)
                die("conf.tsv ligne %d : ref invalide", lineno);
            c->ref_t[c->nref] = t;
            c->ref_lat[c->nref] = a;
            c->ref_lon[c->nref] = b;
            c->nref++;
        } else if (!strcmp(key, "station")) {
            if (c->nst >= MAX_STATIONS || sscanf(line, "%*s %15s %lf %lf", id, &a, &b) != 3)
                die("conf.tsv ligne %d : station invalide", lineno);
            memcpy(c->st_id[c->nst], id, sizeof id);
            c->st_lat[c->nst] = a;
            c->st_lon[c->nst] = b;
            c->nst++;
        } else {
            die("conf.tsv ligne %d : directive inconnue %s", lineno, key);
        }
    }
    free(buf);
    if (!c->ep0 || !c->ep1 || !c->tr0 || !c->tr1 || c->nref == 0)
        die("conf.tsv incomplet");
}
 
void fill_nearest(const float *in, float *out, int nx, int ny)
{
    size_t n = (size_t)nx * (size_t)ny;
    int *queue = xmalloc(n, sizeof(int));
    size_t head = 0, tail = 0;
    for (size_t k = 0; k < n; k++) {
        out[k] = in[k];
        if (!isnan(in[k]))
            queue[tail++] = (int)k;
    }
    if (tail == 0) {
        free(queue);
        return;
    }
    while (head < tail) {
        int k = queue[head++], x = k % nx, y = k / nx;
        for (int dy = -1; dy <= 1; dy++)
            for (int dx = -1; dx <= 1; dx++) {
                int xx = x + dx, yy = y + dy;
                if ((dx || dy) && xx >= 0 && yy >= 0 && xx < nx && yy < ny) {
                    size_t kk = (size_t)yy * (size_t)nx + (size_t)xx;
                    if (isnan(out[kk])) {
                        out[kk] = out[k];
                        queue[tail++] = (int)kk;
                    }
                }
            }
    }
    free(queue);
}
 
static int kernel(double sigma, double **w)
{
    int r = (int)(4.0 * sigma + 0.5);
    *w = xmalloc((size_t)(2 * r + 1), sizeof(double));
    double s = 0;
    for (int i = -r; i <= r; i++) {
        (*w)[i + r] = exp(-0.5 * i * i / (sigma * sigma));
        s += (*w)[i + r];
    }
    for (int i = 0; i <= 2 * r; i++)
        (*w)[i] /= s;
    return r;
}
 
static int clampi(int v, int lo, int hi) { return v < lo ? lo : v > hi ? hi : v; }
 
void gauss_nearest(const float *in, float *out, int nx, int ny, double sigma)
{
    double *w;
    int r = kernel(sigma, &w);
    size_t n = (size_t)nx * (size_t)ny;
    double *tmp = xmalloc(n, sizeof(double));
    /* noyau séparable : passe en latitude puis en longitude */
    for (int y = 0; y < ny; y++)
        for (int x = 0; x < nx; x++) {
            double s = 0;
            for (int k = -r; k <= r; k++)
                s += w[k + r] * in[(size_t)clampi(y + k, 0, ny - 1) * (size_t)nx + (size_t)x];
            tmp[(size_t)y * (size_t)nx + (size_t)x] = s;
        }
    for (int y = 0; y < ny; y++)
        for (int x = 0; x < nx; x++) {
            double s = 0;
            for (int k = -r; k <= r; k++)
                s += w[k + r] * tmp[(size_t)y * (size_t)nx + (size_t)clampi(x + k, 0, nx - 1)];
            out[(size_t)y * (size_t)nx + (size_t)x] = (float)s;
        }
    free(tmp);
    free(w);
}
 
void gauss_masked(const float *in, float *out, int nx, int ny, double sigma)
{
    size_t n = (size_t)nx * (size_t)ny;
    float *v = xmalloc(n, sizeof(float)), *m = xmalloc(n, sizeof(float));
    float *sv = xmalloc(n, sizeof(float)), *sm = xmalloc(n, sizeof(float));
    for (size_t k = 0; k < n; k++) {
        int ok = !isnan(in[k]);
        v[k] = ok ? in[k] : 0.0f;
        m[k] = ok ? 1.0f : 0.0f;
    }
    gauss_nearest(v, sv, nx, ny, sigma);
    gauss_nearest(m, sm, nx, ny, sigma);
    for (size_t k = 0; k < n; k++)
        out[k] = isnan(in[k]) ? NAN : sv[k] / (sm[k] > 1e-6f ? sm[k] : 1e-6f);
    free(v);
    free(m);
    free(sv);
    free(sm);
}
 
void city_block(const uint8_t *mask, int *dist, int nx, int ny)
{
    const int inf = 1 << 28;
    for (int y = 0; y < ny; y++)
        for (int x = 0; x < nx; x++) {
            size_t k = (size_t)y * (size_t)nx + (size_t)x;
            if (!mask[k]) {
                dist[k] = 0;
                continue;
            }
            int up = y > 0 ? dist[k - (size_t)nx] : 0, left = x > 0 ? dist[k - 1] : 0;
            int d = (up < left ? up : left) + 1;
            dist[k] = d < inf ? d : inf;
        }
    for (int y = ny - 1; y >= 0; y--)
        for (int x = nx - 1; x >= 0; x--) {
            size_t k = (size_t)y * (size_t)nx + (size_t)x;
            if (!dist[k])
                continue;
            int dn = y < ny - 1 ? dist[k + (size_t)nx] : 0, rt = x < nx - 1 ? dist[k + 1] : 0;
            int d = (dn < rt ? dn : rt) + 1;
            if (d < dist[k])
                dist[k] = d;
        }
}
 
double bilinear(const Grid *g, const float *f, double lat, double lon)
{
    double fy = (lat - g->lat0) / g->dlat, fx = (lon - g->lon0) / g->dlon;
    if (fx < 0 || fy < 0 || fx > g->nx - 1 || fy > g->ny - 1)
        return NAN;
    int x0 = (int)fx, y0 = (int)fy;
    if (x0 >= g->nx - 1)
        x0 = g->nx - 2;
    if (y0 >= g->ny - 1)
        y0 = g->ny - 2;
    if (x0 < 0 || y0 < 0)
        return NAN;
    double ax = fx - x0, ay = fy - y0;
    size_t k = (size_t)y0 * (size_t)g->nx + (size_t)x0;
    double a = f[k], b = f[k + 1], c = f[k + (size_t)g->nx], d = f[k + (size_t)g->nx + 1];
    return (1 - ay) * ((1 - ax) * a + ax * b) + ay * ((1 - ax) * c + ax * d);
}
 
double dist_km(double lat, double lon, double la0, double lo0)
{
    return KM_PER_DEG * hypot(lat - la0, (lon - lo0) * cos(la0 * DEG));
}
[code 7]

Le suivi du centre et les rafales maxi

analyze.c, 344 lignes

Le script ne cherche pas le point le plus bas de toute la carte : il sauterait d'une dépression à l'autre. Il suit Austen de proche en proche :

  1. Accrochage. Jeudi à 12 UTC, minimum de la pression lissée (environ 45 km) dans un rayon de 350 km autour de la position de référence.
  2. Suivi. D'heure en heure, position prolongée à la vitesse récente, minimum cherché dans 140 km plus 45 km par heure d'écart. Il remonte aussi le temps.
  3. Garde-fous. Arrêt si le minimum touche le bord du cercle (une pente, pas un centre fermé), s'il sort de la zone de données valides, ou si la pression varie de plus de 2 hPa + 1,5 hPa par heure.
  4. Le « D ». Position affinée entre deux points de grille ; pression = vraie valeur la plus basse, non lissée, dans un rayon de 70 km.

Il calcule aussi la rafale maxi de l'épisode en chaque point et les séries horaires aux 38 stations.

analyze.c
Fichier brut
/* analyze.c : « austen analyze <conf.tsv> <run.grd> <sortie.json> <gmax.grd> »
 *
 * 1. Suivi du centre de la dépression :
 *    - accrochage au premier instant >= première position de référence, minimum de la pression lissée
 *      dans 350 km (puis 500 km) autour de la référence ;
 *    - suivi heure par heure, vers l'avant puis vers l'arrière : position prolongée à la vitesse
 *      récente, recherche dans 140 km + 45 km par heure d'écart ;
 *    - arrêt si le minimum touche le bord du cercle (pas de centre fermé), s'il sort de l'intérieur
 *      des données valides, ou si la pression saute de plus de 2 hPa + 1,5 hPa par heure.
 *    La pression retenue est la valeur brute la plus basse à moins de 70 km du minimum lissé.
 * 2. Rafale maxi de l'épisode (ep0, ep1] en chaque point → gmax.grd.
 * 3. Séries de rafales aux stations (interpolation bilinéaire) et maxi sur l'épisode.
 */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
 
#include "field.h"
#include "json.h"
#include "util.h"
 
typedef struct {
    const Grid *g;
    int nt;
    const float *raw;          /* [nt][np] pression brute */
    float *smooth;             /* [nt][np] pression lissée */
    uint8_t *interior;         /* [np] */
    double dkm;
} Tracker;
 
typedef struct {
    int ok;
    double lat, lon, p;
} Hit;
 
static double ref_lat_lon(const Conf *c, int64_t t, double *lon)
{
    if (t <= c->ref_t[0]) {
        *lon = c->ref_lon[0];
        return c->ref_lat[0];
    }
    for (int k = 0; k + 1 < c->nref; k++)
        if (t >= c->ref_t[k] && t <= c->ref_t[k + 1]) {
            double f = (double)(t - c->ref_t[k]) / (double)(c->ref_t[k + 1] - c->ref_t[k]);
            *lon = c->ref_lon[k] + f * (c->ref_lon[k + 1] - c->ref_lon[k]);
            return c->ref_lat[k] + f * (c->ref_lat[k + 1] - c->ref_lat[k]);
        }
    *lon = c->ref_lon[c->nref - 1];
    return c->ref_lat[c->nref - 1];
}
 
static void tracker_init(Tracker *tr, const Grid *g, const float *raw, int nt)
{
    size_t np = grid_points(g);
    tr->g = g;
    tr->nt = nt;
    tr->raw = raw;
    tr->dkm = fabs(g->dlat) * KM_PER_DEG;
    double sigma = fmax(0.8, 45.0 / tr->dkm);
    tr->smooth = xmalloc((size_t)nt * np, sizeof(float));
    float *filled = xmalloc(np, sizeof(float));
    for (int t = 0; t < nt; t++) {
        const float *f = raw + (size_t)t * np;
        /* trous comblés au plus proche avant lissage : une valeur fixe créerait un faux minimum */
        fill_nearest(f, filled, g->nx, g->ny);
        float mx = -INFINITY;
        for (size_t k = 0; k < np; k++)
            if (!isnan(filled[k]) && filled[k] > mx)
                mx = filled[k];
        for (size_t k = 0; k < np; k++)
            if (isnan(filled[k]))
                filled[k] = mx;
        gauss_nearest(filled, tr->smooth + (size_t)t * np, g->nx, g->ny, sigma);
    }
    free(filled);
    /* intérieur : à plus de « edge » points (≈ 60 km) des bords de grille et de toute zone sans donnée */
    int edge = (int)fmax(2.0, round(60.0 / tr->dkm));
    uint8_t *fin = xmalloc(np, 1);
    for (int y = 0; y < g->ny; y++)
        for (int x = 0; x < g->nx; x++) {
            size_t k = (size_t)y * (size_t)g->nx + (size_t)x;
            fin[k] = !isnan(raw[k]) && x >= edge && y >= edge && x < g->nx - edge && y < g->ny - edge;
        }
    int *dist = xmalloc(np, sizeof(int));
    city_block(fin, dist, g->nx, g->ny);
    tr->interior = xmalloc(np, 1);
    for (size_t k = 0; k < np; k++)
        tr->interior[k] = dist[k] > edge;
    free(dist);
    free(fin);
}
 
static Hit find(const Tracker *tr, int t, double la0, double lo0, double R)
{
    const Grid *g = tr->g;
    size_t np = (size_t)g->nx * (size_t)g->ny;
    const float *s = tr->smooth + (size_t)t * np, *raw = tr->raw + (size_t)t * np;
    Hit h = {0, 0, 0, 0};
    /* boîte englobante du disque pour ne pas parcourir toute la grille */
    double dlat = R / KM_PER_DEG + fabs(g->dlat), dlon = R / (KM_PER_DEG * fmax(0.1, cos(la0 * DEG))) + fabs(g->dlon);
    int y0 = (int)floor((la0 - dlat - g->lat0) / g->dlat), y1 = (int)ceil((la0 + dlat - g->lat0) / g->dlat);
    int x0 = (int)floor((lo0 - dlon - g->lon0) / g->dlon), x1 = (int)ceil((lo0 + dlon - g->lon0) / g->dlon);
    y0 = y0 < 0 ? 0 : y0;
    x0 = x0 < 0 ? 0 : x0;
    y1 = y1 > g->ny - 1 ? g->ny - 1 : y1;
    x1 = x1 > g->nx - 1 ? g->nx - 1 : x1;
    double best = INFINITY;
    int bj = -1, bk = -1;
    for (int j = y0; j <= y1; j++)
        for (int k = x0; k <= x1; k++) {
            double la = g->lat0 + j * g->dlat, lo = g->lon0 + k * g->dlon;
            if (dist_km(la, lo, la0, lo0) > R)
                continue;
            double v = s[(size_t)j * (size_t)g->nx + (size_t)k];
            if (v < best) {            /* premier minimum dans l'ordre des lignes, comme numpy.argmin */
                best = v;
                bj = j;
                bk = k;
            }
        }
    if (bj < 0)
        return h;
    if (!tr->interior[(size_t)bj * (size_t)g->nx + (size_t)bk])
        return h;   /* bord du domaine : la dépression n'est pas (ou plus) dans le modèle */
    double la = g->lat0 + bj * g->dlat, lo = g->lon0 + bk * g->dlon;
    if (dist_km(la, lo, la0, lo0) > R - 2.5 * tr->dkm)
        return h;   /* minimum collé au bord du cercle : pente, pas centre fermé */
    /* pression : valeur brute la plus basse à moins de 70 km */
    double p = INFINITY;
    int r = (int)ceil(70.0 / tr->dkm) + 2;
    double coslat = fmax(0.1, cos(la * DEG));
    int rx = (int)ceil(70.0 / (fabs(g->dlon) * KM_PER_DEG * coslat)) + 2;
    for (int j = bj - r; j <= bj + r; j++)
        for (int k = bk - rx; k <= bk + rx; k++) {
            if (j < 0 || k < 0 || j >= g->ny || k >= g->nx)
                continue;
            if (dist_km(g->lat0 + j * g->dlat, g->lon0 + k * g->dlon, la, lo) > 70.0)
                continue;
            double v = raw[(size_t)j * (size_t)g->nx + (size_t)k];
            if (!isnan(v) && v < p)
                p = v;
        }
    if (!isfinite(p))
        return h;
    /* affinage parabolique de la position sur le champ lissé */
    double fla = la, flo = lo;
    if (bj > 0 && bj < g->ny - 1 && bk > 0 && bk < g->nx - 1) {
        size_t c = (size_t)bj * (size_t)g->nx + (size_t)bk;
        double a = s[c - (size_t)g->nx], b = s[c], d = s[c + (size_t)g->nx];
        double den = a - 2 * b + d;
        if (den > 0)
            fla += 0.5 * (a - d) / den * g->dlat;
        a = s[c - 1];
        d = s[c + 1];
        den = a - 2 * b + d;
        if (den > 0)
            flo += 0.5 * (a - d) / den * g->dlon;
    }
    h.ok = 1;
    h.lat = fla;
    h.lon = flo;
    h.p = p;
    return h;
}
 
static void walk(const Tracker *tr, const int64_t *times, Hit *found, int seed, int dir)
{
    int n = seed;
    int have_prev = 0;
    double pla = 0, plo = 0, pdt = 1;
    for (;;) {
        int m = n + dir;
        if (m < 0 || m >= tr->nt)
            return;
        double dth = fabs((double)(times[m] - times[n])) / 3600.0;
        double la = found[n].lat, lo = found[n].lon;
        if (have_prev) {
            la += (la - pla) * dth / pdt;
            lo += (lo - plo) * dth / pdt;
        }
        double R = 140.0 + 45.0 * dth;
        Hit r = find(tr, m, la, lo, R);
        if (!r.ok || fabs(r.p - found[n].p) > 2.0 + 1.5 * dth)
            return;
        pla = found[n].lat;
        plo = found[n].lon;
        pdt = dth;
        have_prev = 1;
        found[m] = r;
        n = m;
    }
}
 
int cmd_analyze(int argc, char **argv)
{
    if (argc != 4)
        die("usage : austen analyze <conf.tsv> <run.grd> <sortie.json> <gmax.grd>");
    Conf c;
    conf_read(argv[0], &c);
    Grd *d = grd_read(argv[1]);
    const Grid *g = &d->g;
    size_t np = grid_points(g);
 
    /* ---- trajectoire, sur les instants de la fenêtre [tr0, tr1] */
    int first = -1, nt = 0;
    for (int k = 0; k < d->nt_msl; k++)
        if (d->t_msl[k] >= c.tr0 && d->t_msl[k] <= c.tr1) {
            if (first < 0)
                first = k;
            nt++;
        }
    Hit *found = xcalloc((size_t)(nt ? nt : 1), sizeof(Hit));
    const int64_t *times = first >= 0 ? d->t_msl + first : NULL;
    Tracker tr = {0};
    if (nt > 0) {
        tracker_init(&tr, g, d->msl + (size_t)first * np, nt);
        int seed = 0;
        for (int k = 0; k < nt; k++)
            if (times[k] >= c.ref_t[0]) {
                seed = k;
                break;
            }
        double lo0, la0 = ref_lat_lon(&c, times[seed], &lo0);
        Hit h = find(&tr, seed, la0, lo0, 350.0);
        if (!h.ok)
            h = find(&tr, seed, la0, lo0, 500.0);
        if (!h.ok)   /* centre pas encore dans le domaine (AROME) : on essaie plus tard */
            for (int k = seed + 1; k < nt && !h.ok; k++) {
                la0 = ref_lat_lon(&c, times[k], &lo0);
                h = find(&tr, k, la0, lo0, 300.0);
                if (h.ok)
                    seed = k;
            }
        if (h.ok) {
            found[seed] = h;
            walk(&tr, times, found, seed, +1);
            walk(&tr, times, found, seed, -1);
        }
    }
 
    Json j;
    js_init(&j);
    js_obj(&j);
    js_key(&j, "track");
    js_arr(&j);
    char iso[20];
    for (int k = 0; k < nt; k++)
        if (found[k].ok) {
            fmt_iso(times[k], iso);
            js_obj(&j);
            js_key(&j, "t");
            js_str(&j, iso);
            js_key(&j, "lat");
            js_num(&j, found[k].lat, 3);
            js_key(&j, "lon");
            js_num(&j, found[k].lon, 3);
            js_key(&j, "p");
            js_num(&j, found[k].p, 1);
            js_end_obj(&j);
        }
    js_end_arr(&j);
 
    /* ---- rafales */
    if (d->nt_gust > 0) {
        float *gmax = xmalloc(np, sizeof(float));
        for (size_t k = 0; k < np; k++)
            gmax[k] = NAN;
        for (int t = 0; t < d->nt_gust; t++) {
            if (!(d->t_gust[t] > c.ep0 && d->t_gust[t] <= c.ep1))
                continue;
            const float *f = d->gust + (size_t)t * np;
            for (size_t k = 0; k < np; k++)
                if (!isnan(f[k]) && (isnan(gmax[k]) || f[k] > gmax[k]))
                    gmax[k] = f[k];
        }
        Grd out = {0};
        out.g = *g;
        out.nt_gust = 1;
        int64_t tend = c.ep1;
        out.t_gust = &tend;
        out.gust = gmax;
        out.t_msl = NULL;
        grd_write(argv[3], &out);
        double dom = -INFINITY;
        for (size_t k = 0; k < np; k++)
            if (!isnan(gmax[k]) && gmax[k] > dom)
                dom = gmax[k];
 
        /* séries aux stations sur (tr0, tr1] */
        js_key(&j, "gust_times");
        js_arr(&j);
        for (int t = 0; t < d->nt_gust; t++)
            if (d->t_gust[t] > c.tr0 && d->t_gust[t] <= c.tr1) {
                fmt_iso(d->t_gust[t], iso);
                js_str(&j, iso);
            }
        js_end_arr(&j);
        js_key(&j, "stations");
        js_obj(&j);
        for (int s = 0; s < c.nst; s++) {
            double mx = NAN;
            int64_t tmx = 0;
            js_key(&j, c.st_id[s]);
            js_obj(&j);
            js_key(&j, "series");
            js_arr(&j);
            for (int t = 0; t < d->nt_gust; t++) {
                if (!(d->t_gust[t] > c.tr0 && d->t_gust[t] <= c.tr1))
                    continue;
                double v = bilinear(g, d->gust + (size_t)t * np, c.st_lat[s], c.st_lon[s]);
                v = isnan(v) ? NAN : round(v * 10.0) / 10.0;
                js_num(&j, v, 1);
                if (!isnan(v) && d->t_gust[t] > c.ep0 && d->t_gust[t] <= c.ep1 && (isnan(mx) || v > mx)) {
                    mx = v;
                    tmx = d->t_gust[t];
                }
            }
            js_end_arr(&j);
            js_key(&j, "max");
            js_num(&j, mx, 1);
            js_key(&j, "tmax");
            if (isnan(mx))
                js_null(&j);
            else {
                fmt_iso(tmx, iso);
                js_str(&j, iso);
            }
            js_end_obj(&j);
        }
        js_end_obj(&j);
        js_key(&j, "gmax_domain");
        js_num(&j, dom, 1);
        free(gmax);
    }
    js_end_obj(&j);
    js_save(&j, argv[2]);
    js_free(&j);
    free(found);
    free(tr.smooth);
    free(tr.interior);
    grd_free(d);
    return 0;
}
[code 8]

La moyenne des modèles

ensemble.c, 287 lignes

Chaque carte de rafales maxi est interpolée sur une grille commune de 0,025°, puis moyennée point par point (minimum, maximum, nombre de modèles à 100 km/h). Par département : masque par test point dans polygone, maximum, médiane et 90e centile.

ensemble.c
Fichier brut
/* ensemble.c : « austen ensemble <departements.bin> <sortie> modèle=gmax.grd ... »
 *
 * Interpole chaque carte de rafales maxi sur une grille commune de 0,025° (celle d'AROME) ; un modèle
 * peut recevoir plusieurs runs (« modèle=f1,f2 »), moyennés d'abord entre eux. Puis
 * calcule point par point la moyenne, le minimum, le maximum et le nombre de modèles >= 100 km/h.
 * Écrit <sortie>.grd (couches : moyenne, mini, maxi, n100, puis un modèle par couche) et
 * <sortie>.json (noms des couches ; par département : maximum, médiane et 90e centile de la
 * moyenne, maximum de chaque modèle).
 *
 * departements.bin (écrit par la partie PHP) : "AUPO", u32 nb ; pour chaque département :
 * u8 longueur du code, code, u32 nb d'anneaux extérieurs ; par anneau : u32 nb de points, f64 lon, lat.
 */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
 
#include "field.h"
#include "json.h"
#include "util.h"
 
#define LON0 -7.5
#define LON1 10.0
#define LAT0 41.2
#define LAT1 52.4
#define STEP 0.025
#define MAX_MODELS 16
 
typedef struct {
    char code[8];
    int nring;
    int *npt;
    double **xy;      /* lon, lat alternés */
    double bb[4];     /* lon mini, lat mini, lon maxi, lat maxi */
} Dept;
 
static uint32_t rd32(const uint8_t *b, size_t len, size_t *p)
{
    if (*p + 4 > len)
        die("departements.bin tronqué");
    uint32_t v;
    memcpy(&v, b + *p, 4);
    *p += 4;
    return v;
}
 
static Dept *read_depts(const char *path, int *n)
{
    size_t len, p = 4;
    uint8_t *b = read_file(path, 64 << 20, &len);
    if (len < 8 || memcmp(b, "AUPO", 4) != 0)
        die("departements.bin : format inconnu");
    uint32_t nd = rd32(b, len, &p);
    if (nd > 1000)
        die("departements.bin : trop de départements");
    Dept *d = xcalloc(nd, sizeof *d);
    for (uint32_t k = 0; k < nd; k++) {
        if (p >= len || b[p] == 0 || b[p] >= sizeof d[k].code || p + 1 + b[p] > len)
            die("departements.bin : code invalide");
        memcpy(d[k].code, b + p + 1, b[p]);
        p += 1 + (size_t)b[p];
        uint32_t nr = rd32(b, len, &p);
        if (nr > 10000)
            die("departements.bin : trop d'anneaux");
        d[k].nring = (int)nr;
        d[k].npt = xmalloc(nr, sizeof(int));
        d[k].xy = xmalloc(nr, sizeof(double *));
        d[k].bb[0] = d[k].bb[1] = INFINITY;
        d[k].bb[2] = d[k].bb[3] = -INFINITY;
        for (uint32_t r = 0; r < nr; r++) {
            uint32_t np = rd32(b, len, &p);
            if (np > 10000000 || p + (size_t)np * 16 > len)
                die("departements.bin : anneau tronqué");
            d[k].npt[r] = (int)np;
            d[k].xy[r] = xmalloc((size_t)np * 2, sizeof(double));
            memcpy(d[k].xy[r], b + p, (size_t)np * 16);
            p += (size_t)np * 16;
            for (uint32_t q = 0; q < np; q++) {
                double x = d[k].xy[r][2 * q], y = d[k].xy[r][2 * q + 1];
                d[k].bb[0] = fmin(d[k].bb[0], x);
                d[k].bb[1] = fmin(d[k].bb[1], y);
                d[k].bb[2] = fmax(d[k].bb[2], x);
                d[k].bb[3] = fmax(d[k].bb[3], y);
            }
        }
    }
    free(b);
    *n = (int)nd;
    return d;
}
 
/* point dans un polygone : règle pair-impair (rayon horizontal) */
static int inside_ring(const double *xy, int n, double x, double y)
{
    int in = 0;
    for (int i = 0, j = n - 1; i < n; j = i++) {
        double xi = xy[2 * i], yi = xy[2 * i + 1], xj = xy[2 * j], yj = xy[2 * j + 1];
        if ((yi > y) != (yj > y) && x < (xj - xi) * (y - yi) / (yj - yi) + xi)
            in = !in;
    }
    return in;
}
 
static int cmp_f(const void *a, const void *b)
{
    float x = *(const float *)a, y = *(const float *)b;
    return x < y ? -1 : x > y;
}
 
/* centile à interpolation linéaire, comme numpy.percentile (méthode par défaut) */
static double percentile(const float *sorted, size_t n, double q)
{
    double pos = q / 100.0 * (double)(n - 1);
    size_t lo = (size_t)floor(pos), hi = lo + 1 < n ? lo + 1 : lo;
    return sorted[lo] + (pos - (double)lo) * (sorted[hi] - sorted[lo]);
}
 
int cmd_ensemble(int argc, char **argv)
{
    if (argc < 3)
        die("usage : austen ensemble <departements.bin> <sortie> modèle=gmax.grd ...");
    int nm = argc - 2;
    if (nm > MAX_MODELS)
        die("trop de modèles");
    Grid cg = {(int)lround((LON1 - LON0) / STEP) + 1, (int)lround((LAT1 - LAT0) / STEP) + 1, LAT0, STEP, LON0, STEP};
    size_t np = grid_points(&cg);
    int nl = 4 + nm;
    float *layers = xmalloc((size_t)nl * np, sizeof(float));
    char names[MAX_MODELS][32];
    for (int m = 0; m < nm; m++) {
        const char *arg = argv[2 + m], *eq = strchr(arg, '=');
        if (!eq || eq == arg || (size_t)(eq - arg) >= sizeof names[m])
            die("argument modèle=fichier attendu : %s", arg);
        memcpy(names[m], arg, (size_t)(eq - arg));
        names[m][eq - arg] = 0;
        /* « modèle=f1,f2,... » : plusieurs runs du même modèle, moyennés point par point avant la
         * moyenne entre modèles, pour qu'un modèle à 8 runs par jour ne pèse pas plus qu'un autre */
        float *dst = layers + (size_t)(4 + m) * np;
        float *sum = xcalloc(np, sizeof(float));
        int *cnt = xcalloc(np, sizeof(int));
        char *list = xstrdup(eq + 1), *cur = list;
        int nfiles = 0;
        while (cur && *cur) {
            char *comma = strchr(cur, ',');
            if (comma)
                *comma = 0;
            Grd *d = grd_read(cur);
            if (d->nt_gust != 1)
                die("%s : une seule couche attendue", cur);
            for (int y = 0; y < cg.ny; y++)
                for (int x = 0; x < cg.nx; x++) {
                    double v = bilinear(&d->g, d->gust, cg.lat0 + y * cg.dlat, cg.lon0 + x * cg.dlon);
                    size_t q = (size_t)y * (size_t)cg.nx + (size_t)x;
                    if (!isnan(v)) {
                        sum[q] += (float)v;
                        cnt[q]++;
                    }
                }
            grd_free(d);
            nfiles++;
            cur = comma ? comma + 1 : NULL;
        }
        if (!nfiles)
            die("aucun fichier pour %s", names[m]);
        for (size_t q = 0; q < np; q++)
            dst[q] = cnt[q] ? sum[q] / (float)cnt[q] : NAN;
        free(list);
        free(sum);
        free(cnt);
    }
    float *mean = layers, *mn = layers + np, *mx = layers + 2 * np, *n100 = layers + 3 * np;
    for (size_t k = 0; k < np; k++) {
        double s = 0;
        int c = 0, over = 0;
        float lo = NAN, hi = NAN;
        for (int m = 0; m < nm; m++) {
            float v = layers[(size_t)(4 + m) * np + k];
            if (isnan(v))
                continue;
            s += v;
            c++;
            over += v >= 100.0f;
            lo = isnan(lo) || v < lo ? v : lo;
            hi = isnan(hi) || v > hi ? v : hi;
        }
        mean[k] = c ? (float)(s / c) : NAN;
        mn[k] = lo;
        mx[k] = hi;
        n100[k] = (float)over;
    }
    char path[4096];
    if (snprintf(path, sizeof path, "%s.grd", argv[1]) >= (int)sizeof path)
        die("chemin trop long");
    int64_t *lt = xmalloc((size_t)nl, sizeof(int64_t));
    for (int l = 0; l < nl; l++)
        lt[l] = l;
    Grd out = {cg, 0, nl, NULL, lt, NULL, layers};
    grd_write(path, &out);
 
    /* départements */
    int nd;
    Dept *dep = read_depts(argv[0], &nd);
    Json j;
    js_init(&j);
    js_obj(&j);
    js_key(&j, "layers");
    js_arr(&j);
    js_str(&j, "mean");
    js_str(&j, "min");
    js_str(&j, "max");
    js_str(&j, "n100");
    for (int m = 0; m < nm; m++)
        js_str(&j, names[m]);
    js_end_arr(&j);
    js_key(&j, "mean_max");
    double mm = -INFINITY;
    for (size_t k = 0; k < np; k++)
        if (!isnan(mean[k]) && mean[k] > mm)
            mm = mean[k];
    js_num(&j, mm, 1);
    js_key(&j, "depts");
    js_arr(&j);
    float *vals = xmalloc(np, sizeof(float));
    for (int k = 0; k < nd; k++) {
        size_t nv = 0;
        double per[MAX_MODELS];
        for (int m = 0; m < nm; m++)
            per[m] = NAN;
        int y0 = (int)floor((dep[k].bb[1] - cg.lat0) / cg.dlat), y1 = (int)ceil((dep[k].bb[3] - cg.lat0) / cg.dlat);
        int x0 = (int)floor((dep[k].bb[0] - cg.lon0) / cg.dlon), x1 = (int)ceil((dep[k].bb[2] - cg.lon0) / cg.dlon);
        for (int y = y0 < 0 ? 0 : y0; y <= y1 && y < cg.ny; y++)
            for (int x = x0 < 0 ? 0 : x0; x <= x1 && x < cg.nx; x++) {
                double lo = cg.lon0 + x * cg.dlon, la = cg.lat0 + y * cg.dlat;
                int in = 0;
                for (int r = 0; r < dep[k].nring && !in; r++)
                    in = inside_ring(dep[k].xy[r], dep[k].npt[r], lo, la);
                if (!in)
                    continue;
                size_t q = (size_t)y * (size_t)cg.nx + (size_t)x;
                if (!isnan(mean[q]))
                    vals[nv++] = mean[q];
                for (int m = 0; m < nm; m++) {
                    float v = layers[(size_t)(4 + m) * np + q];
                    if (!isnan(v) && (isnan(per[m]) || v > per[m]))
                        per[m] = v;
                }
            }
        if (nv < 5)
            continue;
        qsort(vals, nv, sizeof(float), cmp_f);
        js_obj(&j);
        js_key(&j, "code");
        js_str(&j, dep[k].code);
        js_key(&j, "max");
        js_num(&j, round(vals[nv - 1]), 0);
        js_key(&j, "med");
        js_num(&j, round(percentile(vals, nv, 50)), 0);
        js_key(&j, "p90");
        js_num(&j, round(percentile(vals, nv, 90)), 0);
        js_key(&j, "models");
        js_obj(&j);
        for (int m = 0; m < nm; m++)
            if (!isnan(per[m])) {
                js_key(&j, names[m]);
                js_num(&j, round(per[m]), 0);
            }
        js_end_obj(&j);
        js_end_obj(&j);
    }
    js_end_arr(&j);
    js_end_obj(&j);
    if (snprintf(path, sizeof path, "%s.json", argv[1]) >= (int)sizeof path)
        die("chemin trop long");
    js_save(&j, path);
    js_free(&j);
    free(vals);
    for (int k = 0; k < nd; k++) {
        for (int r = 0; r < dep[k].nring; r++)
            free(dep[k].xy[r]);
        free(dep[k].xy);
        free(dep[k].npt);
    }
    free(dep);
    free(lt);
    free(layers);
    return 0;
}
[code 9]

Les cartes et l'animation

render.c, 737 lignes

Tout le dessin est fait ici : couleurs de l'échelle Meteociel, traits anticrénelés calculés par distance au segment, isobares par l'algorithme des « marching squares » (avec le cas du col tranché par la moyenne des quatre coins), textes par FreeType avec un halo blanc, image finale en WebP. Projection longitude-latitude corrigée par le cosinus de la latitude moyenne.

render.c
Fichier brut
/* render.c : « austen render <spec.txt> » : dessine les cartes décrites dans le fichier de
 * description (un bloc par image, terminé par « end »), écrit par la partie PHP.
 *
 * Directives d'un bloc (le texte libre va jusqu'à la fin de la ligne) :
 *   out <fichier.webp>            quality <0..100>          width <px>
 *   box <lon0> <lon1> <lat0> <lat1>                         background <#rrggbb>
 *   geo <rep> <départements 0|1> <échelle des traits>       (coast.bin, border.bin, dept.bin)
 *   field <fichier.grd> <couche> gust|n100 [<nb de modèles>]
 *   isobars <fichier.grd> <indice pmer> <lissage 0|1>
 *   center <lon> <lat> <pression>                           (« D » rouge et pression)
 *   labels <fichier.tsv>   (lon lat texte : point noir et valeur, rafales par station)
 *   dots <fichier.tsv>     (lon lat valeur : pastilles colorées, observations)
 *   title <texte>   sub <texte>   right <texte>   sign <texte>
 *   image <fichier.webp>   (fond : image existante, en projection lon/lat simple, dont la taille fixe celle
 *                           de la carte ; sert aux images satellite)
 *   geocolor <#rrggbb>     (couleur des côtes, frontières et départements)
 *   idw <fichier.tsv> <rayon km>   (champ coloré interpolé entre stations : lon lat valeur ; pondération
 *                                   par l'inverse du carré de la distance, rien au-delà du rayon)
 *
 * Projection : longitude-latitude, échelle horizontale corrigée par le cosinus de la latitude
 * moyenne (suffisant à l'échelle de la France). Traits anticrénelés par distance au segment,
 * isobares par « marching squares », textes par FreeType, image finale en WebP (libwebp).
 */
#include <ft2build.h>
#include FT_FREETYPE_H
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <webp/decode.h>
#include <webp/encode.h>
 
#include "field.h"
#include "util.h"
 
#define PT (100.0 / 72.0)   /* un point typographique à 100 ppp, comme les cartes matplotlib d'origine */
#define FONT_BOLD "/usr/share/fonts/truetype/dejavu/DejaVuSans-Bold.ttf"
#define FONT_REG "/usr/share/fonts/truetype/dejavu/DejaVuSans.ttf"
 
/* échelle Meteociel des rafales : 30 cases, borne basse de chaque case */
static const double GB[30] = {0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90,
                              100, 110, 120, 130, 140, 150, 160, 180, 200, 220, 240};
static const uint32_t GC[30] = {0xffffff, 0x97e6ff, 0x33ccff, 0x0099ff, 0x00ff99, 0x33cc66, 0x66cc33, 0x66ff00,
                                0xa4f22f, 0xb7cf0e, 0xd6f017, 0xcc9900, 0xff9900, 0xff9966, 0xcc9999, 0xcc6633,
                                0xcc3333, 0xff0d0d, 0xc60000, 0x800000, 0x800050, 0xa00077, 0xcc00cc, 0xff00ff,
                                0xff40ff, 0xff80ff, 0xffa0ff, 0xffffff, 0xe3e3e3, 0xc6c6c6};
static const uint32_t N100C[7] = {0xffffff, 0xfff3bf, 0xffd43b, 0xfd7e14, 0xe03131, 0x862e9c, 0x3b0a45};
#define BAD 0xe9ecef
 
typedef struct {
    int w, h;
    uint8_t *rgb;
    float *cov;        /* couverture des traits en cours (0..1), réutilisée par couche */
    double lon0, lon1, lat0, lat1;
} Canvas;
 
static FT_Library ft;
static FT_Face face_bold, face_reg;
 
/* ------------------------------------------------------------ pixels */
static void blend(Canvas *c, int x, int y, uint32_t col, double a)
{
    if (x < 0 || y < 0 || x >= c->w || y >= c->h || a <= 0)
        return;
    if (a > 1)
        a = 1;
    uint8_t *p = c->rgb + 3 * ((size_t)y * (size_t)c->w + (size_t)x);
    double r = (col >> 16) & 255, g = (col >> 8) & 255, b = col & 255;
    p[0] = (uint8_t)lround(p[0] + (r - p[0]) * a);
    p[1] = (uint8_t)lround(p[1] + (g - p[1]) * a);
    p[2] = (uint8_t)lround(p[2] + (b - p[2]) * a);
}
 
static void proj(const Canvas *c, double lon, double lat, double *x, double *y)
{
    *x = (lon - c->lon0) / (c->lon1 - c->lon0) * c->w;
    *y = (c->lat1 - lat) / (c->lat1 - c->lat0) * c->h;
}
 
static uint32_t parse_color(const char *s)
{
    unsigned v;
    if (s[0] != '#' || strlen(s) != 7 || sscanf(s + 1, "%6x", &v) != 1)
        die("couleur invalide : %s", s);
    return v;
}
 
/* ------------------------------------------------------------ traits */
/* Couverture d'un segment épais (demi-largeur hw) : max par pixel, pour qu'une polyligne faite de
 * nombreux petits segments ne s'assombrisse pas à chaque jointure. */
static void seg_cov(Canvas *c, double x0, double y0, double x1, double y1, double hw)
{
    double pad = hw + 1.0;
    int xa = (int)floor(fmin(x0, x1) - pad), xb = (int)ceil(fmax(x0, x1) + pad);
    int ya = (int)floor(fmin(y0, y1) - pad), yb = (int)ceil(fmax(y0, y1) + pad);
    if (xb < 0 || yb < 0 || xa >= c->w || ya >= c->h)
        return;
    xa = xa < 0 ? 0 : xa;
    ya = ya < 0 ? 0 : ya;
    xb = xb >= c->w ? c->w - 1 : xb;
    yb = yb >= c->h ? c->h - 1 : yb;
    double dx = x1 - x0, dy = y1 - y0, l2 = dx * dx + dy * dy;
    for (int y = ya; y <= yb; y++)
        for (int x = xa; x <= xb; x++) {
            double px = x + 0.5 - x0, py = y + 0.5 - y0;
            double t = l2 > 0 ? (px * dx + py * dy) / l2 : 0;
            t = t < 0 ? 0 : t > 1 ? 1 : t;
            double ex = px - t * dx, ey = py - t * dy;
            double cv = hw + 0.5 - sqrt(ex * ex + ey * ey);
            if (cv <= 0)
                continue;
            float *q = &c->cov[(size_t)y * (size_t)c->w + (size_t)x];
            if (cv > *q)
                *q = (float)(cv > 1 ? 1 : cv);
        }
}
 
static void cov_clear(Canvas *c) { memset(c->cov, 0, (size_t)c->w * (size_t)c->h * sizeof(float)); }
 
static void cov_paint(Canvas *c, uint32_t col, double alpha)
{
    for (int y = 0; y < c->h; y++)
        for (int x = 0; x < c->w; x++) {
            float v = c->cov[(size_t)y * (size_t)c->w + (size_t)x];
            if (v > 0)
                blend(c, x, y, col, v * alpha);
        }
}
 
/* fichiers de polylignes : u32 nb ; par ligne : u32 nb de points, f32 lon, lat */
static void draw_polylines(Canvas *c, const char *path, double lw_px, uint32_t col, double alpha)
{
    size_t len, p = 0;
    uint8_t *b = read_file(path, 256 << 20, &len);
    uint32_t n;
    if (len < 4)
        die("%s : tronqué", path);
    memcpy(&n, b, 4);
    p = 4;
    cov_clear(c);
    for (uint32_t k = 0; k < n; k++) {
        uint32_t np;
        if (p + 4 > len)
            die("%s : tronqué", path);
        memcpy(&np, b + p, 4);
        p += 4;
        if ((size_t)np * 8 > len - p)
            die("%s : tronqué", path);
        float prev[2] = {0, 0};
        for (uint32_t q = 0; q < np; q++) {
            float pt[2];
            memcpy(pt, b + p + (size_t)q * 8, 8);
            if (q) {
                double x0, y0, x1, y1;
                proj(c, prev[0], prev[1], &x0, &y0);
                proj(c, pt[0], pt[1], &x1, &y1);
                seg_cov(c, x0, y0, x1, y1, lw_px / 2);
            }
            prev[0] = pt[0];
            prev[1] = pt[1];
        }
        p += (size_t)np * 8;
    }
    cov_paint(c, col, alpha);
    free(b);
}
 
/* ------------------------------------------------------------ texte */
static size_t utf8_next(const char *s, size_t i, uint32_t *cp)
{
    const unsigned char *u = (const unsigned char *)s;
    if (u[i] < 0x80) {
        *cp = u[i];
        return i + 1;
    }
    int n = (u[i] & 0xe0) == 0xc0 ? 1 : (u[i] & 0xf0) == 0xe0 ? 2 : (u[i] & 0xf8) == 0xf0 ? 3 : 0;
    uint32_t v = u[i] & (0x3f >> n);
    size_t j = i + 1;
    for (int k = 0; k < n; k++, j++) {
        if ((u[j] & 0xc0) != 0x80) {
            *cp = 0xfffd;
            return j;
        }
        v = v << 6 | (u[j] & 0x3f);
    }
    *cp = n ? v : 0xfffd;
    return j;
}
 
typedef struct {
    int w, h, base;      /* largeur, hauteur, ligne de base depuis le haut */
    float *a;            /* couverture 0..1 */
} TextMask;
 
static TextMask text_mask(const char *s, double size_px, int bold)
{
    FT_Face f = bold ? face_bold : face_reg;
    if (FT_Set_Pixel_Sizes(f, 0, (FT_UInt)lround(size_px)))
        die("taille de police refusée");
    int asc = (int)(f->size->metrics.ascender >> 6), desc = (int)(-f->size->metrics.descender >> 6);
    int width = 0;
    for (size_t i = 0; s[i];) {
        uint32_t cp;
        i = utf8_next(s, i, &cp);
        if (!FT_Load_Char(f, cp, FT_LOAD_DEFAULT))
            width += (int)(f->glyph->advance.x >> 6);
    }
    TextMask m = {width + 2, asc + desc + 2, asc + 1, NULL};
    m.a = xcalloc((size_t)m.w * (size_t)m.h, sizeof(float));
    int pen = 1;
    for (size_t i = 0; s[i];) {
        uint32_t cp;
        i = utf8_next(s, i, &cp);
        if (FT_Load_Char(f, cp, FT_LOAD_RENDER))
            continue;
        FT_GlyphSlot g = f->glyph;
        for (unsigned r = 0; r < g->bitmap.rows; r++)
            for (unsigned q = 0; q < g->bitmap.width; q++) {
                int x = pen + g->bitmap_left + (int)q, y = m.base - g->bitmap_top + (int)r;
                if (x < 0 || y < 0 || x >= m.w || y >= m.h)
                    continue;
                float v = g->bitmap.buffer[r * (unsigned)g->bitmap.pitch + q] / 255.0f;
                float *d = &m.a[(size_t)y * (size_t)m.w + (size_t)x];
                if (v > *d)
                    *d = v;
            }
        pen += (int)(g->advance.x >> 6);
    }
    return m;
}
 
/* ancre : ax 0 gauche, 1 centre, 2 droite ; ay 0 haut, 1 milieu, 2 bas */
static void draw_text(Canvas *c, double x, double y, int ax, int ay, const char *s, double size_px, int bold,
                      uint32_t col, uint32_t halo, double halo_px, uint32_t *box_bg, double box_pad)
{
    TextMask m = text_mask(s, size_px, bold);
    int ox = (int)lround(x - (ax == 0 ? 0 : ax == 1 ? m.w / 2.0 : m.w));
    int oy = (int)lround(y - (ay == 0 ? 0 : ay == 1 ? m.h / 2.0 : m.h));
    if (box_bg) {
        int pad = (int)lround(box_pad);
        for (int yy = oy - pad; yy < oy + m.h + pad; yy++)
            for (int xx = ox - pad; xx < ox + m.w + pad; xx++)
                blend(c, xx, yy, *box_bg, 1.0);
    }
    if (halo_px > 0) {
        /* halo : dilatation de la couverture sur un disque, comme un contour de texte */
        int r = (int)ceil(halo_px);
        for (int yy = -r; yy < m.h + r; yy++)
            for (int xx = -r; xx < m.w + r; xx++) {
                float best = 0;
                for (int dy = -r; dy <= r; dy++)
                    for (int dx = -r; dx <= r; dx++) {
                        double d = sqrt((double)(dx * dx + dy * dy));
                        if (d > halo_px + 0.5)
                            continue;
                        int sx = xx + dx, sy = yy + dy;
                        if (sx < 0 || sy < 0 || sx >= m.w || sy >= m.h)
                            continue;
                        float v = m.a[(size_t)sy * (size_t)m.w + (size_t)sx] * (float)fmin(1.0, halo_px + 0.5 - d);
                        if (v > best)
                            best = v;
                    }
                blend(c, ox + xx, oy + yy, halo, best);
            }
    }
    for (int yy = 0; yy < m.h; yy++)
        for (int xx = 0; xx < m.w; xx++)
            blend(c, ox + xx, oy + yy, col, m.a[(size_t)yy * (size_t)m.w + (size_t)xx]);
    free(m.a);
}
 
static void disc(Canvas *c, double x, double y, double r, uint32_t fill, double edge_px, uint32_t edge)
{
    double R = r + edge_px;
    for (int yy = (int)floor(y - R - 1); yy <= (int)ceil(y + R + 1); yy++)
        for (int xx = (int)floor(x - R - 1); xx <= (int)ceil(x + R + 1); xx++) {
            double d = hypot(xx + 0.5 - x, yy + 0.5 - y);
            if (edge_px > 0)
                blend(c, xx, yy, edge, R + 0.5 - d);
            blend(c, xx, yy, fill, r + 0.5 - d);
        }
}
 
/* ------------------------------------------------------------ champ coloré */
static uint32_t gust_color(double v)
{
    int i = 0;
    for (int k = 0; k < 30; k++)
        if (v >= GB[k])
            i = k;
    return GC[i];
}
 
static void draw_field(Canvas *c, const Grd *d, int layer, int n100, int nmodels)
{
    if (layer < 0 || layer >= d->nt_gust)
        die("couche %d absente", layer);
    size_t np = grid_points(&d->g);
    const float *f = d->gust + (size_t)layer * np;
    for (int y = 0; y < c->h; y++) {
        double lat = c->lat1 - (y + 0.5) / c->h * (c->lat1 - c->lat0);
        for (int x = 0; x < c->w; x++) {
            double lon = c->lon0 + (x + 0.5) / c->w * (c->lon1 - c->lon0);
            double v = bilinear(&d->g, f, lat, lon);
            uint32_t col;
            if (isnan(v))
                col = BAD;
            else if (n100) {
                /* nombre de modèles : valeur entière la plus proche du point de grille */
                double fy = (lat - d->g.lat0) / d->g.dlat, fx = (lon - d->g.lon0) / d->g.dlon;
                int k = (int)lround(f[(size_t)lround(fy) * (size_t)d->g.nx + (size_t)lround(fx)]);
                k = k < 0 ? 0 : k > nmodels ? nmodels : k > 6 ? 6 : k;
                col = N100C[k];
            } else
                col = gust_color(v);
            uint8_t *p = c->rgb + 3 * ((size_t)y * (size_t)c->w + (size_t)x);
            p[0] = (uint8_t)(col >> 16);
            p[1] = (uint8_t)(col >> 8);
            p[2] = (uint8_t)col;
        }
    }
}
 
/* ------------------------------------------------------------ isobares */
typedef struct {
    double x, y;
    int level;
} LabelPos;
 
static void draw_isobars(Canvas *c, const Grd *d, int idx, int smooth)
{
    if (idx < 0 || idx >= d->nt_msl)
        die("échéance de pression %d absente", idx);
    const Grid *g = &d->g;
    size_t np = grid_points(g);
    float *P = xmalloc(np, sizeof(float));
    if (smooth)   /* AROME : moyenne normalisée, les trous ne comptent pas (voir field.c) */
        gauss_masked(d->msl + (size_t)idx * np, P, g->nx, g->ny, 2.0);
    else
        memcpy(P, d->msl + (size_t)idx * np, np * sizeof(float));
    /* cellules utiles : le cadre plus un degré de marge */
    int x0 = (int)floor((c->lon0 - 1 - g->lon0) / g->dlon), x1 = (int)ceil((c->lon1 + 1 - g->lon0) / g->dlon);
    int y0 = (int)floor((c->lat0 - 1 - g->lat0) / g->dlat), y1 = (int)ceil((c->lat1 + 1 - g->lat0) / g->dlat);
    x0 = x0 < 0 ? 0 : x0;
    y0 = y0 < 0 ? 0 : y0;
    x1 = x1 > g->nx - 1 ? g->nx - 1 : x1;
    y1 = y1 > g->ny - 1 ? g->ny - 1 : y1;
    LabelPos lab[256];
    int nlab = 0, per_level[64] = {0};
    cov_clear(c);
    for (int j = y0; j < y1; j++)
        for (int i = x0; i < x1; i++) {
            size_t k = (size_t)j * (size_t)g->nx + (size_t)i;
            double v[4] = {P[k], P[k + 1], P[k + 1 + (size_t)g->nx], P[k + (size_t)g->nx]};
            if (isnan(v[0]) || isnan(v[1]) || isnan(v[2]) || isnan(v[3]))
                continue;
            double lo = fmin(fmin(v[0], v[1]), fmin(v[2], v[3])), hi = fmax(fmax(v[0], v[1]), fmax(v[2], v[3]));
            /* coins (i,j) (i+1,j) (i+1,j+1) (i,j+1) en coordonnées géographiques */
            double cx[4] = {i, i + 1.0, i + 1.0, i}, cy[4] = {j, j, j + 1.0, j + 1.0};
            for (int L = 940; L <= 1048; L += 2) {
                if (L < lo || L > hi)
                    continue;
                int cs = (v[0] > L) | (v[1] > L) << 1 | (v[2] > L) << 2 | (v[3] > L) << 3;
                if (cs == 0 || cs == 15)
                    continue;
                /* points de passage sur les arêtes 0:(0-1) 1:(1-2) 2:(2-3) 3:(3-0) */
                double ex[4], ey[4];
                for (int e = 0; e < 4; e++) {
                    int a = e, b = (e + 1) & 3;
                    double t = (L - v[a]) / (v[b] - v[a]);
                    t = !isfinite(t) ? 0.5 : t < 0 ? 0 : t > 1 ? 1 : t;
                    double gx = cx[a] + t * (cx[b] - cx[a]), gy = cy[a] + t * (cy[b] - cy[a]);
                    proj(c, g->lon0 + gx * g->dlon, g->lat0 + gy * g->dlat, &ex[e], &ey[e]);
                }
                int segs[2][2], ns = 0;
                switch (cs) {
                case 1: case 14: segs[0][0] = 3; segs[0][1] = 0; ns = 1; break;
                case 2: case 13: segs[0][0] = 0; segs[0][1] = 1; ns = 1; break;
                case 3: case 12: segs[0][0] = 3; segs[0][1] = 1; ns = 1; break;
                case 4: case 11: segs[0][0] = 1; segs[0][1] = 2; ns = 1; break;
                case 6: case 9: segs[0][0] = 0; segs[0][1] = 2; ns = 1; break;
                case 7: case 8: segs[0][0] = 2; segs[0][1] = 3; ns = 1; break;
                case 5: case 10: {
                    /* col : on tranche avec la moyenne des quatre coins */
                    int center_hi = (v[0] + v[1] + v[2] + v[3]) / 4 > L;
                    if ((cs == 5) == center_hi) {
                        segs[0][0] = 3; segs[0][1] = 2; segs[1][0] = 0; segs[1][1] = 1;
                    } else {
                        segs[0][0] = 3; segs[0][1] = 0; segs[1][0] = 1; segs[1][1] = 2;
                    }
                    ns = 2;
                    break;
                }
                }
                for (int s = 0; s < ns; s++) {
                    double ax = ex[segs[s][0]], ay = ey[segs[s][0]], bx = ex[segs[s][1]], by = ey[segs[s][1]];
                    seg_cov(c, ax, ay, bx, by, 0.7 * PT / 2);
                    /* étiquettes : niveaux multiples de 4, espacées d'au moins 180 px, loin des bords */
                    int li = (L - 940) / 2;
                    double mx = (ax + bx) / 2, my = (ay + by) / 2;
                    if (L % 4 || nlab >= 256 || per_level[li] >= 4 || mx < 40 || my < 60 || mx > c->w - 40 ||
                        my > c->h - 30)
                        continue;
                    int okp = 1;
                    for (int q = 0; q < nlab && okp; q++)
                        okp = hypot(lab[q].x - mx, lab[q].y - my) > 180;
                    if (okp) {
                        lab[nlab++] = (LabelPos){mx, my, L};
                        per_level[li]++;
                    }
                }
            }
        }
    cov_paint(c, 0x1a1a1a, 1.0);
    for (int q = 0; q < nlab; q++) {
        char t[8];
        snprintf(t, sizeof t, "%d", lab[q].level);
        draw_text(c, lab[q].x, lab[q].y, 1, 1, t, 8.5 * PT, 0, 0x111111, 0xffffff, 2.0 * PT / 2, NULL, 0);
    }
    free(P);
}
 
/* ------------------------------------------------------------ fichiers tsv */
typedef struct {
    double lon, lat;
    char text[64];
} Row;
 
static Row *read_tsv(const char *path, int *n)
{
    size_t len;
    char *b = (char *)read_file(path, 16 << 20, &len);
    Row *r = NULL;
    int cnt = 0;
    char *cur_ = b;
    for (char *line = next_line(&cur_); line; line = next_line(&cur_)) {
        Row x;
        char *t1 = strchr(line, '\t'), *t2 = t1 ? strchr(t1 + 1, '\t') : NULL;
        if (!t1 || !t2 || strlen(t2 + 1) >= sizeof x.text)
            die("%s : ligne invalide", path);
        char *end;
        x.lon = strtod(line, &end);
        if (end != t1)
            die("%s : longitude invalide", path);
        x.lat = strtod(t1 + 1, &end);
        if (end != t2)
            die("%s : latitude invalide", path);
        strcpy(x.text, t2 + 1);
        r = xrealloc(r, (size_t)cnt + 1, sizeof *r);
        r[cnt++] = x;
    }
    free(b);
    *n = cnt;
    return r;
}
 
/* ------------------------------------------------------------ champ interpolé entre stations */
static void draw_idw(Canvas *c, const char *path, double radius_km)
{
    int n;
    Row *r = read_tsv(path, &n);
    if (n == 0) {
        free(r);
        return;
    }
    /* grille de recherche : cellules de 0,25° pour ne regarder que les stations voisines */
    const double cell = 0.25;
    int gx = (int)ceil((c->lon1 - c->lon0) / cell) + 1, gy = (int)ceil((c->lat1 - c->lat0) / cell) + 1;
    int *head = xmalloc((size_t)gx * (size_t)gy, sizeof(int)), *next = xmalloc((size_t)n, sizeof(int));
    double *val = xmalloc((size_t)n, sizeof(double));
    for (int k = 0; k < gx * gy; k++)
        head[k] = -1;
    for (int q = 0; q < n; q++) {
        val[q] = strtod(r[q].text, NULL);
        int cx = (int)floor((r[q].lon - c->lon0) / cell), cy = (int)floor((r[q].lat - c->lat0) / cell);
        if (cx < 0 || cy < 0 || cx >= gx || cy >= gy || !isfinite(val[q])) {
            next[q] = -2;
            continue;
        }
        next[q] = head[cy * gx + cx];
        head[cy * gx + cx] = q;
    }
    double coslat = cos((c->lat0 + c->lat1) / 2 * DEG);
    int reach = (int)ceil(radius_km / (KM_PER_DEG * coslat) / cell) + 1;
    for (int y = 0; y < c->h; y++) {
        double lat = c->lat1 - (y + 0.5) / c->h * (c->lat1 - c->lat0);
        for (int x = 0; x < c->w; x++) {
            double lon = c->lon0 + (x + 0.5) / c->w * (c->lon1 - c->lon0);
            int cx = (int)floor((lon - c->lon0) / cell), cy = (int)floor((lat - c->lat0) / cell);
            double sw = 0, sv = 0, dmin = INFINITY;
            for (int yy = cy - reach; yy <= cy + reach; yy++)
                for (int xx = cx - reach; xx <= cx + reach; xx++) {
                    if (xx < 0 || yy < 0 || xx >= gx || yy >= gy)
                        continue;
                    for (int q = head[yy * gx + xx]; q >= 0; q = next[q]) {
                        double d = dist_km(lat, lon, r[q].lat, r[q].lon);
                        if (d > radius_km)
                            continue;
                        double w = 1.0 / (d * d + 25.0);   /* 5 km de lissage au plus près d'une station */
                        sw += w;
                        sv += w * val[q];
                        if (d < dmin)
                            dmin = d;
                    }
                }
            if (sw <= 0)
                continue;
            /* fondu sur le dernier tiers du rayon, pour ne pas finir en bord franc */
            double a = (radius_km - dmin) / (radius_km * 0.35);
            blend(c, x, y, gust_color(sv / sw), a > 1 ? 1 : a);
        }
    }
    free(head);
    free(next);
    free(val);
    free(r);
}
 
/* ------------------------------------------------------------ bloc */
typedef struct {
    char out[1024], geo[1024], field[1024], iso[1024], labels[1024], dots[1024], image[1024], idw[1024];
    char title[256], sub[256], right[256], sign[128];
    int quality, width, layer, n100, nmodels, geo_depts, iso_idx, iso_smooth, has_center;
    double box[4], geo_lw, clon, clat, cp, idw_km;
    uint32_t bg, geocolor;
    int has_geocolor;
} Spec;
 
static Grd *cache_grd(const char *path)
{
    static char cpath[2][1024];
    static Grd *cg[2];
    static int next;
    for (int k = 0; k < 2; k++)
        if (cg[k] && !strcmp(cpath[k], path))
            return cg[k];
    int k = next;
    next = 1 - next;
    grd_free(cg[k]);
    cg[k] = grd_read(path);
    snprintf(cpath[k], sizeof cpath[k], "%s", path);
    return cg[k];
}
 
static void render_one(const Spec *s)
{
    if (s->width < 100 || s->width > 4000 || !(s->box[1] > s->box[0]) || !(s->box[3] > s->box[2]))
        die("bloc invalide pour %s", s->out);
    Canvas c;
    c.lon0 = s->box[0];
    c.lon1 = s->box[1];
    c.lat0 = s->box[2];
    c.lat1 = s->box[3];
    if (s->image[0]) {
        /* fond imposé (satellite) : sa taille fixe celle de la carte, projection lon/lat simple */
        size_t len;
        uint8_t *raw = read_file(s->image, 64 << 20, &len);
        int iw = 0, ih = 0;
        if (!WebPGetInfo(raw, len, &iw, &ih) || iw < 100 || ih < 50 || iw > 4000 || ih > 4000)
            die("image de fond illisible : %s", s->image);
        uint8_t *px = WebPDecodeRGB(raw, len, &iw, &ih);
        if (!px)
            die("image de fond illisible : %s", s->image);
        c.w = iw;
        c.h = ih;
        c.rgb = xmalloc((size_t)c.w * (size_t)c.h, 3);
        memcpy(c.rgb, px, (size_t)c.w * (size_t)c.h * 3);
        WebPFree(px);
        free(raw);
    } else {
        double k = cos((c.lat0 + c.lat1) / 2 * DEG);
        c.w = s->width;
        c.h = (int)lround(c.w * (c.lat1 - c.lat0) / ((c.lon1 - c.lon0) * k));
        if (c.h < 50 || c.h > 4000)
            die("hauteur d'image invalide");
        c.rgb = xmalloc((size_t)c.w * (size_t)c.h, 3);
        for (size_t q = 0; q < (size_t)c.w * (size_t)c.h; q++) {
            c.rgb[3 * q] = (uint8_t)(s->bg >> 16);
            c.rgb[3 * q + 1] = (uint8_t)(s->bg >> 8);
            c.rgb[3 * q + 2] = (uint8_t)s->bg;
        }
    }
    c.cov = xmalloc((size_t)c.w * (size_t)c.h, sizeof(float));
    if (s->field[0])
        draw_field(&c, cache_grd(s->field), s->layer, s->n100, s->nmodels);
    if (s->idw[0])
        draw_idw(&c, s->idw, s->idw_km);
    if (s->geo[0]) {
        char p[1100];
        uint32_t line = s->has_geocolor ? s->geocolor : 0x111111, dept = s->has_geocolor ? s->geocolor : 0x3a3a3a;
        if (s->geo_depts) {
            snprintf(p, sizeof p, "%s/dept.bin", s->geo);
            draw_polylines(&c, p, 0.35 * PT * s->geo_lw, dept, s->has_geocolor ? 0.55 : 0.75);
        }
        snprintf(p, sizeof p, "%s/border.bin", s->geo);
        draw_polylines(&c, p, 0.9 * PT * s->geo_lw, line, 1.0);
        snprintf(p, sizeof p, "%s/coast.bin", s->geo);
        draw_polylines(&c, p, 0.8 * PT * s->geo_lw, line, 1.0);
    }
    if (s->iso[0])
        draw_isobars(&c, cache_grd(s->iso), s->iso_idx, s->iso_smooth);
    if (s->has_center) {
        double x, y;
        proj(&c, s->clon, s->clat, &x, &y);
        draw_text(&c, x, y, 1, 1, "D", 19 * PT, 1, 0xcc0000, 0xffffff, 3 * PT / 2, NULL, 0);
        double x2, y2;
        proj(&c, s->clon, s->clat - 0.42, &x2, &y2);
        char t[16];
        snprintf(t, sizeof t, "%d", (int)lround(s->cp));
        draw_text(&c, x2, y2, 1, 0, t, 11 * PT, 1, 0xcc0000, 0xffffff, 3 * PT / 2, NULL, 0);
    }
    if (s->dots[0]) {
        int n;
        Row *r = read_tsv(s->dots, &n);
        for (int q = 0; q < n; q++) {   /* fichier déjà trié du plus faible au plus fort */
            double x, y;
            proj(&c, r[q].lon, r[q].lat, &x, &y);
            disc(&c, x, y, sqrt(58.0) * PT / 2 - 0.6, gust_color(strtod(r[q].text, NULL)), 0.6 * PT, 0x111111);
        }
        free(r);
    }
    if (s->labels[0]) {
        int n;
        Row *r = read_tsv(s->labels, &n);
        for (int q = 0; q < n; q++) {
            double x, y;
            proj(&c, r[q].lon, r[q].lat, &x, &y);
            if (strncmp(r[q].text, "*", 1) != 0)   /* « *texte » : étiquette sans point */
                disc(&c, x, y, 3.2 * PT / 2, 0x111111, 0.8 * PT / 2, 0xffffff);
            double tx, ty;
            proj(&c, r[q].lon + 0.07, r[q].lat + 0.05, &tx, &ty);
            draw_text(&c, tx, ty, 0, 2, r[q].text[0] == '*' ? r[q].text + 1 : r[q].text, 10.5 * PT, 1, 0x111111,
                      0xffffff, 2.6 * PT / 2, NULL, 0);
        }
        free(r);
    }
    if (s->title[0])
        draw_text(&c, 0.012 * c.w, 0.012 * c.h, 0, 0, s->title, 13 * PT, 1, 0x0b2a4a, 0xffffff, 3 * PT / 2, NULL, 0);
    if (s->sub[0])
        draw_text(&c, 0.012 * c.w, 0.048 * c.h, 0, 0, s->sub, 10.5 * PT, 0, 0x0b2a4a, 0xffffff, 3 * PT / 2, NULL, 0);
    if (s->right[0]) {
        uint32_t bgc = 0x111111;
        draw_text(&c, 0.988 * c.w - 5, 0.012 * c.h + 5, 2, 0, s->right, 10.5 * PT, 0, 0xffffff, 0, 0, &bgc, 5);
    }
    if (s->sign[0])
        draw_text(&c, 0.988 * c.w, 0.988 * c.h, 2, 2, s->sign, 9 * PT, 0, 0x111111, 0xffffff, 3 * PT / 2, NULL, 0);
    uint8_t *webp = NULL;
    size_t sz = WebPEncodeRGB(c.rgb, c.w, c.h, c.w * 3, (float)s->quality, &webp);
    if (!sz)
        die("encodage WebP en échec : %s", s->out);
    write_file_atomic(s->out, webp, sz);
    WebPFree(webp);
    free(c.rgb);
    free(c.cov);
}
 
static void copy_rest(char *dst, size_t n, const char *line, const char *key)
{
    const char *v = line + strlen(key);
    while (*v == ' ')
        v++;
    if (strlen(v) >= n)
        die("valeur trop longue pour %s", key);
    strcpy(dst, v);
}
 
int cmd_render(int argc, char **argv)
{
    if (argc != 1)
        die("usage : austen render <spec.txt>");
    if (FT_Init_FreeType(&ft) || FT_New_Face(ft, FONT_BOLD, 0, &face_bold) || FT_New_Face(ft, FONT_REG, 0, &face_reg))
        die("polices DejaVu introuvables");
    size_t len;
    char *buf = (char *)read_file(argv[0], 16 << 20, &len);
    Spec s;
    memset(&s, 0, sizeof s);
    s.quality = 85;
    s.width = 960;
    s.bg = 0xffffff;
    s.geo_lw = 1.0;
    int count = 0, lineno = 0;
    char *cur_ = buf;
    for (char *line = next_line(&cur_); line; line = next_line(&cur_)) {
        lineno++;
        char key[16] = "";
        sscanf(line, "%15s", key);
        int ok = 1;
        if (!strcmp(key, "out")) copy_rest(s.out, sizeof s.out, line, key);
        else if (!strcmp(key, "quality")) ok = sscanf(line, "%*s %d", &s.quality) == 1;
        else if (!strcmp(key, "width")) ok = sscanf(line, "%*s %d", &s.width) == 1;
        else if (!strcmp(key, "box")) ok = sscanf(line, "%*s %lf %lf %lf %lf", &s.box[0], &s.box[1], &s.box[2], &s.box[3]) == 4;
        else if (!strcmp(key, "background")) { char col[16]; ok = sscanf(line, "%*s %15s", col) == 1; if (ok) s.bg = parse_color(col); }
        else if (!strcmp(key, "geo")) ok = sscanf(line, "%*s %1023s %d %lf", s.geo, &s.geo_depts, &s.geo_lw) == 3;
        else if (!strcmp(key, "field")) {
            char pal[16];
            s.nmodels = 0;
            ok = sscanf(line, "%*s %1023s %d %15s %d", s.field, &s.layer, pal, &s.nmodels) >= 3;
            s.n100 = !strcmp(pal, "n100");
        }
        else if (!strcmp(key, "isobars")) ok = sscanf(line, "%*s %1023s %d %d", s.iso, &s.iso_idx, &s.iso_smooth) == 3;
        else if (!strcmp(key, "center")) ok = (s.has_center = sscanf(line, "%*s %lf %lf %lf", &s.clon, &s.clat, &s.cp) == 3);
        else if (!strcmp(key, "labels")) copy_rest(s.labels, sizeof s.labels, line, key);
        else if (!strcmp(key, "image")) copy_rest(s.image, sizeof s.image, line, key);
        else if (!strcmp(key, "idw")) ok = sscanf(line, "%*s %1023s %lf", s.idw, &s.idw_km) == 2 && s.idw_km > 0 && s.idw_km < 500;
        else if (!strcmp(key, "geocolor")) { char col[16]; ok = sscanf(line, "%*s %15s", col) == 1; if (ok) { s.geocolor = parse_color(col); s.has_geocolor = 1; } }
        else if (!strcmp(key, "dots")) copy_rest(s.dots, sizeof s.dots, line, key);
        else if (!strcmp(key, "title")) copy_rest(s.title, sizeof s.title, line, key);
        else if (!strcmp(key, "sub")) copy_rest(s.sub, sizeof s.sub, line, key);
        else if (!strcmp(key, "right")) copy_rest(s.right, sizeof s.right, line, key);
        else if (!strcmp(key, "sign")) copy_rest(s.sign, sizeof s.sign, line, key);
        else if (!strcmp(key, "end")) {
            if (!s.out[0])
                die("spec ligne %d : bloc sans sortie", lineno);
            render_one(&s);
            count++;
            int q = s.quality, w = s.width;
            uint32_t bg = s.bg;
            double lw = s.geo_lw;
            char geo[1024];
            int gd = s.geo_depts;
            memcpy(geo, s.geo, sizeof geo);
            memset(&s, 0, sizeof s);
            s.quality = q;      /* réglages d'aspect conservés d'un bloc au suivant */
            s.width = w;
            s.bg = bg;
            s.geo_lw = lw;
            s.geo_depts = gd;
            memcpy(s.geo, geo, sizeof geo);
        }
        else ok = 0;
        if (!ok)
            die("spec ligne %d invalide : %.60s", lineno, line);
    }
    free(buf);
    printf("{\"images\":%d}\n", count);
    return 0;
}
[code 10]

Le chef d'orchestre

pipeline.php, 566 lignes

Toutes les 20 minutes : téléchargement, décodage, analyse, puis moyenne et cartes seulement si un nouveau run est arrivé. Les cartes sont rendues en parallèle, un processus par modèle (moins de 20 secondes en tout). Il assemble data.json pour la page et publie fichier par fichier.

pipeline.php
Fichier brut
<?php
/**
 * Chaîne complète : pipeline.php [--runs modèle:run ...] [--no-fetch] [--force]
 * téléchargement (fetch.php) → décodage, analyse, moyenne, cartes (binaire C) → data.json → publication.
 * Les cartes ne sont refaites que si un nouveau run est arrivé (--force pour les refaire quand même).
 */
declare(strict_types=1);
require __DIR__ . '/lib.php';
 
const FR_BOX = [-7.5, 10.0, 41.2, 52.4];      /* cartes de rafales */
const AN_BOX = [-12.0, 11.0, 41.5, 55.0];     /* animation : on voit arriver la dépression */
const SHOW = ['Ouessant', 'Brest', 'Lorient', 'Belle-Île', 'Saint-Brieuc', 'Rennes', 'Nantes', 'La Hague', 'Caen',
              'Granville', 'Cap de la Hève', 'Dieppe', 'Rouen', 'Boulogne', 'Lille', 'Paris', 'Le Mans', 'Tours',
              'Orléans', 'Reims', 'Chassiron', 'Pointe du Raz', 'Dunkerque', 'Beauvais', 'Alençon'];
const JOURS = ['dimanche', 'lundi', 'mardi', 'mercredi', 'jeudi', 'vendredi', 'samedi'];
 
/* heure de Paris (UTC+2 jusqu'au 25/10) : « jeudi 22 à 20 h » */
function fr_time(int $t): string
{
    $d = (new DateTimeImmutable('@' . $t))->setTimezone(new DateTimeZone('Europe/Paris'));
    return JOURS[(int) $d->format('w')] . ' ' . $d->format('j') . ' à ' . $d->format('G') . ' h';
}
 
function ucfirst_fr(string $s): string { return mb_strtoupper(mb_substr($s, 0, 1)) . mb_substr($s, 1); }
 
/* ---------------------------------------------------------------- étapes C */
function grd_path(string $m, string $r): string { return data_dir() . "/grd/$m/$r.grd"; }
function derived(string $m, string $r): string { return data_dir() . "/derived/$m/$r.json"; }
 
function process(string $model, string $run, string $conf): bool
{
    $raw = data_dir() . "/raw/$model/$run";
    $grd = grd_path($model, $run);
    $js = derived($model, $run);
    ensure_dir(dirname($grd));
    ensure_dir(dirname($js));
    if (!is_file($grd) && !is_file($js)) {
        if (!is_file("$raw/.complete")) {
            return false;
        }
        [$code, , $err] = run([BIN, 'decode', $raw, $grd], 900);
        if ($code !== 0) {
            logm('ÉCHEC décodage', $model, $run, trim($err));
            return false;
        }
    }
    if (is_file($grd) && (!is_file($js) || filemtime($js) < filemtime($grd) || filemtime($js) < filemtime(BIN))) {
        [$code, , $err] = run([BIN, 'analyze', $conf, $grd, $js, data_dir() . "/derived/$model/{$run}_gmax.grd"], 900);
        if ($code !== 0) {
            logm('ÉCHEC analyse', $model, $run, trim($err));
            return false;
        }
    }
    foreach (glob("$raw/*.grib2") ?: [] as $f) {   /* les GRIB bruts ne servent plus */
        unlink($f);
    }
    if (is_dir($raw)) {
        touch("$raw/.decoded");
    }
    return true;
}
 
/** Runs analysés (le .json suffit : la grille d'un ancien run est effacée, sa trajectoire reste). */
function decoded_runs(string $m): array
{
    $out = [];
    foreach (glob(data_dir() . "/derived/$m/*.json") ?: [] as $p) {
        $r = basename($p, '.json');
        if (valid_run($r)) {
            $out[] = $r;
        }
    }
    sort($out);
    return $out;
}
 
/** Garde la grille des seuls runs encore affichés (dernier run, run de prévision) : 300 Mo par run AROME. */
function prune_grids(array $keep): void
{
    foreach (glob(data_dir() . '/grd/*/*.grd') ?: [] as $p) {
        $m = basename(dirname($p));
        $r = basename($p, '.grd');
        if (!in_array($r, $keep[$m] ?? [], true)) {
            unlink($p);
        }
    }
}
 
/**
 * Runs retenus pour la « moyenne de tous les runs » : partis avant le début de l'épisode et allant au moins
 * jusqu'à 6 h avant sa fin (les plus anciens s'arrêtent parfois un peu avant).
 */
function allruns_selection(array $all): array
{
    $out = [];
    foreach ($all as $m => $runs) {
        if (!MODELS[$m]['gust']) {
            continue;
        }
        foreach ($runs as $r) {
            if (run_time($r) > ep0() || !is_file(data_dir() . "/derived/$m/{$r}_gmax.grd")) {
                continue;
            }
            $d = json_in(derived($m, $r));
            $last = $d['gust_times'] ? utc(end($d['gust_times'])) : 0;
            if ($last >= ep1() - 6 * 3600) {
                $out[$m][] = $r;
            }
        }
    }
    return $out;
}
 
/* ---------------------------------------------------------------- rendu */
/** Lance plusieurs « austen render » en parallèle, un par fichier de description. */
function render_parallel(array $specs): void
{
    $procs = [];
    foreach ($specs as $spec) {
        $p = proc_open([BIN, 'render', $spec], [0 => ['file', '/dev/null', 'r'], 1 => ['file', '/dev/null', 'w'], 2 => ['pipe', 'w']], $pipes);
        if (!is_resource($p)) {
            throw new RuntimeException('rendu impossible');
        }
        $procs[] = [$p, $pipes[2], $spec];
    }
    foreach ($procs as [$p, $err, $spec]) {
        $e = stream_get_contents($err);
        fclose($err);
        if (proc_close($p) !== 0) {
            logm('ÉCHEC rendu', basename($spec), trim((string) $e));
        }
    }
}
 
function spec_block(array $kv): string
{
    $s = '';
    foreach ($kv as $k => $v) {
        $v = str_replace(["\n", "\r"], ' ', (string) $v);
        $s .= "$k $v\n";
    }
    return $s . "end\n";
}
 
function grd_times(string $path): array
{
    $h = fopen($path, 'rb');
    $head = fread($h, 56);
    ['nt_msl' => $nm, 'nt_gust' => $ng] = unpack('a4magic/Vver/Vnx/Vny/Vnt_msl/Vnt_gust', $head);
    $tm = $nm ? array_values(unpack("q$nm", fread($h, 8 * $nm))) : [];
    $tg = $ng ? array_values(unpack("q$ng", fread($h, 8 * $ng))) : [];
    fclose($h);
    return [$tm, $tg];
}
 
function runs_label(array $runs): string
{
    $runs = array_values(array_unique($runs));
    sort($runs);
    $f = fn($r) => substr($r, 6, 2) . '/' . substr($r, 4, 2) . ' ' . substr($r, 8, 2) . ' UTC';
    $label = 'Runs du ' . implode(', ', array_map($f, $runs));
    return mb_strlen($label) > 48 ? 'Runs du ' . $f($runs[0]) . ' au ' . $f(end($runs)) : $label;
}
 
function render_all(array $latest, array $forecast, array $stations, array $allSel): array
{
    $tmp = data_dir() . '/render';
    ensure_dir($tmp);
    $site = site_dir();
    ensure_dir("$site/maps");
    $geo = data_dir() . '/geo';
    $ens = data_dir() . '/derived/ensemble';
    $gustFc = array_filter($forecast, fn($m) => MODELS[$m]['gust'], ARRAY_FILTER_USE_KEY);
    $label = runs_label(array_values($gustFc));
    $period = 'Du ' . fr_time(ep0()) . ' au ' . fr_time(ep1()) . ' (heure de Paris)';
    $common = ['box' => implode(' ', FR_BOX), 'geo' => "$geo 1 1", 'sign' => 'Meteociel.fr'];
 
    /* étiquettes : moyenne des modèles aux stations retenues, puis chaque modèle */
    $byName = array_column($stations, null, 'name');
    $lab = fn(callable $val) => implode('', array_map(fn($n) => $val($byName[$n]) === null ? ''
        : sprintf("%.6f\t%.6f\t%d\n", $byName[$n]['lon'], $byName[$n]['lat'], $val($byName[$n])), array_filter(SHOW, fn($n) => isset($byName[$n]))));
    atomic_write("$tmp/lab_mean.tsv", $lab(fn($s) => $s['mean']));
 
    $maps = spec_block(['out' => "$site/maps/gust-mean.webp", 'width' => 1200, 'quality' => 88] + $common + [
        'field' => "$ens.grd 0 gust", 'labels' => "$tmp/lab_mean.tsv",
        'title' => 'Rafales maxi de l\'épisode, moyenne de ' . count($gustFc) . ' modèles (km/h)', 'sub' => $period, 'right' => $label]);
    $maps .= spec_block(['out' => "$site/maps/gust-n100.webp", 'field' => "$ens.grd 3 n100 " . count($gustFc),
        'title' => 'Combien de modèles voient des rafales de 100 km/h ou plus ?', 'right' => $label,
        'sub' => 'Sur ' . count($gustFc) . ' modèles, du ' . fr_time(ep0()) . ' au ' . fr_time(ep1())] + $common);
    if ($allSel && is_file(data_dir() . '/derived/ensemble_all.grd')) {
        $nr = array_sum(array_map('count', $allSel));
        atomic_write("$tmp/lab_all.tsv", $lab(fn($s) => $s['all'] ?? null));
        $maps .= spec_block(['out' => "$site/maps/gust-allruns.webp", 'width' => 1200, 'field' => data_dir() . '/derived/ensemble_all.grd 0 gust',
            'labels' => "$tmp/lab_all.tsv", 'title' => "Rafales maxi de l'épisode, moyenne de $nr runs de " . count($allSel) . ' modèles (km/h)',
            'sub' => $period, 'right' => 'Tous les runs partis avant jeudi 20 h'] + $common);
    }
    $layer = 4;
    foreach ($gustFc as $m => $r) {
        atomic_write("$tmp/lab_$m.tsv", $lab(fn($s) => isset($s['models'][$m]) ? $s['models'][$m] : null));
        $maps .= spec_block(['out' => "$site/maps/gust-$m.webp", 'width' => 1000, 'field' => "$ens.grd $layer gust",
            'labels' => "$tmp/lab_$m.tsv", 'title' => MODELS[$m]['name'] . ' : rafales maxi de l\'épisode (km/h)',
            'sub' => MODELS[$m]['gnote'], 'right' => 'Run du ' . substr($r, 6, 2) . '/' . substr($r, 4, 2) . ' ' . substr($r, 8, 2) . ' UTC'] + $common);
        $layer++;
    }
    atomic_write("$tmp/maps.txt", $maps);
    $specs = ["$tmp/maps.txt"];
 
    /* animation : une image par heure de 12 UTC jeudi à 0 UTC samedi, un processus par modèle */
    $frames = [];
    foreach ($latest as $m => $r) {
        if (!MODELS[$m]['gust']) {
            continue;
        }
        $grd = grd_path($m, $r);
        [$tm, $tg] = grd_times($grd);
        $track = [];
        foreach (json_in(derived($m, $r))['track'] as $p) {
            $track[utc($p['t'])] = $p;
        }
        $dir = "$site/frames/$m";
        ensure_dir($dir);
        $txt = '';
        $times = [];
        $rt = run_time($r);
        for ($t = utc('2026-10-22T12:00Z'); $t <= utc('2026-10-24T00:00Z'); $t += 3600) {
            $im = array_search($t, $tm, true);
            $ig = array_search($t, $tg, true);
            if ($im === false || $ig === false) {
                continue;
            }
            $b = ['out' => "$dir/" . gmdate('YmdH', $t) . '.webp', 'width' => 960, 'quality' => 82, 'box' => implode(' ', AN_BOX),
                  'geo' => "$geo 1 0.9", 'field' => "$grd $ig gust", 'isobars' => "$grd $im " . ($m === 'arome' ? 1 : 0)];
            if (isset($track[$t])) {
                $b['center'] = sprintf('%.4f %.4f %.1f', $track[$t]['lon'], $track[$t]['lat'], $track[$t]['p']);
            }
            $b += ['title' => MODELS[$m]['name'] . ' : rafales et pression',
                   'sub' => ucfirst_fr(fr_time($t)) . ' (heure de Paris)',
                   'right' => sprintf('Run %s UTC du %s/%s, +%dh', substr($r, 8, 2), substr($r, 6, 2), substr($r, 4, 2), intdiv($t - $rt, 3600)),
                   'sign' => 'Meteociel.fr'];
            $txt .= spec_block($b);
            $times[] = gmdate('YmdH', $t);
        }
        atomic_write("$tmp/frames_$m.txt", $txt);
        $specs[] = "$tmp/frames_$m.txt";
        $frames[$m] = ['run' => $r, 'times' => $times];
    }
    render_parallel($specs);
    json_out("$site/frames/frames.json", $frames);
    return $frames;
}
 
/* ---------------------------------------------------------------- data.json */
function station_table(array $forecast, array $allSel = []): array
{
    /* moyenne de tous les runs : maxi de chaque run moyenné par modèle, puis entre modèles */
    $allMax = [];
    foreach ($allSel as $m => $runs) {
        foreach ($runs as $r) {
            foreach (json_in(derived($m, $r))['stations'] ?? [] as $sid => $v) {
                if ($v['max'] !== null) {
                    $allMax[(string) $sid][$m][] = $v['max'];
                }
            }
        }
    }
    $per = [];
    foreach ($forecast as $m => $r) {
        if (MODELS[$m]['gust']) {
            $per[$m] = json_in(derived($m, $r))['stations'] ?? [];
        }
    }
    $out = [];
    foreach (STATIONS as [$id, $name, $lat, $lon]) {
        $models = $tmax = [];
        foreach ($per as $m => $st) {
            $v = $st[$id] ?? null;
            if ($v !== null && $v['max'] !== null) {
                $models[$m] = (int) round($v['max']);
                $tmax[$m] = $v['tmax'];
            }
        }
        $raw = array_values(array_filter(array_map(fn($m) => $per[$m][$id]['max'] ?? null, array_keys($per)), fn($v) => $v !== null));
        $perModel = array_map(fn($v) => array_sum($v) / count($v), $allMax[$id] ?? []);
        $out[] = ['id' => $id, 'name' => $name, 'lat' => $lat, 'lon' => $lon, 'models' => $models, 'tmax' => $tmax,
                  'all' => $perModel ? (int) round(array_sum($perModel) / count($perModel)) : null,
                  'mean' => $raw ? (int) round(array_sum($raw) / count($raw)) : null,
                  'lo' => $raw ? (int) round(min($raw)) : null, 'hi' => $raw ? (int) round(max($raw)) : null];
    }
    return $out;
}
 
/** Positions horaires : les modèles à pas de 3 ou 6 h sont interpolés entre deux points suivis. */
function hourly(array $points): array
{
    $out = [];
    for ($k = 0; $k + 1 < count($points); $k++) {
        $a = $points[$k];
        $b = $points[$k + 1];
        $ta = utc($a['t']);
        $gap = (utc($b['t']) - $ta) / 3600;
        if ($gap > 6.01) {
            continue;
        }
        for ($h = 0; $h <= (int) round($gap); $h++) {
            $f = $gap ? $h / $gap : 0;
            $out[iso($ta + $h * 3600)] = [$a['lat'] + $f * ($b['lat'] - $a['lat']), $a['lon'] + $f * ($b['lon'] - $a['lon']), $a['p'] + $f * ($b['p'] - $a['p'])];
        }
    }
    if (count($points) === 1) {
        $out[$points[0]['t']] = [$points[0]['lat'], $points[0]['lon'], $points[0]['p']];
    }
    return $out;
}
 
/**
 * Trajectoire « observée » : position du centre au départ de chaque run (l'analyse du modèle, qui
 * assimile les observations), moyennée entre modèles à chaque heure d'analyse. Après la dernière
 * analyse, et sur 6 h au plus, la prévision à courte échéance des derniers runs prend le relais
 * (marquée « f ») pour placer le centre sur les images satellite les plus récentes.
 */
function analysis_track(array $sel, array $tracks): array
{
    $byTime = [];
    foreach ($sel['all'] as $m => $runs) {
        foreach ($runs as $r) {
            $t0 = iso(run_time($r));
            foreach (json_in(derived($m, $r))['track'] as $p) {
                if ($p['t'] === $t0) {
                    $byTime[$t0][$m] = [$p['lat'], $p['lon'], $p['p']];
                }
            }
        }
    }
    ksort($byTime);
    $out = [];
    foreach ($byTime as $t => $ms) {
        $v = array_values($ms);
        $out[] = ['t' => $t, 'lat' => round(array_sum(array_column($v, 0)) / count($v), 3),
                  'lon' => round(array_sum(array_column($v, 1)) / count($v), 3),
                  'p' => round(array_sum(array_column($v, 2)) / count($v), 1), 'n' => count($v), 'models' => $ms];
    }
    if ($out) {
        $last = utc(end($out)['t']);
        $hly = array_map(fn($tr) => hourly($tr['points']), $tracks);
        for ($t = $last + 3600; $t <= $last + 6 * 3600; $t += 3600) {
            $k = iso($t);
            $pts = array_values(array_filter(array_map(fn($h) => $h[$k] ?? null, $hly)));
            if (count($pts) >= 2) {
                $out[] = ['t' => $k, 'lat' => round(array_sum(array_column($pts, 0)) / count($pts), 3),
                          'lon' => round(array_sum(array_column($pts, 1)) / count($pts), 3),
                          'p' => round(array_sum(array_column($pts, 2)) / count($pts), 1), 'n' => count($pts), 'f' => true];
            }
        }
    }
    return $out;
}
 
function build(array $sel, array $frames): void
{
    $tracks = $history = $series = [];
    foreach ($sel['all'] as $m => $runs) {
        $history[$m] = [];
        foreach ($runs as $r) {
            $d = json_in(derived($m, $r));
            $tr = $d['track'];
            if (!$tr) {
                continue;
            }
            $pm = array_reduce($tr, fn($c, $p) => $c === null || $p['p'] < $c['p'] ? $p : $c);
            $at6 = null;
            foreach ($tr as $p) {
                if ($p['t'] === '2026-10-23T06:00Z') {
                    $at6 = $p;
                }
            }
            $history[$m][] = ['run' => $r, 'pmin' => $pm['p'], 'tmin' => $pm['t'], 'lat' => $pm['lat'], 'lon' => $pm['lon'], 'at6' => $at6,
                              'track' => array_values(array_filter($tr, fn($p) => substr($p['t'], 14, 2) === '00' && (int) substr($p['t'], 11, 2) % 3 === 0))];
            if ($r === ($sel['latest'][$m] ?? null)) {
                $tracks[$m] = ['run' => $r, 'points' => $tr];
                if (isset($d['stations'])) {
                    $series[$m] = ['run' => $r, 'times' => $d['gust_times'],
                                   'values' => array_map(fn($v) => $v['series'], $d['stations'])];
                }
            }
        }
    }
    $hly = array_map(fn($tr) => hourly($tr['points']), $tracks);
    $cons = [];
    for ($t = utc('2026-10-22T06:00Z'); $t <= ep1(); $t += 3600) {
        $k = iso($t);
        $pts = array_values(array_filter(array_map(fn($h) => $h[$k] ?? null, $hly)));
        if (count($pts) >= 3) {
            $ps = array_column($pts, 2);
            $cons[] = ['t' => $k, 'lat' => round(array_sum(array_column($pts, 0)) / count($pts), 3),
                       'lon' => round(array_sum(array_column($pts, 1)) / count($pts), 3),
                       'p' => round(array_sum($ps) / count($ps), 1), 'pmin' => round(min($ps), 1), 'pmax' => round(max($ps), 1), 'n' => count($pts)];
        }
    }
    $analysis = analysis_track($sel, $tracks);
    /* moyenne de tous les runs : pour chaque heure, moyenne des runs de chaque modèle, puis entre modèles */
    $perModel = [];
    foreach ($sel['all'] as $m => $runs) {
        foreach ($runs as $r) {
            foreach (hourly(json_in(derived($m, $r))['track']) as $k => $p) {
                $perModel[$k][$m][] = $p;
            }
        }
    }
    ksort($perModel);
    $consAll = [];
    foreach ($perModel as $k => $ms) {
        if ($k < '2026-10-22T06:00Z' || utc($k) > ep1() || count($ms) < 3) {
            continue;
        }
        $avg = array_map(fn($ps) => [array_sum(array_column($ps, 0)) / count($ps), array_sum(array_column($ps, 1)) / count($ps), array_sum(array_column($ps, 2)) / count($ps)], $ms);
        $v = array_values($avg);
        $consAll[] = ['t' => $k, 'lat' => round(array_sum(array_column($v, 0)) / count($v), 3), 'lon' => round(array_sum(array_column($v, 1)) / count($v), 3),
                      'p' => round(array_sum(array_column($v, 2)) / count($v), 1), 'n' => array_sum(array_map('count', $ms))];
    }
    $ens = json_in(data_dir() . '/derived/ensemble.json');
    $names = json_in(data_dir() . '/geo/departements.json');
    $depts = array_map(fn($d) => ['code' => $d['code'], 'nom' => $names[$d['code']] ?? $d['code']] + $d, $ens['depts']);
    usort($depts, fn($a, $b) => $b['max'] <=> $a['max'] ?: strcmp($a['code'], $b['code']));
    $gustFc = array_filter($sel['forecast'], fn($m) => MODELS[$m]['gust'], ARRAY_FILTER_USE_KEY);
    $data = [
        'updated' => gmdate('Y-m-d\TH:i:s\Z'),
        'episode' => [iso(ep0()), iso(ep1())],
        'models' => array_map(fn($v) => array_intersect_key($v, array_flip(['name', 'res', 'color', 'gust', 'gnote', 'org'])), MODELS),
        'selection' => ['latest' => $sel['latest'], 'forecast' => $sel['forecast']],
        'tracks' => $tracks, 'consensus' => $cons, 'history' => $history, 'series' => $series,
        'stations' => station_table($sel['forecast'], allruns_selection($sel['all'])), 'depts' => $depts,
        'ensemble_models' => array_map(fn($m, $r) => [$m, $r], array_keys($gustFc), array_values($gustFc)),
        'frames' => $frames, 'track_box' => [-21.0, 11.0, 45.0, 55.8], 'analysis' => $analysis, 'consensus_all' => $consAll,
        'allruns' => allruns_selection($sel['all']),
        'legend' => ['bounds' => GUST_BOUNDS, 'colors' => GUST_COLORS],
    ];
    json_out(site_dir() . '/data.json', $data);
    if (!is_file(site_dir() . '/basemap.json')) {
        copy(data_dir() . '/geo/basemap.json', site_dir() . '/basemap.json');
    }
}
 
/* ---------------------------------------------------------------- publication */
/** Copie le dossier de travail vers le dossier publié, fichier par fichier (temporaire caché puis renommage). */
function publish(): void
{
    $src = site_dir();
    $dst = pub_dir();
    $n = 0;
    $it = new RecursiveIteratorIterator(new RecursiveDirectoryIterator($src, FilesystemIterator::SKIP_DOTS));
    foreach ($it as $f) {
        if ($f->isLink() || !$f->isFile() || str_starts_with($f->getFilename(), '.')) {
            continue;
        }
        $rel = substr($f->getPathname(), strlen($src) + 1);
        $to = "$dst/$rel";
        if (is_file($to) && filesize($to) === $f->getSize() && filemtime($to) >= $f->getMTime()) {
            continue;
        }
        ensure_dir(dirname($to));
        $tmp = dirname($to) . '/.' . basename($to) . '.new';
        if (!copy($f->getPathname(), $tmp) || !chmod($tmp, 0644) || !rename($tmp, $to)) {
            throw new RuntimeException("publication impossible : $rel");
        }
        $n++;
    }
    logm('publié', (string) $n, 'fichiers');
}
 
/* ---------------------------------------------------------------- principal */
function main(array $argv): int
{
    $runsArg = [];
    $noFetch = $force = false;
    for ($k = 1; $k < count($argv); $k++) {
        match ($argv[$k]) {
            '--no-fetch' => $noFetch = true,
            '--force' => $force = true,
            '--runs' => null,
            default => $runsArg[] = $argv[$k],
        };
    }
    if (time() > utc('2026-10-24T12:00Z') && !$runsArg) {
        logm('épisode terminé, plus de mise à jour automatique');
        return 0;
    }
    ensure_dir(data_dir());
    $lock = fopen(data_dir() . '/.pipeline.lock', 'c');
    if (!$lock || !flock($lock, LOCK_EX | LOCK_NB)) {
        logm('déjà en cours');
        return 0;
    }
    $conf = write_conf();
    $todo = $runsArg ? array_map(fn($a) => explode(':', $a, 2), $runsArg) : array_map(fn($m) => [$m, null], array_keys(MODELS));
    foreach ($todo as [$model, $run]) {
        if (!isset(MODELS[$model]) || ($run !== null && !valid_run($run))) {
            logm('argument ignoré', $model);
            continue;
        }
        if (!$noFetch) {
            [$code, $out, $err] = run(array_merge([PHP_BINARY, __DIR__ . '/fetch.php', $model], $run ? [$run] : []), 1800);
            file_put_contents(data_dir() . '/fetch.log', $out . $err, FILE_APPEND);
            if (preg_match('/^RUN (\d{10})$/m', $out, $mm)) {
                $run = $mm[1];
            }
        }
        if ($run) {
            logm($model, $run, process($model, $run, $conf) ? 'ok' : 'non traité');
        }
    }
    $latest = $forecast = $all = [];
    foreach (array_keys(MODELS) as $m) {
        $runs = decoded_runs($m);
        $all[$m] = $runs;
        if (!$runs) {
            continue;
        }
        $latest[$m] = end($runs);
        $before = array_values(array_filter($runs, fn($r) => run_time($r) <= ep0()));
        if ($before) {
            $forecast[$m] = end($before);
        }
    }
    $sel = ['latest' => $latest, 'forecast' => $forecast, 'all' => $all];
    $sp = data_dir() . '/selection.json';
    if (!$force && is_file($sp) && json_in($sp) == $sel) {
        logm('aucun nouveau run');
        return 0;
    }
    $args = [BIN, 'ensemble', data_dir() . '/geo/departements.bin', data_dir() . '/derived/ensemble'];
    foreach ($forecast as $m => $r) {
        if (MODELS[$m]['gust']) {
            $args[] = "$m=" . data_dir() . "/derived/$m/{$r}_gmax.grd";
        }
    }
    [$code, , $err] = run($args, 600);
    if ($code !== 0) {
        logm('ÉCHEC moyenne', trim($err));
        return 1;
    }
    $allSel = allruns_selection($all);
    if ($allSel) {
        $a2 = [BIN, 'ensemble', data_dir() . '/geo/departements.bin', data_dir() . '/derived/ensemble_all'];
        foreach ($allSel as $m => $runs) {
            $a2[] = "$m=" . implode(',', array_map(fn($r) => data_dir() . "/derived/$m/{$r}_gmax.grd", $runs));
        }
        [$code, , $err] = run($a2, 900);
        if ($code !== 0) {
            logm('ÉCHEC moyenne de tous les runs', trim($err));
        }
    }
    $stations = station_table($forecast, $allSel);
    $frames = render_all($latest, $forecast, $stations, $allSel);
    build($sel, $frames);
    run([PHP_BINARY, __DIR__ . '/export_code.php'], 300);
    publish();
    json_out($sp, $sel);
    $keep = [];
    foreach (array_keys(MODELS) as $m) {
        $keep[$m] = array_values(array_unique(array_filter([$latest[$m] ?? null, $forecast[$m] ?? null])));
    }
    prune_grids($keep);
    return 0;
}
 
exit(main($argv));
[code 11]

Le suivi des stations

obs.php, 343 lignes

Toutes les 6 minutes, les paquets d'observations de Météo-France : rafale, vent moyen, pression. Les fichiers viennent d'un dossier inscriptible par le site web : ils sont lus par une petite fonction C qui refuse liens et FIFO, et toute valeur non numérique est écartée.

obs.php
Fichier brut
<?php
/**
 * Suivi observé : obs.php — rafales et pression des stations Météo-France (paquets 6 minutes).
 * Écrit obs.json, maps/obs-gust.webp et vigilance-austen.json directement dans le dossier publié.
 * Le dossier des observations et la vigilance sont inscriptibles par le site : lecture par
 * « austen catsafe » (ni lien, ni FIFO, taille plafonnée), valeurs non numériques écartées.
 */
declare(strict_types=1);
require __DIR__ . '/lib.php';
 
const T0 = '2026-10-22T12:00Z';                 /* début des séries (avant l'épisode) */
const T1 = '2026-10-24T06:00Z';
const OBS_BOX = [-6.0, 9.0, 45.5, 51.3];        /* stations retenues pour la carte observée */
const FILE_RE = '/^paquetobs-stations-infrahoraire-6m_(\d{4}-\d\d-\d\dT\d\d_\d\d)_00Z\.csv$/';
 
/** Nombre fini ou null (« nan », « inf » et le texte sont écartés). */
function num(?string $x, float $scale = 1.0): ?float
{
    if ($x === null || $x === '' || !is_numeric($x)) {
        return null;
    }
    $v = (float) $x * $scale;
    return is_finite($v) ? round($v, 1) : null;
}
 
/** Une ligne par station retenue : [lat, lon, rafale, vent moyen, direction, pmer, direction de la rafale]. */
function parse_csv(string $text, array $follow): array
{
    $rows = [];
    $lines = explode("\n", $text);
    $head = str_getcsv(array_shift($lines), ';', '"', '');
    $col = array_flip($head);
    foreach (['lat', 'lon', 'geo_id_insee', 'raf10', 'ff', 'dd', 'pmer', 'ddraf10'] as $k) {
        if (!isset($col[$k])) {
            return [];
        }
    }
    foreach ($lines as $line) {
        if ($line === '') {
            continue;
        }
        $r = str_getcsv($line, ';', '"', '');
        if (count($r) !== count($head)) {
            continue;
        }
        $sid = $r[$col['geo_id_insee']];
        if (!is_numeric($r[$col['lat']]) || !is_numeric($r[$col['lon']]) || !preg_match('/^\d{8}$/', $sid)) {
            continue;
        }
        $la = (float) $r[$col['lat']];
        $lo = (float) $r[$col['lon']];
        if (!is_finite($la) || !is_finite($lo)) {
            continue;
        }
        if (!($lo >= OBS_BOX[0] && $lo <= OBS_BOX[1] && $la >= OBS_BOX[2] && $la <= OBS_BOX[3]) && !isset($follow[$sid])) {
            continue;
        }
        $raf = num($r[$col['raf10']], 3.6);
        if ($raf !== null && $raf > 250) {         /* valeur aberrante (> 69 m/s) */
            $raf = null;
        }
        $pm = num($r[$col['pmer']], 0.01);
        if ($pm !== null && !($pm > 930 && $pm < 1060)) {
            $pm = null;
        }
        $rows[$sid] = [$la, $lo, $raf, num($r[$col['ff']], 3.6), num($r[$col['dd']]), $pm, num($r[$col['ddraf10']])];
    }
    return $rows;
}
 
/**
 * Records de rafale d'octobre (fxi, km/h) par station, tirés des climatologies de collecteur (21 Mo) :
 * extraits une fois dans un petit cache, refait seulement si la source change.
 */
function october_records(): array
{
    $src = '/chemin/vers/climatologie/monthly_records_10.json';
    $cache = data_dir() . '/records_oct.json';
    $st = @lstat($src);
    if ($st === false || ($st['mode'] & 0170000) !== 0100000) {
        return is_file($cache) ? json_in($cache)['rec'] : [];
    }
    if (is_file($cache)) {
        $c = json_in($cache);
        if (($c['mtime'] ?? 0) === $st['mtime']) {
            return $c['rec'];
        }
    }
    $raw = read_untrusted($src, 64_000_000);
    $d = $raw !== null ? json_decode($raw, true, 64) : null;
    if (!is_array($d)) {
        return [];
    }
    $rec = [];
    $walk = function ($o) use (&$walk, &$rec): void {
        if (!is_array($o)) {
            return;
        }
        if (isset($o['id'], $o['records']) && is_string($o['id']) && preg_match('/^\d{8}$/', $o['id'])) {
            $f = $o['records']['fxi'] ?? null;
            if (is_array($f) && is_numeric($f['value'] ?? null) && is_string($f['date'] ?? null)) {
                $rec['s' . $o['id']] = [round((float) $f['value'], 1), substr($f['date'], 0, 10), (int) ($o['history_year_count'] ?? 0)];
            }
            return;
        }
        foreach ($o as $v) {
            $walk($v);
        }
    };
    $walk($d);
    json_out($cache, ['mtime' => $st['mtime'], 'rec' => $rec]);
    return $rec;
}
 
function merge_vigilance(string $outdir): void
{
    $raw = read_untrusted(VIGILANCE, 5_000_000);
    $v = $raw !== null ? json_decode($raw, true) : null;
    if (!is_array($v)) {
        return;
    }
    $path = "$outdir/vigilance-austen.json";
    $acc = is_file($path) ? json_in($path) : ['update' => null, 'comment' => '', 'depts' => []];
    foreach ((is_array($v['periods'] ?? null) ? $v['periods'] : []) as $p) {
        if (!is_array($p)) {
            continue;
        }
        if (is_string($p['comment'] ?? null) && stripos($p['comment'], 'AUSTEN') !== false) {
            $acc['comment'] = $p['comment'];
        }
        foreach ((is_array($p['depts'] ?? null) ? $p['depts'] : []) as $code => $dv) {
            $w = $dv['pheno']['1'] ?? null;
            if (!is_array($w) || !preg_match('/^(\d{2}|2[AB])$/', (string) $code)) {
                continue;
            }
            $sub = is_array($dv['pheno']['9'] ?? null) ? $dv['pheno']['9'] : [];
            $cur = $acc['depts'][$code] ?? ['c' => 1, 'orange' => [], 'yellow' => [], 'vs' => 1];
            $cur['c'] = max($cur['c'], (int) ($w['max'] ?? 1));
            $cur['vs'] = max($cur['vs'], (int) ($sub['max'] ?? 1));
            foreach ((is_array($w['frise'] ?? null) ? $w['frise'] : []) as $f) {
                if (!is_array($f) || !is_string($f['b'] ?? null) || !is_string($f['e'] ?? null)) {
                    continue;
                }
                $key = (int) ($f['c'] ?? 0) >= 3 ? 'orange' : ((int) ($f['c'] ?? 0) === 2 ? 'yellow' : null);
                if ($key && !in_array([$f['b'], $f['e']], $cur[$key], true)) {
                    $cur[$key][] = [$f['b'], $f['e']];
                    sort($cur[$key]);
                }
            }
            $acc['depts'][$code] = $cur;
        }
    }
    if (is_string($v['update_time'] ?? null)) {
        $acc['update'] = $v['update_time'];
    }
    json_out($path, $acc);
}
 
function main(): int
{
    if (time() > utc('2026-10-24T12:00Z')) {
        return 0;   /* épisode terminé : obs.json reste figé */
    }
    $out = pub_dir();
    ensure_dir($out);
    $follow = array_flip(array_column(STATIONS, 0));
    $cachePath = data_dir() . '/obs_cache.json';
    $cache = is_file($cachePath) ? json_in($cachePath) : [];
    $t0 = utc(T0);
    $t1 = utc(T1);
    foreach (scandir(OBS_DIR) ?: [] as $fn) {
        if (!preg_match(FILE_RE, $fn, $m)) {
            continue;
        }
        $t = utc(str_replace('_', ':', $m[1]) . 'Z');
        if ($t < $t0 || $t > $t1) {
            continue;
        }
        $st = @lstat(OBS_DIR . "/$fn");
        if ($st === false || ($st['mode'] & 0170000) !== 0100000) {
            continue;   /* lien, FIFO ou autre : ignoré */
        }
        $key = gmdate('YmdHi', $t);
        if (isset($cache[$key]) && $cache[$key][0] === $st['mtime']) {
            continue;
        }
        $text = read_untrusted(OBS_DIR . "/$fn", 5_000_000);
        if ($text !== null) {
            $cache[$key] = [$st['mtime'], parse_csv($text, $follow)];
        }
    }
    ksort($cache);
    json_out($cachePath, $cache);
 
    $names = [];
    $csv = read_untrusted(STATIONS_CSV, 5_000_000) ?? '';
    foreach (explode("\n", $csv) as $line) {
        $r = explode(';', $line);
        if (count($r) >= 6 && preg_match('/^\d{8}$/', $r[0])) {
            $names[$r[0]] = [$r[2], $r[5]];
        }
    }
    $ep0 = ep0();
    $iso = fn(string $k) => substr($k, 0, 4) . '-' . substr($k, 4, 2) . '-' . substr($k, 6, 2) . 'T' . substr($k, 8, 2) . ':' . substr($k, 10, 2) . 'Z';
 
    $records = october_records();
    /* stations suivies : séries complètes */
    $followOut = [];
    foreach (STATIONS as [$sid, $short]) {
        $ser = [];
        foreach ($cache as $k => [, $rows]) {
            if (isset($rows[$sid])) {
                $r = $rows[$sid];
                $ser[] = [$iso((string) $k), $r[2], $r[3], $r[5]];
            }
        }
        $ep = array_values(array_filter($ser, fn($s) => $s[1] !== null && utc($s[0]) >= $ep0));
        $all = array_values(array_filter($ser, fn($s) => $s[1] !== null));
        $pr = array_values(array_filter($ser, fn($s) => $s[3] !== null));
        $maxOf = fn(array $a, int $i) => $a ? array_reduce($a, fn($c, $s) => $c === null || $s[$i] > $c[$i] ? $s : $c) : null;
        $minOf = fn(array $a, int $i) => $a ? array_reduce($a, fn($c, $s) => $c === null || $s[$i] < $c[$i] ? $s : $c) : null;
        $mx = $maxOf($ep, 1);
        $mxa = $maxOf($all, 1);
        $pmin = $minOf($pr, 3);
        $dmax = null;
        if ($mx) {
            $k = str_replace(['-', 'T', ':', 'Z'], '', $mx[0]);
            $dmax = $cache[$k][1][$sid][6] ?? null;
        }
        $ro = $records['s' . $sid] ?? null;
        $followOut[$sid] = ['name' => $short, 'series' => $ser, 'max' => $mx[1] ?? null, 'tmax' => $mx[0] ?? null,
            'rec_oct' => $ro ? ['v' => $ro[0], 'd' => $ro[1], 'y' => $ro[2]] : null,
            'max_all' => $mxa[1] ?? null, 'tmax_all' => $mxa[0] ?? null, 'dmax' => $dmax,
            'pmin' => $pmin[3] ?? null, 'tpmin' => $pmin[0] ?? null];
    }
 
    /* toutes les stations du cadre : maxima depuis le début de l'épisode */
    $allst = [];
    $last = null;
    foreach ($cache as $k => [, $rows]) {
        $t = utc($iso((string) $k));
        $last = $iso((string) $k);
        if ($t < $ep0) {
            continue;
        }
        foreach ($rows as $sid => $r) {
            $a = $allst[$sid] ?? ['lat' => $r[0], 'lon' => $r[1], 'max' => null, 'tmax' => null, 'pmin' => null, 'tpmin' => null];
            if ($r[2] !== null && ($a['max'] === null || $r[2] > $a['max'])) {
                $a['max'] = $r[2];
                $a['tmax'] = $iso((string) $k);
            }
            if ($r[5] !== null && ($a['pmin'] === null || $r[5] < $a['pmin'])) {
                $a['pmin'] = $r[5];
                $a['tpmin'] = $iso((string) $k);
            }
            $allst[$sid] = $a;
        }
    }
    $row = function (string $sid, array $a) use ($names, $records): array {
        $ro = $records['s' . $sid] ?? null;
        return ['id' => $sid, 'name' => $names[$sid][0] ?? $sid, 'alt' => $names[$sid][1] ?? '',
            'lat' => $a['lat'], 'lon' => $a['lon'], 'max' => $a['max'], 'tmax' => $a['tmax'], 'pmin' => $a['pmin'], 'tpmin' => $a['tpmin'],
            'rec_oct' => $ro ? ['v' => $ro[0], 'd' => $ro[1], 'y' => $ro[2]] : null];
    };
    $top = array_filter($allst, fn($a) => $a['max'] !== null);
    uasort($top, fn($a, $b) => $b['max'] <=> $a['max']);
    $low = array_filter($allst, fn($a) => $a['pmin'] !== null);
    uasort($low, fn($a, $b) => $a['pmin'] <=> $b['pmin']);
    $started = $last !== null && utc($last) > $ep0;
    /* les identifiants à 8 chiffres deviennent des clés entières en PHP : on les rend en texte */
    $list = function (array $a, int $n) use ($row): array {
        $out = [];
        foreach ($a as $sid => $v) {
            if (count($out) >= $n) {
                break;
            }
            $out[] = $row(sprintf('%08d', (int) $sid), $v);
        }
        return $out;
    };
    $res = ['updated' => gmdate('Y-m-d\TH:i:s\Z'), 'last' => $last, 'started' => $started, 'follow' => $followOut,
        'top' => $list($top, 25), 'lowest' => $list($low, 8),
        /* records d'octobre battus ou égalés, stations d'au moins 10 ans d'historique */
        'records' => $list(array_filter($top, function (array $a, $sid) use ($records): bool {
            $ro = $records['s' . sprintf('%08d', (int) $sid)] ?? null;
            return $ro && $ro[2] >= 10 && $a['max'] >= $ro[0];
        }, ARRAY_FILTER_USE_BOTH), 40),
        'map' => array_values(array_map(fn($a) => [round($a['lon'], 3), round($a['lat'], 3), (int) round($a['max'])], $top))];
 
    /* carte d'abord : la page recharge l'image dès qu'obs.json change */
    if ($started && $top) {
        render_obs_map($out, $top, $last);
    }
    json_out("$out/obs.json", $res);
    if (time() <= ep1()) {
        merge_vigilance($out);
    }
    run([PHP_BINARY, __DIR__ . '/verify.php'], 120);   /* classement des modèles, à chaque nouvelle observation */
    logm('obs', (string) $last, 'stations', (string) count($allst));
    return 0;
}
 
function render_obs_map(string $out, array $top, string $last): void
{
    $tmp = data_dir() . '/render';
    ensure_dir($tmp);
    $asc = $top;
    uasort($asc, fn($a, $b) => $a['max'] <=> $b['max']);   /* les plus fortes dessinées en dernier */
    atomic_write("$tmp/obs_dots.tsv", implode('', array_map(fn($a) => sprintf("%.5f\t%.5f\t%.1f\n", $a['lon'], $a['lat'], $a['max']), $asc)));
    $placed = [];
    $lab = '';
    foreach ($top as $a) {   /* étiquettes des plus fortes, sans chevauchement */
        if ($a['max'] < 70 || count($placed) > 45) {
            continue;
        }
        foreach ($placed as [$x, $y]) {
            if (abs($a['lon'] - $x) < 0.42 && abs($a['lat'] - $y) < 0.16) {
                continue 2;
            }
        }
        $lab .= sprintf("%.5f\t%.5f\t*%d\n", $a['lon'], $a['lat'], (int) round($a['max']));
        $placed[] = [$a['lon'], $a['lat']];
    }
    atomic_write("$tmp/obs_labels.tsv", $lab);
    $d = (new DateTimeImmutable($last))->setTimezone(new DateTimeZone('Europe/Paris'));
    $jours = ['dimanche', 'lundi', 'mardi', 'mercredi', 'jeudi', 'vendredi', 'samedi'];
    $when = $jours[(int) $d->format('w')] . ' ' . $d->format('j') . ' à ' . $d->format('G') . ' h ' . $d->format('i');
    $geo = data_dir() . '/geo';
    $spec = "out $out/maps/.obs-gust.webp\nwidth 1200\nquality 88\nbackground #f4f1ea\nbox " . implode(' ', [-6.0, 9.0, 45.6, 51.3]) .
        "\nidw $tmp/obs_dots.tsv 45\ngeo $geo 1 1\ndots $tmp/obs_dots.tsv\nlabels $tmp/obs_labels.tsv\n" .
        "title Rafales maxi observées depuis jeudi 20 h (km/h)\nright Stations Météo-France\n" .
        "sub Données jusqu'au $when (heure de Paris)\nsign Meteociel.fr\nend\n";
    atomic_write("$tmp/obs_spec.txt", $spec);
    ensure_dir("$out/maps");
    [$code, , $err] = run([BIN, 'render', "$tmp/obs_spec.txt"], 120);
    if ($code === 0) {
        rename("$out/maps/.obs-gust.webp", "$out/maps/obs-gust.webp");
    } else {
        logm('ÉCHEC carte observée', trim($err));
    }
}
 
exit(main());
[code 12]

Le classement des modèles

verify.php, 217 lignes

Chaque run est comparé à ce qui a été mesuré après son départ, aux 38 stations : rafale maxi de chaque station sur les heures que le run couvre, et rafale heure par heure (la plus forte des dix mesures de 6 minutes de l'heure, huit mesures au moins). Le centre de la dépression est comparé à la trajectoire analysée, position et pression, à partir de 3 h après le départ du run.

Les écarts de tous les runs d'un même modèle sont moyennés ; le vainqueur est le modèle dont l'écart moyen sur la rafale maxi est le plus faible. Le calcul est refait à chaque nouvelle observation.

verify.php
Fichier brut
<?php
/**
 * Classement des modèles : verify.php — chaque run est comparé à ce qui a été observé, puis les runs
 * d'un même modèle sont moyennés ; le vainqueur est le modèle dont l'écart moyen est le plus faible.
 *
 * Rafales (38 stations suivies) : pour chaque heure H de l'épisode déjà observée, rafale maxi prévue sur
 * l'heure (H-1, H] contre la plus forte rafale relevée sur la même heure (paquets 6 minutes, 8 sur 10
 * au moins). Deux écarts par run : sur la rafale maxi de chaque station (sur les heures que le run
 * couvre), et heure par heure (qui juge aussi le moment). Un run n'est jugé que sur ce qui suit son départ.
 * Centre de la dépression : position et pression prévues contre la trajectoire analysée (analyses des
 * modèles, à partir de 3 h après le départ du run).
 */
declare(strict_types=1);
require __DIR__ . '/lib.php';
 
function dist_km(float $la, float $lo, float $la0, float $lo0): float
{
    return 111.19 * hypot($la - $la0, ($lo - $lo0) * cos(deg2rad($la0)));
}
 
/** Rafale observée maxi par station et par heure (fin d'heure, epoch) : [sid][H] = km/h. */
function observed_hourly(array $cache): array
{
    $acc = [];
    foreach ($cache as $k => [, $rows]) {
        $k = (string) $k;
        $t = gmmktime((int) substr($k, 8, 2), (int) substr($k, 10, 2), 0, (int) substr($k, 4, 2), (int) substr($k, 6, 2), (int) substr($k, 0, 4));
        $h = (int) (ceil($t / 3600) * 3600);          /* 6 min jusqu'à H inclus → heure H */
        foreach ($rows as $sid => $r) {
            if ($r[2] === null) {
                continue;
            }
            $sid = sprintf('%08d', (int) $sid);
            $a = $acc[$sid][$h] ?? [0, -1.0];
            $acc[$sid][$h] = [$a[0] + 1, max($a[1], (float) $r[2])];
        }
    }
    $out = [];
    foreach ($acc as $sid => $hours) {
        foreach ($hours as $h => [$n, $mx]) {
            if ($n >= 8) {
                $out[$sid][$h] = $mx;
            }
        }
    }
    return $out;
}
 
function mean(array $v): ?float { return $v ? array_sum($v) / count($v) : null; }
 
/**
 * Écarts d'une prévision de rafales sur les heures (max(h0, départ + 1 h) .. h1] : $fc(id, heure) donne la
 * rafale prévue ou null. Retourne null si moins de 10 stations sont comparables.
 */
function gust_score(callable $fc, int $rt, array $obs, array $ids, int $h0, int $h1): ?array
{
    $errH = $errMax = [];
    foreach ($ids as $sid) {
        if (!isset($obs[$sid])) {
            continue;
        }
        $fmax = $omax = null;
        for ($h = max($h0, $rt + 3600); $h <= $h1; $h += 3600) {
            $o = $obs[$sid][$h] ?? null;
            if ($o !== null) {
                $omax = max($omax ?? $o, $o);
            }
            $f = $fc($sid, $h);
            if ($f === null) {
                continue;
            }
            $fmax = max($fmax ?? $f, $f);
            if ($o !== null) {
                $errH[] = $f - $o;
            }
        }
        if ($fmax !== null && $omax !== null) {
            $errMax[] = $fmax - $omax;
        }
    }
    if (count($errMax) < 10) {
        return null;
    }
    return ['n_st' => count($errMax), 'n_h' => count($errH),
            'mae_max' => round(mean(array_map('abs', $errMax)), 1), 'bias_max' => round(mean($errMax), 1),
            'mae_h' => $errH ? round(mean(array_map('abs', $errH)), 1) : null, 'bias_h' => $errH ? round(mean($errH), 1) : null];
}
 
/** Série horaire d'un run : [id][heure] = rafale prévue. */
function run_series(string $m, string $r): array
{
    $d = json_in(data_dir() . "/derived/$m/$r.json");
    $out = [];
    foreach ($d['stations'] ?? [] as $sid => $v) {
        foreach ($d['gust_times'] as $i => $t) {
            if (($v['series'][$i] ?? null) !== null) {
                $out[sprintf('%08d', (int) $sid)][utc($t)] = $v['series'][$i];
            }
        }
    }
    return $out;
}
 
/** Moyenne de plusieurs séries de runs : d'abord par modèle, puis entre modèles. */
function blend_series(array $byModel): array
{
    $perModel = [];
    foreach ($byModel as $m => $runs) {
        $acc = [];
        foreach ($runs as $r) {
            foreach (run_series($m, $r) as $sid => $hours) {
                foreach ($hours as $h => $v) {
                    $acc[$sid][$h][] = $v;
                }
            }
        }
        foreach ($acc as $sid => $hours) {
            foreach ($hours as $h => $vs) {
                $perModel[$sid][$h][] = array_sum($vs) / count($vs);
            }
        }
    }
    return array_map(fn($hours) => array_map(fn($vs) => array_sum($vs) / count($vs), $hours), $perModel);
}
 
function main(): int
{
    $pub = pub_dir();
    $cachePath = data_dir() . '/obs_cache.json';
    if (!is_file($cachePath) || !is_file("$pub/data.json")) {
        return 0;
    }
    $obs = observed_hourly(json_in($cachePath));
    $data = json_in("$pub/data.json");
    $analysis = array_values(array_filter($data['analysis'] ?? [], fn($a) => empty($a['f'])));
    $anaAt = array_column($analysis, null, 't');
    $lastObs = 0;
    foreach ($obs as $hours) {
        $lastObs = max($lastObs, max(array_keys($hours)));
    }
    $h0 = ep0() + 3600;
    $h1 = min(ep1(), $lastObs);
    $ids = array_column(STATIONS, 0);
    $models = [];
    foreach (array_keys(MODELS) as $m) {
        $runs = [];
        foreach (glob(data_dir() . "/derived/$m/*.json") ?: [] as $p) {
            $r = basename($p, '.json');
            if (!valid_run($r)) {
                continue;
            }
            $rt = run_time($r);
            $d = json_in($p);
            $row = ['run' => $r];
            /* ---- rafales */
            if (MODELS[$m]['gust'] && isset($d['stations']) && $h1 >= $h0) {
                $ser = run_series($m, $r);
                $sc = gust_score(fn($sid, $h) => $ser[$sid][$h] ?? null, $rt, $obs, $ids, $h0, $h1);
                if ($sc) {
                    $row += $sc;
                }
            }
            /* ---- centre de la dépression */
            $dk = $dp = [];
            foreach ($d['track'] as $pt) {
                $t = utc($pt['t']);
                $a = $anaAt[$pt['t']] ?? null;
                if ($a && $t >= $rt + 3 * 3600) {
                    $dk[] = dist_km($pt['lat'], $pt['lon'], $a['lat'], $a['lon']);
                    $dp[] = $pt['p'] - $a['p'];
                }
            }
            if ($dk) {
                $row += ['n_c' => count($dk), 'km' => round(mean($dk)), 'dp' => round(mean($dp), 1), 'adp' => round(mean(array_map('abs', $dp)), 1)];
            }
            if (count($row) > 1) {
                $runs[] = $row;
            }
        }
        usort($runs, fn($a, $b) => strcmp($a['run'], $b['run']));
        $g = array_values(array_filter($runs, fn($x) => isset($x['mae_max'])));
        $c = array_values(array_filter($runs, fn($x) => isset($x['km'])));
        $models[$m] = ['runs' => $runs, 'n_runs' => count($g),
            'mae_max' => $g ? round(mean(array_column($g, 'mae_max')), 1) : null,
            'bias_max' => $g ? round(mean(array_column($g, 'bias_max')), 1) : null,
            'mae_h' => $g ? round(mean(array_filter(array_column($g, 'mae_h'), fn($v) => $v !== null)) ?? 0, 1) : null,
            'n_track' => count($c), 'km' => $c ? round(mean(array_column($c, 'km'))) : null,
            'adp' => $c ? round(mean(array_column($c, 'adp')), 1) : null, 'dp' => $c ? round(mean(array_column($c, 'dp')), 1) : null];
    }
    /* deux concurrents de plus : la moyenne des modèles (derniers runs partis avant l'épisode, comme la carte)
       et la moyenne de tous les runs partis avant l'épisode, chaque modèle moyenné sur ses runs d'abord */
    if ($h1 >= $h0) {
        $fc = array_filter($data['selection']['forecast'] ?? [], fn($m) => MODELS[$m]['gust'] ?? false, ARRAY_FILTER_USE_KEY);
        foreach (['mean' => array_map(fn($r) => [$r], $fc), 'allruns' => $data['allruns'] ?? []] as $key => $byModel) {
            if (!$byModel) {
                continue;
            }
            $ser = blend_series($byModel);
            $sc = gust_score(fn($sid, $h) => $ser[$sid][$h] ?? null, ep0(), $obs, $ids, $h0, $h1);
            if ($sc) {
                $models[$key] = $sc + ['runs' => [], 'n_runs' => array_sum(array_map('count', $byModel)), 'n_track' => 0,
                                       'km' => null, 'adp' => null, 'dp' => null, 'pseudo' => true];
            }
        }
    }
    $ranked = array_keys(array_filter($models, fn($v) => $v['mae_max'] !== null));
    usort($ranked, fn($a, $b) => $models[$a]['mae_max'] <=> $models[$b]['mae_max']);
    $rankedC = array_keys(array_filter($models, fn($v) => $v['km'] !== null));
    usort($rankedC, fn($a, $b) => ($models[$a]['km'] + 20 * $models[$a]['adp']) <=> ($models[$b]['km'] + 20 * $models[$b]['adp']));
    json_out("$pub/verif.json", ['updated' => gmdate('Y-m-d\TH:i:s\Z'),
        'window' => $h1 >= $h0 ? [iso($h0 - 3600), iso($h1)] : null, 'hours' => $h1 >= $h0 ? intdiv($h1 - $h0, 3600) + 1 : 0,
        'analyses' => count($analysis), 'models' => $models, 'ranking' => $ranked, 'ranking_center' => $rankedC,
        'winner' => $ranked[0] ?? null]);
    return 0;
}
 
exit(main());
[code 13]

Les images satellite

sat.php, 160 lignes

L'archive satellite du site ne garde que 24 heures : les images de la tempête sont copiées pour le bilan, avec la même lecture sûre et un contrôle de l'en-tête WebP. Une image infrarouge sur trois est ensuite composée par le binaire C : côtes et départements, et le « D » du centre analysé, interpolé entre deux analyses. Une image est refaite quand un nouveau run déplace ce centre.

sat.php
Fichier brut
<?php
/**
 * Images satellite de la tempête : sat.php
 * 1. Sauvegarde : l'archive du site ne garde que 24 heures ; IR 10,8, sandwich, visible couleur et FCI de
 *    la fenêtre sont copiés. Le dossier source est inscriptible par le site : lecture par « austen catsafe »
 *    (ni lien, ni FIFO, 30 Mo maxi), noms vérifiés par expression régulière, en-tête WebP contrôlé.
 * 2. Vue France : chaque image infrarouge (toutes les 5 minutes), composée par le binaire C avec les côtes,
 *    les départements et le « D » de la position analysée du centre (data.json, clé « analysis »).
 *    Une image est refaite quand la position du centre à son heure change (nouveau run).
 * 3. Vue Europe : l'infrarouge Europe de Meteociel (Meteosat 11 en balayage rapide, toutes les 5 minutes,
 *    côtes déjà tracées), téléchargé en HTTPS, contrôlé (signature PNG, taille) et réenregistré en JPEG.
 */
declare(strict_types=1);
require __DIR__ . '/lib.php';
 
const SAT_SRC = '/chemin/vers/archive-satellite';
const SAT_RE = '/^sat_(ir108|sandwich|visc|fcisand|fciir)_(2026102[234]T\d{4})(_m)?\.webp$/';
const SAT_BOX = [-8.0, 13.0, 39.5, 52.5];      /* emprise des images (projection lon/lat simple) */
const JOURS = ['dimanche', 'lundi', 'mardi', 'mercredi', 'jeudi', 'vendredi', 'samedi'];
 
function archive(string $dst): int
{
    $n = 0;
    foreach (scandir(SAT_SRC) ?: [] as $f) {
        if (!preg_match(SAT_RE, $f, $m) || $m[2] < '20261022T1200' || $m[2] > '20261024T0300' || is_file("$dst/$f")) {
            continue;
        }
        $data = read_untrusted(SAT_SRC . "/$f", 30_000_000);
        if ($data !== null && str_starts_with($data, 'RIFF') && substr($data, 8, 4) === 'WEBP') {
            atomic_write("$dst/$f", $data);
            $n++;
        }
    }
    return $n;
}
 
/** Position du centre à l'instant $t, interpolée entre deux points distants de 6 h au plus. */
function center_at(array $analysis, int $t): ?array
{
    for ($k = 0; $k + 1 < count($analysis); $k++) {
        $a = $analysis[$k];
        $b = $analysis[$k + 1];
        $ta = utc($a['t']);
        $tb = utc($b['t']);
        if ($t >= $ta && $t <= $tb && $tb - $ta <= 6 * 3600) {
            $f = $tb > $ta ? ($t - $ta) / ($tb - $ta) : 0;
            return [round($a['lat'] + $f * ($b['lat'] - $a['lat']), 2), round($a['lon'] + $f * ($b['lon'] - $a['lon']), 2),
                    round($a['p'] + $f * ($b['p'] - $a['p']))];
        }
    }
    return null;
}
 
function animate(string $src, string $pub): int
{
    $dir = "$pub/sat/ir";
    ensure_dir($dir);
    $indexPath = "$pub/sat/index.json";
    $index = is_file($indexPath) ? json_in($indexPath) : ['frames' => []];
    $known = array_column($index['frames'], null, 'k');
    $data = is_file("$pub/data.json") ? json_in("$pub/data.json") : [];
    $analysis = $data['analysis'] ?? [];
    $geo = data_dir() . '/geo';
    $spec = '';
    $frames = [];
    foreach (glob("$src/sat_ir108_*_m.webp") ?: [] as $f) {
        if (!preg_match('/_(\d{8})T(\d{2})(\d{2})_m\.webp$/', $f, $m)) {
            continue;
        }
        $t = gmmktime((int) $m[2], (int) $m[3], 0, (int) substr($m[1], 4, 2), (int) substr($m[1], 6, 2), (int) substr($m[1], 0, 4));
        $k = gmdate('YmdHi', $t);
        $c = center_at($analysis, $t);
        $frames[] = ['k' => $k, 'd' => $c];
        if (isset($known[$k]) && $known[$k]['d'] === $c && is_file("$dir/$k.webp")) {
            continue;
        }
        $d = (new DateTimeImmutable('@' . $t))->setTimezone(new DateTimeZone('Europe/Paris'));
        $when = JOURS[(int) $d->format('w')] . ' ' . $d->format('j') . ' à ' . $d->format('G') . ' h ' . $d->format('i');
        $b = ['out' => "$dir/$k.webp", 'quality' => 80, 'image' => $f, 'box' => implode(' ', [SAT_BOX[0], SAT_BOX[1], SAT_BOX[2], SAT_BOX[3]]),
              'geo' => "$geo 1 0.8", 'geocolor' => '#f2c94c', 'title' => 'Austen vue par Meteosat (infrarouge 10,8 µm)',
              'sub' => mb_strtoupper(mb_substr($when, 0, 1)) . mb_substr($when, 1) . ' (heure de Paris)', 'right' => 'MSG, EUMETSAT', 'sign' => 'Meteociel.fr'];
        if ($c) {
            $b['center'] = sprintf('%.2f %.2f %d', $c[1], $c[0], $c[2]);
        }
        foreach ($b as $kk => $v) {
            $spec .= $kk . ' ' . str_replace(["\n", "\r"], ' ', (string) $v) . "\n";
        }
        $spec .= "end\n";
    }
    if ($spec !== '') {
        $tmp = data_dir() . '/render/sat_spec.txt';
        atomic_write($tmp, $spec);
        [$code, , $err] = run([BIN, 'render', $tmp], 600);
        if ($code !== 0) {
            logm('ÉCHEC images satellite', trim($err));
            return 0;
        }
    }
    usort($frames, fn($a, $b) => strcmp($a['k'], $b['k']));
    $index['updated'] = gmdate('Y-m-d\TH:i:s\Z');
    $index['box'] = SAT_BOX;
    $index['frames'] = $frames;
    json_out($indexPath, $index);
    return substr_count($spec, "end\n");
}
 
/** Vue Europe : images de Meteociel, une toutes les 5 minutes, réenregistrées en JPEG (environ 100 Ko). */
function europe(string $pub): int
{
    $dir = "$pub/sat/eu";
    ensure_dir($dir);
    $indexPath = "$pub/sat/index.json";
    $n = 0;
    $end = min(time() - 600, utc('2026-10-24T03:00Z'));
    for ($t = utc('2026-10-22T12:00Z'); $t <= $end && $n < 80; $t += 300) {
        $k = gmdate('YmdHi', $t);
        if (is_file("$dir/$k.jpg")) {
            continue;
        }
        $url = sprintf('https://modeles20.meteociel.fr/satellite/archives/msgrss_%s/satirmsgeu-%s-%s.png', gmdate('Ymd', $t), gmdate('H', $t), gmdate('i', $t));
        [$code, $png] = http_get($url, [], null, 3 << 20);
        if ($code !== 200 || !str_starts_with($png, "\x89PNG\r\n\x1a\n")) {
            continue;   /* image absente chez Meteociel : on réessaiera au passage suivant */
        }
        $im = @imagecreatefromstring($png);
        if ($im === false || imagesx($im) > 2000 || imagesy($im) > 2000) {
            continue;
        }
        $tmp = "$dir/.$k.jpg.new";
        imagejpeg($im, $tmp, 75);
        chmod($tmp, 0644);
        rename($tmp, "$dir/$k.jpg");
        $n++;
        usleep(300_000);   /* téléchargements espacés */
    }
    $list = array_map(fn($f) => basename($f, '.jpg'), glob("$dir/*.jpg") ?: []);
    sort($list);
    $index = is_file($indexPath) ? json_in($indexPath) : ['frames' => []];
    $index['eu'] = $list;
    json_out($indexPath, $index);
    return $n;
}
 
function main(): int
{
    if (time() > utc('2026-10-24T12:00Z')) {
        return 0;
    }
    $dst = data_dir() . '/sat_archive';
    ensure_dir($dst);
    $n = archive($dst);
    $r = animate($dst, pub_dir());
    $e = europe(pub_dir());
    if ($n || $r || $e) {
        logm('satellite :', (string) $n, 'images copiées,', (string) $r, 'images composées,', (string) $e, 'images Europe');
    }
    return 0;
}
 
exit(main());
[code 14]

Outils C communs

util.c, 290 lignes

Allocations contrôlées (débordement de multiplication vérifié), lecture de fichier sans suivre de lien, écriture atomique, format de grille, dates UTC sans dépendre du fuseau du système, et la lecture sûre utilisée par la partie PHP.

util.c
Fichier brut
#define _GNU_SOURCE
#include "util.h"
 
#include <errno.h>
#include <fcntl.h>
#include <stdarg.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <sys/stat.h>
#include <time.h>
#include <unistd.h>
 
void die(const char *fmt, ...)
{
    va_list ap;
    va_start(ap, fmt);
    fputs("austen: ", stderr);
    vfprintf(stderr, fmt, ap);
    fputc('\n', stderr);
    va_end(ap);
    exit(1);
}
 
void *xmalloc(size_t n, size_t size)
{
    if (size && n > SIZE_MAX / size)
        die("allocation trop grande");
    size_t total = n * size;
    void *p = malloc(total ? total : 1);
    if (!p)
        die("mémoire épuisée");
    return p;
}
 
void *xcalloc(size_t n, size_t size)
{
    void *p = calloc(n ? n : 1, size ? size : 1);
    if (!p)
        die("mémoire épuisée");
    return p;
}
 
void *xrealloc(void *p, size_t n, size_t size)
{
    if (size && n > SIZE_MAX / size)
        die("allocation trop grande");
    size_t total = n * size;
    p = realloc(p, total ? total : 1);
    if (!p)
        die("mémoire épuisée");
    return p;
}
 
char *xstrdup(const char *s)
{
    size_t n = strlen(s) + 1;
    char *d = xmalloc(n, 1);
    memcpy(d, s, n);
    return d;
}
 
uint8_t *read_file(const char *path, size_t max, size_t *len)
{
    int fd = open(path, O_RDONLY | O_NOFOLLOW | O_CLOEXEC);
    if (fd < 0)
        die("ouverture impossible : %s (%s)", path, strerror(errno));
    struct stat st;
    if (fstat(fd, &st) != 0 || !S_ISREG(st.st_mode))
        die("pas un fichier ordinaire : %s", path);
    if ((uintmax_t)st.st_size > max)
        die("fichier trop gros : %s", path);
    size_t n = (size_t)st.st_size;
    uint8_t *buf = xmalloc(n + 1, 1);
    size_t got = 0;
    while (got < n) {
        ssize_t r = read(fd, buf + got, n - got);
        if (r < 0 && errno == EINTR)
            continue;
        if (r <= 0)
            die("lecture incomplète : %s", path);
        got += (size_t)r;
    }
    buf[n] = 0;
    close(fd);
    *len = n;
    return buf;
}
 
void write_file_atomic(const char *path, const void *data, size_t len)
{
    size_t pl = strlen(path);
    char *tmp = xmalloc(pl + 5, 1);
    memcpy(tmp, path, pl);
    memcpy(tmp + pl, ".new", 5);
    int fd = open(tmp, O_WRONLY | O_CREAT | O_TRUNC | O_NOFOLLOW | O_CLOEXEC, 0644);
    if (fd < 0)
        die("écriture impossible : %s (%s)", tmp, strerror(errno));
    const uint8_t *p = data;
    size_t done = 0;
    while (done < len) {
        ssize_t w = write(fd, p + done, len - done);
        if (w < 0 && errno == EINTR)
            continue;
        if (w <= 0)
            die("écriture incomplète : %s", tmp);
        done += (size_t)w;
    }
    if (fchmod(fd, 0644) != 0 || close(fd) != 0)
        die("fermeture impossible : %s", tmp);
    if (rename(tmp, path) != 0)
        die("renommage impossible : %s", path);
    free(tmp);
}
 
size_t grid_points(const Grid *g)
{
    if (g->nx <= 0 || g->ny <= 0 || (size_t)g->nx * (size_t)g->ny > MAX_POINTS)
        die("grille invalide %dx%d", g->nx, g->ny);
    return (size_t)g->nx * (size_t)g->ny;
}
 
/* ------------------------------------------------------------------ .grd */
typedef struct {
    char magic[4];
    uint32_t version, nx, ny, nt_msl, nt_gust;
    double lat0, dlat, lon0, dlon;
} GrdHeader;
 
Grd *grd_read(const char *path)
{
    size_t len;
    uint8_t *buf = read_file(path, (size_t)8 << 30, &len);
    if (len < sizeof(GrdHeader))
        die("%s : en-tête tronqué", path);
    GrdHeader h;
    memcpy(&h, buf, sizeof h);
    if (memcmp(h.magic, "AUGR", 4) != 0 || h.version != 1)
        die("%s : format inconnu", path);
    if (h.nt_msl > MAX_TIMES || h.nt_gust > MAX_TIMES)
        die("%s : trop d'échéances", path);
    Grd *d = xcalloc(1, sizeof *d);
    d->g = (Grid){(int)h.nx, (int)h.ny, h.lat0, h.dlat, h.lon0, h.dlon};
    size_t np = grid_points(&d->g);
    d->nt_msl = (int)h.nt_msl;
    d->nt_gust = (int)h.nt_gust;
    size_t need = sizeof h + (h.nt_msl + h.nt_gust) * sizeof(int64_t) + (h.nt_msl + h.nt_gust) * np * sizeof(float);
    if (len != need)
        die("%s : taille incohérente", path);
    const uint8_t *p = buf + sizeof h;
    d->t_msl = xmalloc(h.nt_msl, sizeof(int64_t));
    d->t_gust = xmalloc(h.nt_gust, sizeof(int64_t));
    memcpy(d->t_msl, p, h.nt_msl * sizeof(int64_t));
    p += h.nt_msl * sizeof(int64_t);
    memcpy(d->t_gust, p, h.nt_gust * sizeof(int64_t));
    p += h.nt_gust * sizeof(int64_t);
    d->msl = xmalloc((size_t)h.nt_msl * np, sizeof(float));
    d->gust = xmalloc((size_t)h.nt_gust * np, sizeof(float));
    memcpy(d->msl, p, (size_t)h.nt_msl * np * sizeof(float));
    p += (size_t)h.nt_msl * np * sizeof(float);
    memcpy(d->gust, p, (size_t)h.nt_gust * np * sizeof(float));
    free(buf);
    return d;
}
 
void grd_write(const char *path, const Grd *d)
{
    size_t np = grid_points(&d->g);
    GrdHeader h = {{'A', 'U', 'G', 'R'}, 1, (uint32_t)d->g.nx, (uint32_t)d->g.ny, (uint32_t)d->nt_msl,
                   (uint32_t)d->nt_gust, d->g.lat0, d->g.dlat, d->g.lon0, d->g.dlon};
    size_t len = sizeof h + (size_t)(d->nt_msl + d->nt_gust) * (sizeof(int64_t) + np * sizeof(float));
    uint8_t *buf = xmalloc(len, 1), *p = buf;
    memcpy(p, &h, sizeof h);
    p += sizeof h;
    memcpy(p, d->t_msl, (size_t)d->nt_msl * sizeof(int64_t));
    p += (size_t)d->nt_msl * sizeof(int64_t);
    memcpy(p, d->t_gust, (size_t)d->nt_gust * sizeof(int64_t));
    p += (size_t)d->nt_gust * sizeof(int64_t);
    memcpy(p, d->msl, (size_t)d->nt_msl * np * sizeof(float));
    p += (size_t)d->nt_msl * np * sizeof(float);
    memcpy(p, d->gust, (size_t)d->nt_gust * np * sizeof(float));
    write_file_atomic(path, buf, len);
    free(buf);
}
 
void grd_free(Grd *d)
{
    if (!d)
        return;
    free(d->t_msl);
    free(d->t_gust);
    free(d->msl);
    free(d->gust);
    free(d);
}
 
/* ------------------------------------------------------------------ dates */
int64_t epoch_utc(int y, int mo, int d, int h, int mi, int s)
{
    /* jours depuis 1970 (algorithme de Howard Hinnant), sans dépendre du fuseau du système */
    y -= mo <= 2;
    int64_t era = (y >= 0 ? y : y - 399) / 400;
    int64_t yoe = y - era * 400;
    int64_t doy = (153 * (mo + (mo > 2 ? -3 : 9)) + 2) / 5 + d - 1;
    int64_t doe = yoe * 365 + yoe / 4 - yoe / 100 + doy;
    int64_t days = era * 146097 + doe - 719468;
    return days * 86400 + h * 3600 + mi * 60 + s;
}
 
int64_t parse_iso(const char *s)
{
    int y, mo, d, h, mi;
    char z;
    if (sscanf(s, "%4d-%2d-%2dT%2d:%2d%c", &y, &mo, &d, &h, &mi, &z) != 6 || z != 'Z')
        return -1;
    if (mo < 1 || mo > 12 || d < 1 || d > 31 || h > 23 || mi > 59)
        return -1;
    return epoch_utc(y, mo, d, h, mi, 0);
}
 
void fmt_iso(int64_t t, char out[20])
{
    time_t tt = (time_t)t;
    struct tm tm;
    gmtime_r(&tt, &tm);
    snprintf(out, 20, "%04d-%02d-%02dT%02d:%02dZ", tm.tm_year + 1900, tm.tm_mon + 1, tm.tm_mday, tm.tm_hour, tm.tm_min);
}
 
char *next_line(char **cur)
{
    char *s = *cur;
    while (*s == '\n')
        s++;
    if (!*s) {
        *cur = s;
        return NULL;
    }
    char *e = strchr(s, '\n');
    if (e) {
        *e = 0;
        *cur = e + 1;
    } else {
        *cur = s + strlen(s);
    }
    return s;
}
 
double now_epoch(void)
{
    struct timespec ts;
    clock_gettime(CLOCK_REALTIME, &ts);
    return (double)ts.tv_sec + (double)ts.tv_nsec * 1e-9;
}
 
int cmd_catsafe(int argc, char **argv)
{
    /* « austen catsafe <fichier> <octets maxi> » : copie le fichier sur la sortie standard sans suivre
     * de lien (O_NOFOLLOW), sans se bloquer sur une FIFO (O_NONBLOCK) et seulement s'il est ordinaire.
     * Sert à la partie PHP, qui ne sait pas ouvrir un fichier avec ces drapeaux. */
    if (argc != 2)
        die("usage : austen catsafe <fichier> <octets maxi>");
    char *end;
    unsigned long long max = strtoull(argv[1], &end, 10);
    if (*end || max == 0 || max > (1ull << 30))
        die("taille maxi invalide");
    int fd = open(argv[0], O_RDONLY | O_NOFOLLOW | O_NONBLOCK | O_CLOEXEC);
    if (fd < 0)
        die("ouverture refusée : %s", argv[0]);
    struct stat st;
    if (fstat(fd, &st) != 0 || !S_ISREG(st.st_mode) || (unsigned long long)st.st_size > max)
        die("fichier refusé : %s", argv[0]);
    char buf[65536];
    unsigned long long total = 0;
    for (;;) {
        ssize_t r = read(fd, buf, sizeof buf);
        if (r < 0 && errno == EINTR)
            continue;
        if (r < 0)
            die("lecture impossible : %s", argv[0]);
        if (r == 0)
            break;
        total += (unsigned long long)r;
        if (total > max)
            die("fichier trop gros : %s", argv[0]);
        if (fwrite(buf, 1, (size_t)r, stdout) != (size_t)r)
            die("écriture impossible");
    }
    close(fd);
    return fflush(stdout) ? 1 : 0;
}
[code 15]

L'écriture JSON

json.c, 114 lignes

Le C écrit du JSON mais n'en lit jamais : ses réglages lui arrivent par un fichier texte strictement vérifié.

json.c
Fichier brut
#include "json.h"
 
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
 
#include "util.h"
 
static void put(Json *j, const char *s, size_t n)
{
    if (j->len + n + 1 > j->cap) {
        size_t c = j->cap ? j->cap : 4096;
        while (c < j->len + n + 1)
            c *= 2;
        j->buf = xrealloc(j->buf, c, 1);
        j->cap = c;
    }
    memcpy(j->buf + j->len, s, n);
    j->len += n;
    j->buf[j->len] = 0;
}
 
static void puts_(Json *j, const char *s) { put(j, s, strlen(s)); }
 
/* virgule avant un élément si ce n'est pas le premier du conteneur courant */
static void sep(Json *j)
{
    if (j->depth > 0 && j->need_comma[j->depth])
        put(j, ",", 1);
    if (j->depth > 0)
        j->need_comma[j->depth] = 1;
}
 
void js_init(Json *j) { memset(j, 0, sizeof *j); }
void js_free(Json *j) { free(j->buf); memset(j, 0, sizeof *j); }
 
static void open_(Json *j, const char *c)
{
    sep(j);
    puts_(j, c);
    if (j->depth >= 63)
        die("JSON trop profond");
    j->depth++;
    j->need_comma[j->depth] = 0;
}
 
void js_obj(Json *j) { open_(j, "{"); }
void js_arr(Json *j) { open_(j, "["); }
void js_end_obj(Json *j) { j->depth--; puts_(j, "}"); }
void js_end_arr(Json *j) { j->depth--; puts_(j, "]"); }
 
static void quoted(Json *j, const char *s)
{
    put(j, "\"", 1);
    for (const unsigned char *p = (const unsigned char *)s; *p; p++) {
        char e[8];
        if (*p == '"' || *p == '\\') {
            e[0] = '\\';
            e[1] = (char)*p;
            put(j, e, 2);
        } else if (*p < 0x20) {
            snprintf(e, sizeof e, "\\u%04x", *p);
            puts_(j, e);
        } else {
            put(j, (const char *)p, 1);
        }
    }
    put(j, "\"", 1);
}
 
void js_key(Json *j, const char *k)
{
    sep(j);
    quoted(j, k);
    put(j, ":", 1);
    j->need_comma[j->depth] = 0;   /* la valeur qui suit ne prend pas de virgule */
}
 
void js_str(Json *j, const char *s) { sep(j); quoted(j, s); }
 
void js_num(Json *j, double v, int decimals)
{
    sep(j);
    if (!isfinite(v)) {
        puts_(j, "null");
        return;
    }
    char b[64];
    snprintf(b, sizeof b, "%.*f", decimals, v);
    /* pas de zéros inutiles : 12.50 → 12.5, 3.0 → 3 */
    if (decimals > 0) {
        char *e = b + strlen(b) - 1;
        while (*e == '0')
            *e-- = 0;
        if (*e == '.')
            *e = 0;
    }
    if (strcmp(b, "-0") == 0)
        strcpy(b, "0");
    puts_(j, b);
}
 
void js_int(Json *j, long long v)
{
    sep(j);
    char b[32];
    snprintf(b, sizeof b, "%lld", v);
    puts_(j, b);
}
 
void js_null(Json *j) { sep(j); puts_(j, "null"); }
 
void js_save(Json *j, const char *path) { write_file_atomic(path, j->buf ? j->buf : "", j->len); }
[code 16]

Le point d'entrée

main.c, 30 lignes

Quatre sous-commandes : decode, analyze, ensemble, render, plus la lecture sûre.

main.c
Fichier brut
/* austen : chaîne C de la page tempête Austen (décodage GRIB2, analyse, moyenne, cartes).
 * Sous-commandes : decode, analyze, ensemble, render. Les réglages viennent de la partie PHP. */
#include <stdio.h>
#include <string.h>
 
#include "util.h"
 
int cmd_decode(int argc, char **argv);
int cmd_analyze(int argc, char **argv);
int cmd_ensemble(int argc, char **argv);
int cmd_render(int argc, char **argv);
int cmd_catsafe(int argc, char **argv);
 
int main(int argc, char **argv)
{
    if (argc < 2)
        die("usage : austen decode|analyze|ensemble|render|catsafe ...");
    const char *c = argv[1];
    if (!strcmp(c, "decode"))
        return cmd_decode(argc - 2, argv + 2);
    if (!strcmp(c, "analyze"))
        return cmd_analyze(argc - 2, argv + 2);
    if (!strcmp(c, "ensemble"))
        return cmd_ensemble(argc - 2, argv + 2);
    if (!strcmp(c, "render"))
        return cmd_render(argc - 2, argv + 2);
    if (!strcmp(c, "catsafe"))
        return cmd_catsafe(argc - 2, argv + 2);
    die("sous-commande inconnue : %s", c);
}
[code 17]

La compilation

Makefile, 29 lignes

C11, tous les avertissements traités comme des erreurs, protections de pile et de tas, binaire à adresse aléatoire et réadressage en lecture seule. Trois bibliothèques seulement : libaec (CCSDS), FreeType (texte), libwebp (images).

Makefile
Fichier brut
# austen : chaîne C de la page tempête Austen.
# Dépendances : libaec (décompression CCSDS des GRIB2), freetype (texte des cartes), libwebp (images).
CC      = gcc
WARN    = -Wall -Wextra -Wformat=2 -Wshadow -Wconversion -Wno-sign-conversion -Wvla -Werror
HARDEN  = -fstack-protector-strong -D_FORTIFY_SOURCE=2 -fPIE -fno-strict-overflow -fstack-clash-protection
CFLAGS  = -std=c11 -O2 -g $(WARN) $(HARDEN) $(shell pkg-config --cflags freetype2)
LDFLAGS = -pie -Wl,-z,relro,-z,now,-z,noexecstack
LDLIBS  = -laec $(shell pkg-config --libs freetype2) -lwebp -lm
 
SRC  = util.c grib2.c json.c field.c decode.c analyze.c ensemble.c render.c main.c
OBJ  = $(SRC:.c=.o)
 
austen: $(OBJ)
	$(CC) $(LDFLAGS) -o $@ $^ $(LDLIBS)
 
%.o: %.c *.h
	$(CC) $(CFLAGS) -c -o $@ $<
 
# Banc d'essai du décodeur (comparaison à eccodes) et version instrumentée ASan + UBSan.
test/check_grib: test/check_grib.c grib2.c util.c grib2.h util.h
	$(CC) -std=c11 -O1 -g $(WARN) -fsanitize=address,undefined -fno-sanitize-recover=all -I. -o $@ test/check_grib.c grib2.c util.c -laec -lm
 
test/fuzz_grib: test/fuzz_grib.c grib2.c util.c grib2.h util.h
	$(CC) -std=c11 -O1 -g $(WARN) -fsanitize=address,undefined -fno-sanitize-recover=all -I. -o $@ test/fuzz_grib.c grib2.c util.c -laec -lm
 
clean:
	rm -f $(OBJ) austen test/check_grib test/fuzz_grib
 
.PHONY: clean
[code 18]

Les fonds de carte

geo_prep.php, 183 lignes

Préparation unique : lecture des fichiers Natural Earth (format shapefile décodé à la main) et des contours des départements, simplification de Douglas-Peucker pour la carte de trajectoire.

geo_prep.php
Fichier brut
<?php
/**
 * Préparation unique des fonds de carte : geo_prep.php <rép. Natural Earth> <departements.geojson> <sortie>
 *  - coast.bin, border.bin, dept.bin : polylignes pour le rendu C (u32 nb ; u32 nb points, f32 lon, lat) ;
 *  - departements.bin : anneaux extérieurs pour les masques de « austen ensemble » ;
 *  - departements.json : code → nom ;
 *  - basemap.json : côtes et frontières simplifiées (Douglas-Peucker) pour la carte de trajectoire.
 */
declare(strict_types=1);
require __DIR__ . '/lib.php';
 
const CLIP = [-30.0, 20.0, 38.0, 62.0];          /* cadre utile des cartes, plus un degré de marge */
const TRACK_BOX = [-21.0, 11.0, 45.0, 55.8];     /* carte de trajectoire de la page */
 
/** Lecture d'un .shp de polylignes ou polygones (format ESRI) : liste de lignes [[lon, lat], ...]. */
function read_shp(string $path): array
{
    $b = file_get_contents($path);
    if ($b === false || strlen($b) < 100 || unpack('N', $b)[1] !== 9994) {
        throw new RuntimeException("shapefile invalide : $path");
    }
    $lines = [];
    $p = 100;
    $n = strlen($b);
    while ($p + 8 <= $n) {
        $len = unpack('N', $b, $p + 4)[1] * 2;          /* longueur en mots de 16 bits, gros-boutiste */
        $c = $p + 8;
        if ($len < 4 || $c + $len > $n) {
            throw new RuntimeException("enregistrement tronqué : $path");
        }
        $type = unpack('V', $b, $c)[1];
        if ($type === 3 || $type === 5) {                /* polyligne ou polygone */
            ['np' => $np, 'nq' => $nq] = unpack('Vnp/Vnq', $b, $c + 36);
            if ($c + 44 + 4 * $np + 16 * $nq > $c + $len) {
                throw new RuntimeException("géométrie tronquée : $path");
            }
            $parts = array_values(unpack("V$np", $b, $c + 44));
            $parts[] = $nq;
            $pts = array_values(unpack('e' . (2 * $nq), $b, $c + 44 + 4 * $np));
            for ($k = 0; $k < $np; $k++) {
                $line = [];
                for ($q = $parts[$k]; $q < $parts[$k + 1]; $q++) {
                    $line[] = [$pts[2 * $q], $pts[2 * $q + 1]];
                }
                $lines[] = $line;
            }
        }
        $p = $c + $len;
    }
    return $lines;
}
 
function clip(array $lines): array
{
    return array_values(array_filter($lines, function (array $l): bool {
        foreach ($l as [$x, $y]) {
            if ($x > CLIP[0] - 1 && $x < CLIP[1] + 1 && $y > CLIP[2] - 1 && $y < CLIP[3] + 1) {
                return count($l) >= 2;
            }
        }
        return false;
    }));
}
 
function write_lines(string $path, array $lines): void
{
    $out = pack('V', count($lines));
    foreach ($lines as $l) {
        $out .= pack('V', count($l));
        foreach ($l as [$x, $y]) {
            $out .= pack('gg', $x, $y);
        }
    }
    atomic_write($path, $out);
}
 
/** Douglas-Peucker itératif. */
function simplify(array $xy, float $tol): array
{
    $n = count($xy);
    if ($n < 3) {
        return $xy;
    }
    $keep = array_fill(0, $n, false);
    $keep[0] = $keep[$n - 1] = true;
    $stack = [[0, $n - 1]];
    while ($stack) {
        [$a, $b] = array_pop($stack);
        if ($b <= $a + 1) {
            continue;
        }
        [$px, $py] = $xy[$a];
        $dx = $xy[$b][0] - $px;
        $dy = $xy[$b][1] - $py;
        $norm = hypot($dx, $dy);
        $best = -1.0;
        $bi = -1;
        for ($i = $a + 1; $i < $b; $i++) {
            $d = $norm > 0 ? abs($dx * ($xy[$i][1] - $py) - $dy * ($xy[$i][0] - $px)) / $norm
                           : hypot($xy[$i][0] - $px, $xy[$i][1] - $py);
            if ($d > $best) {
                $best = $d;
                $bi = $i;
            }
        }
        if ($best > $tol) {
            $keep[$bi] = true;
            $stack[] = [$a, $bi];
            $stack[] = [$bi, $b];
        }
    }
    return array_values(array_filter($xy, fn($k) => $keep[$k], ARRAY_FILTER_USE_KEY));
}
 
function basemap_layer(array $lines, float $tol, bool $coast): array
{
    $out = [];
    foreach ($lines as $l) {
        $inside = 0;
        foreach ($l as [$x, $y]) {
            $inside += ($x > TRACK_BOX[0] - 2 && $x < TRACK_BOX[1] + 2 && $y > TRACK_BOX[2] - 2 && $y < TRACK_BOX[3] + 2);
        }
        if ($inside < 2) {
            continue;
        }
        $s = simplify($l, $tol);
        if (count($s) < 2) {
            continue;
        }
        if ($coast) {   /* îlots minuscules écartés */
            $xs = array_column($s, 0);
            $ys = array_column($s, 1);
            if (max(max($xs) - min($xs), max($ys) - min($ys)) < 0.12) {
                continue;
            }
        }
        $out[] = array_map(fn($p) => [round($p[0], 2), round($p[1], 2)], $s);
    }
    return $out;
}
 
function main(array $argv): int
{
    if (count($argv) !== 4) {
        fwrite(STDERR, "usage : geo_prep.php <rép. Natural Earth> <departements.geojson> <sortie>\n");
        return 2;
    }
    [, $ne, $geojson, $out] = $argv;
    ensure_dir($out);
    $coast = clip(read_shp("$ne/physical/ne_10m_coastline.shp"));
    $border = clip(read_shp("$ne/cultural/ne_10m_admin_0_boundary_lines_land.shp"));
    write_lines("$out/coast.bin", $coast);
    write_lines("$out/border.bin", $border);
 
    $g = json_in($geojson);
    $deptLines = [];
    $names = [];
    $bin = 'AUPO' . pack('V', count($g['features']));
    foreach ($g['features'] as $f) {
        $code = (string) $f['properties']['code'];
        $names[$code] = (string) $f['properties']['nom'];
        $polys = $f['geometry']['type'] === 'Polygon' ? [$f['geometry']['coordinates']] : $f['geometry']['coordinates'];
        $bin .= chr(strlen($code)) . $code . pack('V', count($polys));
        foreach ($polys as $poly) {
            foreach ($poly as $ring) {
                $deptLines[] = $ring;
            }
            $outer = $poly[0];                      /* masques : anneau extérieur seulement */
            $bin .= pack('V', count($outer));
            foreach ($outer as [$x, $y]) {
                $bin .= pack('ee', $x, $y);
            }
        }
    }
    write_lines("$out/dept.bin", clip($deptLines));
    atomic_write("$out/departements.bin", $bin);
    json_out("$out/departements.json", $names);
    json_out("$out/basemap.json", ['coast' => basemap_layer($coast, 0.025, true), 'border' => basemap_layer($border, 0.03, false)]);
    logm('fonds de carte :', (string) count($coast), 'côtes,', (string) count($border), 'frontières,', (string) count($names), 'départements');
    return 0;
}
 
exit(main($argv));