Arrêtez de saturer le disque avec des GeoTIFF annuels : comment SunCast streame du Parquet sans corrompre stdout

Réduisez la pression disque et simplifiez l’aval analytique en branchant un binaire GDAL/OpenMP sur un pipe binaire strict, consommé par PyArrow en un RecordBatch Parquet par jour.
Auteur·rice

Nicolas Decoopman

Date de publication

17 avril 2026

Mots clés

Géospatial, Production, HPC, C++17, GDAL, OpenMP, Parquet, PyArrow, Streaming, subprocess, Slurm, MNT

Imaginez que vous deviez publier, pour chaque département français, une grille raster complète des heures de lever et de coucher du soleil sur 365 ou 366 jours, à la résolution d’un MNT.

L’approche « fichier unique » mène vite à un GeoTIFF annuel avec des centaines de bandes : techniquement viable (GDAL, tuiles, BIGTIFF=IF_NEEDED), mais lourd pour le stockage, peu pratique pour l’analyse tabulaire, et pénible à versionner ou à pousser sur un cluster.

SunCast propose une autre voie : garder le moteur numérique en C++/GDAL/OpenMP, mais éjecter les résultats jour par jour vers Python via un flux binaire sur stdout, puis matérialiser du Parquet compressé (Snappy) avec PyArrow — sans écrire d’énorme raster intermédiaire.

Voici pourquoi stdout devient un contrat, et comment le code évite de le corrompre.

1 Le problème : le disque et le format ne suffisent pas

Pour beaucoup d’équipes géo ou climat, le réflexe est d’écrire un raster géoréférencé par bande ou par paire lever/coucher. Chaque pixel porte une information fine (minutes depuis minuit, souvent en int16), mais le produit dérivé est un cube énorme : largeur × hauteur × nombre de bandes × précision.

Même bien compressé (LZW, tuiles 512×512), cela reste un artefact monolithique. Pour l’aval data (jointures, filtres par jour, agrégations hors outils raster), on se retrouve à ré-encoder ou à multiplier les copies sur disque partagé HPC.

Le goulot n’est pas seulement la CPU : c’est le cycle de vie des fichiers entre le binaire natif et l’orchestration Python.

2 L’approche : un protocole binaire maison sur stdout

L’idée est de séparer les canaux :

  • stdout : uniquement des octets structurés (magic, dimensions, géotransform, puis blocs jour par jour).
  • stderr : messages humains (démarrage, progression, erreurs).

Le binaire solar_calculator expose --stream : sans --output, il ne produit pas de GeoTIFF annuel ; il appelle DemProcessor::streamBinaryOutput et pousse le flux vers stdout. Côté Python, run_solar_parquet.py ouvre un subprocess.Popen avec stdout=subprocess.PIPE, lit exactement le nombre d’octets attendu à chaque étape, et écrit un RecordBatch Parquet par jour.

2.1 Pourquoi c’est pertinent en production ?

  1. Mémoire Python bornée : on ne matérialise jamais l’année entière en tableaux Python ; une journée à la fois suffit pour write_batch.
  2. Échec rapide : si le process C++ meurt au milieu du flux, la lecture stricte lève une erreur au lieu d’écrire un Parquet silencieusement incohérent.
  3. HPC : le script Slurm du dépôt lance python -m src.solar.run_solar_parquet --index … : chaque tâche array consomme le même pattern (binaire + pipe + partition dept=).

3 Implémentation technique : du producteur au consommateur

3.1 Point d’entrée C++ : --stream et discipline stdout

La validation impose --output sauf en mode stream ; en stream, la configuration verbeuse évite std::cout pour ne pas polluer le flux binaire — seul std::cerr annonce le démarrage.

if (!streamMode && outputPath.empty()) {
    std::cerr << "Error: Output file is required (--output) unless in --stream mode" << std::endl;
    printUsage(argv[0]);
    return 1;
}
// ...
if (streamMode) {
    std::cerr << "Starting binary stream for " << inputPath << " (Year " << year << ")" << std::endl;
    success = processor.streamBinaryOutput(inputPath, year, timezoneOffset);
}

(Fichier src/main.cpp.)

3.2 En-tête puis blocs journaliers côté ProcessDEM

Après lecture du MNT en mémoire et calcul du nombre de jours (bissextile), le producteur écrit :

  1. Cinq octets magiques SOLAR.
  2. width, height, daysInYear en int32.
  3. Six double de géotransform GDAL.
  4. Pour chaque jour calendaire : currentDayOfYear (int32), puis deux tableaux int16 de width × height minutes (sentinel -1 pour nodata ou cas non calculables).

OpenMP parallélise le calcul par jour ; après chaque jour, std::cout.flush() vide le buffer vers le pipe. La progression (Processed day …) part sur stderr tous les dix jours.

const char magic[] = "SOLAR";
std::cout.write(magic, 5);
std::cout.write(reinterpret_cast<const char*>(&width), sizeof(int32_t));
std::cout.write(reinterpret_cast<const char*>(&height), sizeof(int32_t));
std::cout.write(reinterpret_cast<const char*>(&daysInYear), sizeof(int32_t));
std::cout.write(reinterpret_cast<const char*>(geoTransform), 6 * sizeof(double));
std::cout.flush();
// … pour chaque jour …
std::cout.write(reinterpret_cast<const char*>(&currentDayOfYear), sizeof(int32_t));
std::cout.write(reinterpret_cast<const char*>(sunriseBuffer), totalPixels * sizeof(int16_t));
std::cout.write(reinterpret_cast<const char*>(sunsetBuffer), totalPixels * sizeof(int16_t));
std::cout.flush();

(Fichier src/ProcessDEM.cpp, fonction streamBinaryOutput.)

Le compromis explicite : un flush par jour augmente le coût en syscalls, mais le consommateur avance de façon synchronisée et une panne laisse une trace plus localisable qu’un énorme buffer noyé.

3.3 Consommateur Python : read_exact, schéma, une ligne Parquet par jour

stderr du sous-processus est relié à sys.stderr : les logs C++ passent au terminal sans être mélangés au parseur binaire qui ne lit que proc.stdout.

def read_exact(proc, size):
    data = proc.stdout.read(size)
    if len(data) < size:
        raise EOFError(f"Expected {size} bytes, got {len(data)}")
    return data

proc = subprocess.Popen(
    cmd,
    stdout=subprocess.PIPE,
    stderr=sys.stderr,
    bufsize=1024 * 1024,
)

Le schéma Parquet encode chaque jour comme une ligne avec deux colonnes list<int16> de longueur width*height. Chaque itération construit un RecordBatch unique et appelle writer.write_batch(batch) avec compression Snappy.

schema = pa.schema([
    ('day', pa.int32()),
    ('sunrise', pa.list_(pa.int16())),
    ('sunset', pa.list_(pa.int16())),
])
writer = pq.ParquetWriter(parquet_file, schema, compression='snappy')
for _ in range(days_in_year):
    day_id = struct.unpack('i', read_exact(proc, 4))[0]
    sunrise_bytes = read_exact(proc, array_bytes)
    sunrise_array = np.frombuffer(sunrise_bytes, dtype=np.int16)
    # … idem sunset …
    batch = pa.RecordBatch.from_arrays([…], schema=schema)
    writer.write_batch(batch)

(Fichier src/solar/run_solar_parquet.py.)

Les partitions Hive-style (data/parquet/dept=<code>/data.parquet) et un metadata.json (dimensions, géotransform) complètent le livrable pour la viz ou les jointures spatiales ultérieures.

3.4 Une leçon de « contrat vivant »

Tout message de diagnostic en mode stream doit vivre sur stderr. Le constructeur DemProcessor logue l’état OpenMP sur std::cerr dans tous les cas : sans cela, un simple std::cout émis avant le flux aurait cassé le magic SOLAR — exactement le genre de bug que l’on voit quand la toolchain cluster diffère de la machine de développement.

Le commentaire de ProcessDEM.h décrit désormais le même en-tête et la même structure par jour que le code : la source de vérité reste l’implémentation et le lecteur Python, mais la doc d’en-tête suit.

4 Stratégie d’adoption : ne pas improviser le binaire

  • Spécifier le protocole : endianness (implicitement celle de la plateforme pour struct Python et écritures natives C++), taille des blocs, comportement sur nodata (-1).
  • Tests d’intégration légers : golden file sur les N premiers jours ou contrôle des longueurs de listes (le dépôt fournit src/utils/inspect_parquet.py pour un smoke test rapide).
  • Règle d’équipe : interdiction de std::cout dans le chemin partagé constructeur / bibliothèques quand --stream existe ; revue de code ou grep en CI.
  • Slurm : aligner --array sur la liste des départements cibles dans src/config.py, comme le script submit_job.slurm.

5 Conclusion : gains concrets et compromis

Gains : pas de GeoTIFF annuel obligatoire sur le chemin Parquet ; Parquet columnar + Snappy pour l’outil analytique ; charge CPU portée par OpenMP dans le natif ; empreinte mémoire Python stable.

Compromis : protocole ad hoc (pas Arrow IPC ni gRPC) ; flush par jour ; lignes Parquet très « larges » (listes de millions d’entiers) peu adaptées à du SQL naïf sans transformation ; couplage version binaire / lecteur Python.

Le dépôt complet, scripts et calculateur, est sur GitHub : https://github.com/NCSdecoopman/SunCast. Une présentation plus large du projet est aussi disponible dans la section projets.

6 Références et lectures complémentaires

  1. GDAL/OGR project. Raster I/O tutorial — lecture MNT, géotransform, bonnes pratiques d’écriture GeoTIFF (compression, tuiles, BigTIFF).
  2. Apache Arrow / PyArrow. Parquet file format — schémas, ParquetWriter, RecordBatch, compression Snappy.
  3. Python Software Foundation. subprocess — Popen, pipes, gestion de stdout/stderr en production.
  4. OpenMP ARB. OpenMP Application Programming Interface — parallélisation de boucles sur données raster en mémoire.