1. Introduction
1.1 Fondements biophysiques des indices spectraux
Depuis le lancement du satellite Earth Resources Technology Satellite (ERTS-1 / Landsat 1) en 1972, la télédétection optique spatiale constitue un outil central pour observer et quantifier les systèmes terrestres et aquatiques (Rouse et al., 1974 ; Huete, 1988 ; Gao, 1996). Les capteurs optiques mesurent le rayonnement solaire réfléchi par la surface terrestre dans des bandes de longueurs d’onde discrètes du spectre électromagnétique. Les matériaux de surface présentent des signatures diagnostiques d’absorption et de diffusion, déterminées par leurs propriétés biophysiques et biochimiques (Gitelson et al., 1996 ; McFeeters, 1996 ; Xu, 2006).
Par exemple, la végétation photosynthétique verte absorbe fortement la lumière rouge (640–670 nm) grâce aux pigments chlorophylliens, tandis que la structure spongieuse du mésophylle des feuilles diffuse intensément le proche infrarouge (PIR/NIR, 850–880 nm) (Tucker, 1979). À l’inverse, l’eau libre absorbe presque entièrement le rayonnement PIR tout en réfléchissant modérément les longueurs d’onde vertes (McFeeters, 1996). Les surfaces bâties imperméables présentent, quant à elles, une réflectance dans l’infrarouge à ondes courtes (SWIR) relativement élevée par rapport au PIR (Zha et al., 2003).
Pour condenser les données multispectrales de grande dimension en indicateurs biophysiques sensibles, les spécialistes de la télédétection élaborent des indices spectraux, c’est-à-dire des combinaisons algébriques de bandes conçues pour mettre en évidence certaines caractéristiques tout en atténuant les effets de l’angle d’éclairement solaire, des ombres topographiques et du bruit atmosphérique (Liu et Huete, 1995 ; Rondeaux et al., 1996). L’indice de végétation par différence normalisée (NDVI) (Rouse et al., 1974) a établi le format standard du rapport normalisé. Par la suite, des indices spécialisés ont été développés : suppression de l’influence du sol (SAVI, MSAVI, OSAVI), résistance aux effets atmosphériques (ARVI), limitation de la saturation en présence d’une biomasse élevée (EVI, GNDVI), délimitation de l’eau de surface (NDWI, MNDWI, AWEI), cartographie des zones bâties (NDBI, IBI), humidité du couvert végétal (NDMI, MSI), détection des sols nus (BSI) et couverture neigeuse (NDSI).
1.2 Difficultés techniques des flux de travail opérationnels
Même si le calcul algébrique d’un indice spectral est élémentaire, son intégration dans des chaînes scientifiques reproductibles soulève plusieurs difficultés récurrentes :
- Hétérogénéité des bandes selon les capteurs : les différentes missions satellitaires attribuent des désignations différentes à des canaux spectraux équivalents (par exemple, le PIR correspond à B08 sur Sentinel-2 MSI, à B5 sur Landsat 8/9 OLI et à B4 sur Landsat 4–7 TM). Un indexage des couches codé en dur peut introduire des erreurs difficiles à détecter.
- Incompatibilité spatiale : la combinaison de bandes de résolutions natives différentes (par exemple, bandes visibles/PIR à 10 m et SWIR à 20 m sur Sentinel-2), sans contrôles explicites, peut provoquer un décalage spatial silencieux ou un échec de l’exécution.
- Piège des nombres numériques (DN) : les rapports normalisés sont invariants à un facteur d’échelle commun, mais les formules sensibles à l’échelle et comportant des constantes additives ou des racines non linéaires (SAVI, EVI, MSAVI, AWEI, IBI, par exemple) sont fortement déformées lorsqu’elles sont calculées sur des DN entiers non convertis en réflectance physique de surface.
- Singularités et division par zéro : les rapports normalisés peuvent avoir des dénominateurs nuls ou presque nuls (|A + B| < 10⁻⁶), notamment dans les ombres sombres, les eaux profondes ou les pixels de bord. Ils produisent alors des valeurs non finies IEEE 754 (Inf, NaN), susceptibles de fausser les analyses zonales ultérieures.
- Sécurité du prototypage de formules : l’analyse dynamique naïve de chaînes au moyen de eval(parse(text = ...)) introduit des vulnérabilités critiques, notamment l’exécution de commandes arbitraires, ainsi que des plantages non gérés en environnement de production.
- Dépendances de packages hérités : l’écosystème spatial de R passe du package historique raster au moteur terra, plus performant et fondé sur C++ (Hijmans, 2023), alors que plusieurs packages existants conservent des dépendances anciennes et lourdes.
1.3 Objectifs et contributions de GeoIndexR
Pour résoudre ces difficultés opérationnelles, nous avons développé GeoIndexR. Le package fournit une interface de haut niveau unifiée (geo_index(), geo_indices(), geo_index_all()), des noyaux algébriques vectorisés (calc_*()), un moteur intelligent de résolution des bandes doté de préréglages pour les capteurs Sentinel-2 et Landsat 4–9, une arithmétique tolérante aux singularités, un compilateur de formules personnalisées fondé sur l’arbre syntaxique abstrait (AST) et un registre centralisé de métadonnées couvrant 16 indices. Ses seules dépendances sont R de base et terra.
2. Contexte et travaux connexes
Le tableau 1 compare les fonctionnalités techniques de GeoIndexR à celles d’autres logiciels libres de télédétection en R et en Python.
Tableau 1. Matrice comparative des fonctionnalités de GeoIndexR et des logiciels existants de télédétection.
| Fonctionnalité | GeoIndexR | RStoolbox | spyndex (Python) |
|---|---|---|---|
| Moteur principal | terra (C++) | raster / terra | xarray / rasterio |
| Interface unifiée | Oui (geo_index) | Oui (spectralIndices) | Oui (compute) |
| Noyaux vectorisés | Oui (calc_*) | Non | Non |
| Préréglages de capteurs | Sentinel-2, Landsat 4–9, générique | Sentinel-2, Landsat 5, 7, 8 | Catalogue étendu |
| Heuristiques par expressions régulières | Oui | Partielles | Aucune |
| Alerte sur la mise à l’échelle | Oui, automatique | Non | Partielle |
| Division sécurisée | Oui, tolérance 10⁻⁶ | Non (Inf/NaN) | Non |
| Moteur AST sécurisé | Oui, liste blanche | Non | Non |
| Catalogue de métadonnées | Oui | Limité | Oui (JSON) |
| Dépendances | R de base + terra | Importantes (caret, etc.) | NumPy, Pandas |
RStoolbox (Leutner et al., 2023) a contribué au développement de l’analyse spectrale dans R, mais son architecture reste liée à l’ancien package raster et comporte des dépendances graphiques importantes. Le package Python spyndex (Montero et al., 2023) propose un catalogue étendu, mais n’est pas disponible nativement dans les chaînes de traitement R et ne dispose pas d’une désinfection des formules au niveau AST. GeoIndexR comble cette lacune en proposant une solution moderne, sécurisée et native à terra.
3. Architecture et conception de GeoIndexR
3.1 Objectifs de conception
GeoIndexR privilégie la programmation défensive, une modularité sans surcharge et une gestion de la mémoire évolutive. Les opérations manipulent directement les pointeurs C++ sous-jacents de terra sans créer de copies mémoire inutiles, afin d’assurer de bonnes performances sur de grands rasters.
3.2 Architecture logicielle
L’architecture de GeoIndexR comprend cinq couches modulaires : le registre des indices, le moteur de résolution des bandes, le noyau arithmétique numérique, le moteur de formules personnalisées fondé sur AST et la suite de diagnostic de post-traitement.
Lors de l’appel à geo_index(), le pipeline coordonne la validation des rasters, la résolution des bandes, l’application du facteur d’échelle et le calcul par le noyau.
3.3 Registre des indices spectraux
La fonction index_registry() constitue la base de connaissances du package et maintient des métadonnées structurées pour 16 indices évalués par les pairs. Les utilisateurs peuvent consulter les indices pris en charge avec list_indices() ou récupérer des profils techniques détaillés au moyen de get_index_meta().
3.4 Correspondance des bandes et préréglages des capteurs
GeoIndexR prend en charge trois mécanismes de résolution des bandes via resolve_bands() : l’affectation explicite par l’utilisateur, les préréglages prédéfinis des capteurs (Sentinel-2, Landsat 8/9 et Landsat 4–7) et les heuristiques automatiques fondées sur des expressions régulières.
3.5 Mise à l’échelle de la réflectance
Le paramètre scale_factor = 10000 convertit les nombres numériques entiers en réflectance physique de surface [0, 1]. Pour les indices sensibles à l’échelle, tels que EVI ou SAVI, GeoIndexR examine automatiquement les plages de valeurs du raster. Si des valeurs non mises à l’échelle supérieures à 10 sont détectées en l’absence de facteur d’échelle, un avertissement de diagnostic explicite est émis.
3.6 Formules d’indices personnalisées
Avant l’évaluation d’une formule personnalisée, l’arbre d’expression est inspecté récursivement et comparé à des listes strictes d’opérateurs autorisés (+, −, *, /, ^), de fonctions (sqrt, log, exp, abs, etc.) et de symboles déclarés. Ce mécanisme, implémenté dans formula_engine.R, neutralise les risques d’injection de code.
3.7 Validation et gestion des erreurs
Avant le calcul, check_spatial_compatibility() vérifie que les couches ont des dimensions et des emprises spatiales identiques. La routine safe_divide() applique une tolérance aux singularités :
« safe_divide(x, y, tol = 1e-6) = NA if |y| < tol or is.na(x, y) (1) »Équation (1) — division sécurisée
Sur les objets SpatRaster, safe_divide() exécute des appels compilés à terra::ifel(), évitant les copies de tableaux en mémoire et convertissant proprement les singularités en valeurs NA standard.
4. Indices spectraux
4.1 Catalogue des indices Implémentés
Le tableau 2 présente les formules mathématiques, les paramètres et les références bibliographiques fondamentales des 16 indices pris en charge par GeoIndexR.
Tableau 2. Formulation mathématique et paramètres définis par défaut des indices spectraux pris en charge.
| Indice (domaine) | Formule mathématique | Paramètres | Référence |
|---|---|---|---|
| NDVI (végétation) | (PIR − ROUGE) / (PIR + ROUGE) | — | Rouse, 1974 |
| SAVI (végétation) | ((PIR − ROUGE) / (PIR + ROUGE + L)) × (1 + L) | L = 0,5 | Huete, 1988 |
| EVI (végétation) | G × (PIR − ROUGE) / (PIR + C1×ROUGE − C2×BLEU + L) | G = 2,5 | Liu, 1995 |
| MSAVI (végétation) | (2×PIR + 1 − √((2×PIR + 1)² − 8×(PIR − ROUGE))) / 2 | — | Qi, 1994 |
| OSAVI (végétation) | ((PIR − ROUGE) / (PIR + ROUGE + θ)) × (1 + θ) | θ = 0,16 | Rondeaux, 1996 |
| ARVI (végétation) | (PIR − rb) / (PIR + rb), rb = ROUGE − γ×(BLEU − ROUGE) | γ = 1,0 | Kaufman, 1992 |
| GNDVI (végétation) | (PIR − VERT) / (PIR + VERT) | — | Gitelson, 1996 |
| NDWI (eau) | (VERT − PIR) / (VERT + PIR) | — | McFeeters, 1996 |
| MNDWI (eau) | (VERT − SWIR1) / (VERT + SWIR1) | — | Xu, 2006 |
| AWEI (eau) | 4×(VERT − SWIR1) − (0,25×PIR + 2,75×SWIR2) | — | Feyisa, 2014 |
| NDBI (urbain) | (SWIR1 − PIR) / (SWIR1 + PIR) | — | Zha, 2003 |
| IBI (urbain) | (NDBI − (SAVI + MNDWI)/2) / (NDBI + (SAVI + MNDWI)/2) | L = 0,5 | Xu, 2007 |
| NDMI (humidité) | (PIR − SWIR1) / (PIR + SWIR1) | — | Gao, 1996 |
| MSI (humidité) | SWIR1 / PIR | — | Rock, 1986 |
| BSI (sol) | ((SWIR1 + ROUGE) − (PIR + BLEU)) / ((SWIR1 + ROUGE) + (PIR + BLEU)) | — | Rikimaru, 2002 |
| NDSI (neige) | (VERT − SWIR1) / (VERT + SWIR1) | — | Hall, 1995 |
4.2 Fondements biophysiques et formulations
Le catalogue couvre plusieurs domaines de suivi environnemental : vigueur de la végétation et absorption de la chlorophylle (NDVI, GNDVI), réduction de l’influence du sol dans les couverts végétaux clairsemés (SAVI, MSAVI, OSAVI), résistance aux aérosols atmosphériques (ARVI), limitation de la saturation en biomasse élevée (EVI), délimitation des eaux libres (NDWI, MNDWI, AWEI), cartographie des surfaces urbaines imperméables (NDBI, IBI), teneur en eau du couvert végétal (NDMI, MSI), détection des sols nus (BSI) et couverture neigeuse (NDSI).
5. Tests logiciels
5.1 Cadre de tests unitaires
GeoIndexR comprend une suite complète de tests reposant sur testthat (Wickham, 2023), composée de 47 blocs de tests unitaires répartis dans 9 fichiers. Elle vérifie la précision des calculs d’indices, la résolution des préréglages de capteurs, l’analyse des formules AST, les cas limites, la propagation des valeurs NA et la gestion des facteurs d’échelle. Le taux de réussite rapporté est de 100 %, en moins de deux secondes.
5.2 Intégrité numérique et protection contre les singularités
Les tests vérifient que les dénominateurs nuls (0/0), les dénominateurs proches de zéro (|den| < 10⁻⁶) et les nombres non finis (Inf, NaN) sont convertis en NA sans déclencher d’exceptions non interceptées, ce qui préserve l’intégrité des données pour les modèles spatiaux ultérieurs.
5.3 Génération de données de référence synthétiques
Pour permettre des tests reproductibles sans télécharger des scènes satellitaires de plusieurs gigaoctets, get_example_data() génère un objet SpatRaster à 6 bandes (10 × 10 pixels, EPSG:32631). Il simule des réflectances de surface réalistes pour une canopée forestière (PIR élevé), des plans d’eau (forte absorption dans le vert selon le scénario simulé) et des zones urbaines ou des sols nus (SWIR1 élevé), avec un pixel NA en bordure pour vérifier la propagation des valeurs manquantes.
6. Exemples reproductibles
Indice personnalisé (AST)
ci <- geo_index_custom(
img, formula = '(nir / green) - 1',
bands = c(nir = 'nir', green = 'green')
)6.1 Flux de travail à un ou plusieurs indices
Le calcul d’indices spectraux ne nécessite que quelques lignes de code:
Indices NDVI / NDWI / NDBI / NDMI
library(GeoIndexR); library(terra)
img <- get_example_data()
ndvi <- geo_index(img, 'NDVI')
stack <- geo_indices(img, c('NDVI', 'NDWI', 'NDBI', 'NDMI'))6.2 Extraction exhaustive par lots
La fonction geo_index_all() calcule automatiquement tous les indices compatibles avec les bandes disponibles:
Extraction exhaustive
all_idx <- geo_index_all(img)
names(all_idx) #16 couches calculées automatiquement6.3 Profilage statistique des distributions
La fonction index_summary() calcule des statistiques descriptives pour les couches d’indices, comme illustré dans le tableau 3.
Tableau 3. Statistiques descriptives récapitulatives calculées sur l’ensemble des couches d’indices à l’aide de la fonction index_summary().
| Indice | Min | Max | Moyenne | Médiane | Écart-type | NA (%) |
|---|---|---|---|---|---|---|
| NDVI | -0,360 | 0,711 | 0,114 | 0,1095 | 0,177 | 2,06 |
| NDWI | -0,666 | 0,366 | -0,129 | -0,1925 | 0,191 | 2,06 |
| NDBI | -0,543 | 0,449 | -0,010 | 0,0050 | 0,152 | 2,06 |
| NDMI | -0,449 | 0,543 | 0,010 | -0,0050 | 0,152 | 2,06 |
6.4 Visualisation géospatiale thématique
La fonction plot_index() produit des cartes thématiques des indices spectraux, en appliquant des palettes adaptées aux catégories : végétation, eau, zones urbaines et humidité.
6.5 Prototypage sécurisé d’indices personnalisés
Le moteur AST permet de prototyper un indice de chlorophylle personnalisé (CIgreen = PIR/Vert − 1) tout en validant la formule avant son évaluation:
7. Evaluation de performances
7.1 Évolutivité mémoire et traitement hors mémoire
Grâce à l’architecture de traitement par blocs de terra, GeoIndexR traite les grands rasters en flux, à l’aide de caches temporaires sur disque lorsque les scènes dépassent la mémoire vive disponible. Cela permet de traiter de manière fiable des tuiles Sentinel-2 complètes (10 980 × 10 980 pixels) sans débordement mémoire.
7.2 Vitesse de calcul et vectorisation
Les noyaux vectorisés calc_*() utilisent des opérations arithmétiques compilées en C++. Le calcul du NDVI sur une grille de 10 000 × 10 000 cellules (100 millions de cellules) prend moins de 1,8 seconde sur un matériel standard, selon les résultats rapportés par l’article.
8. Discussion
GeoIndexR traite plusieurs causes récurrentes d’échec dans les calculs de télédétection en combinant la résolution automatique des bandes selon les capteurs, la gestion défensive des singularités, l’analyse sécurisée des formules au niveau AST et une interopérabilité directe avec terra. Le cadre vise à garantir la reproductibilité mathématique tout en conservant une empreinte logicielle minimale.
9. Limites
GeoIndexR n’effectue pas de modélisation du transfert radiatif atmosphérique (par exemple 6S ou Sen2Cor). Les bandes d’entrée doivent avoir la même résolution spatiale et la même emprise avant le calcul des indices. Les utilisateurs qui travaillent avec des bandes de résolutions différentes doivent les rééchantillonner, par exemple avec terra::resample().
10. Travaux futures
Les prochaines versions devraient étendre les préréglages de capteurs aux constellations commerciales (PlanetScope SuperDove, WorldView-3), ajouter la lecture par plages des GeoTIFF optimisés pour le cloud (COG) et mettre en œuvre un traitement parallèle optionnel par blocs pour les cubes raster multitemporels.
11. Conclusion
GeoIndexR fournit un cadre moderne, robuste et extensible pour le calcul d’indices spectraux dans R. En reliant les cubes raster bruts à un traitement défensif des indices biophysiques, il favorise des flux de travail reproductibles et standardisés en observation de la Terre et en science des données géospatiales.
Disponibilité des données et du code
GeoIndexR est un package R libre distribué sous licence GPL-3. Le code source, la documentation, les vignettes et les tests sont disponibles à l’adresse : https://github.com/sowsalim01/GeoIndexR.
Remerciements
Le développement de GeoIndexR a été inspiré par la communauté de l’observation de la Terre en accès libre et par les auteurs du package terra.
Aucun commentaire publié pour le moment.