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.