<figure><img src="image_1.png" /></figure>

# Traitement d'images satellites avec Python

Première édition

Samuel Foucher

Philippe Apparicio

Yacine Bouroubi

Mickaël Germain

Étienne Clabaut

2025-03-24
# Table des matières

## Préface 1

Un manuel sous la forme d'une ressource éducative libre 2

Comment lire ce manuel? 4

Comment utiliser les données du livre pour reproduire les exemples? 4

Structure du livre 5

Remerciements 5

Introduction aux images de télédétection 6

Ressources en ligne 6

Listes des *librairies* utilisés 6

## À propos des auteurs 7

## Partie 1. Importation, manipulation et visualisation de données satellitaires 8

### 1 Introduction au langage Python 9

1.1 Les distributions 9

1.2 Les styles de programmation en Python 10

1.2.1 Les outils de programmation 10

1.3 Bonnes pratiques 11

1.3.1 Création d'un environnement virtuel 11

1.3.2 Création d'un environnement de travail local (avancé) 11

1.4 Les structures de base en Python 13

1.4.1 Les listes 13

1.4.2 Les tuples 13

1.4.3 Les ensembles (Sets) 13

1.5 Dictionnaires 13

1.6 Programmation objet 14

### 2 Importation et manipulation de données spatiales 15

2.1 Préambule 15

2.1.1 Objectifs 15

2.1.2 Bibliothèques 15

2.1.3 Données 16

2.2 Importation d'images 16

2.2.1 Formats des images 17

2.2.2 Métadonnées des images 21

2.3 Manipulation des images 22

2.3.1 Manipulation de la matrice de pixels 22
<table><tbody><tr><td>2.3.2</td><td>Information de base</td><td>22</td></tr><tr><td>2.3.3</td><td>Découpage et indexation de la matrice</td><td>24</td></tr><tr><td>2.3.4</td><td>Changement de projection cartographique (à venir)</td><td>26</td></tr><tr><td>2.4</td><td>Données en géoscience</td><td>26</td></tr><tr><td>2.4.1</td><td>xarray</td><td>26</td></tr><tr><td><b>3</b></td><td><b>Réhaussement et visualisation d'images</b></td><td><b>28</b></td></tr><tr><td>3.1</td><td>Préambule</td><td>28</td></tr><tr><td>3.1.1</td><td>Objectifs</td><td>28</td></tr><tr><td>3.1.3</td><td>Bibliothèques</td><td>28</td></tr><tr><td>3.1.4</td><td>Données</td><td>29</td></tr><tr><td>3.2</td><td>Visualisation en Python</td><td>30</td></tr><tr><td>3.3</td><td>Réhaussements visuels</td><td>33</td></tr><tr><td>3.3.1</td><td>Statistiques d'une image</td><td>33</td></tr><tr><td>3.3.2</td><td>Réhaussements linéaires</td><td>37</td></tr><tr><td>3.3.3</td><td>Réhaussements non linéaires</td><td>41</td></tr><tr><td>3.3.4</td><td>Composés colorés</td><td>46</td></tr><tr><td><b>4</b></td><td><b>Transformations spectrales</b></td><td><b>49</b></td></tr><tr><td>4.1</td><td>Préambule</td><td>49</td></tr><tr><td>4.1.1</td><td>Objectifs</td><td>49</td></tr><tr><td>4.1.2</td><td>Librairies</td><td>49</td></tr><tr><td>4.1.3</td><td>Images utilisées</td><td>50</td></tr><tr><td>4.2</td><td>Qu'est ce que l'information spectrale?</td><td>51</td></tr><tr><td>4.3</td><td>Indices spectraux</td><td>52</td></tr><tr><td>4.4</td><td>Quiz</td><td>57</td></tr><tr><td><b>5</b></td><td><b>Transformations spatiales</b></td><td><b>58</b></td></tr><tr><td>5.1</td><td>Préambule</td><td>58</td></tr><tr><td>5.1.1</td><td>Objectifs</td><td>58</td></tr><tr><td>5.1.2</td><td>Librairies</td><td>58</td></tr><tr><td>5.1.3</td><td>Images utilisées</td><td>59</td></tr><tr><td>5.2</td><td>Analyse fréquentielle</td><td>60</td></tr><tr><td>5.2.1</td><td>La transformée de Fourier</td><td>60</td></tr><tr><td>5.2.2</td><td>Filtrage fréquentiel</td><td>61</td></tr><tr><td>5.2.3</td><td>L'aliasing</td><td>62</td></tr><tr><td>5.3</td><td>Filtrage d'image</td><td>65</td></tr><tr><td>5.3.1</td><td>Filtrage linéaire stationnaire</td><td>66</td></tr><tr><td>5.4</td><td>Gestion des bordures</td><td>69</td></tr><tr><td>5.4.1</td><td>Filtrage adaptatif</td><td>71</td></tr><tr><td>5.5</td><td>Segmentation</td><td>74</td></tr><tr><td>5.5.1</td><td>Super-pixel</td><td>75</td></tr><tr><td>5.5.2</td><td>Fusion des segments par graphe de proximité</td><td>78</td></tr><tr><td>5.5.3</td><td>Approche objet</td><td>81</td></tr></tbody></table>
## Partie 3. Classifications d'images

83

## 6 Classifications d'images supervisées

84

<table><tr><td>6.1 Préambule</td><td>84</td></tr><tr><td>6.1.1 Objectifs</td><td>84</td></tr><tr><td>6.1.2 Librairies</td><td>84</td></tr><tr><td>6.1.3 Images utilisées</td><td>85</td></tr><tr><td>6.2 Principes généraux</td><td>86</td></tr><tr><td>6.2.1 Comportement d'un modèle</td><td>86</td></tr><tr><td>6.2.2 Pipeline</td><td>87</td></tr><tr><td>6.2.3 Construction d'un ensemble d'entraînement</td><td>87</td></tr><tr><td>6.3 Analyse préliminaire des données</td><td>94</td></tr><tr><td>6.4 Mesures de performance d'une méthode de classification</td><td>98</td></tr><tr><td>6.5 Méthodes non paramétriques</td><td>99</td></tr><tr><td>6.5.1 Méthode des parallélépipèdes</td><td>99</td></tr><tr><td>6.5.2 Plus proches voisins</td><td>102</td></tr><tr><td>6.5.3 Méthodes par arbre de décision</td><td>107</td></tr><tr><td>6.6 Méthodes paramétriques</td><td>110</td></tr><tr><td>6.6.1 Méthode Bayésienne naïve</td><td>111</td></tr><tr><td>6.6.2 Analyse discriminante quadratique (ADQ)</td><td>112</td></tr></table>

## Bibliographie

115
# Liste des Figures

<table><tr><td>1</td><td>Licence Creative Commons du livre</td><td>3</td></tr><tr><td>2</td><td>Téléchargement de l'intégralité du livre</td><td>5</td></tr><tr><td>1.1</td><td>Client Jupyter Lab</td><td>12</td></tr><tr><td>2.1</td><td>La librairie NumPy est le fondement de nombreuses bibliothèques scientifiques (d'après (Harris 2020)).</td><td>23</td></tr><tr><td>2.2</td><td>Vue d'ensemble des opérations de base des matrices avec NumPy</td><td>24</td></tr><tr><td>2.3</td><td>Organisation d'un Dataset dans xarray</td><td>27</td></tr><tr><td>4.1</td><td>Positions des bandes spectrales pour quelques capteurs (source)</td><td>51</td></tr><tr><td>4.2</td><td>Visualisation des points d'une image Sentinel-2 pour trois classes</td><td>57</td></tr><tr><td>5.1</td><td>Graphe d'adjacence de régions, d'après (Jaworek-Korjakowska (2018)). Chaque nœud est un segment, un lien est formé uniquement si les segments se touchent (par exemple le segment 6 ne touche que la région 5). La fonction graph.rag_mean_color produit un graphe à partir d'une segmentation et de l'image originale. Chaque nœud tient la couleur de chaque segment dans un attribut appelé 'mean color'.</td><td>79</td></tr><tr><td>6.1</td><td>Exemples de sur et sous-apprentissage.</td><td>87</td></tr><tr><td>6.2</td><td>Étapes standards dans un entraînement.</td><td>88</td></tr><tr><td>6.3</td><td>Frontières de décision pour le classificateur K-NN</td><td>104</td></tr><tr><td>6.4</td><td>Frontières de décision pour des arbres de décision de différente profondeur</td><td>108</td></tr><tr><td>6.5</td><td>Frontières de décision pour un classificateur Bayésien naif</td><td>111</td></tr></table>
# Liste des Tables

<table><tr><td>2.1</td><td>Type de données de NumPy</td><td>24</td></tr><tr><td>4.1</td><td>Noms des bandes Sentinel-2</td><td>53</td></tr></table>
# Préface

**Résumé** : Ce livre vise à décrire une panoplie de méthodes de traitement d'images satellites avec le langage Python. Celles et ceux souhaitant migrer progressivement d'un autre logiciel d'imagerie et de télédétection vers Python trouveront dans cet ouvrage les éléments pour une transition en douceur. La philosophie de ce livre est de donner toutes les clefs de compréhension et de mise en œuvre des méthodes abordées dans Python. La présentation des méthodes est basée sur une approche compréhensive et intuitive plutôt que mathématique, sans pour autant négliger la rigueur mathématique ou statistique. Des rappels sur les fondements en télédétection pourront apparaître au besoin afin d'éclairer les approches techniques. Plusieurs éditions régulières sont prévues sachant que ce domaine évolue constamment.

Ce projet est en cours d'écriture et le contenu n'est pas complet

<figure><img src="image_2.png" /></figure>

**Remerciements** : Ce manuel a été réalisé avec le soutien de la fabriqueREL. Fondée en 2019, la fabriqueREL est portée par divers établissements d'enseignement supérieur du Québec et agit en collaboration avec les services de soutien pédagogique et les bibliothèques. Son but est de faire des ressources éducatives libres (REL) le matériel privilégié en enseignement supérieur au Québec.

**Mise en page** : Samuel Foucher, Philippe Apparicio et Marie-Hélène Gadbois Del Carpio.

© Samuel Foucher, Philippe Apparicio, Yacine Bouroubi, Mickaël Germain et Étienne Clabaut.

**Pour citer cet ouvrage** : Foucher S., Apparicio P., Bouroubi Y., Germain M. et Clabaut, E. (2025). *Traitement d'images satellites avec Python*. Université de Sherbrooke, Département de géomatique appliquée. fabriqueREL.
Licence CC BY-SA.
<figure><img src="image_3.png" /></figure>

Sauf indications contraires, le contenu de ce manuel électronique
est disponible en vertu des termes de la [Licence Creative Commons](https://creativecommons.org/licenses/by-sa/4.0/) Attribution - Partage dans les mêmes conditions 4.0 International.

Vous êtes autorisé-e à :

**Partager** – copier, distribuer et communiquer le matériel par tous moyens et sous tous formats.

**Adapter** – remixer, transformer et créer à partir du matériel pour toute utilisation, y compris commerciale.

Selon les conditions suivantes :

**Paternité** – Vous devez citer le nom des auteurs originaux.

**Mêmes conditions** – Si vous remixez, transformez, ou créez à partir du matériel composant l'Œuvre originale, vous devez diffuser l'Œuvre modifiée avec la même licence.

UDS

Université de
Sherbrooke

<figure><img src="image_4.png" /></figure>

fabrique REL
RESSOURCES ÉDUCATIVES LIBRES

# Un manuel sous la forme d'une ressource éducative libre

## Pourquoi un manuel sous licence libre?

Les logiciels libres sont aujourd'hui très répandus. Comparativement aux logiciels propriétaires, l'accès au code source permet à quiconque de l'utiliser, de le modifier, de le dupliquer et de le partager. Le logiciel Python, dans lequel sont mises en œuvre les méthodes de traitement d'images satellites décrites dans ce livre, est d'ailleurs à la fois un langage de programmation et un logiciel libre (sous la licence publique générale GNU GPL2). Par analogie aux logiciels libres, il existe aussi des **ressources éducatives libres** (REL) « dont la licence accorde les permissions désignées par les 5R (**Retenir** — **Réutiliser** — **Réviser** — **Remixer** — **Redistribuer**) et donc permet nécessairement la modification » (*fabriqueREL*). La licence de ce livre, CC BY-SA (figure 1), permet donc de :

- **Retenir**, c'est-à-dire télécharger et imprimer gratuitement le livre. Notez qu'il aurait été plutôt surprenant d'écrire un livre payant sur un logiciel libre et donc gratuit. Aussi, nous aurions été très embarrassés que des personnes étudiantes avec des ressources financières limitées doivent payer pour avoir accès au livre, sans pour autant savoir préalablement si le contenu est réellement adapté à leurs besoins.

- **Réutiliser**, c'est-à-dire utiliser la totalité ou une section du livre sans limitation et sans compensation financière. Cela permet ainsi à d'autres personnes enseignantes de l'utiliser dans le cadre d'activités pédagogiques.

- **Réviser**, c'est-à-dire modifier, adapter et traduire le contenu en fonction d'un besoin pédagogique précis puisqu'aucun manuel n'est parfait, tant s'en faut! Le livre a d'ailleurs été écrit intégralement dans Python avec Quatro. Quiconque peut ainsi télécharger gratuitement le code source du livre sur GitHub et le modifier à sa guise (voir l'encadré intitulé *Suggestions d'adaptation du manuel*).

- **Remixer**, c'est-à-dire « combiner la ressource avec d'autres ressources dont la licence le permet aussi pour créer une nouvelle ressource intégrée » (*fabriqueREL*).

- **Redistribuer**, c'est-à-dire distribuer, en totalité ou en partie le manuel ou une version révisée sur d'autres canaux que le site Web du livre (par exemple, sur le site Moodle de votre université ou en faire une version imprimée).
La licence de ce livre, CC BY-SA (figure 1), oblige donc à :

- Attribuer la paternité de l'auteur dans vos versions dérivées, ainsi qu'une mention concernant les grandes modifications apportées, en utilisant la formulation suivante :

Samuel Foucher, Apparicio Philippe, Mickaël Germain, Yacine Bouroubi et Étienne Clabaut (2025). *Traitement d'images satellites avec Python*. Université de Sherbrooke, Département de géomatique appliquée. fabriqueREL.
Licence CC BY-SA.

- Utiliser la même licence ou une licence similaire à toutes versions dérivées.

<figure><img src="image_5.png" /><figcaption>FIGURE 1 - Licence Creative Commons du livre</figcaption></figure>

<figure><img src="image_6.png" /></figure>

## Astuce

### Suggestions d'adaptation du manuel

Pour chaque méthode de traitement d'image abordée dans le livre, une description détaillée et une mise en œuvre dans Python sont disponibles. Par conséquent, plusieurs adaptations du manuel sont possibles :

- Conserver uniquement les chapitres sur les méthodes ciblées dans votre cours.

- En faire une version imprimée et la distribuer aux personnes étudiantes.

- Modifier la description d'une ou de plusieurs méthodes en effectuant les mises à jour directement dans les chapitres.

- Insérer ses propres jeux de données dans les sections intitulées *Mise en œuvre dans Python*.

- Modifier les tableaux et figures.

- Ajouter une série d'exercices.

- Modifier les quiz de révision.

- Rédiger un nouveau chapitre.

- Modifier des syntaxes en Python. Plusieurs *librairies* Python peuvent être utilisées pour mettre en œuvre telle ou telle méthode. Ces derniers évoluent aussi très vite et de nouvelles *librairies* sont proposées fréquemment! Par conséquent, il peut être judicieux de modifier une syntaxe Python du livre en fonction de ses habitudes de programmation en Python (utilisation d'autres *librairies* que ceux utilisés dans le manuel par exemple) ou de bien mettre à jour une syntaxe à la suite de la parution d'une nouvelle *librairie* plus performante ou intéressante.

- Toute autre adaptation qui permet de répondre au mieux à un besoin pédagogique.
# Comment lire ce manuel?

Le livre comprend plusieurs types de blocs de texte qui en facilitent la lecture.

<figure><img src="image_7.png" /></figure>

## Package

### Bloc packages

Habituellement localisé au début d'un chapitre, il comprend la liste des packages Python utilisés pour un chapitre.

<figure><img src="image_8.png" /></figure>

## Objectif

### Bloc objectif

Il comprend une description des objectifs d'un chapitre ou d'une section.

<figure><img src="image_9.png" /></figure>

## Note

### Bloc notes

Il comprend une information secondaire sur une notion, une idée abordée dans une section.

<figure><img src="image_10.png" /></figure>

## Aller plus loin

### Bloc pour aller plus loin

Il comprend des références ou des extensions d'une méthode abordée dans une section.

<figure><img src="image_11.png" /></figure>

## Astuce

### Bloc astuce

Il décrit un élément qui vous facilitera la vie : une propriété statistique, un package, une fonction, une syntaxe Python.

<figure><img src="image_12.png" /></figure>

## Attention

### Bloc attention

Il comprend une notion ou un élément important à bien maîtriser.

<figure><img src="image_13.png" /></figure>

## Exercice

### Bloc exercice

Il comprend un court exercice de révision à la fin de chaque chapitre.

# Comment utiliser les données du livre pour reproduire les exemples?

Ce livre comprend des exemples détaillés et appliqués en Python pour chacune des méthodes abordées. Ces exemples se basent sur des jeux de données ouverts et mis à disposition avec le livre. Ils sont disponibles sur le repo GitHub dans le sous-dossier data, à l'adresse https://github.com/serie-tele-pyton/TraitementImagesVol1/tree/main/data.
Une autre option est de télécharger le *repo* complet du livre directement sur *GitHub* (https://github.com/serie-tele-pyton/TraitementImagesVol1) en cliquant sur le bouton **Code**, puis le bouton **Download ZIP** (figure 2). Les données se trouvent alors dans le sous-dossier nommé **data**.

<figure><img src="image_14.png" /><figcaption>FIGURE 2 – Téléchargement de l’intégralité du livre</figcaption></figure>

## Structure du livre

Le livre est organisé autour de quatre grandes parties.

**Partie 1. Importation et manipulation de données spatiales.** Dans cette première partie, nous voyons comment importer, manipuler, visualiser et exporter des données spatiales de type image (ou de type matriciel) avec Python, principalement avec les *packages rasterio*, *xarray* et * numpy* (chapitre 2). Ce chapitre vous permettra de maîtriser la manipulation à bas niveau de différents types d'imagerie. Différents exemples et exercises sont disponibles avec différents capteurs satellites (multi-spectral, RGB-NIR, SAR, etc.)

**Partie 2. Transformations des données spatiales.** Cette deuxième partie comprend deux chapitres : les transformations spectrales (chapitre 4) et les transformations spatiales (chapitre 5).

**Partie 3. Classifications d'images.** Cette troisième partie comprend deux chapitres : les classifications supervisées (chapitre 6) et non supervisées (à venir dans une prochaine édition).

**Partie 4. Données massives.** (à venir dans une édition future). Cette quatrième et dernière partie comprend un seul chapitre qui est dédié aux plateformes de mégadonnées, notamment Google Earth Engine.

## Remerciements

De nombreuses personnes ont contribué à l'élaboration de ce manuel.

Ce projet a bénéficié du soutien pédagogique et financier de la *fabrique REL* (ressources éducatives libres). Les différentes rencontres avec le comité de suivi nous ont permis de comprendre l'univers des ressources éducatives
libres (REL) et notamment leurs fameux 5R (Retenir — Réutiliser — Réviser — Remixer — Redistribuer), de
mieux définir le besoin pédagogique visé par ce manuel, d'identifier des ressources pédagogiques et des outils
pertinents pour son élaboration. Ainsi, nous remercions chaleureusement les membres de la fabrique REL pour
leur soutien inconditionnel :

— José-Miguel Escobar-Zuniga, bibliothécaire à l'Université de Sherbrooke.

— Marianne Dubé, coordonnatrice de la fabriqueREL, Université de Sherbrooke.

— Claude Potvin, conseiller en formation, Service de soutien à l'enseignement, Université Laval.

Nous remercions chaleureusement les personnes étudiantes du Baccalauréat en géomatique appliquée à l'environnement et du Microprogramme de 1er cycle en géomatique appliquée du Département de géomatique appliquée de l'Université de Sherbrooke qui utilisent le manuel dans le cadre de cours.

Introduction aux images de télédétection

L'imagerie numérique a pris une place importante dans notre vie de tous les jours depuis une quinzaine d'années.
Ces images sont prises généralement au niveau du sol (imagerie proximale) avec seulement trois couleurs dans le
domaine de la vision humaine (rouge, vert et bleu). Dans la suite du manuel, on parlera d'images du domaine de
la vision par ordinateur ou images en vision pour faire plus court.

Les images de télédétection ont des particularités et des propriétés qui les différencient des images “classiques”,
dont les cinq principales sont:

1. Les images sont géoréférencées. Cela signifie que pour chaque pixel nous pouvons y associer une position géographique ou cartographique.

2. Le point de vue est très différent. Ces images sont prises avec une vue d'en haut (Nadir) ou oblique avec une distance qui peut être très grande; on parle d'images distales.

3. Elles possèdent plus que 3 bandes. Contrairement aux images en vision, les images de télédétection possèdent bien souvent plus que trois bandes. Il n'est pas rare de trouver quatre bandes (Pléiade), 13 bandes (Sentinel-2, Landsat) et même 200 bandes pour des capteurs hyperspectraux.

4. Elles peuvent être calibrées. Les valeurs numérique de l'image peuvent être converties en quantités physiques (luminance, réflectance, section efficace, etc.) via une fonction de calibration.

5. Elles sont de grande taille. Il n'est pas rare de manipuler des images qui font plusieurs dizaines de milliers de pixels en dimension.

Ressources en ligne

Listes des librairies utilisés

Dans ce livre, nous utilisons de nombreux packages Python que vous pouvez installer en une seule fois (voir section 1.3.1) ou chapitre par chapitre.
# À propos des auteurs

**Samuel Foucher** est professeur au **Département de géomatique appliquée** de l'Université de Sherbrooke. Il y enseigne aux **programmes de 1er et 2e cycles de géomatique** les cours *Traitement numérique des images de télédétection*, *Base de données géospatiales* et *Apprentissage profond appliqué à l'observation de la Terre*. Ses intérêts de recherche portent sur le traitement d'images et l'application de l'IA aux données géospatiales.

**Philippe Apparicio** est professeur titulaire au **Département de géomatique appliquée** de l'Université de Sherbrooke. Il y enseigne aux **programmes de 1er et 2e cycles de géomatique** les cours *Transport et mobilité durable*, *Modélisation et analyse spatiale* et *Géomatique appliquée à la gestion urbaine*. Durant les dernières années, il a offert plusieurs formations aux Écoles d'été du Centre interuniversitaire québécois de statistiques sociales (CIQSS). Géographe de formation, ses intérêts de recherche incluent la justice et l'équité environnementale, la mobilité durable, les pollutions atmosphérique et sonore, et le vélo en ville. Il a publié une centaine d'articles scientifiques dans différents domaines des études urbaines et de la géographie mobilisant la géomatique et l'analyse spatiale.

**Mickaël Germain** est professeur agrégé au **Département de géomatique appliquée** de l'Université de Sherbrooke. Ses expertises comprennent les systèmes d'information géographique (SIG) et leur utilisation sur Internet (SIG web), les systèmes décisionnels, les bases de données à référence spatiale et le forage de données massives. Il s'intéresse également à la cartographie web, à la mobilité, à l'Internet des objets (IoT), à l'apprentissage machine et à la fusion des données multisources pour l'aide à la décision. Ses systèmes sont utilisés dans divers domaines de la géomatique tels que les ressources naturelles, la géologie, les transports et les écosystèmes. De plus, il se passionne pour le traitement, la gestion et l'analyse des données de télédétection pour extraire des paramètres descriptifs facilitant l'identification, la segmentation ou la classification de l'information.

**Yacine Bouroubi** est professeur agrégé au **Département de géomatique appliquée** de l'Université de Sherbrooke. Il possède plus de vingt ans d'expérience dans le domaine du traitement des images numériques appliqué à diverses thématiques des sciences6. Ses travaux de recherche se concentrent principalement sur l'utilisation de techniques de traitement d'images pour résoudre des problèmes géographiques. En tant que professeur, il partage son expertise avec les étudiants et contribue à l'avancement des connaissances dans son domaine. Bien que les détails spécifiques de sa carrière ne soient pas fournis, son expérience étendue et son poste actuel témoignent de sa contribution significative à la géographie et au traitement d'images numériques.

**Étienne Clabaut** est professeur associé au **Département de géomatique appliquée** de l'Université de Sherbrooke. Étienne détient un doctorat en télédétection de l'Université de Sherbrooke. Il se spécialise dans l'utilisation de l'apprentissage automatique pour l'exploitation de l'imagerie satellite.
# Partie 1. Importation, manipulation et visualisation de données satellitaires
# 1 Introduction au langage Python

Dans ce chapitre, nous présentons quelques éléments essentiels du langage Python qui nous seront utiles dans ce manuel. Python est un langage très riche et peut aboutir à des projets logiciels très sophistiqués. Il est important de comprendre que la programmation Python n'est pas ici une fin en soi, mais plutôt un outil de scriptage et de manipulation des données satellitaires.

<figure><img src="image_15.png" /></figure>

## Objectif

### Objectifs d'apprentissage visés dans ce chapitre

À la fin de ce chapitre, vous devriez être en mesure de :

- connaître les principales distributions de Python;

- installer un environnement d'exécution du code de cet ouvrage;

- comprendre les structures de base du langage Python;

Python, créé par Guido van Rossum en 1991, est un langage de programmation polyvalent et facile à apprendre, souvent comparé à un couteau suisse numérique pour sa simplicité et sa polyvalence. Comme un outil multifonction, Python peut être utilisé pour une variété de tâches, du développement web à l'analyse de données, en passant par l'intelligence artificielle.

## 1.1 Les distributions

Il existe plusieurs distributions du langage Python, ces distributions sont comme différentes saveurs de votre glace préférée - chacune a ses propres caractéristiques uniques, mais elles sont toutes fondamentalement Python. Voici un aperçu des principales distributions :

- CPython est la distribution "vanille" officielle, comme la recette originale de Python. Elle est ainsi le choix idéal pour la compatibilité et la conformité aux standards.

- Anaconda. Pensez-y comme à un sundae tout garni. Elle vient avec de nombreuses bibliothèques scientifiques préinstallées, ce qui est idéal pour l'analyse de données et l'apprentissage automatique (machine learning).

- Miniconda est une distribution légère de Python qui vous permet d'ajouter au besoin d'autres bibliothèques.

- PyPy : est une version turbo de Python, optimisée pour la vitesse.

Chaque distribution a ses forces, que ce soit la simplicité, la vitesse ou des fonctionnalités spécifiques. Le choix dépend donc de vos besoins, comme choisir entre une glace simple ou une glace royal (banana split) élaboré.
## 1.2 Les styles de programmation en Python

Il existe plusieurs approches pour programmer en Python. La plus directe est en version interactive en tapant `python` et de rentrer des commandes ligne par ligne.

### 1.2.1 Les outils de programmation

Un code python prend la forme d'un simple fichier texte avec l'extension `.py` et peut être modifié avec un simple éditeur de texte. Cependant, il n'y aura pas de rétroactions immédiates de l'interpréteur Python, ce qui rend la correction d'erreurs (débogage) beaucoup plus laborieux.

Un IDE (*Integrated Developement Environnement*) est comme une boîte à outils complète pour les programmeurs, vous trouverez :

- Un éditeur de texte amélioré pour écrire votre code, avec des fonctionnalités comme la coloration syntaxique qui rend le code plus lisible.

- Un compilateur qui transforme votre code en instructions que l'ordinateur peut comprendre.

- Un débogueur pour trouver et corriger les erreurs, tel un détective numérique.

- Des outils d'automatisation qui effectuent des tâches répétitives, comme un assistant virtuel pour le codage.

- L'accès à la documentation des différentes librairies.

Ces outils intégrés permettent aux développeurs de travailler plus efficacement, en passant moins de temps à jongler entre différentes applications et plus de temps à produire du code.

Voici quelques options populaires :

- **PyCharm**: est un des outils les plus utilisés dans l'industrie. Il offre une multitude de fonctionnalités comme l'autocomplétion intelligente et le débogage intégré, idéal pour les grands projets. Cependant, cet outil peut être assez gourmand en mémoire et en CPU. Notez qu'il existe une version gratuite et une version professionnelle.

- **Visual Studio Code**: gratuit, léger mais puissant, il est personnalisable avec des extensions pour Python.

- **Spyder**: logiciel libre et gratuit, orienté vers les applications scientifiques.

- **Jupyter Notebooks**: imaginez un cahier interactif pour le code. Idéal pour l'analyse de données et l'apprentissage, il permet de mélanger code, texte et visualisations. Des services gratuits dans le **cloud** sont disponibles comme Google Colab et Kaggle. Ces environnements sont néanmoins moins appropriés pour des grands projets et le débogage. Jupyter peut souffrir de problèmes de reproductibilité dus à l'exécution arbitraire des cellules.

- **Marimo**: se veut une version plus moderne des notebooks jupyter. Contrairement à Jupyter, Marimo garantit la cohérence entre le code, les sorties et l'état du programme. Il analyse intelligemment les relations entre les cellules et réexécute automatiquement celles qui sont affectées par des changements, éliminant ainsi les problèmes d'état caché.
## 1.3 Bonnes pratiques

Python est un langage très dynamique, qui évolue constamment. Cela pose certains défis pour la gestion du code à long terme. Il est fortement conseillé d'utiliser des environnements virtuels pour gérer vos différentes bibliothèques (*libraries*). Voici quelques bonnes pratiques à suivre :

1. **N'installez par la toute dernière version de Python :** Il est recommandé d'installer 1 ou 2 version antérieure, par exemple si 3.13 est la version plus récente, installer plutôt la version 3.11. Les versions trop récentes peuvent être instables. La version de python désirée peut être spécifiée au moment de la création d'un environnement virtuel (voir plus bas). Vous pouvez afficher la liste des versions de python avec la commande `conda search --full-name python`.

2. **N'utilisez pas de version obsolète de Python.** Cela peut sembler contradictoire avec le point précédent mais c'est l'excès inverse. Si vous utilisez une version trop ancienne alors toutes vos libraries cesseront d'évoluer et peuvent devenir obsolètes.

3. **Utilisez des environnements virtuels.** Pensez-y comme à des compartiments séparées pour chaque projet. Cela évite les conflits entre les différentes versions de bibliothèques (*libraries*)et garde votre système propre. Par exemple, si vous souhaitez vérifier une nouvelle version de Python, utilisez un environnement : `conda create --name test python=3.11`

4. **Vérifiez l'installation.** Après l'installation, ouvrez un terminal et tapez `python --version` pour vous assurer que tout fonctionne correctement.

### 1.3.1 Création d'un environnement virtuel

Il y a deux façons d'installer un environnement virtuel selon votre distribution de Python:

1. **Option 1.** Vous utilisez Anaconda ou Miniconda. La commande `conda` est utilisée pour créer un environnement test avec Python 3.10:

```
1 conda env -n test python=3.10
```

```
2 conda activate test
```

2. **Option 2.** Vous utilisez CPython

```
1 conda env -n test python=3.10
```

```
2 conda activate test
```

### 1.3.2 Création d'un environnement de travail local (avancé)

**Note:** les notebooks peuvent fonctionner localement uniquement sous Linux ou avec WSL2.

Les notebooks Python fonctionnent par défaut dans l'environnement Google Colab. Si vous souhaitez faire fonctionner ces notebook localement, vous pouvez installer un environnement local avec un serveur Jupyter. Il suffit de suivre les étapes suivantes:

1. Installer WSL2 sous Windows
1 Introduction au langage Python

2. Installer vscode

3. Installer Miniconda

4. Faire une installation du contenu du livre soit en utilisant une commande `git clone` ou en récupérant le `.zip` du livre

5. Ouvrir WSL2 et placer vous dans le répertoire du livre `TraitementImagesPythonVol1`. Assurez vous que vous avez accès à conda en tapant `conda --version`

6. Lancer la commande `conda env create -f jupyter_env.yaml`

7. Activer le nouvel environnement: `conda activate jupyter_env`

8. Le serveur jupyter peut ensuite être lancé avec la commande suivante: `jupyter lab --ip='*' -- NotebookApp. token=''' -- NotebookApp. password='''`

Une fenêtre devrait alors apparaître dans votre fureteur. Dans le menu de gauche vous pouvez accéder aux notebooks dans le répertoire `notebooks`:

<figure><img src="image_16.png" /><figcaption>FIGURE 1.1 – Fenêtre principale du serveur Jupyter Lab.</figcaption></figure>
## 1.4 Les structures de base en Python

Il y a essentiellement deux structures de données que Python manipule : les listes et les dictionnaires.

### 1.4.1 Les listes

Les listes sont comme des boîtes extensibles où vous pouvez ranger différents types d'objets:

- Représentées par des crochets : `[1, 2, 3, "python"]`.

- Ordonnées et modifiables (*mutable*), vous pouvez récupérer une valeur par sa position avec `[]`.

- Permettent les doublons (deux fois la même valeur).

- Idéales pour stocker des collections d'éléments que vous voulez modifier

### 1.4.2 Les tuples

Les tuples sont similaires aux listes, mais les boîtes sont scellées :

- Représentés par des parenthèses : `(1, 2, 3, "python")`.

- Ordonnés mais non modifiables (*immutable*).

- Permettent les doublons.

- Souvent utilisés pour stocker des données qui ne doivent pas changer (comme des paramètres).

### 1.4.3 Les ensembles (Sets)

Les ensembles sont comme des boîtes magiques qui ne gardent qu'un exemplaire de chaque objet :

- Représentés par des accolades : `{1, 2, 3}`.

- Non ordonnés et modifiables.

- N'autorisent pas les doublons.

- Utiles pour éliminer les doublons et effectuer des opérations mathématiques sur des ensembles.

## 1.5 Dictionnaires

Les dictionnaires sont comme des boîtes avec des étiquettes sur chacune d'elle :

- Représentés par des accolades avec des paires clé-valeur : `{"nom": "Python", "année": 1991}`.

- Non ordonnés et modifiables.

- Les clés doivent être uniques, mais les valeurs peuvent être dupliquées

- Utiles pour stocker des données associatives ou pour créer des tables de recherche rapide
## 1.6 Programmation objet

La programmation orientée objet (POO) en Python est comme construire avec des blocs LEGO. Chaque objet est un bloc LEGO avec ses propres caractéristiques (attributs) et capacités (méthodes). Les classes sont les plans pour créer ces blocs. Par exemple, une classe "Voiture" pourrait avoir des attributs comme "couleur" et "vitesse", et des méthodes comme "démarrer" et "accélérer".

Python rend la POO accessible avec des fonctionnalités conviviales:

1. **Encapsulation**: comme emballer un cadeau, elle cache les détails internes d'un objet.

2. **Héritage**: permet de créer de nouvelles classes basées sur des classes existantes, comme un enfant héritant des traits de ses parents.

3. **Polymorphisme**: permet à différents objets de répondre au même message de manière unique, comme si différents animaux répondaient différemment à "fais du bruit".

Ces caractéristiques font de Python un excellent choix pour apprendre et appliquer les concepts de la POO, rendant le code plus organisé et réutilisable

<figure><img src="image_17.png" /></figure>

### Package

#### Liste des packages utilisés dans ce chapitre

- Pour importer et manipuler des fichiers géographiques :

- ` numpy ` pour manipuler des données matricielles.

- ` rasterio ` pour importer et manipuler des données matricielles.

- Pour construire des cartes et des graphiques :

- ` matplotlib ` est certainement le *package* le plus complet pour l'affichage général.

- ` seaborn ` pour construire des graphiques plus détaillés en particulier pour les statistiques.
# 2 Importation et manipulation de données spatiales

## 2.1 Préambule

Assurez-vous de lire ce préambule avant d'exécutez le reste du notebook.

### 2.1.1 Objectifs

Dans ce chapitre, nous abordons quelques formats d'images ainsi que leur lecture. Ce chapitre est aussi disponible sous la forme d'un notebook Python:

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /></figure>

#### Objectif

**Objectifs d'apprentissage visés dans ce chapitre**

À la fin de ce chapitre, vous devriez être en mesure de :

- connaître les principales bibliothèques Python pour lire une image;

- accéder à l'information d'une image avant de la lire;

- comprendre les principaux formats pour une image

- manipuler la matrice de la donnée d'une image avec numpy

### 2.1.2 Bibliothèques

Les bibliothèques qui vont être explorées dans ce chapitre sont les suivantes:

- SciPy

- NumPy

- opencv-python · PyPI

- scikit-image

- Rasterio

- Xarray

- rioxarray

Dans l'environnement Google Colab, seul `rioxarray` et `gdal` doivent être installés:

1. `!apt-get update`

2. `!apt-get install gdal-bin libgdal-dev`

3. `!pip install -q rioxarray`
Vérifier les importations:

1 import numpy as np

2 import rioxarray as rxr

3 from scipy import signal

4 import xarray as xr

5 import xrscipy

6 import matplotlib.pyplot as plt

### 2.1.3 Données

Nous utilisons ces images dans ce chapitre:

1 import gdown

2 gdown/download()

→ 'https://drive.google.com/uc?export=download&confirm=pbef&id=1a6Ypg0g10y4AJt9XWKWfnR12NW1XhNg_',
→ output= 'RGBNIR_of_S2A.tif')

3 gdown/download()

→ 'https://drive.google.com/uc?export=download&confirm=pbef&id=1a4PQ68Ru8zBphbQ22j0sgJ4D2quw-Wo6',
→ output= 'landsat7.tif')

4 gdown/download()

→ 'https://drive.google.com/uc?export=download&confirm=pbef&id=1_zwCLN-x7XJcNHJCH6Z8upEdUxtVtvs1',
→ output= 'berkeley.jpg')

5 !wget https://raw.github.com/sfoucher/TraitementImagesPythonVol1/refs/heads/main/images/ gdown/download()

Vérifiez que vous êtes capable de les lire:

1 with rxr.open_rasterio('berkeley.jpg', mask_and_scale= True ) as img rgb:

2 print(img rgb)

3 with rxr.open_rasterio('RGBNIR_of_S2A.tif', mask_and_scale= True ) as img rgnir:

4 print(img rgnir)

## 2.2 Importation d'images

La première étape avant tout traitement est d'accéder à la donnée image pour qu'elle soit manipulée par le langage Python. L'imagerie satellite présente certains défis notamment en raison de la taille parfois très importante des images. Il existe maintenant certaines bibliothèques, comme `Xarray`, qui visent à optimiser la lecture et l'écriture de grandes images. Il est donc conseillé de toujours garder un oeil sur l'espace mémoire occupé par les variables Python représentant les images. La librairie principale en géomatique qui permettre d'importer (et d'exporter) de l'imagerie est la librairie `GDAL` qui rassemble la plupart des formats sous forme de `driver` (ou pilote en français).
Dans le domaine de la géomatique, il faut prêter attention à trois caractéristiques principales des images:

1. La **matrice des données** elle-même qui contient les valeurs brutes des pixels. Cette matrice sera souvent un cube à trois dimensions. En Python, ce cube sera le plus souvent un objet de la librairie `NumPy` (voir section).

2. La **dynamique des images** c.-à-d le format de stockage des valeurs individuelles (octet, entier, double, etc.). Ce format décide principalement de la résolution radiométrique et des valeurs minimales et maximales supportées.

3. Le **nombre de bandes** spectrales de l'image qui est souvent supérieur à trois et peut atteindre plusieurs centaines de bandes pour certains capteurs (notamment hyperspectraux).

4. La **métadonnée** qui va transporter l'information auxiliaire de l'image comme les dimensions et la position de l'image, la date, etc. Cette donnée auxiliaire prendra souvent la forme d'un dictionnaire Python. Elle contiendra aussi l'information de géoréférence.

Les différents formats se distinguent principalement sur la manière dont ces trois caractéristiques sont gérées.

### 2.2.1 Formats des images

Il existe de nombreux formats numériques pour la donnée de type image parfois appelé donnée matricielle ou donnée raster. La librairie GDAL rassemble la plupart des formats matriciels rencontrés en géomatique (voir Raster drivers — GDAL documentation pour une liste complète).

On peut distinguer deux grandes familles de format:

1. Les formats de type **RVB** issus de l'imagerie numérique grand public comme **JPEG**, **png**, etc. Ces formats ne supportent généralement que trois bandes au maximum (rouge, vert et bleu) et des valeurs de niveaux de gris entre 0 et 255 (format dit 8 bits ou `uint8`).

2. Les **géo-formats** issus des domaines scientifiques ou techniques comme GeoTIFF, HDF5, NetCDF, etc. qui peuvent inclure plus que trois bandes et des dynamiques plus élevées (16 bits ou même float).

Les formats RVB restent très utilisés en Python notamment par les bibliothèques dites de vision par ordinateur (*Computer Vision*) comme OpenCV et sickit-image ainsi que les grandes bibliothèques en apprentissage profond (PyTorch, Tensorflow).

<figure><img src="image_20.png" /></figure>

**Package**

**Installation de gdal dans un système Linux**

— Pour installer GDAL :

```
!apt-get update
!apt-get install gdal-bin libgdal-dev
```
2.2.1.1 Formats de type RVB

Les premiers formats pour de l'imagerie à une bande (monochrome) et à trois bandes (image couleur rouge-vert-
bleu) sont issus du domaine des sciences de l'ordinateur. On trouvera, entre autres, les formats pbm, png et jpeg.
Ces formats supportent peu de métadonnées et sont placées dans un entête (header) très limité. Cependant, ils
restent très populaires dans le domaine de la vision par ordinateur et sont très utilisés en apprentissage profond
en particulier. Pour la lecture des images RVB, on peut utiliser les bibliothèques Rasterio, PIL ou OpenCV.

Lecture avec la librairie PIL

La librairie PIL retourne un objet de type `PngImageFile`, l'affichage de l'image se fait directement dans la cellule de sortie.

1 from PIL import Image

2 img = Image.open('modis-aqua.PNG')

3 img
<figure><img src="image_21.png" /></figure>

### Lecture avec la librairie OpenCV

La librairie OpenCV est aussi très populaire en vision par ordinateur. La fonction `imread` donne directement un objet de type `NumPy` en sortie.

1. `import cv2`

2. `img = cv2.imread('modis-aqua.PNG')`

3. `img`

```python
array([[[17, 50, 33],
     [15, 49, 31],
     [14, 48, 30],
     ...,
     [23, 56, 36],
     [23, 55, 36],
     [22, 55, 36]]],
    [[18, 51, 34],
     [16, 50, 32]]
```
```
[15, 49, 32],
[27, 59, 40],
[28, 60, 41],
[27, 60, 41]]

[[18, 53, 35],
[18, 52, 34],
[18, 51, 34],
[31, 64, 44],
[34, 66, 47],
[33, 65, 46]]



[[34, 74, 48],
[35, 73, 48],
[34, 70, 46],
[41, 74, 54],
[41, 73, 54],
[41, 73, 54]]

[[36, 76, 50],
[36, 74, 49],
[35, 71, 47],
[37, 70, 51],
[38, 71, 51],
[38, 71, 51]]

[[36, 76, 50],
[35, 73, 48],
[33, 69, 45],
[31, 63, 44],
[33, 65, 46],
[33, 66, 46]]], dtype=uint8)
```

### Lecture avec la librairie RasterIO

Rien ne nous empêche de lire une image de format RVB avec `RasterIO` comme décrit dans ci-dessous. Vous noterez cependant les avertissements concernant l'absence de géoréférence pour ce type d'image.

1. `import rasterio`

2. `img = rasterio.open('modis-aqua.PNG')`

3. `img`

```python
<open DatasetReader name='modis-aqua.PNG' mode='r'>
```

### 2.2.1.2 Le format GeoTiff

Le format GeoTIFF est une extension du format TIFF (Tagged Image File Format) qui permet d'incorporer des métadonnées géospatiales directement dans un fichier image. Développé initialement par Dr. Niles Ritter au Jet Propulsion Laboratory de la NASA dans les années 1990, GeoTIFF est devenu un standard de facto pour le stockage et l'échange d'images géoréférencées dans les domaines de la télédétection et des systèmes d'information géographique (SIG). Ce format supporte plus que trois bandes aussi longtemps que ces bandes sont de même dimension.

Le format GeoTIFF est très utilisé et est largement supporté par les bibliothèques et logiciels géospatiaux, notamment GDAL (Geospatial Data Abstraction Library), qui offre des capacités de lecture et d'écriture pour ce format. Cette compatibilité étendue a contribué à son adoption généralisée dans la communauté géospatiale.
Standardisation par l'OGC

Le standard GeoTIFF proposé par l'Open Geospatial Consortium (OGC) en 2019 formalise et étend les spécifications originales du format GeoTIFF, offrant une norme robuste pour l'échange d'images géoréférencées. Cette standardisation, connue sous le nom d'OGC GeoTIFF 1.1 (2019), apporte plusieurs améliorations et clarifications importantes.

2.2.1.3 Le format COG

Une innovation récente dans l'écosystème GeoTIFF est le format *Cloud Optimized GeoTIFF (COG)*, conçu pour faciliter l'utilisation de fichiers GeoTIFF hébergés sur des serveurs web HTTP. Le COG permet aux utilisateurs et aux logiciels d'accéder à des parties spécifiques du fichier sans avoir à le télécharger entièrement, ce qui est particulièrement utile pour les applications basées sur l'infonuagique.

2.2.2 Métadonnées des images

La manière la plus directe d'accéder à la métadonnée d'une image est d'utiliser les commandes rio info de la
librairie Rasterio ou gdalinfo de la librairie gdal. Le résultat est imprimé dans la sortie standard ou sous forme
d'un dictionnaire Python.

1

!gdalinfo RGBNIR_of_S2A.tif

Warning 1: TIFFReadDirectory:Sum of Photometric type-related color channels and ExtraSamples doesn't match SamplesPerPixel. Defining non-color
-> channels as ExtraSamples.
Driver: GTiff/GeoTIFF
Files: RGBNIR_of_S2A.tif
RGBNIR_of_S2A.tif.aux.xml
Size is 2074, 1926
Coordinate System is:
PROJCS["WGS 84 / UTM zone 18N",
    GEOGCS["WGS 84",
        DATUM["WGS 1984",
            SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
            AUTHORITY["EPSG","4326"]],
        PROJECTION["Transverse_Mercator"],
        PARAMETER["latitude_of_origin",0],
        PARAMETER["central_meridian",-75],
        PARAMETER["scale_factor",0.9996],
        PARAMETER["false_easting",50000],
        PARAMETER["false_northing",0],
        UNIT["metre",1,
            AUTHORITY["EPSG","9001"]],
        AXIS["Easting",EAST],
        AXIS["Northing",NORTH],
        AUTHORITY["EPSG","32618"]]
    Origin = (731780.000000000000000,5040800.00000000000000)
    Pixel Size = (10.000000000000000,-10.000000000000000)
Metadata:
    AREA_OR_POINT=Area
    TIFFTAG_IMAGEDESCRIPTION=subset_RGBNIR_of_S2A_MSIL2A_20240625T153941_N0510_R011_T18TYR_20240625T221903
    TIFFTAG RESOLUTIONUNIT=1 (unitless)
    TIFFTAG_XRESOLUTION=1
    TIFFTAG_YRESOLUTION=1
Image Structure Metadata:
    INTERLEAVE=BAND
Corner Coordinates:
    Upper Left ( 731780.000, 5040800.000) ( 72d 2' 3.11"W, 45d28'55.98"N)
    Lower Left ( 731780.000, 5021540.000) ( 72d 2'35.69"W, 45d18'32.70"N)
    Upper Right ( 752520.000, 5040800.000) ( 71d46' 9.19"W, 45d28'30.08"N)
    Lower Right ( 752520.000, 5021540.000) ( 71d46'44.67"W, 45d18' 6.95"N)
    Center ( 742150.000, 5031170.000) ( 71d54'23.16"W, 45d23'31.71"N)
Band 1 Block=2074x1926 Type=UINT16, ColorInterp=Gray
    Min=86.000 Max=15104.000
    Minimum=86.000, Maximum=15104.000, Mean=1426.625, StdDev=306.564
```
Metadata:
    STATISTICS_MAXIMUM=15104
    STATISTICS_MEAN=1426.6252674912
    STATISTICS_MINIMUM=86
    STATISTICS_STDDEV=306.56427126942
    STATISTICS_VALID_percent=100
Band 2 Block=2074x1926 Type=UINT16, ColorInterp=Undefined
    Min=1139.000 Max=14352.000
    Minimum=1139.000, Maximum=14352.000, Mean=1669.605, StdDev=310.919
Metadata:
    STATISTICS_MAXIMUM=14352
    STATISTICS_MEAN=1669.6050060032
    STATISTICS_MINIMUM=1139
    STATISTICS_STDDEV=310.91935787639
    STATISTICS_VALID_percent=100
Band 3 Block=2074x1926 Type=UINT16, ColorInterp=Undefined
    Min=706.000 Max=15280.000
    Minimum=706.000, Maximum=15280.000, Mean=1471.392, StdDev=385.447
Metadata:
    STATISTICS_MAXIMUM=15280
    STATISTICS_MEAN=1471.3923473736
    STATISTICS_MINIMUM=706
    STATISTICS_STDDEV=385.44654593014
    STATISTICS_VALID_percent=100
Band 4 Block=2074x1926 Type=UINT16, ColorInterp=Undefined
    Min=1067.000 Max=15642.000
    Minimum=1067.000, Maximum=15642.000, Mean=4393.945, StdDev=1037.934
Metadata:
    STATISTICS_MAXIMUM=15642
    STATISTICS_MEAN=4393.94485025
    STATISTICS_MINIMUM=1067
    STATISTICS_STDDEV=1037.933939728
    STATISTICS_VALID_percent=100
```

Le plus simple est d'utiliser la fonction rio info :

```
1 !rio info RGBNIR_of_S2A.tif --indent 2 --verbose
```

## 2.3 Manipulation des images

### 2.3.1 Manipulation de la matrice de pixels

La donnée brute de l'image est généralement contenue dans un cube matricielle à trois dimensions (deux dimensions spatiales et une dimension spectrale). Comme exposé précédemment, la librairie dite "fondationnelle" pour la manipulation de matrices en Python est NumPy. Cette librairie contient un nombre très important de fonctionnalités couvrant l'algèbre linéaire, les statistiques, etc.; elle constitue la fondation de nombreuses bibliothèques en traitement numérique (voir (figure 2.1))

### 2.3.2 Information de base

Les deux informations de base à afficher sur une matrice sont 1) les dimensions de la matrice et 2) le format de stockage (le type). Pour cela, on peut utiliser le code ci-dessous, dont le résultat nous informe que la matrice a trois dimensions et une taille de (442, 553, 3) et un type `uint8` qui représente 1 octet (8 bit). Par conséquent, la matrice a 442 lignes, 553 colonnes et 3 canaux ou bandes. Il faut prêter une attention particulière aux valeurs minimales et maximales tolérées par le type de la donnée comme indiqué dans le (tableau 2.1) (voir aussi Data types — NumPy v2.1 Manual).

```
1 import cv2
```

```
2 img = cv2.imread('modis-aqua.PNG')
```

```
3 print('Nombre de dimensions: ', img Dim)
```
<figure><img src="image_22.png" /><figcaption>FIGURE 2.1 – La librairie NumPy est le fondement de nombreuses bibliothèques scientifiques (d’après (Harris 2020)).</figcaption></figure>
```python
print('Dimensions de la matrice: ', img.shape)
```

```python
print('Type de la donnée: ', img.dtype)
```

```
Nombre de dimensions: 3
Dimensions de la matrice: (442, 553, 3)
Type de la donnée: uint8
```

<table><caption>TABLEAU 2.1 – Type de données de NumPy</caption><thead><tr><th>datatype</th><th>Nom</th><th>Taille (bits)</th><th>Min</th><th>Max</th></tr></thead><tbody><tr><td>uint8</td><td>char</td><td>8</td><td>0</td><td>255</td></tr><tr><td>int8</td><td>signed char</td><td>8</td><td>-127</td><td>128</td></tr><tr><td>uint16</td><td>unsigned short</td><td>16</td><td>0</td><td>-32768</td></tr><tr><td>int16</td><td>short</td><td>16</td><td>0</td><td>655355</td></tr></tbody></table>

### 2.3.3 Découpage et indexation de la matrice

L'indexation et le découpage (*slicing*) des matrices dans NumPy sont des techniques essentielles pour manipuler efficacement les données multidimensionnelles en Python, offrant une syntaxe puissante et flexible pour accéder et modifier des sous-ensembles spécifiques d'éléments dans les tableaux (voir figure 2.2). Indexer une matrice consiste à accéder à une valeur dans la matrice pour une position particulière, la syntaxe générale est **matrice[ligne, colonne, bande]** et est similaire à la manipulation des listes en Python. Les indices commencent à 0 et se termine à la **taille-1** de l'axe considéré.

**a Data structure**

$$x = \begin{bmatrix} 0 & 1 & 2 \\ 3 & 4 & 5 \\ 6 & 7 & 8 \\ 9 & 10 & 11 \end{bmatrix}$$

<figure><img src="image_23.png" /></figure>

**b Indexing (view)**

$$x[:,1:] \rightarrow \begin{bmatrix} 0 & 1 & 2 \\ 3 & 4 & 5 \\ 6 & 7 & 8 \\ 9 & 10 & 11 \end{bmatrix}$$

with slices

Slices are start:end:step,
any of which can be left blank

$$xOODS$$

$$\begin{bmatrix} 0 & 1 & 2 \\ 3 & 4 & 5 \\ 6 & 7 & 8 \\ 9 & 10 & 11 \end{bmatrix}$$

with slices
with steps

**c Indexing (copy)**

$$x[1,2] \rightarrow 5 \text{ with scalars}$$

$$x[x > 9] \rightarrow \begin{bmatrix} 10 & 11 \end{bmatrix} \text{ with masks}$$

$$x[\begin{bmatrix} 0 & 1 \\ 1 & 2 \end{bmatrix}] \rightarrow [x[0,1], x[1,2]] \rightarrow \begin{bmatrix} 1 & 5 \end{bmatrix} \text{ with arrays}$$

$$x[\begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix}] \rightarrow x[\begin{bmatrix} 1 & 1 \\ 1 & 2 \end{bmatrix}, \begin{bmatrix} 1 & 0 \\ 1 & 1 \end{bmatrix}] \rightarrow \begin{bmatrix} 4 & 3 \\ 7 & 6 \end{bmatrix} \text{ with arrays with broadcasting}$$

**d Vectorization**

$$\begin{bmatrix} 0 & 1 & 2 \\ 3 & 4 & 5 \\ 6 & 7 & 8 \\ 9 & 10 & 11 \end{bmatrix} + \begin{bmatrix} 1 & 1 & 2 \\ 1 & 1 & 3 \\ 1 & 1 & 4 \\ 1 & 1 & 5 \end{bmatrix} \rightarrow \begin{bmatrix} 1 & 2 \\ 4 & 5 \\ 7 & 8 \\ 10 & 11 \end{bmatrix}$$

**e Broadcasting**

$$\begin{bmatrix} 0 \\ 3 \\ 6 \\ 9 \end{bmatrix} \times \begin{bmatrix} 1 & 2 \\ 1 & 3 \\ 1 & 4 \\ 1 & 5 \end{bmatrix} \rightarrow \begin{bmatrix} 0 & 0 \\ 3 & 6 \\ 6 & 12 \\ 9 & 18 \end{bmatrix}$$

**f Reduction**

$$\begin{bmatrix} 0 & 1 & 2 \\ 3 & 4 & 5 \\ 6 & 7 & 8 \\ 9 & 10 & 11 \end{bmatrix} \xrightarrow{\text{sum axis 1}} \begin{bmatrix} 3 \\ 12 \\ 21 \\ 30 \end{bmatrix}$$

**g Example**

In [1]: import numpy as np

In [2]: x = np.arange(12)

In [3]: x = x.reshape(4, 3)

In [4]: x

Out [4]:

array([0, 1, 2],
[3, 4, 5],
[6, 7, 8],
[9, 10, 11]])

In [5]: np.mean(x, axis=0)

Out [5]: array([4.5, 5.5, 6.5])

In [6]: x = x - np.mean(x, axis=0)

In [7]: x

Out [7]:

array([-4.5, -4.5, -4.5],
[-1.5, -1.5, -1.5],
[ 1.5,  1.5,  1.5],
[ 4.5,  4.5,  4.5]])

<caption>FIGURE 2.2 – Vue d'ensemble des opérations de base des matrices avec NumPy</caption>

Le découpage (ou *slicing* en anglais) consiste à produire une nouvelle matrice qui est un sous-ensemble de la matrice d'origine. Un découpage se fait avec le symbole ':', la syntaxe générale pour définir un découpage
est `[début:fin:pas]`. Si on ne spécifie pas début ou fin alors les valeurs 0 ou `dimension-1` sont considérées implicitement. Quelques exemples: * choisir un pixel en particulier avec toutes les bandes: `matrice[1,1,:] *` choisir la colonne 2: `matrice[:,2:]`

La syntaxe de base pour le découpage (*slicing*) des tableaux NumPy repose sur l'utilisation des deux-points (:) à l'intérieur des crochets d'indexation. Cette notation permet de sélectionner des plages d'éléments de manière concise et intuitive. La structure générale du découpage est `matrice[start:stop:step]`, où : 1. `start` représente l'index de départ (inclus) 2. `stop` indique l'index de fin (exclu) 3. `step` définit l'intervalle entre chaque élément sélectionné

Si l'un de ces paramètres est omis, NumPy utilise des valeurs par défaut : 0 pour `start`, la taille du tableau pour `stop`, et 1 pour `step`. Par exemple, pour un tableau unidimensionnel `array`, on peut extraire les éléments du deuxième au quatrième avec `array[1:4]`. Pour sélectionner tous les éléments à partir du troisième, on utiliserait `array[2:]`. Cette syntaxe s'applique également aux tableaux multidimensionnels, où chaque dimension est séparée par une virgule. Ainsi, pour une matrice 2D m, `m[0:2, 1:3]` sélectionnerait une sous-matrice 2x2 composée des deux premières lignes et des deuxième et troisième colonnes. L'indexation négative est également supportée, permettant de compter à partir de la fin du tableau. Par exemple, `a[-3:]` sélectionnerait les trois derniers éléments d'un tableau.

```python
import cv2
img = cv2.imread('modis-aqua.PNG')
img_col = img[:1,:]
print('Nombre de dimensions: ', img_col Dim)
print('Dimensions de la matrice: ', img_col.shape)
```

```
Nombre de dimensions: 2
Dimensions de la matrice: (442, 3)
```

<figure><img src="image_24.png" /></figure>

### Aller plus loin

#### Une vue versus une copie

Avec NumPy, les manipulations peuvent créer des vues ou des copies. Une vue est une simple représentation de la même donnée originale alors qu'une copie est un nouvel espace mémoire.

Par défaut, un découpage créé une vue.

On peut vérifier si l'espace mémoire est partagé avec `np.shares_memory(arr, slice_arr)`.

On peut toujours forcer une copie avec la méthode `copy()`.

### 2.3.3.1 Masquage

L'utilisation d'un masque est un outil important en traitement d'image car la plupart des images de télédétection contiennent des pixels non valides qu'il faut exclure des traitements (ce que l'on appelle le *no data* en Anglais). Il y a plusieurs raisons possibles pour la présence de pixels non valides:

1. L'image est projetée dans une grille cartographique et certaines zones, généralement situées en dehors de l'empreinte au sol du capteur, sont à exclure.

2. La présence de nuages que l'on veut exclure.

3. La présence de pixels erronés dus à des problèmes de capteurs.

4. La présence de valeurs non numériques (*not a number* ou *nan*)
La librairie NumPy fournit des mécanismes pour exclure automatiquement certaines valeurs.

### 2.3.4 Changement de projection cartographique (à venir)

## 2.4 Données en géoscience

Les données en géoscience contiennent beaucoup de métadonnées et peuvent être composées de différentes variables avec différentes unités, résolution, etc. Ces données sont aussi souvent étiquetées avec des dates sur certains axes, des coordonnées géographiques, des identifiants d'expériences, etc. Par conséquent, utiliser seulement des matrices est souvent incomplet (Hoyer et Hamman 2017).

Calibration, unités, données manquantes, données éparses.

### 2.4.1 xarray

Xarray est une puissante bibliothèque Python qui améliore les matrices multidimensionnelles de type numpy en y ajoutant des étiquettes, des dimensions, des coordonnées et des attributs. Elle fournit deux structures de données principales : `DataArray` (un tableau étiqueté à n dimensions) et `Dataset` (une base de données de tableaux multidimensionnels en mémoire).

Les caractéristiques principales sont les suivantes:

- Opérations sur les dimensions nommées au lieu des numéros d'axe

- Sélection et opérations basées sur les étiquettes

- Diffusion automatique de tableaux basée sur les noms de dimensions

- Alignement de type base de données avec des étiquettes de coordonnées

- Suivi des métadonnées grâce à des dictionnaires Python

### 2.4.1.1 Avantages

La bibliothèque réduit considérablement la complexité du code et améliore la lisibilité du code pour les applications de calcul scientifique dans divers domaines, notamment la physique, l'astronomie, les géosciences, la bio-informatique, l'ingénierie, la finance et l'apprentissage profond. Elle s'intègre de manière transparente avec NumPy et pandas tout en restant compatible avec l'écosystème Python au sens large.

### 2.4.1.2 DataArray

Un tableau multidimensionnel étiqueté avec des propriétés clées :

- valeurs : Les données réelles du tableau

- dims : Dimensions nommées (par exemple, « x », « y », « z »)

- coords : Dictionnaire de tableaux étiquetant chaque point

- attrs : Stockage de métadonnées arbitraires

- name : Identifiant facultatif
### 2.4.1.3 Dataset

Un conteneur de type dictionnaire de **DataArrays** avec des dimensions alignées, contenant :

- **dims**: Dictionnaire de correspondance entre les noms des dimensions et les longueurs

- **data_vars**: Dictionnaire des variables du DataArray

- **coords**: Dictionnaire des variables de coordonnées

- **attrs**: Stockage des métadonnées

Les principales différences sont les suivantes :

- **DataArray** contient un seul tableau avec des étiquettes

- Le **Dataset** contient plusieurs DataArrays alignés.

Ces trois structures prennent en charge les opérations de type dictionnaire et les calculs de coordination tout en conservant les métadonnées.

<figure><img src="image_25.png" /><figcaption>FIGURE 2.3 – Organisation d'un Dataset dans xarray</figcaption></figure>
# 3 Réhaussement et visualisation d'images

Assurez-vous de lire ce préambule avant d'exécutez le reste du notebook.

## 3.1 Préambule

### 3.1.1 Objectifs

Dans ce chapitre, nous abordons quelques techniques de réhaussement et de visualisation d'images. Ce chapitre est aussi disponible sous la forme d'un notebook Python:

<figure><img src="image_26.png" /></figure>

<figure><img src="image_27.png" /></figure>

**Objectif**

**Objectifs d'apprentissage visés dans ce chapitre**

À la fin de ce chapitre, vous devriez être en mesure de :

- exploiter les statistiques d'une image pour améliorer la visualisation;

- calculer les histogrammes de valeurs;

- appliquer une transformation linéaire ou non linéaire pour améliorer une visualisation;

- comprendre le principe des composés colorés;

### 3.1.2

### 3.1.3 Bibliothèques

Les bibliothèques qui vont être explorées dans ce chapitre sont les suivantes:

- SciPy

- NumPy

- opencv-python · PyPI

- scikit-image

- Rasterio

- Xarray

- rioxarray

Dans l'environnement Google Colab, seul rioxarray et GDAL doivent être installés:
1 **%capture**

2 !apt-get update

3 !apt-get install gdal-bin libgdal-dev

Dans l'environnement Google Colab, il convient de s'assurer que les librairies sont installées:

1 **%capture**

2 !pip install -qU matplotlib rioxarray xrscipy scikit-image

Vérifier les importations:

1 import numpy as np

2 import rioxarray as rxr

3 from scipy import signal

4 import xarray as xr

5 import xrscipy

6 import matplotlib.pyplot as plt

### 3.1.4 Données

Nous utiliserons les images suivantes dans ce chapitre:

1 **%capture**

2 import gdown

3

4 gdown.download() 
→ 'https://drive.google.com/uc?export=download&confirm=pbef&id=1a6Ypg0g1Oy4AJt9XWKWfnR12NW1XhNg_',
→ output= 'RGBNIR_of_S2A.tif')

5 gdown.download() 
→ 'https://drive.google.com/uc?export=download&confirm=pbef&id=1a603L_ab0fU7h94K22At8qtBuLMGErwo',
→ output= 'sentinel2.tif')

6 gdown.download() 
→ 'https://drive.google.com/uc?export=download&confirm=pbef&id=1_zwCLN-x7XJcNHJCH6Z8upEdUxtVtvs1',
→ output= 'berkeley.jpg')

7 gdown.download() 
→ 'https://drive.google.com/uc?export=download&confirm=pbef&id=1dM6IVqjba6GHwTLmI7CpX8GP2z5txUq6',
→ output= 'SAR.tif')

8 gdown.download() 
→ 'https://drive.google.com/uc?export=download&confirm=pbef&id=1a4PQ68Ru8zBphbQ22j0sgJ4D2quw-Wo6',
→ output= 'landsat7.tif')

Vérifiez que vous êtes capable de les lire :
1

with rxr.open_rasterio('berkeley.jpg', mask_and_scale= True) as img rgb:

2

print(img rgb)

3

with rxr.open_rasterio('sentinel2.tif', mask_and_scale= True) as img s2:

4

print(img s2)

5

with rxr.open_rasterio('RGBNIR_of_S2A.tif', mask_and_scale= True) as img rgbnir:

6

print(img rgbnir)

7

with rxr.open_rasterio('SAR.tif', mask_and_scale= True) as img SAR:

8

print(img SAR)

## 3.2 Visualisation en Python

ID'emblée, il faut mentionner que Python n'est pas vraiment fait pour visualiser de la donnée de grande taille, le niveau d'interactivité est aussi assez limité. Pour une visualisation interactives, il est plutôt conseillé d'utiliser un outil comme QGIS. Néanmoins, il est possible de visualiser de petites images avec la librairie ` matplotlib ` qui est la librairie principale de visualisation en Python. Cette librairie est extrêmement riche et versatile, nous ne présenterons que les bases nécessaires pour démarrer. Le lecteur désirant aller plus loin pourra consulter les nombreux tutoriels disponibles comme celui-ci.

La fonction de base pour créer une figure est `subplots`, la largeur et la hauteur en pouces de la figure peuvent être contrôlées via le paramètre `figsize`:

1

import matplotlib.pyplot as plt

2

fig, ax= plt.subplots(figsize=(5, 4))

3

plt.show()

<figure><img src="image_28.png" /></figure>
3 Réhaussement et visualisation d'images

Pour l'affichage des images, la fonction **imshow** permet d'afficher une matrice 2D à une dimension en format *float* ou une matrice RVB avec 3 bandes. Il est important que les dimensions de la matrice soient dans l'ordre hauteur, largeur et bande.

1
import matplotlib.pyplot as plt
2
fig, ax= plt.subplots(figsize=(6, 5))
3
plt.imshow(img_rgbnir[0].data)
4
plt.show()

<figure><img src="image_29.png" /></figure>

Pour un affichage à trois bandes, les valeurs seront ramenées sur une échelle de 0 à 1, il est donc nécessaire de normaliser les valeurs avant l'affichage:

1
import matplotlib.pyplot as plt
2
fig, ax= plt.subplots(figsize=(6, 5))
3
plt.imshow(img_rgbnir.data transpose(1,2,0)/2500.0)
4
plt.show()
3 Réhaussement et visualisation d'images

<figure><img src="image_30.png" /></figure>

On remarquera les valeurs des axes x et y avec une origine en haut à gauche. Ceci est un référentiel purement matriciel (lignes et colonnes); autrement dit, il n'y a pas ici de géoréférence. Pour pallier à cette limitation, les libraries rasterio et xarray proposent une extension de la fonction imshow permettant d'afficher les coordonnées cartographiques ainsi qu'un contrôle la dynamique de l'image:

1 import rioxarray as rxr

2 fig, ax= plt.subplots(figure=(6, 5))

3 img_rgbnir sel(band=[1,2,3]). plot .imshow( vmin=86, vmax=5000)

4 ax.set_title('Imshow avec rioxarray')

5 plt.show()
<figure><img src="image_31.png" /></figure>

## 3.3 Réhaussements visuels

Le réhaussement visuel d'une image vise principalement à améliorer la qualité visuelle d'une image en améliorant le contraste, la dynamique ou la texture d'une image. De manière générale, ce réhaussement ne modifie pas la donnée d'origine mais il est appliquée dynamiquement à l'affichage pour des fins d'inspection visuelle. Le réhaussement nécessite généralement une connaissance des caractéristiques statistiques d'une image. Ces statistiques sont ensuite exploitées pour appliquer diverses transformations linéaires ou non linéaires.

### 3.3.1 Statistiques d'une image

On peut considérer un ensemble de statistique pour chacune des bandes d'une image:

- valeurs minimales et maximales

- valeurs moyennes,

- Quartiles (1er quartile, médiane et 3ième quartile), quantiles et percentiles.

- écart-type, et coefficients d'asymétrie (*skewness*) et d'applatissement (*kurtosis*)

Ces statistiques doivent être calculées pour chaque bande d'une image multispectrale.
3 Réhaussement et visualisation d'images

En ligne de commande, `gdalinfo` permet d'interroger rapidement un fichier image pour connaitre ces statistiques univariées de base:

1

```
!gdalinfo -stats landsat7.tif
```

```
Driver: GTiff/GeoTIFF
Files: landsat7.tif
        landsat7.tif.aux.xml
Size is 2181, 1917
Coordinate System is:
PROJCS["WGS 84 / Pseudo-Mercator",
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,
AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4326"]],
PROJECTION["Mercator_1SP",
PARAMETER["central_meridian",0],
PARAMETER["scale_factor",1],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]],
AXIS["X",EAST],
AXIS["Y",NORTH],
EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"],
AUTHORITY["EPSG","3857"]]
Origin = (-13651650.00000000000000,4576290.00000000000000)
Pixel Size = (30.0000000000000000,-30.0000000000000000)
Metadata:
AREA_OR=Area
OVR_RE SAMPLING_ALG=NEAREST
TIFFTAG_RE SOLUTIONUNIT=1 (unitless)
TIFFTAG_XRESOLUTION=1
TIFFTAG_YRESOLUTION=1
Image Structure Metadata:
COMPRESSION=DEFLATE
INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left (-13651650.000, 4576290.000) (122d43' 5.49"W, 37d58'40.08"N)
Lower Left (-13651650.000, 4518780.000) (122d43' 5.49"W, 37d34'10.00"N)
Upper Right (-13586220.000, 4576290.000) (122d 2'49.53"W, 37d58'40.08"N)
Lower Right (-13586220.000, 4518780.000) (122d 2'49.53"W, 37d34'10.00"N)
Center (-13618935.000, 4547535.000) (122d20'27.51"W, 37d46'26.05"N)
Band 1 Block=512x512 Type=Byte, ColorInterp=Red
Min=19.000 Max=233.000
Minimum=19.000, Maximum=233.000, Mean=98.433, StdDev=21.164
NoData Value=0
Overviews: 1091x959, 546x480
Metadata:
STATISTICS_MAXIMUM=233
STATISTICS_MEAN=98.433096940153
STATISTICS_MINIMUM=19
STATISTICS_STDDEV=21.164021026458
Band 2 Block=512x512 Type=Byte, ColorInterp=Green
Min=19.000 Max=178.000
Minimum=19.000, Maximum=178.000, Mean=55.068, StdDev=22.204
NoData Value=0
Overviews: 1091x959, 546x480
Metadata:
STATISTICS_MAXIMUM=178
STATISTICS_MEAN=55.067787534804
STATISTICS_MINIMUM=19
STATISTICS_STDDEV=22.203571974581
Band 3 Block=512x512 Type=Byte, ColorInterp=Blue
Min=19.000 Max=187.000
Minimum=19.000, Maximum=187.000, Mean=43.341, StdDev=20.330
NoData Value=0
Overviews: 1091x959, 546x480
Metadata:
STATISTICS_MAXIMUM=187
STATISTICS_MEAN=43.340507443056
STATISTICS_MINIMUM=19
STATISTICS_STDDEV=20.32987736339
```

Les librairies de base comme `rasterio` et `xarray` produisent facilement un sommaire des statistiques de base avec la fonction `stats`:
1 import rasterio as rio

2 import numpy as np

3 with rio.open('landsat7.tif') as src:

4 stats= src.stats()

5 print.stats()

La librairie `xarray` donne accès à des fonctionnalités plus sophistiquées comme le calcul des quantiles:

1 import rioxarray as riox

2 with riox.open_rasterio('landsat7.tif', masked= True ) as src:

3 print(src)

4 quantiles = src.quantile(dim=[ 'x','y'], q=[.025,.25,.5,.75,.975])

5 quantiles

<xarray.DataArray (band: 3, y: 1917, x: 2181)> Size: 50MB
[12542931 values with dtype浮点型]
Coordinates:
  * band         (band) int64 24B 1 2 3
  * x            (x) float64 17kB -1.365e+07 -1.365e+07 ... -1.359e+07
  * y            (y) float64 15kB 4.576e+06 4.576e+06 ... 4.519e+06 4.519e+06
  spatial_ref   int64 8B 0

Attributes:

<table>
  <tr>
    <td>AREA_OR_Point:</td>
    <td>Area</td>
  </tr>
  <tr>
    <td>OVR_RESAMPLING_ALG:</td>
    <td>NEAREST</td>
  </tr>
  <tr>
    <td>TIFFTAG RESOLUTIONUNIT:</td>
    <td>1 (unitless)</td>
  </tr>
  <tr>
    <td>TIFFTAG_XRESOLUTION:</td>
    <td>1</td>
  </tr>
  <tr>
    <td>TIFFTAG_YRESOLUTION:</td>
    <td>1</td>
  </tr>
  <tr>
    <td>STATISTICS_MAXIMUM:</td>
    <td>233</td>
  </tr>
  <tr>
    <td>STATISTICS_MEAN:</td>
    <td>98.433096940153</td>
  </tr>
  <tr>
    <td>STATISTICS_MINIMUM:</td>
    <td>19</td>
  </tr>
  <tr>
    <td>STATISTICS_STDDEV:</td>
    <td>21.164021026458</td>
  </tr>
  <tr>
    <td>scale_factor:</td>
    <td>1.0</td>
  </tr>
  <tr>
    <td>add_offset:</td>
    <td>0.0</td>
  </tr>
</table>

<xarray.DataArray (quantile: 5, band: 3)> Size: 120B
array([ 54., 19., 19.],
    [ 85., 38., 27.],
    [ 99., 54., 38.],
    [111., 69., 57.],
    [140., 102., 89.]])

Coordinates:

* band (band) int64 24B 1 2 3

* quantile (quantile) float64 40B 0.025 0.25 0.5 0.75 0.975

3.3.1.1 Calcul de l'histogramme

Le calcul d'un histogramme pour une image (une bande) permet d'avoir une vue plus détaillée de la répartition
des valeurs radiométriques. Le calcul d'un histogramme nécessite minimalement de faire le choix du nombre de
barre ( bins ou de la largeur ). Un bin est un intervalle de valeurs pour lequel on peut calculer le nombre de
valeurs observées dans l'image. La fonction de base pour ce type de calcul est la fonction numpy.histogram() :

1 import numpy as np

2 array = np.randomonde, 10, 100) # 100 valeurs aléatoires entre 0 et 10

3 hist, bin_limites = np.histogram( array, density = True )

4 print('valeurs : ', hist)

5 print(' ; imites : ', bin_limites)

valeurs : [0.1 0.1333333 0.0777778 0.1111111 0.0777778 0.0777778
0.1666667 0.1111111 0.1333333 0.1222222]
; imites : [0. 0.9 1.8 2.7 3.6 4.5 5.4 6.3 7.2 8.1 9.]
3 Réhaussement et visualisation d'images

Le calcul se fait avec 10 intervalles par défaut.

1 fig, ax= plt.subplots(figsize=(5, 4))

2 plt.barbin_limites[:-1], hist)

3 plt.show()

<figure><img src="image_32.png" /></figure>

Pour des besoins de visualisation, le calcul des valeurs extrêmes de l'histogramme peut aussi se faire via les
quantiles comme discutés auparavant.

Visualisation des histogrammes

La librairie rasterio est probablement l'outil le plus simples pour visualiser rapidement des histogrammes sur
une image multi-spectrale:

1 import rasterio as rio

2 from rasterio .plot import show_histogram

3 with rio.open('RGBNIR_of_S2A.tif') as src:

4 show_histogram(src, bins=50, lw=0.0, stacked=False, alpha=0.3, histtype='stepfilled', title="Histogram")
<figure><img src="image_33.png" /></figure>

### 3.3.2 Réhaussements linéaires

Le réhaussement linéaire (linear stretch) d'une image est la forme la plus simple de réhaussement, elle consiste à 1) optimiser les valeurs des pixels d'une image afin de maximiser la dynamique disponibles à l'affichage, ou 2) à changer le format de stockage des valeurs (de 8 bits à 16 bits):

$$ \text{nouvelle valeur d'un pixel} = \frac{\text{valeur d'un pixel} - min_0}{max_0 - min_0} \times (max_1 - min_1) + min_1 \quad (3.1) $$

Par cette opération, on passe de la dynamique de départ ($max_0 - min_0$) vers la dynamique cible ($max_1 - min_1$). Bien que cette opération semble triviale, il est important d'être conscient des trois contraintes suivantes:

1. Faire attention à la dynamique cible, ainsi, pour sauvegarder une image en format 8 bit, on utilisera alors $max_1 = 255$ et $min_1 = 0$.

2. Préservation de la valeur de no data : il faut faire attention à la valeur $min_1$ dans le cas d'une valeur présente pour no_data. Par exemple, si no_data=0 alors il faut s'assurer que $min_1 > 0$.

3. Précision du calcul : si possible réaliser la division ci-dessus en format float

#### 3.3.2.1 Cas des histogrammes asymétriques

Dans certains cas, la distribution de valeurs est très asymétrique et présente une longue queue avec des valeurs extrêmes élevées (à droite ou à gauche de l'histogramme). Le cas des images SAR est particulièrement représentatif de ce type de données. En effet, celles-ci peuvent présenter une distribution de valeurs de type exponentiel. Il est alors préférable d'utiliser des percentiles au préalable afin d'explorer la forme de l'histogramme et la distribution des valeurs:
NO_DATA_FLOAT=-999.0

# on prend tous les pixels de la première bande

values = img_SAR[0].values.flatten().astype(float)

# on exclut les valeurs invalides

values = values[~np.isnan(values)]

# on exclut le no data

values = values[values!=NO_DATA_FLOAT]

# calcul des percentiles

percentiles_position= (0,0.1,1,2,50,98,99,99.9,100)

percentiles= dict(zip(percentiles_position, np))

print(percentiles)

{0: np.float64(8.172580237442162e-06), 0.1: np.float64(1.588739885482937e-05), 1: np.float64(8.657756850880105e-05), 2:

→ np.float64(0.00018846066552214325), 50: np.float64(0.012372820172458887), 98: np.float64(0.1719470709562302), 99:

→ np.float64(0.27963151514529694), 99.9: np.float64(1.5235805057287233), 100: np.float64(483.223876953125)

On constate que la valeur médiane (0.012) est très faible, ce qui signifie que 50% des valeurs sont inférieures
à cette valeur alors que la valeur maximale (483) est 10 000 fois plus élevée! Une manière de visualiser cette
distribution de valeurs est d'utiliser boxplot et violinplot de la librairie matplotlib:

1 fig, ax = plt.subplots(nrows=2, ncols=1, figsize=(6, 4), sharex=True)

2 ax[0].set_title('Distribution de la bande 0 de img_SAR', fontsize='small')

3 ax[0].grid(True)

4 ax[0].violinplot(values, orientation ='horizontal',
                  quantiles =(0.01,0.02,0.50,0.98,0.99),

6 showmeans=False,

7 showmedians=True)

8 ax[1].set_xlabel('Valeur des pixels')

9 ax[1].grid(True)

10 bplot = ax[1].boxplot(values, notch = True, orientation ='horizontal')

11 plt.tight_layout()

12 plt.show()
3 Réhaussement et visualisation d'images

<figure><img src="image_34.png" /></figure>

Afin de visualiser correctement l'histogramme, il faut se limiter à un intervalle de valeurs plus réduit. Dans le code ci-dessous, on impose à la fonction np.histogramme de compter les valeurs de pixels dans des intervalles de valeurs fixés par la fonction np.linspace(percentiles[0.1], percentiles[99.9], 50) où percentiles[0.1] et percentiles[99.9] sont les 0.1% et 99.9% percentiles respectivement:

```python
1 hist, bin_edges = np.histogram(values,
2                                bins= np.linspace(percentiles[0.1],
3                                percentiles[99.9], 50),
4                                density = True)

6 fig, ax = plt.subplots(nrows = 2, ncols = 1, figsize = (6, 5), sharex = True)
7 ax[0].bar bin_edges [:-1],
8         hist *(bin_edges [1] -bin_edges [0]),
9         width = (bin_edges [1] -bin_edges [0]),
10        edgecolor = 'w')
11 ax[0].set_title("Distribution de probabilité (PDF)")
12 ax[0].set_ylabel("Densité de probabilité")
13 ax[0].grid(True)

15 ax[1].plot (bin_edges [:-1],
16           hist .cumsum()*(bin_edges [1]-bin_edges [0]))
17 ax[1].set_title("Distribution de probabilité cumulée (CDF)")
18 ax[1].set_xlabel("Valeur du pixel")
19 ax[1].set_ylabel("Probabilité cumulée")
20 ax[1].grid(True)
```
21 plt.tight_layout()

22 plt.show()

<figure><img src="image_35.png" /></figure>

Au niveau de l'affichage avec `matplotlib`, la dynamique peut être contrôlée directement avec les paramètres `vmin` et `vmax` comme ceci:

1 fig, ax = plt.subplots(nrows=2, ncols=2, figsize=(6, 5), sharex=True, sharey=True)

2 [ax axis('off') for a in ax.flatten()]]

3 ax[0,0].imshow(img_SAR[0].values, vmin=percentiles[0], vmax=percentiles[100])

4 ax[0,0].set_title("f"0% - 100%={percentiles[0]:2.1f} - {percentiles[100]:2.1f}")

5 ax[0,1].imshow(img_SAR[0].values, vmin=percentiles[0.1], vmax=percentiles[99.9])

6 ax[0,1].set_title("f"0.1% - 99.9%={percentiles[0.1]:2.1f} - {percentiles[99.9]:2.1f}")

7 ax[1,0].imshow(img_SAR[0].values, vmin=percentiles[1], vmax=percentiles[99])

8 ax[1,0].set_title("f"1% - 99%={percentiles[1]:2.1f} - {percentiles[99]:2.1f}")

9 ax[1,1].imshow(img_SAR[0].values, vmin=percentiles[2], vmax=percentiles[98])

10 ax[1,1].set_title("f"2% - 98%={percentiles[2]:2.1f} - {percentiles[98]:2.1f}")

11 plt.tight_layout()
0% - 100%=0.0 - 483.2

0.1% - 99.9%=0.0 - 1.5

<figure><img src="image_36.png" /></figure>

<figure><img src="image_37.png" /></figure>

1% - 99%=0.0 - 0.3

2% - 98%=0.0 - 0.2

<figure><img src="image_38.png" /></figure>

<figure><img src="image_39.png" /></figure>

### 3.3.3 Réhaussements non linéaires

#### 3.3.3.1 Réhaussement par fonctions

Le réhaussenent par fonction consiste à appliquer une fonction non linéaire afin de modifier la dynamique de l'image. Par exemple, pour une image radar, une transformation populaire est d'afficher les valeurs de rétrodiffusion en décibel (dB) avec la fonction `log10()`.

1. `percentiles_position`= (0,0.1,1,2,50,98,99,99.9,100)

2. `values`= `xr_apply_ufunc(lambda x: 10 * np.log10(x), img_SAR[0]).data`

3. `percentiles_db`= `dict(zip(percentiles_position, np))`

4. `print(percentiles_db)`

(0: np.float64(-50.87641143798828), 0.1: np.float64(-47.98947440338135), 1: np.float64(-40.625951309204105), 2: np.float64(-37.24779266357422),
    ↪→ 50: np.float64(-19.075313568115234), 98: np.float64(-7.6460519409179675), 99: np.float64(-5.534139237403945), 99.9:
    ↪→ np.float64(1.8286541925668764), 100: np.float64(26.84148406982422))

Les boites à moustache (boxplots) ont une bien meilleure distribution qui est en effet très proche d'une distribution normale gaussienne:
3 Réhaussement et visualisation d'images

<figure><img src="image_40.png" /></figure>

On obtient ainsi les images suivantes:
<figure><img src="image_41.png" /></figure>

### 3.3.3.2 Égalisation d'histogramme

L'égalisation d'histogramme consiste à modifier les valeurs des pixels d'une image source afin que la distribution cumulée des valeurs (CDF) devienne similaire à celle d'une image cible. La CDF (Cumulative Distribution Function) est simplement la somme cumulée des valeurs de l'histogramme:

$$CDF_{source}(i) = \frac{1}{K} \sum_{j=0}^{j \le i} hist_{source}(j)$$

avec K choisit de façon à ce que la dernière valeur soit égale à 1 ($CDF_{source}(i_{max}) = 1$). De la même manière, $CDF_{cible}$ est la CDF d'une image cible. La formule générale pour l'égalisation d'histogramme est la suivante:

$$j = CDF_{cible}^{-1}(CDF_{source}(i))$$

On peut choisir $CDF_{cible}$ comme correspondant à une image où chaque valeur de pixel est équiprobable (d'où le terme égalisation), ce qui veut dire $hist_{cible}(j) = 1/L$ avec L égale au nombre de valeurs possibles dans l'image (par exemple $L = 256$).

$$j = L \times CDF_{source}(i)$$

On peut appliquer cette procédure sur l'image SAR en dB de la façon suivante:
1 values = np.sort(10 * np.log10(img_SAR[0].data .flatten()))

2 cdf_x = np.linspace(values[0], values[-1], 1000)

3 cdf_source = np.interp(cdf_x, values, np.arange(len(values))/len(values)*255)

4 values_eq =np.interp(np.log10(img_SAR[0].data), cdf_x, cdf_source).astype('uint8')

5 plt.imshow(values_eq)

6 plt.axis('off')

<figure><img src="image_42.png" /></figure>

3.3.3.3 Palettes de couleur

Les palettes de couleurs sont appliquées dynamiquement à l'affichage sur une image à une seule bande. La librairie `matplotlib` contient un nombre considérable de palettes.

1 from matplotlib import colormaps

2 list(colormaps)

Voici quelques exemples ci-dessous, les valeurs de l'image doivent être normalisées entre 0 et 1 ou entre 0 et 255 sinon les paramètres `vmin` et `vmax` doivent être spécifiés. On peut observer comment ces palettes révèlent les détails de l'image malgré une image originalement très sombre.

1 fig, ax = plt.subplots(nrows=2, ncols=2, figsize=(6, 5), sharex=True, sharey=True)

2 [a axis('off') for a in ax .flatten()]]

3 ax[0,0].imshow(img_SAR[0].data, vmin=percentiles[2], vmax=percentiles[98], cmap='jet')

4 ax[0,0].set_title('f"jet")

5 ax[0,1].imshow(img_SAR[0].data, vmin=percentiles[2], vmax=percentiles[98], cmap='hot')

6 ax[0,1].set_title('f"hot")

7 ax[1,0].imshow(img_SAR[0].data, vmin=percentiles[2], vmax=percentiles[98], cmap='hsv')

8 ax[1,0].set_title('f"hsv")

9 ax[1,1].imshow(img_SAR[0].data, vmin=percentiles[2], vmax=percentiles[98], cmap='terrain')
10 ax[1,1].set_title(f"terrain")

11 plt.tight_layout()

<figure><img src="image_43.png" /><figcaption>hsv</figcaption></figure>

<figure><img src="image_44.png" /><figcaption>hot</figcaption></figure>

<figure><img src="image_45.png" /><figcaption>terrain</figcaption></figure>

<figure><img src="image_46.png" /></figure>

Il peut être utile d'ajouter une barre de couleurs afin d'indiquer la correspondance entre les couleurs et les valeurs numériques:

1 import matplotlib as mpl

2 fig, ax = plt.subplots(figsize=(6, 6))

3 cmap=mpl.colormaps.get_cmap('jet').with_extremes(under='white', over='magenta')

4 h= plt.imshow(img_SAR[0].data, norm=mplcolors.LogNorm(vmin=percentiles[2], vmax=percentiles[98]),
                 cmap=cmap)

5 fig.colorbar(h, ax=ax, orientation='horizontal', label="Intensité", extend='both')

6 axaxis('off')
<figure><img src="image_47.png" /></figure>

<figure><img src="image_48.png" /></figure>

### 3.3.4 Composés colorés

Le système visuel humain est sensible seulement à la partie visible du spectre électromagnétique qui compose les couleurs de l'arc-en-ciel du bleu au rouge. L'ensemble des couleurs du spectre visible peut être obtenu à partir du mélange de trois couleurs primaires (rouge, vert et bleu). Ce système de décomposition à trois couleurs est à la base de la plupart des systèmes de visualisation ou de représentation de l'information de couleur. Si on prend le cas des images Sentinel-2, 12 bandes sont disponibles, plusieurs composés couleurs sont donc possibles (voir le site de Copernicus). Voici quelques exemples possibles, chaque composé mettant en valeur des propriétés différentes de la surface.

1 import rioxarray as rxr

2 fig, ax= plt.subplots(nrows=2, ncols=2, figsize=(8, 6), sharex=True, sharey=True)

3 img_s2 sel(band=[4,3,2]).plot.imshow(vmin=86, vmax=4000, ax=ax[0,0])

4 ax[0,0].set_title('RVB')

5 img_s2 sel(band=[8,3,2]).plot.imshow(vmin=86, vmax=4000, ax=ax[0,1])

6 ax[0,1].set_title('NIR,V,B')

7 img_s2 sel(band=[12,8,4]).plot.imshow(vmin=86, vmax=4000, ax=ax[1,0])

8 ax[1,0].set_title('SWIR2,NIR,R')

9 img_s2 sel(band=[12,11,4]).plot.imshow(vmin=86, vmax=4000, ax=ax[1,1])
3 Réhaussement et visualisation d'images

10

ax[1,1].set_title('SWIR2,SWIR1,NIR')

11

plt.tight_layout()

12

plt.show()

<figure><img src="image_49.png" /></figure>
# Partie 2. Transformations des données satellitaires
# 4 Transformations spectrales

## 4.1 Préambule

Assurez-vous de lire ce préambule avant d'exécuter le reste du notebook.

### 4.1.1 Objectifs

Dans ce chapitre, nous abordons l'exploitation de la dimension spectrale des images satellites. Ce chapitre est aussi disponible sous la forme d'un notebook Python:

<figure><img src="image_50.png" /></figure>

<figure><img src="image_51.png" /></figure>

#### Objectif

#### Objectifs d'apprentissage visés dans ce chapitre

À la fin de ce chapitre, vous devriez être en mesure de :

- comprendre le principe des indices spectraux;

- calculer différents indices avec spyndex;

- analyser le gain en information des indices;

### 4.1.2 Librairies

Les librairies qui vont être explorées dans ce chapitre sont les suivantes:

- SciPy

- NumPy

- spyindex

- Rasterio

- Xarray

- rioxarray

Dans l'environnement Google Colab, seul rioxarray doit être installés

Vérifiez les importations:

<footnote n="1">1. `%%capture`</footnote>

<footnote n="2">2. `!pip install -qU matplotlib rioxarray xrscipy scikit-image pyarrow spyndex`</footnote>
1 import numpy as np

2 import rioxarray as rxr

3 from Kapropy import signal

4 import xarray as xr

5 import xrscipy

6 import matplotlib.pyplot as plt

7 import spyndex

8 import rasterio as rio

4.1.3 Images utilisées

Nous utilisons les images suivantes dans ce chapitre:

1 %%capture

2 import gdown

3

4 gdownerna(↓
  ↪→ 'https://drive.google.com/uc?export=download&confirm=pbef&id=1a6Ypg0g1Oy4AJt9XWKWfnR12NW1XhNg_',
  ↪→ output= 'RGBNIR_of_S2A.tif')

5 gdown. download()
  ↪→ 'https://drive.google.com/uc?export=download&confirm=pbef&id=1a603L_ab0fU7h94K22At8qtBuLMGErwo',
  ↪→ output= 'sentinel2.tif')

6 gdown. download()
  ↪→ 'https://drive.google.com/uc?export=download&confirm=pbef&id=1_zwCLN-x7XJcNHJCH6Z8upEdUxtVtvs1',
  ↪→ output= 'berkeley.jpg')

7 gdown. download()
  ↪→ 'https://drive.google.com/uc?export=download&confirm=pbef&id=1dM6IVqjba6GHwTLmI7CpX8GP2z5txUq6',
  ↪→ output= 'SAR.tif')

8 gdown. download()
  ↪→ 'https://drive.google.com/uc?export=download&confirm=pbef&id=1aAq7crc_LoaLC3kG3HkQ6Fv5JfG0mswg',
  ↪→ output= 'carte.tif')

9 gdown. download()
  ↪→ 'https://drive.google.com/uc?export=download&confirm=pbef&id=1iCZNYTv0qEZRzPhe22nPdpV4Ks7NsY3b',
  ↪→ output= 'ASCIIdata_splib07b_rsSentinel2.zip')

10 !unzip -q ASCIIdata_splib07b_rsSentinel2.zip

Vérifiez que vous êtes capable de les lire :

1 with rxr.open_rasterio('berkeley.jpg', mask_and_scale= True ) as img rgb:
2 print(img rgb)

3 with rxr.open_rasterio('RGBNIR_of_S2A.tif', mask_and_scale= True ) as img_rgnir:
4 print(img_rgnir)

5 with rxr.open_rasterio('sentinel2.tif', mask_and_scale= True ) as img_s2:
6 print(img_s2)
7

8

```python
with rxr.open_rasterio('carte.tif', mask_and_scale= True) as imgCarte:
    print(imgCarte)
```

## 4.2 Qu'est ce que l'information spectrale?

L'information spectrale touche à l'exploitation de la dimension spectrale des images (c.à.d le long des bandes spectrales de l'image). La taille de cette dimension spectrale dépend du type de capteurs considéré. Un capteur à très haute résolution spectrale par exemple aura très peu de bandes (4 ou 5). Un capteur multispectral pourra contenir une quinzaine de bande. À l'autre extrême, on trouvera les capteurs hyperspectraux qui peuvent contenir des centaines de bandes spectrales.

<figure><img src="image_52.png" /><figcaption>FIGURE 4.1 – Positions des bandes spectrales pour quelques capteurs (source)</figcaption></figure>

Pour une surface donnée, la forme des valeurs le long de l'axe spectrale caractérise le type de matériau observé ainsi que son état. On parle souvent alors de signature spectrale. On peut voir celle-ci comme une généralisation de la couleur d'un matériau au delà des bandes visibles du spectre. L'exploitation de ces signatures spectrales est probablement un des principes les plus importants en télédétection qui le distingue de la vison par ordinateur. L'USGS maintient une base de données spectrales acquises en laboratoire (Kokaly et Klein 2017). On peut observer sur la figure ci-dessous comment la forme et l'amplitude de trois signatures différentes peut changer en fonction du type de surface.

```python
Text(0, 0.5, 'Réflectance')
```

Exemples de signatures spectrales pour trois surfaces différentes pour les bandes spectrales de Sentinel-2
<figure><img src="image_53.png" /></figure>

## 4.3 Indices spectraux

Il existe une vaste littérature sur les indices spectraux, le choix d'un indice plutôt qu'un autre dépend fortement de l'application visée, nous allons simplement couvrir les principes de base ici. Le principe d'un indice spectral consiste à mettre en valeur certaines caractéristiques saillantes du spectre comme des pentes, des gradients, etc.

La librairie Python *Awesome Spectral Indices* maintient une liste de plus de 200 indices spectraux (radar et optiques). La liste complète est affichable avec la commande suivante:

1

### spyndex.indices

```
SpectralIndices([('AFRI1600', 'AFRI2100', 'ANDWI', 'ARI', 'ARI2', 'ARVI', 'ATSAVI', 'AVI', 'AWEInsh', 'AWEIsh', 'BAI', 'BAIM', 'BAIS2', 'BCC',
'BI', 'BITM', 'BIXS', 'BLFEI', 'BNDVI', 'BRBA', 'BWDRI', 'BaI', 'CCI', 'CIG', 'CIRE', 'CRI550', 'CRI700', 'CSI', 'CSIT', 'CVI', 'DBI',
'DBSI', 'DPDD', 'DSI', 'DSWI1', 'DSWI2', 'DSWI3', 'DSWI4', 'DSWI5', 'DVI', 'DVIplus', 'DpRVHH', 'DpRVVV', 'EBBI', 'EBI', 'EMBI', 'ENDVI',
'EVI', 'EVI2', 'EVIv', 'ExG', 'ExGR', 'ExR', 'FAI', 'FCVI', 'GARI', 'GBNDVI', 'GCC', 'GDVI', 'GEMI', 'GLI', 'GM1', 'GM2', 'GNDVI', 'GOSAVI',
'GRNDVI', 'GRVI', 'GSAVI', 'GVMI', 'IAVI', 'IBI', 'IKAW', 'IPVI', 'IRECI', 'LSWI', 'MBI', 'MBWI', 'MCARI', 'MCARI1', 'MCARI2', 'MCARI705',
'MCARIOSAVI', 'MCARIOSAVI705', 'MGRVI', 'MIRBI', 'MLSWI26', 'MLSWI27', 'MNDVI', 'MNDWI', 'MNLI', 'MRBVI', 'MSAVI', 'MSI', 'MSR', 'MSR705',
'MTCI', 'MTVI1', 'MTVI2', 'MuWIR', 'NBAI', 'NBLI', 'NBLIOLI', 'NBR', 'NBR2', 'NBRSWIR', 'NBRTI', 'NBRT2', 'NBRT3', 'NBRplus', 'NBSIMS',
'NBUI', 'ND705', 'NDBI', 'NDBaI', 'NDCI', 'NDDI', 'NDGI', 'NDGlaI', 'NDII', 'NDISIb', 'NDISIg', 'NDISImndwi', 'NDISIndwi', 'NDISIr', 'NDMI',
'NDPI', 'NDPoli', 'NDPonI', 'NDREI', 'NDSI', 'NDSII', 'NDSIWV', 'NDSInw', 'NDSWIR', 'NDSAII', 'NDSol', 'NDTI', 'NDVI', 'NDVI705', 'NDVIMNDWI',
'NDVIT', 'NDWI', 'NDWIns', 'NDYI', 'NGRDI', 'NHFD', 'NIRv', 'NIRvH2', 'NIRvP', 'NLI', 'NMDI', 'NRFIg', 'NRFIr', 'NSDS', 'NSDSII', 'NSDSI2',
'NSDSI3', 'NSTv1', 'NSTv2', 'NWI', 'NormG', 'NormNIR', 'NormR', 'OCVI', 'OSAVI', 'OSI', 'PI', 'PISI', 'PSRI', 'QpRVI', 'RCC', 'RDVI', 'REDSI',
'RENDVI', 'RFDI', 'RGBVI', 'RGRI', 'RI', 'RI4XS', 'RNDVI', 'RVI', 'S2REP', 'S2WI', 'S3', 'SARVI', 'SAVI', 'SAVI2', 'SAVIT', 'SEVI', 'SI',
'SIPI', 'SLAVI', 'SR', 'SR2', 'SR3', 'SR555', 'SR705', 'SWI', 'SWM', 'SeLI', 'TCARI', 'TCARIOSAVI', 'TCARIOSAVI705', 'TCI', 'TDVI', 'TGI',
'TRRVI', 'TSAVI', 'TIVI', 'TVI', 'TWI', 'TriVI', 'UI', 'VARI', 'VARI700', 'VDDPI', 'VHVVD', 'VHVVP', 'VHVVR', 'VI6T', 'VI700', 'VIBI', 'VIG',
'VVVHD', 'VVVHR', 'VVVHS', 'VgNIRBI', 'VrNIRBI', 'WDRVI', 'WDVI', 'WII', 'WI2', 'WI2015', 'WRI', 'bNIRv', 'kEVI', 'kIPVI', 'kNDVI', 'kRVI',
'kVARI', 'mND705', 'mSR705', 'sNIRvLSWI', 'sNIRvNDPI', 'sNIRvNDVILSWIP', 'sNIRvNDVILSWIS', 'sNIRvSWIR'])
```
Le détail d'un indice particulier, par exemple le 'NDVI', est aussi affichable:

1 `spyndex.indices["NDVI"]`

SpectralIndex(NDVI: Normalized Difference Vegetation Index)

* Application Domain: vegetation

* Bands/Parameters: ['N', 'R']

* Formula: $(N-R)/(N+R)$

* Reference: https://ntrs.nasa.gov/citations/19740022614

`spyndex` pré-suppose une nomenclature prédéfinie des bandes, on peut voir la correspondance sur le tableau ci-dessous:

1 `spyndex.bands`

Bands(['A', 'B', 'G', 'G1', 'N', 'N2', 'R', 'RE1', 'RE2', 'RE3', 'S1', 'S2', 'T', 'T1', 'T2', 'WV', 'Y'])

<table><caption>TABLEAU 4.1 – Noms des bandes Sentinel-2</caption><thead><tr><th>Index</th><th>Noms</th><th>Spyndex</th><th>Noms</th></tr></thead><tbody><tr><td>1</td><td>B01</td><td>A</td><td>Aérosol</td></tr><tr><td>2</td><td>B02</td><td>B</td><td>Bleu</td></tr><tr><td>3</td><td>B03</td><td>G</td><td>Vert</td></tr><tr><td>4</td><td>B04</td><td>R</td><td>Rouge</td></tr><tr><td>5</td><td>B05</td><td>RE1</td><td>Red edge 1</td></tr><tr><td>6</td><td>B06</td><td>RE1</td><td>Red edge 2</td></tr><tr><td>7</td><td>B07</td><td>RE2</td><td>Red edge 3</td></tr><tr><td>8</td><td>B08</td><td>N</td><td>Proche-infrarouge 1</td></tr><tr><td>9</td><td>B08A</td><td>N2</td><td>Proche-infrarouge 2</td></tr><tr><td>10</td><td>B09</td><td>-</td><td>Vapeur d'eau</td></tr><tr><td>11</td><td>B11</td><td>S1</td><td>Infra-rouge onde courte 1</td></tr><tr><td>12</td><td>B12</td><td>S2</td><td>Infra-rouge onde courte 1</td></tr></tbody></table>

Deux options sont possibles, on peut soit renommer les noms des bandes avec `xarray` ou "mapper" les noms vers les noms appropriés. Regardons les dimensions de notre jeux de données:

1 `img_s2.dims`

{'band', 'y', 'x'

On peut simplement changer les index (`coordinates`) de la dimension `band`:

1 `img_s2 = img_s2 sel(band = list(range(1,13))).assign�ers({'band': ['A', 'B', 'G', 'R', 'RE1', 'RE2', 'RE3', 'N', 'N2', 'WV', 'S1', 'S2']})`

2 `img_s2 = img_s2/10000 # normalisation en réflectance`
```python
from rasterio import plot
idx = spyndex.computeIndex(
    index = ["NDVI", "GNDVI", "SAVI"],
    params = {
        "N": img_s2 sel(band = "N"),
        "R": img_s2 sel(band = "R"),
        "G": img_s2 sel(band = "G"),
        "L": 0.5
    }
)

# Plot the indices (and the RGB image for comparison)
fig, ax = plt.subplots(2, 2, figsize = (9, 9))
ax.axis('off') for a in ax.flatten())
plot.show(img_s2 sel(band = ["R", "G", "B"]).data / 0.3, ax=ax[0,0], title = "RGB")
plot.show(idx sel(index = "NDVI"), ax=ax[0,1], title = "NDVI")
plot.show(idx sel(index = "GNDVI"), ax=ax[1,0], title = "GNDVI")
plot.show(idx sel(index = "SAVI"), ax=ax[1,1], title = "SAVI")
```
RGB

<figure><img src="image_54.png" /></figure>

NDVI

<figure><img src="image_55.png" /></figure>

GNDVI

<figure><img src="image_56.png" /></figure>

SAVI

<figure><img src="image_57.png" /></figure>

On peut vérifier l'utilité des indices en vérifiant leur séparabilité pour certaines classes d'intérêts. Nous reprenons ici l'exemple de la section 6.2.3 pour vérifier l'utilité des indices NDVI, NDWI et NDBI:

```python
import pandas as pd
```

```python
import seaborn as sns
```

```python
# On sélectionne trois classes
```

```python
class selected= [1,3,9]
```
```python
df = pd.concat([gdf[gdf['class'] == c] for c in classSelected], ignore_index=True)
idx['Land Cover'] = [nom_classes[l] for l in df['class'].list()]

# Compute the desired spectral indices
idx = spyndex.computeIndex(
    index = ['NDVI', 'NDWI', 'NDBI'],
    params = {
        "N": df['SR_B8'],
        "R": df['SR_B4'],
        "G": df['SR_B3'],
        "S1": df['SR_B11']
    }
)

colors = [couleurs_classes[c] for c in classSelected]

# Plot a pairplot to check the indices behaviour
plt.figure(figsize=(15, 15))
g = sns.PairGrid(idx, hue = "Land Cover", palette = sns.color palette(colors))
g.map_lower(sns.scatterplot)
g.map_upper(sns.kdeplot, fill = True, alpha = .5)
g.map_diag(sns.kdeplot, fill = True)
g.add legend()
plt.show()
```
4 Transformations spectrales

<figure><img src="image_58.png" /><figcaption>FIGURE 4.2 – Visualisation des points d'une image Sentinel-2 pour trois classes</figcaption></figure>

4.4 Quiz

Utilisez la version html.
# 5 Transformations spatiales

## 5.1 Préambule

Assurez-vous de lire ce préambule avant d'exécuter le reste du notebook.

### 5.1.1 Objectifs

Dans ce chapitre, nous abordons quelques techniques de traitement d'images dans le domaine spatial uniquement.
Ce chapitre est aussi disponible sous la forme d'un notebook Python sur Google Colab:

<figure><img src="image_59.png" /></figure>

<figure><img src="image_60.png" /></figure>

#### Objectif

#### Objectifs d'apprentissage visés dans ce chapitre

À la fin de ce chapitre, vous devriez être en mesure de :

- comprendre le principe de la décomposition de Fourier;

- comprendre le principe de la convolution;

- appliquer un filtrage local à l'aide d'une fenêtre;

- segmenter une image en super-pixels et calculer leurs propriétés

### 5.1.2 Librairies

Les librairies utilisées dans ce chapitre sont les suivantes:

- SciPy

- NumPy

- opencv-python · PyPI

- scikit-image

- Rasterio

- Xarray

- rioxarray

Dans l'environnement Google Colab, seul rioxarray doit être installé:

1. `%%capture`

2. `!pip install -qU matplotlib rioxarray xrscipy scikit-image`
Vérifiez les importations:

1 import numpy as np

2 import numpy.fft

3 import rioxarray as rxr

4 from scipy import signal, ndimage

5 import xarray as xr

6 import xrscipy

7 import matplotlib.pyplot as plt

8 from skimage import data, measure, graph, segmentation, color

9 from skimage.color import rgb2gray

10 from skimage import slic, mark_boundaries

11 import pandas as pd

5.1.3 Images utilisées

Nous utilisons les images suivantes dans ce chapitre:

1 %%capture

2 import gdown

3

4 gdown.download()  
→ 'https://drive.google.com/uc?export=download&confirm=pbef&id=1a6Ypg0g1Oy4AJt9XWKWfnR12NW1XhNg_',  
→ output= 'RGBNIR_of_S2A.tif')

5 gdown.download()  
→ 'https://drive.google.com/uc?export=download&confirm=pbef&id=1a4PQ68Ru8zBphbQ22j0sgJ4D2quw-Wo6',  
→ output= 'landsat7.tif')

6 gdown.download()  
→ 'https://drive.google.com/uc?export=download&confirm=pbef&id=1_zwCLN-x7XJcNHJCH6Z8upEdUxtVtvs1',  
→ output= 'berkeley.jpg')

7 gdown.download()  
→ 'https://drive.google.com/uc?export=download&confirm=pbef&id=1dM6IVqjba6GHwTLmI7CpX8GP2z5txUq6',  
→ output= 'SAR.tif')

Vérifiez que vous êtes capable de les lire :

1 with rxr.open_rasterio('berkeley.jpg', mask_and_scale= True ) as imgRGB:  
2 print(imgRGB)

3 with rxr.open_rasterio('RGBNIR_of_S2A.tif', mask_and_scale= True ) as imgrgbnir:  
4 print(imgrgbnir)

5 with rxr.open_rasterio('SAR.tif', mask_and_scale= True ) as imgSAR:  
6 print(imgSAR)
## 5.2 Analyse fréquentielle

L'analyse fréquentielle, issue du traitement du signal, permet d'avoir un autre point de vue sur les données à partir de ses composantes harmoniques. La modification de ces composantes de Fourier modifie l'ensemble de l'image et permet de corriger des problèmes systématiques comme des artefacts ou du bruit de capteur. Bien que ce domaine soit un peu éloigné de la télédétection, les images issues des capteurs sont toutes sujettes à des étapes de traitement du signal et il faut donc en connaître les grands principes afin de pouvoir comprendre certains enjeux lors des traitements.

### 5.2.1 La transformée de Fourier

La transformée de Fourier permet de transformer une image dans un espace fréquentielle. Cette transformée est complètement réversible. Dans le cas des images numériques, on parle de 2D-DFT (2D-Discrete Fourier Transform) qui est un algorithme optimisé pour le calcul fréquentiel (Cooley et Tukey 1965). La 1D-DFT peu s'écrire simplement comme une projection sur une série d'exponentielles complexes:

$$X[k] = \sum_{n=0...N-1} x[n] \times \exp(-j \times 2\pi \times k \times n/N) \quad (5.1)$$

La transformée inverse prend une forme similaire:

$$x[k] = \frac{1}{N} \sum_{n=0...N-1} X[n] \times \exp(j \times 2\pi \times k \times n/N) \quad (5.2)$$

Le signal d'origine est donc reconstruit à partir d'une somme de sinusoides complexes $exp(j2\pi \frac{k}{N}n)$ de fréquence $k/N$. Noter qu'à partir de $k = N/2$, les sinusoides se répètent à un signe près et forme un miroir des composantes, la convention est alors de mettre ces composantes dans une espace négatif $[-N/2, ..., -1]$.

Dans le cas d'un simple signal périodique à une dimension avec une fréquence de $4/16$ (donc 4 périodes sur 16) on obtient deux pics de fréquence à la position de 4 cycles observés sur $N = 16$ observations. Les puissances de Fourier sont affichées dans un espace fréquentiel en cycles par unité d'espacement de l'échantillon (avec zéro au début) variant entre -1 et +1. Par exemple, si l'espacement des échantillons est en secondes, l'unité de fréquence est cycles/seconde (ou Hz). Dans le cas de N échantillons, le pic sera observé à la fréquence $+/- 4/16 = 0.25$ cycles/secondes. La fréquence d'échantillonnage $F_s$ du signal a aussi beaucoup d'importance aussi et doit être au moins a deux fois la plus haute fréquence observée (ici $F_s > 0.5$) sinon un phénomène de repliement appelé aliasing sera observé.

import math

Fs = 2.0

Ts = 1/Fs

N = 16

arr = xr.DataArray(np.sin(2*math.pi*np.arange(0,N,Ts)*4/16),
                    dims=('x'),
                    coord={ 'x': np.arange(0,N,Ts)})

fourier = np Bridge.fltcbridg (arr)

freq = np Bridge.fltcbridg(fourier.size, d = Ts)
7

axes[1].set_title("Décimée par un facteur 4")

8

plt.show()

<figure><img src="image_61.png" /></figure>

Une façon de réduire le contenu fréquentiel est de filtrer par un filtre passe-bas pour réduire les hautes fréquences
par exemple avec un filtre Gaussien:

1

from scipy.ndimage import gaussian_filter

2

q=4

3

sigma= q*1.1774/math.pi

4

arr = xr.DataArray(gaussian_filter(imgRGB.to numpy(), sigma= (0, sigma, sigma)), dims=( 'band','y','x'),
  → coordinates= {'x': imgRGB coors[ 'x'], 'y': imgRGB coors[ 'y'], 'spatial_ref': 0})

5

6

fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(10, 4))

7

plt.subplots(1, 2, 1)

8

img RGB.astype('int').plot.imshow

9

axes[0].set_title("Originale")

10

plt.subplots(1, 2, 2)

11

arr[: ,::q ,::q].astype('int').plot.imshow(

12

axes[1].set_title("Décimée par un facteur 4")

13

plt.show()
5 Transformations spatiales

<figure><img src="image_62.png" /></figure>

La fonction `decimate` dans `scipy.signal` réalise l'opération de décimation (downsampling) en une seule étape:

```python
import xrscipy.signal as dsp

fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(10, 4))
axes[0].set_title("Originale")
axes[1].set_title("Décimée par un facteur 4")

img RGB. astype('int').plot.imshow rhg="band")

axes[0].set_title("Originale")
axes[1].set_title("Décimée par un facteur 4")

dsp.decimate(imgوص, q=4, dim='x').astype('int').plot.imshow rhg="band")
```

Text (0.5, 1.0, 'Décimée par un facteur 4')
<figure><img src="image_63.png" /></figure>

## 5.3 Filtrage d'image

Le filtrage d'image a plusieurs objectifs en télédétection:

1. La réduction du bruit afin d'améliorer la résolution radiométrique et améliorer la lisibilité de l'image.

2. Le réhaussement de l'image afin d'améliorer le contraste ou faire ressortir les contours.

3. La production de nouvelles caractéristiques, c.-à.-d dérivées de nouvelles images mettant en valeur certaines informations dans l'image comme la texture, les contours, etc.

Il existe de nombreuses méthodes de filtrage dans la littérature qui sont habituellement regroupées en quatre catégories:

1. Le filtrage peut-être global ou local, c.-à.-d qu'il prend en compte soit toute l'image pour filtrer (ex: filtrage par Fourier), soit uniquement avec une fenêtre ou un voisinage local.

2. La fonction de filtrage peut-être linéaire ou non linéaire.

3. La fonction de filtrage peut être stationnaire ou adaptative.

4. Le filtrage peut-être mono-échelle ou multi-échelle.

La librairie Scipy (Multidimensional image processing (scipy.ndimage)) contient une panoplie complète de filtres.
### 5.3.1 Filtrage linéaire stationnaire

Un filtrage linéaire stationnaire consiste à appliquer une même pondération locale des valeurs des pixels dans une fenêtre glissante. La taille de cette fenêtre est généralement un chiffre impair (3,5, etc.) afin de définir une position centrale et une fenêtre symétrique. La valeur calculée à partir de tous les pixels dans la fenêtre est alors attribuée au pixel central.

Le filtre le plus simple est certainement le filtre moyen qui consiste à appliquer le même poids uniforme dans la fenêtre glissante. Par exemple pour un filtre 5x5:

$$F = \frac{1}{25} \begin{bmatrix} 1 & 1 & 1 & 1 & 1 \\ 1 & 1 & 1 & 1 & 1 \\ 1 & 1 & 1 & 1 & 1 \\ 1 & 1 & 1 & 1 & 1 \\ 1 & 1 & 1 & 1 & 1 \end{bmatrix} \quad (5.5)$$

En python, on dispose des fonctions `rolling` et `sliding_windows` définis dans la librairie ` numpy`. Par exemple pour le cas du filtre moyen, on construit une nouvelle vue de l'image avec deux nouvelles dimensions `x_win` et `y_win`:

```python
1 rolling_win = img rgb.rolling(x=5, y=5, min periods=3, center= True).construct(x="x_win", y="y_win",
    → keep_attrs= True)
```

```python
2 print(rolling_win[0,0,1,...])
```

```python
3 print(rolling_win.shape)
```

```python
<xarray.DataArray (x_win: 5, y_win: 5)> Size: 100B
array([[ nan, nan, nan, nan, nan],
    [ nan, nan, 209., 210., 209.],
    [ nan, nan, 213., 214., 212.],
    [ nan, nan, 213., 212., 210.],
    [ nan, nan, 210., 209., 206.]], dtype= float32)
```

Coordinates:

<table><tr><td>band</td><td>int64 8B 1</td></tr><tr><td>x</td><td>float64 8B 1.5</td></tr><tr><td>y</td><td>float64 8B 0.5</td></tr><tr><td>spatial_ref</td><td>int64 8B 0</td></tr></table>

Dimensions without coordinates: x_win, y_win
(3, 771, 1311, 5, 5)

L'avantage de cette approche est qu'il n'y a pas d'utilisation inutile de la mémoire. Noter les `nan` sur les bords de l'image car la fenêtre déborde sur les bordures de l'image. Par la suite un opérateur de moyenne peut être appliqué sur les axes `x_win` et `y_win` correspondant aux fenêtres glissantes.

```python
1 filtre_moyen= rolling_win.mean(dim= ['x_win', 'y_win'], skipna= True)
```

```python
2 fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(8, 4))
```

```python
3 filtre_moyen.astype('int').plot.imshow rhg="band")
```

```python
4 ax.set_title("Filtre moyen 5x5")
```

```python
Text(0.5, 1.0, 'Filtre moyen 5x5')
```
<figure><img src="image_64.png" /></figure>

Lorsque la taille *W* de la fenêtre devient trop grande, il est préférable d'utiliser une convolution dans le domaine fréquentielle. La fonction `fftconvolve` de la librairie `scipy.signal` offre cette possibilité:

1. `kernel` = np_outer(signal.windows.gaussian(70, 8), signal.windows.gaussian(70, 8))
2. `blurred` = signal.fftconvolve(img imp, kernel, mode='same')

### 5.3.1.1 Filtrage par convolution

La façon la plus efficace d'appliquer un filtre linéaire est d'appliquer une convolution. La convolution est généralement très efficace car elle est peut être calculée dans le domaine fréquentiel. Prenons l'exemple du filtre de Scharr (Jahne et S. 1999), qui permet de détecter les contours horizontaux et verticaux:

$$F = \begin{bmatrix} -3 - 3j & 0 - 10j & +3 - 3j \\ -10 + 0j & 0 + 0j & +10 + 0j \\ -3 + 3j & 0 + 10j & +3 + 3j \end{bmatrix} \qquad (5.6)$$

Remarquez l'utilisation de chiffres complexes afin de passer deux filtres différents sur la partie réelle et imaginaire.

1. `scharr` = np_array([[-3-3j, 0-10j, +3-3j], [-10+0j, 0+0j, +10+0j], [-3+3j, 0+10j, +3+3j]]) # Gx + j*Gy
```python
print(img.isel(band=0).shape)
```

```python
grad = signal(convolve2d(img rgb.isel(band=0), scharr, boundary='symm', mode='same'))
```

```python
# on reconstruit un xarray à partir du résultat:
```

```python
arr = xr.DataArray(np.abs在这里), dims=("y", "x"), coord = {'x': img rgb. coordinates["x"], 'y': img rgb. coordinates["y"], 'spatial_ref': 0})
```

```python
print(arr)
```

```python
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(8, 4))
```

```python
arr. plot. imshow()
```

```python
ax. set_title("Amplitude du filtre de Scharr")
```

```python
(771, 1311)
<xarray.DataArray (y: 771, x: 1311)> Size: 8MB
array([[ 65.96969001, 58.8557588, 54.91812087, ..., 1474. ,
     1037.01205393, 389.99487176],
     [ 61.07372594, 39.8246155, 89.18520057, ..., 1763.79647352,
     864.92543031, 270.20362692],
     [ 98.48857802, 112.44554237, 168.10710871, ..., 2110.61365484,
     870.36658943, 204.40156555],
     ...,
     [ 143.17821063, 597.00753764, 2479.42977315, ..., 216.00925906,
     248.33847869, 200.89798406],
     [ 106.07544485, 393.67245268, 2188.78824924, ..., 124.96399481,
     159.90622252, 346.34087255],
     [ 41.59326869, 229.05894438, 1845.1216762 , ..., 175.16278143,
     33.37663854, 414.3911196 ]])

Coordinates:
* x       (x) float64 10kB 0.5 1.5 2.5 ... 1.308e+03 1.31e+03 1.31e+03
* y       (y) float64 6kB 0.5 1.5 2.5 3.5 4.5 ... 767.5 768.5 769.5 770.5
spatial_ref int64 8B 0
```

```python
Text(0.5, 1.0, 'Amplitude du filtre de Scharr')
```

<figure><img src="image_65.png" /></figure>
## 5.4 Gestion des bordures

L'application de filtres à l'intérieur de fenêtres glissantes implique de gérer les bords de l'image, car la fenêtre de traitement va nécessairement déborder de quelques pixels en dehors de l'image (généralement la moitié de la fenêtre déborde). On peut soit décider d'ignorer les valeurs en dehors de l'image en imposant une valeur nan, soit prolonger l'image de quelques lignes et colonnes avec des valeurs miroirs ou constantes.

### 5.4.0.1 Filtrage par une couche convolutionnelle

#### Installation de Pytorch

Cette section nécessite la librairie Pytorch avec un GPU et ne fonctionnera que sur Colab. On peut quand même installer une version locale CPU de pytorch: `pip install -qU torch==2.4.0+cpu`

Une couche convolutionnelle est simplement un ensemble de filtres appliqués sur la donnée d'entrée. Ce type de filtrage est à la base des réseaux dits convolutionnels qui seront abordés dans le tome 2. On peut ici imposer les mêmes filtres de gradient dans la couche convolutionnelle :

```python
import torch
import torch.nn as nn
import numpy as np
import matplotlib.pyplot as plt
normalizedImg = torch.tensor(imgRGB.to numpy())
nchannels = normalizedImg.size()[0] # nombre de canaux de l'image

# On forme une couche convolutionnelle
convLayer = nn.Conv2d(inChannels=ncannels, outChannels=2, kernelSize=3, padding=1, stride=1,
    -> dilation=1)

# Filtre de Sobel
sobel_x = np.array([-3, 0, 3], [-10, 0, 10], [-3, 0, 3]))
sobel_y = np.array([-3, -10, -3], [0, 0, 0], [3, 10, 3]))

# Le filtre (kernel) est formé de deux filtres
kernel = npstack([sobel_x, sobel_y])
kernel = kernel.reshape(2, 1, 3, 3)

# On répète le filtre pour chaque bande
kernel = np tile kernel, (1, nchannels, 1, 1))
print (kernel.shape)
kernel = torch.as tensor (kernel, type = torch.float32)
convLayer weight = nn.Parameter (kernel)
convLayer.bias = nn.Parameter (torch zeros (2,))

input = normalizedImg unsqueeze (0) # il faut ajouter une dimension pour le nombre d'échantillons
print (input.shape)
# Visualize the filters
fig, axes = plt.subplots (1, 2, figsize = (8, 5))
for i in range (2):
```
```python
axs[i].imshow(conv_layers()[i, 0])
axs[i].set_title(f'Filtre {i+1}')
plt.show()
```

```
(2, 3, 3, 3)
torch Size([1, 3, 771, 1311])
```

<figure><img src="image_66.png" /></figure>

Le résultat est alors calculé sur GPU (si disponible):

```python
import torch
import matplotlib.pyplot as plt

output = conv_layers(input)
print(f'Image (BxCxHxW): {input.shape}')
print(f'Sortie (BxFxHxW): {output.shape}')

fig, axes = plt.subplots(1, 2, figsize=(20, 5))
for i in range(2):
    axes[i].imshow(output Detach().data numpy()[0,i], vmin=-5000, vmax=5000, cmap='gray')
    axes[i].set_title(f'Filtrage {i+1}')

plt.show()
```

```
Image (BxCxHxW): torch Size([1, 3, 771, 1311])
Sortie (BxFxHxW): torch Size([1, 2, 771, 1311])
```
<figure><img src="image_67.png" /></figure>

### 5.4.1 Filtrage adaptatif

Les filtrages adaptatifs consistent à appliquer un traitement en fonction du contenu local d'une image. Le filtre n'est alors plus stationnaire et sa réponse peut varier en fonction du contenu local. Ce type de filtre est très utilisé pour filtrer les images SAR (Synthetic Aperture Radar) qui sont dégradées par un bruit multiplicatif que l'on appelle *speckle*. On peut voir un exemple d'une image Sentinel-1 (bande HH) sur la région de Montréal, remarquez que l'image est affichée en dB en appliquant la fonction `log10`.

1. `print(img_SAR.rio resolutions())`

2. `print(img_SAR.rio.crs)`

3. `fig, ax = plt.subplots(1, 1, figsize=(6, 4))`

4. `xr.ufuncs.log10(img_SAR sel(band=1).drop("band")).plot()`

5. `axs.set_title("Image SAR Sentinel-1 (dB)")`

```
(0.00029254428869762705, -0.000287092818453516)
EPSG:4326
```

```
Text(0.5, 1.0, 'Image SAR Sentinel-1 (dB)')
```
5 Transformations spatiales

<figure><img src="image_68.png" /></figure>

Un des filtres les plus simples pour réduire le bruit est d'appliquer un filtre moyen, par exemple un 5 × 5
ci-dessous:

1
`rolling_win = img_SAR sel(band=2).rolling(x=5, y=5, min periods=3, center= True).construct(x="x_win",
  → y="y_win", keep_attrs= True)`

2
`filtre_moyen= rolling_win .mean(dim= ['x_win', 'y_win'], skipna= True)`

3
`fig, axes = plt.subplots(1, 1, figsize=(6, 4))`

4
`xr.ufuncs.log10(filtre_moyen).plot.imshow()`

5
`axes.set_title("Filtrage moyen 5x5 (dB)")`

Text(0.5, 1.0, 'Filtrage moyen 5x5 (dB)')
<figure><img src="image_69.png" /></figure>

Au lieu d'appliquer un filtre moyen de manière indiscriminée, le filtre de Lee (Lee 1986) applique une pondération en fonction du contenu local de l'image *I* dans sa forme la plus simple :

$$ 
\begin{aligned}
I_F &= I_M + K \times (I - I_M) \\
K &= \frac{\sigma_I^2}{\sigma_I^2 + \sigma_{bruit}^2}
\end{aligned}
\qquad (5.7) 
$$

De la sorte, si la variance locale est élevée *K* s'approche de 1 préservant ainsi les détails de l'image *I* sinon l'image moyenne *I*<sub>M</sub> est appliquée.

rolling_win = img_SAR sel(band=2).rolling(x=5, y=5, min periods=3, center= True).construct(x="x_win",
  y="y_win", keep_attrs= True)

filtre_moyen = rolling_win .mean(dim= ['x_win', 'y_win'], skipna= True)

ecart_type= rolling_win .std(dim= ['x_win', 'y_win'], skipna= True)

cv=ecart_type/filtre_moyen

ponderation = (cv - 0.25) / cv

fig, axes = plt .subplots(nrows=1, ncols=2, figsize=(10, 4), sharex = True, sharey = True)

plt .subplot(1, 2, 1)

cv .plot .imshow( vmin=0, vmax=2)

axes [0].set_title("CV")

plt .subplot(1, 2, 2)

ponderation .plot .imshow( vmin=0, vmax=1)

axes [1].set_title("Pondération")

plt .tight_layout()
# 5 Transformations spatiales

<figure><img src="image_70.png" /></figure>

On zoomant sur l'image, on voit clairement que les détails de l'image sont mieux préservés :

<figure><img src="image_71.png" /></figure>

## 5.5 Segmentation

La segmentation d'image consiste à séparer une image en régions homogènes spatialement connexes (segments) où les valeurs sont uniformes selon un certain critère (couleurs, texture, etc.). Une image présente généralement beaucoup de pixels redondants, l'intérêt de ce type de méthode est essentiellement de réduire la quantité de pixels nécessaire. En télédétection, on parle souvent d'approche objet. En vision par ordinateur, on parle parfois de super-pixel. Il existe de nombreuses méthodes de segmentation, la librairie **sickit-image** rend disponible plusieurs implémentations sur des images RVB (Comparison of segmentation and superpixel algorithms — skimage 0.25.0 documentation).
5.5.1 Super-pixel

Ce type de méthode cherche à former des régions homogènes et compactes dans l'image (Achanta et Süsstrunk 2012). Une des méthodes les plus simples est la méthode SLIC (Simple Linear Iterative Clustering), elle combine un regroupement de type k-moyennes avec une distance hybride qui prend en compte les différences de couleur entre pixels mais aussi leur distance par rapport centre du super-pixel:

1. Décomposer l'image en N régions régulières de taille S × S

2. Initialiser les centres Ck de chaque segment k

3. Rechercher les pixels ayant la distance la plus petite dans une région 2S × 2S :

$$
D_{SLIC} = d_{couleur} + \frac{m}{S} d_{xy}
$$

4. Mettre à jour les centre Ck de chaque segment k et réitérer à l'étape 3.

Les régions évoluent rapidement avec les itérations, plus le poids m est élevé, plus la forme du super-pixel est contrainte et ne suivra pas vraiment le contenu de l'image:

img = imgRGB.to numpy().astype('uint8').transpose(1,2,0)

segments_slic1 = slic(img, n_segments=250, compactness=10, sigma=1, start_label=1, max_num_iter=1)

segments_slic2 = slic(img, n_segments=250, compactness=10, sigma=1, start_label=1, max_num_iter=2)

segments_slic100 = slic(img, n_segments=250, compactness=100, sigma=1, start_label=1, max_num_iter=10)

segments_slic100b = slic(img, n_segments=250, compactness=10, sigma=1, start_label=1, max_num_iter=10)

print(f'SLIC nombre de segments: {len(np.unique(segments_slic1))}')

fig, ax = plt.subplots(2, 2, figsize=(10, 6), sharex=True, sharey=True)

ax[0, 0].imshow/mark_boundaries(img, segments_slic1))

ax[0, 0].set_title("Initialisation")

ax[0, 1].imshow/mark_boundaries(img, segments_slic2))

ax[0, 1].set_title('2 itérations')

ax[1, 0].imshow/mark_boundaries(img, segments_slic100))

ax[1, 0].set_title('10 itérations avec m=100')

ax[1, 1].imshow/mark_boundaries(img, segments_slic100b))

ax[1, 1].set_title('10 itérations avec m=10')

for a in aximsel():

a.set_axis_off()

plt.tight_layout()

pht.show()

SLIC nombre de segments: 240
5 Transformations spatiales

Initialisation

<figure><img src="image_72.png" /></figure>

10 itérations avec m=100

2 itérations

<figure><img src="image_73.png" /></figure>

10 itérations avec m=10

<figure><img src="image_74.png" /></figure>

Le nombre de segments initial est probablement le paramètre le plus important. Une manière de l'estimer est d'évaluer l'échelle moyenne des segments homogènes dans l'image à analyser. On observe ci-dessous l'impact de passer d'une échelle 40 x 40 à 20 x 20. En prenant la moyenne de chaque segment, on constate que l'échelle 40 x 40 génère des segments trop grands mélangeant plusieurs classes.

```python
1 from skimage import color, segmentation
2 nRegions = int((img.shape[0] * img.shape[1])/(40*40))
3 print('Nb segments: ', nRegions)
4 segments_slic_40 = slic(img, n_segments=nRegions, compactness=10, sigma=1, start_label=1,
5                       max_num_iter=10)
5 print(f'SLIC nombre de segments: {len(np.unique(segments_slic_40))}')
6 out = color.label2rgb(segments_slic_40, img, kind='avg', bg_label=0)
7 out_40 = segmentation.mark_boundaries(out, segments_slic_40, (0, 0, 0))
8
9 nRegions = int((img.shape[0] * img.shape[1])/(20*20))
10 print('Nb segments: ', nRegions)
11 segments_slic_20 = slic(img, n_segments=nRegions, compactness=10, sigma=1, start_label=1,
12                       max_num_iter=10)
13 print(f'SLIC nombre de segments: {len(np.unique(segments_slic_20))}')
14 out = color.label2rgb(segments_slic_20, img, kind='avg', bg_label=0)
15 out_20 = segmentation.mark_boundaries(out, segments_slic_20, (0, 0, 0))
```
5 Transformations spatiales

16 fig, ax = plt.subplots(2, 1, figsize=(6, 8), sharex=True, sharey=True)
17
18 ax[0]..imshow(out_40)
19 ax[0].set_title("Initialisation avec 631 segments")
20 ax[1]..imshow(out_20)
21 ax[1].set_title('Initialisation avec 2526 segments')
22
23 for a in ax.ravel():
24     a.set_axis_off()
25
26 plt.tight_layout()
27
28 plt.show()

Nb segments: 631
SLIC nombre de segments: 459
Nb segments: 2526
SLIC nombre de segments: 2201
<figure><img src="image_75.png" /><figcaption>## Initialisation avec 2526 segments</figcaption></figure>

<figure><img src="image_76.png" /></figure>

### 5.5.2 Fusion des segments par graphe de proximité

Une segmentation peut produire beaucoup trop de segments. On parle alors de sur-segmentation. Ceci est recherché dans certains cas pour permettre de bien capturer les détails fins de l'image. Cependant, afin de réduire
le nombre de segments, un post-traitement possible est de fusionner les segments similaires selon certaines règles ou distances. Un graphe d'adjacence de régions (figure 5.1) est formé à partir des segments connectés où chaque nœud représente un segment et un lien de proximité (Jaworek-Korjakowska (2018)). À partir de ce graphe, on peut fusionner les nœuds similaires à partir de leur distance radiométrique.

<figure><img src="image_77.png" /><figcaption>FIGURE 5.1 – Graphe d'adjacence de régions, d'après (Jaworek-Korjakowska (2018)). Chaque nœud est un segment, un lien est formé uniquement si les segments se touchent (par exemple le segment 6 ne touche que la région 5). La fonction `graph.rag_mean_color` produit un graphe à partir d'une segmentation et de l'image originale. Chaque nœud tient la couleur de chaque segment dans un attribut appelé `'mean color'`.</figcaption></figure>

def _weight_mean_color(graph, src, dst, n):

```python
        """Fonction pour gérer la fusion des nœuds en recalculant la couleur moyenne.
```

```python
        La méthode suppose que la couleur moyenne de `dst` est déjà calculée.
```

```python
        """
```

```python
diff = graph.nodes[dsthv['mean color'] - graph.nodes[n][['mean color']]
```

```python
diff = np laidg .norm(diff)
```

```python
print(diff)
```

```python
return {'weight': diff}
```

def merge_mean_color(graph, src, dst):

```python
        """Fonction appelée avant la fusion de deux nœuds d'un graphe de distance de couleur moyenne.
```

```python
        Cette méthode calcule la couleur moyenne de `dst`.
```

```python
        """
```

```python
graph.nodes[dsthv['total color'] += graph.nodes(src)[['total color']]
```

```python
graph.nodes[dsthv['pixel count'] += graph.nodes(src)[['pixel count']]
```

```python
graph.nodes[dsthv['mean color']] = (
    graph.nodes[dsthv['total color'] / graph.nodes[dsthv['pixel count']])
```

```python
)
```
```python
g = graph.rag_mean_color(img, segments_slic_20)
print('Nombre de segments:', len(g))
labels2 = graph.merge_hierarchical(
    segments_slic_20,
    g,
    thresh=20,
    rag_copy=False,
    in_place_merge=True,
    merge_func=merge_mean_color,
    weight_func=_weight_mean_color,
)
print('Nombre de segments:', len(g))
out1 = color.label2rgb(segments_slic_20, img, kind='avg', bg_label=0)
out1 = segmentation.mark_boundaries(out1, segments_slic_20, (0, 0, 0))
out2 = color.label2rgb labels2, img, kind='avg', bg_label=0)
out2 = segmentation.mark_boundaries(out2, labels2, (0, 0, 0))
fig, ax = plt.subplots(nrows=2, sharex=True, sharey=True, figsize=(6, 8))
ax[0]..imshow(out1)
ax[0].set_title("Avant fusion")
ax[1]..imshow(out2)
ax[1].set_title("Après fusion")
for a in ax:
    a.axis('off')
plt.tight_layout()
```

Nombre de segments: 2201

Nombre de segments: 1187
<figure><img src="image_78.png" /><figcaption>Après fusion</figcaption></figure>

<figure><img src="image_79.png" /></figure>

### 5.5.3 Approche objet

L'approche objet consiste à traiter chaque segment comme un objet avec un ensemble de propriétés. La librairie `skimage` offre la possibilité d'enrichir chaque segment avec des **propriétés** et de former un tableau:
5 Transformations spatiales

1 properties = ['label', 'area', 'centroid', 'num pixels', 'intensity_mean', 'intensity_std']

2
3 table= measure.regionprops_table labels2, intensity_image= img rgb.to numpy().transpose 1,2,0),
→ properties=properties)

5 table = pd.DataFrame(table)

6 table head(10)

<table>
  <thead>
    <tr>
      <th>label</th>
      <th>area</th>
      <th>centroid-0</th>
      <th>centroid-1</th>
      <th>num pixels</th>
      <th>intensity_mean-0</th>
      <th>intensity_mean-1</th>
      <th>intensity_mean-2</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td>0</td>
      <td>1</td>
      <td>641.0</td>
      <td>15.466459</td>
      <td>69.489860</td>
      <td>641</td>
      <td>136.730103</td>
      <td>132.851791</td>
      <td>117.126366</td>
    </tr>
    <tr>
      <td>1</td>
      <td>2</td>
      <td>480.0</td>
      <td>10.997917</td>
      <td>92.614583</td>
      <td>480</td>
      <td>201.208328</td>
      <td>198.262497</td>
      <td>188.483337</td>
    </tr>
    <tr>
      <td>2</td>
      <td>3</td>
      <td>712.0</td>
      <td>16.683989</td>
      <td>114.776685</td>
      <td>712</td>
      <td>185.349716</td>
      <td>183.113770</td>
      <td>170.994385</td>
    </tr>
    <tr>
      <td>3</td>
      <td>4</td>
      <td>1803.0</td>
      <td>31.974487</td>
      <td>139.379368</td>
      <td>1803</td>
      <td>117.897392</td>
      <td>108.367722</td>
      <td>97.769829</td>
    </tr>
    <tr>
      <td>4</td>
      <td>5</td>
      <td>448.0</td>
      <td>5.004464</td>
      <td>166.542411</td>
      <td>448</td>
      <td>183.511154</td>
      <td>181.276779</td>
      <td>167.720978</td>
    </tr>
    <tr>
      <td>5</td>
      <td>6</td>
      <td>459.0</td>
      <td>9.934641</td>
      <td>191.668845</td>
      <td>459</td>
      <td>133.557739</td>
      <td>133.821350</td>
      <td>129.697174</td>
    </tr>
    <tr>
      <td>6</td>
      <td>7</td>
      <td>355.0</td>
      <td>5.160563</td>
      <td>222.895775</td>
      <td>355</td>
      <td>148.574646</td>
      <td>148.802811</td>
      <td>142.580276</td>
    </tr>
    <tr>
      <td>7</td>
      <td>8</td>
      <td>334.0</td>
      <td>4.904192</td>
      <td>255.904192</td>
      <td>334</td>
      <td>125.973053</td>
      <td>121.197601</td>
      <td>108.973053</td>
    </tr>
    <tr>
      <td>8</td>
      <td>9</td>
      <td>1481.0</td>
      <td>32.279541</td>
      <td>292.865631</td>
      <td>1481</td>
      <td>204.102631</td>
      <td>172.359894</td>
      <td>137.501007</td>
    </tr>
    <tr>
      <td>9</td>
      <td>10</td>
      <td>445.0</td>
      <td>8.013483</td>
      <td>308.053933</td>
      <td>445</td>
      <td>145.373032</td>
      <td>138.182022</td>
      <td>121.402245</td>
    </tr>
  </tbody>
</table>

Ce tableau pourra être exploiter pour une tâche de classification par la suite (on parle alors de classification
objet).
# Partie 3. Classifications d'images
# 6 Classifications d'images supervisées

## 6.1 Préambule

Assurez-vous de lire ce préambule avant d'exécuter le reste du notebook.

### 6.1.1 Objectifs

Dans ce chapitre, nous ferons une introduction générale à l'apprentissage automatique et abordons quelques techniques fondamentales. La librairie centrale utilisée dans ce chapitre sera `sickit-learn`. Ce chapitre est aussi disponible sous la forme d'un notebook Python sur Google Colab:

<figure><img src="image_80.png" /></figure>

<figure><img src="image_81.png" /></figure>

**Objectif**

**Objectifs d'apprentissage visés dans ce chapitre**

À la fin de ce chapitre, vous devriez être en mesure de :

- comprendre les principes de l'apprentissage automatique supervisé;

- mettre en place un pipeline d'entraînement;

- savoir comment évaluer les résultats d'un classificateur;

- visualiser les frontières de décision;

- mettre en place des techniques de classifications comme K-NN et les arbres de décision;

### 6.1.2 Librairies

Les librairies utilisées dans ce chapitre sont les suivantes :

- SciPy

- NumPy

- opencv-python · PyPI

- scikit-image

- Rasterio

- xarray

- rioxarray

- geopandas

- scikit-learn
Dans l'environnement Google Colab, seul rioxarray et xrscipy sont installés:

1
```python
%%capture
2 !pip install -qU matplotlib rioxarray xrscipy supertree
```

Vérifiez les importations nécessaires en premier:

1
```python
import numpy as np
```

### 6.1.3 Images utilisées

Nous utilisons les images suivantes dans ce chapitre:

1
```python
%%capture
```

2
```python
import gdown
```

3
```python
gdown.download()
```

4
→ 'https://drive.google.com/uc?export=download&confirm=pbef&id=1a6Ypg0g10y4AJt9XWKWfnR12NW1XhNg_',
→ output= 'RGBNIR_of_S2A.tif')

5
```python
gdown.download()
```

6
→ 'https://drive.google.com/uc?export=download&confirm=pbef&id=1a4PQ68Ru8zBphbQ22j0sgJ4D2quw-Wo6',
→ output= 'landsat7.tif')

6
```python
gdown.download()
```

7
→ 'https://drive.google.com/uc?export=download&confirm=pbef&id=1_zwCLN-x7XJcNHJCH6Z8upEdUxtVtvs1',
→ output= 'berkeley.jpg')

7
```python
gdown download()
```

→ 'https://drive.google.com/uc?export=download&confirm=pbef&id=1dM6IVqjba6GHwTLmI7CpX8GP2z5txUq6',
→ output= 'SAR.tif')
8

gdown.download()
→ 'https://drive.google.com/uc?export=download&confirm=pbef&id=1aAq7crc_LoaLC3kG3HkQ6Fv5JfG0mswg',
→ output= 'carte.tif')

Vérifiez que vous êtes capable de les lire :

1

2

3

4

5

6

7

8

with rxr.open_rasterio('berkeley.jpg', mask_and_scale= True) as img rgb:
    print(img rgb)
with rxr.open_rasterio('RGBNIR_of_S2A.tif', mask_and_scale= True) as img rgbnir:
    print(img rgbnir)
with rxr.open_rasterio('SAR.tif', mask_and_scale= True) as img SAR:
    print(img SAR)
with rxr.open_rasterio('carte.tif', mask_and_scale= True) as img carte:
    print(img carte)

## 6.2 Principes généraux

Une classification supervisée ou dirigée consiste à attribuer une étiquette (une classe) de manière automatique à chaque point d'un jeu de données. Cette classification peut se faire à l'aide d'une cascade de règles pré-établies (arbre de décision) ou à l'aide de techniques d'apprentissage automatique (machine learning). L'utilisation de règles pré-établies atteint vite une limite car ces règles doivent être fournies manuellement par un expert. Ainsi, l'avantage de l'apprentissage automatique est que les règles de décision sont dérivées automatiquement du jeu de données via une phase dite d'entraînement. On parle souvent de solutions générées par les données (Data Driven Solutions). Cet ensemble de règles est souvent appelé **modèle**. On visualise souvent ces règles sous la forme de frontières de décisions dans l'espace des données. Cependant, un des défis majeurs de ce type de technique est d'être capable de produire des règles qui soient généralisables au-delà du jeu d'entraînement.

Les classifications supervisées ou dirigées présupposent donc que nous avons à disposition un jeu d'entraînement déjà étiqueté. Celui-ci va nous permettre de construire un modèle. Afin que ce modèle soit représentatif et robuste, il nous faut assez de données d'entraînement. Les algorithmes d'apprentissage automatique sont très nombreux et plus ou moins complexes pouvant produire des frontières de décision très complexes et non linéaires.

### 6.2.1 Comportement d'un modèle

Cet exemple tiré de `sickit-learn` illustre les problèmes d'ajustement insuffisant ou **sous-apprentissage** (*underfitting*) et d'ajustement excessif ou **sur-apprentissage** (*overfitting*) et montre comment nous pouvons utiliser la régression linéaire avec un modèle polynomiale pour approximer des fonctions non linéaires. La figure 6.1 montre la fonction que nous voulons approximer, qui est une partie de la fonction cosinus (couleur orange). En outre, les échantillons de la fonction réelle et les approximations de différents modèles sont affichés en bleu. Les modèles ont des caractéristiques polynomiales de différents degrés. Nous pouvons constater qu'une fonction linéaire (polynôme de degré 1) n'est pas suffisante pour s'adapter aux échantillons d'apprentissage. C'est ce qu'on appelle un sous-ajustement (*underfitting*) qui produit un biais systématique quels que soient les points d'entraînement. Un polynôme de degré 4 se rapproche presque parfaitement de la fonction réelle. Cependant, pour des degrés plus élevés, le modèle s'adaptera trop aux données d'apprentissage, c'est-à-dire
qu'il apprendra le bruit des données d'apprentissage. Nous évaluons quantitativement le sur-apprentissage et le sous-apprentissage à l'aide de la validation croisée. Nous calculons l'erreur quadratique moyenne (EQM) sur l'ensemble de validation. Plus elle est élevée, moins le modèle est susceptible de se généraliser correctement à partir des données d'apprentissage.

<figure><img src="image_82.png" /><figcaption>FIGURE 6.1 – Exemples de sur et sous-apprentissage.</figcaption></figure>

On constate aussi que sans les échantillons de validation, nous serions incapables de déterminer la situation de sur-apprentissage, l'erreur sur les points d'entraînement seuls étant excellente pour un degré 15.

### 6.2.2 Pipeline

La construction d'un modèle implique généralement toujours les mêmes étapes illustrées sur la figure 6.2:

1. La préparation des données implique parfois un pré-traitement afin de normaliser les données.

2. Partage des données en trois groupes: entraînement, validation et test.

3. L'apprentissage du modèle sur l'ensemble d'entraînement. Cet apprentissage nécessite de déterminer les valeurs des hyper-paramètres du modèle par l'usager.

4. La validation du modèle sur l'ensemble de validation. Cette étape vise à vérifier que les hyper-paramètres du modèle sont adéquats.

5. Enfin le test du modèle sur un ensemble de données indépendant.

### 6.2.3 Construction d'un ensemble d'entraînement

Les données d'entraînement permettent de construire un modèle. Elles peuvent prendre des formes très variées mais on peut voir cela sous la forme d'un tableau $N \times D$:

1. La taille *N* du jeu de données.
<figure><img src="image_83.png" /><figcaption>FIGURE 6.2 – Étapes standards dans un entraînement.</figcaption></figure>
# 6 Classifications d'images supervisées

2. Chaque entrée définit un échantillon ou un point dans un espace à plusieurs dimensions.

3. Chaque échantillon est décrit par *D* dimensions ou caractéristiques (*features*).

Une façon simple de construire un ensemble d'entraînement est d'échantillonner un produit existant. Nous
utilisons une carte d'occupation des sols qui contient 12 classes différentes.

```python
couleurs_classes = {'NoData': 'black', 'Commercial': 'yellow', 'Nuages': 'lightgrey',
                  'Foret': 'darkgreen', 'Faible_végétation': 'green', 'Sol_nu': 'saddlebrown',
                  'Roche': 'dimgray', 'Route': 'red', 'Urbain': 'orange', 'Eau': 'blue', 'Tourbe':
                  'salmon', 'Végétation éparse': 'darkgoldenrod', 'Roche avec végétation':
                  'darkseagreen'}

nom_classes = [*couleurs_classes.keys()]
couleurs_classes = [*couleurs_classes.values()]
```

On peut visualiser la carte de la façon suivante :

```python
import matplotlib.pyplot as plt
import rioxarray as rxr
cmap_classes = ListedColormap(couleurs_classes)

fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(8, 6))
imgCarte.squeeze().plot.imshow(cmap=cmap_classes, vmin=0, vmax=12)
ax.set_title("Carte d'occupation des sols", fontsize="small")
```

Text(0.5, 1.0, "Carte d'occupation des sols")
## 6 Classifications d'images supervisées

<figure><img src="image_84.png" /></figure>

On peut facilement calculer la fréquence d'occurrences des 12 classes dans l'image à l'aide de numpy:

1. `imgCarte = imgCarte.squeeze()` # nécessaire pour ignorer la dimension du canal

2. `compteClasse = np.unique(imgCarte.data, return_counts=True)`

3. `print(compteClasse)`

```
(array([ 1.,  3.,  4.,  5.,  6.,  7.,  8.,  9., 11., 12., nan],
       dtype浮点型), array([193558, 2104777, 670158, 29523, 14624, 94751, 750046,
       123671,   9079,   4327,    10]))
```

La fréquence d'apparition de chaque classe varie grandement, on parle alors d'un **ensemble déséquilibré**. Ceci est très commun dans la plupart des ensembles d'entraînement, puisque les classes ont très rarement la même fréquence. Par exemple, à la lecture du graphique en barres verticales, on constate que la classe forêt est très présentes contrairement à plusieurs autres classes (notamment, tourbe, végétation éparse, roche, sol nu et nuages).

1. `valeurs, comptes = compteClasse`

2.
```python
3 # Create the histogram
4 plt.figure(figsize=(5, 3))
5 plt.bar(valeurs, comptes/comptes.sum()*100)
6 plt.xlabel("Classes")
7 plt.ylabel("%")
8 plt.title("Fréquences des classes", fontsize="small")
9 plt.xticks(range(len(nom_classes)), nom_classes, rotation=45, ha='right')
10 plt.show()
```

<figure><img src="image_85.png" /></figure>

On peut échantillonner aléatoirement 100 points pour chaque classe:

```python
1 imgCarte = imgCarte.squeeze()
2 class_counts = np.unique(imgCarte.data, return_counts=True)
3
4 # Liste vide des points échantillonnées
5 sampled_points = []
6 class labels= [] # contient les étiquettes des classes
7 for class_label in range(1, 13): # pour chacune des 12 classes
8     # On cherche tous les pixels pour cette étiquette
9     class pixels = np.argwhere(imgCarte.data == class_label)
10
11 # On se limite à 100 pixels par classe
12 nSamples = min(100, len(class pixels))
```
13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

```python
# On les choisit les positions aléatoirement
np.random.seed(0)  # ceci permet de répliquer le tirage aléatoire
sampled_indices = np.random.choice(len(class pixels), n Samples, replace=False)
```

```python
# On prends les positions en lignes, colonnes
sampled pixels = class pixels[sampled indices]
```

```python
# On ajoute les points à la liste
sampled_points.extend(sampled pixels)
class labels.extend(np.array([class label]*n samples)[: ,np newaxis])
```

```python
# Conversion en NumPy array
sampled_points = np.array(sampled points)
class labels = np.array(class labels)
```

```python
# On peut naviguer les points à l'aide de la géoréférence
transformer = rasterio transform AffineTransformer(img carte rio transform())
transform_sampled_points= transformer xy(sampled points [: ,0], sampled points [: ,1])
```

```python
fig, ax = plt subplots(nrows=1, ncols=1, figsize=(8, 6))
img carte squeeze(). plot show(cmap=cmap_classes, vmin=0, vmax=12)
ax scatter transform_sampled_points[0], transform_sampled_points[1], c='w', s=1)  # Plot sampled points
ax set title("Carte d'occupation des sols avec les points échantillonnés", font size="small")
plt show()
```
## 6 Classifications d'images supervisées

<figure><img src="image_86.png" /></figure>

Une fois les points sélectionnés, on ajoute les valeurs des bandes provenant d'une image satellite. Pour cela, on utilise la méthode `sample()` de `rasterio`. Éventuellement, la librairie `geopandas` permet de gérer les données d'entraînement sous la forme d'un tableau transportant aussi l'information de géoréférence. Afin de pouvoir classifier ces points, on ajoute les valeurs radiométriques provenant de l'image Sentinel-2 à 4 bandes `RGBNIR_of_S2A.tif`. Ces valeurs seront stockées dans la colonne `value` sous la forme d'un vecteur en format `string`:

1. `points` = [Point(xy) for xy in zipvere样p点s_ sampled_ points[0], transform_sampled_points[1])]

2. `gdf` = geopandas.GeoDataFrame(range(1, len(points)+1), geometry=points, crs = img_carte.rio.crs)

3. `coord_list` = [(x, y) for x, y in zip(gdf["geometry").x, gdf["geometry").y])]

4. with rasterio.open('RGBNIR_of_S2A.tif') as src:

5. `gdf["value"]` = [x for x in src.sample coord_list])]

6. `gdf['class']` = class labels

7. `gdf.to_csv('sampling_points.csv')` # sauvegarde sous forme csv

8. `gdf head()`
<table><thead><tr><th></th><th></th><th>geometry</th><th>value</th><th>class</th></tr></thead><tbody><tr><td>0</td><td>1</td><td>POINT (740369.77 5032078.683)</td><td>[1894, 1994, 2112, 2318]</td><td>1</td></tr><tr><td>1</td><td>2</td><td>POINT (737542.924 5031770.119)</td><td>[1440, 1650, 1449, 5021]</td><td>1</td></tr><tr><td>2</td><td>3</td><td>POINT (736726.722 5031411.786)</td><td>[1666, 1972, 1819, 3437]</td><td>1</td></tr><tr><td>3</td><td>4</td><td>POINT (736816.305 5027470.128)</td><td>[1858, 2078, 2190, 2436]</td><td>1</td></tr><tr><td>4</td><td>5</td><td>POINT (736746.629 5031362.018)</td><td>[2194, 2304, 2268, 3075]</td><td>1</td></tr></tbody></table>

## 6.3 Analyse préliminaire des données

Une bonne pratique avant d'appliquer une technique d'apprentissage automatique est de regarder les caractéristiques de vos données:

1. Le nombre de dimensions (*features*).

2. Certaines dimensions sont informatives (discriminantes) et d'autres ne sont pas.

3. Le nombre de classes.

4. Le nombre de modes (*clusters*) par classes.

5. Le nombre d'échantillons par classe.

6. La forme des groupes.

7. La séparabilité des classes ou des groupes.

Une manière d'évaluer la séparabilité de vos classes est d'appliquer des modèles Gaussiens sur chacune des classes. Le modèle Gaussien multivarié suppose que les données sont distribuées comme un nuage de points symétrique et unimodale. La distribution d'un point $x$ appartenant à la classe $i$ est la suivante:

$$P(x|Classe = i) = \frac{1}{(2\pi)^{D/2} |\Sigma_i|^{1/2}} \exp \left( -\frac{1}{2} (x - m_i)^T \Sigma_k^{-1} (x - m_i) \right)$$

La méthode **QuadraticDiscriminantAnalysis** permet de calculer les paramètres des Gaussiennes multivariées pour chacune des classes.

On peut calculer une distance entre deux nuages Gaussiens avec la distance dites de Jeffries-Matusita (JM) basée sur la distance de Bhattacharyya $B$ (Jensen 2016):

$$JM_{ij} = 2(1 - e^{-B_{ij}}) \\ B_{ij} = \frac{1}{8}(m_i - m_j)^t \left(\frac{\Sigma_i + \Sigma_j}{2}\right)^{-1} (m_i - m_j) + \frac{1}{2} \ln \left( \frac{|(\Sigma_i + \Sigma_j)/2|}{|\Sigma_i|^{1/2} |\Sigma_j|^{1/2}} \right)$$

Cette distance présuppose que chaque classe $i$ est décrite par son centre $m_i$ et de sa dispersion dans l'espace à $D$ dimensions mesurée par la matrice de covariance $\Sigma_i$. On peut en faire facilement une fonction Python à l'aide de ` numpy`:
```python
def bhattacharyya_distance(m1, s1, m2, s2):
    # Calcul de la covariance moyenne
    s = (s1 + s2) / 2

    # Calcul du premier terme (différence des moyennes)
    m_diff = m1 - m2
    term1 = np.dot(np.dot(m_diff.T, np.lINALG.INV(s)), m_diff) / 8

    # Calcul du second terme (différence de covariances)
    term2 = 0.5 * np.log(np.LINALG.DET(s) / np.SQRT(np.LINALG.DET(s1) * np.LINALG.DET(s2)))

    return term1 + term2

def jeffries_matusita_distance(m1, s1, m2, s2):
    B = bhattacharyya_distance(m1, s1, m2, s2)
    return 2 * (1 - np.exp(-B))
```

La figure ci-dessous illustre différentes situations avec des données simulées ainsi que les distances JM correspondantes :
6 Classifications d'images supervisées

```python
X= X[idx,...]

y = df['class'].to numpy()

y= y[idx]

class labels = np.unique(y). tolist() # on cherche à savoir combien de classes uniques

n_classes = len(class labels)

if max(class labels) > n_classes: # il se peut que certaines classes soit absentes
    y_new= []
    for i,l in enumerate(class labels):
        y_new.extend([i]*sum(y==l))
    y_new = np.array(y_new)

couleurs_classes2= [couleurs_classes[c] for c in np.unique(y). tolist.] # couleurs des classes
nom_classes2= [nom_classes[c] for c in np.unique(y). tolist()
cmap_classes2 = ListedColormap(couleurs_classes2)
```

On peut faire une analyse de séparabilité sur notre ensemble d'entrainement de 10 classes. On obtient un tableau symmétrique de 10x10 valeurs. On observe des valeurs inférieures à 1, indiquant des séparabilités faibles entre ces classes sous l'hypothèse du modèle Gaussien:

```python
qda= QuadraticDiscriminantAnalysis(store covarianc e=True)
qda.fit(X, y_new) # calcul des paramètres des distributions Gaussiennes
JM= []
classes= np.unique(y_new). tolist() # étiquettes uniques des classes
for cl1 in classes:
    for cl2 in classes:
        JM.append(jeffries matusita distance(qda.means_[cl1], qda.covariance_[cl1], qda.means_[cl2],
                                qda.covariance_[cl2]))
JM= np.array(JM). reshape(len classes), len classes)
JM= pd.DataFrame(JM, index=classes, columns=classes)
JM head(10)
```

<table><thead><tr><th></th><th>0</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th></tr></thead><tbody><tr><th>0</th><td>0.000000</td><td>1.931891</td><td>1.809590</td><td>1.761760</td><td>1.156486</td><td>1.326107</td><td>1.319344</td><td>1.830671</td><td>1.873676</td><td>1.700417</td></tr><tr><th>1</th><td>1.931891</td><td>0.000000</td><td>1.082978</td><td>0.918865</td><td>1.788737</td><td>1.527192</td><td>1.331400</td><td>1.901749</td><td>0.854802</td><td>1.133180</td></tr><tr><th>2</th><td>1.809590</td><td>1.082978</td><td>0.000000</td><td>0.266647</td><td>1.428062</td><td>1.255001</td><td>1.198888</td><td>1.947302</td><td>0.193032</td><td>0.782982</td></tr><tr><th>3</th><td>1.761760</td><td>0.918865</td><td>0.266647</td><td>0.000000</td><td>1.413401</td><td>1.219793</td><td>1.127950</td><td>1.929637</td><td>0.377379</td><td>0.840250</td></tr><tr><th>4</th><td>1.156486</td><td>1.788737</td><td>1.428062</td><td>1.413401</td><td>0.000000</td><td>0.397103</td><td>0.596618</td><td>1.956182</td><td>1.517926</td><td>1.036828</td></tr><tr><th>5</th><td>1.326107</td><td>1.527192</td><td>1.255001</td><td>1.219793</td><td>0.397103</td><td>0.000000</td><td>0.167221</td><td>1.976696</td><td>1.248383</td><td>0.660213</td></tr><tr><th>6</th><td>1.319344</td><td>1.331400</td><td>1.198888</td><td>1.127950</td><td>0.596618</td><td>0.167221</td><td>0.000000</td><td>1.956804</td><td>1.207618</td><td>0.660589</td></tr><tr><th>7</th><td>1.830671</td><td>1.901749</td><td>1.947302</td><td>1.929637</td><td>1.956182</td><td>1.976696</td><td>1.956804</td><td>0.000000</td><td>1.966022</td><td>1.886064</td></tr><tr><th>8</th><td>1.873676</td><td>0.854802</td><td>0.193032</td><td>0.377379</td><td>1.517926</td><td>1.248383</td><td>1.207618</td><td>1.966022</td><td>0.000000</td><td>0.741273</td></tr><tr><th>9</th><td>1.700417</td><td>1.133180</td><td>0.782982</td><td>0.840250</td><td>1.036828</td><td>0.660213</td><td>0.660589</td><td>1.886064</td><td>0.741273</td><td>0.000000</td></tr></tbody></table>
## 6 Classifications d'images supervisées

Afin d'évaluer chaque classe, on peut calculer la séparabilité minimale, ce qui nous permet de constater que la classe eau a le maximum de séparabilité avec les autres classes.

<figure><img src="image_87.png" /></figure>

## 6.4 Mesures de performance d'une méthode de classification

Lorsque que l'on cherche à établir la performance d'un modèle, il convient de mesurer la performance du classificateur utilisé. Il existe de nombreuses mesures de performance qui sont toutes dérivées de la matrice de confusion. Cette matrice compare les étiquettes provenant de l'annotation (la vérité terrain) et les étiquettes prédites par un modèle. On peut définir $C(i, j)$ comme étant le nombre de prédicitions dont la vérité terrain indique la classe i qui sont prédites dans la classe j. La fonction `confusion_matrix` permet de faire ce calcul, voici un exemple très simple:

$$ 
\begin{align*}
1 \quad y\_true &= [\text{"cat"}, \text{"ant"}, \text{"cat"}, \text{"cat"}, \text{"ant"}, \text{"bird"}, \text{"bird"}] \\
2 \quad y\_pred &= [\text{"ant"}, \text{"ant"}, \text{"cat"}, \text{"cat"}, \text{"ant"}, \text{"cat"}, \text{"bird"}] \\
3 \quad \text{confusion\_matrix}(y\_true, y\_pred, \text{labels}=[\text{"ant"}, \text{"bird"}, \text{"cat"}])
\end{align*}
$$

```
array([2, 0, 0],
      [0, 1, 1],
      [1, 0, 2]])
```

La fonction `classification_report` permet de générer quelques métriques:
6 Classifications d'images supervisées

et à D dimensions, D hyperplans à D-1 dimensions par classe. Le modèle ici est donc une suite de valeurs min et max pour chacune des bandes et des classes:

```python
def parrallepiped_train(X_train, y_train):
    classes = np.unique(y_train).list()
    clf = []
    for cl in classes:
        data_cl = X_train[y_train == cl]  # on cherche les données pour la classe courante
        limits = []
        for b in range(data_cl.shape[1]):
            limits.append([data_cl[:, b].min(), data_cl[:, b].max()])  # on calcule le min et max pour chaque
        → bande
        clf.append(np.array(limits))
    return clf
```

```python
12 clf = parrallepiped_train(X, y_new)
```

La prédiction consiste à trouver pour chaque point la première limite qui est satisfaite. Notez qu'il n'y a aucun moyen de décider quelle est la meilleure classe si le point appartient à plusieurs classes.

```python
@jit(nopython=True)
def parrallepipedPredict(clf, X_test):
    y_pred = []
    for data in X_test:
        y_pred.append(np.nan)
    for cl, limits in enumerate(clf):
        inside = True
        for b, limit in enumerate(limits):
            inside = inside and (data[b] >= limit[0]) & (data[b] <= limit[1])
            if inside:
                break
        if inside:
            y_pred[-1] = cl
    return np.array(y_pred)
```

On peut appliquer ensuite le modèle sur l'image au complet. Les résultats sont assez mauvais, puisque seule la classe eau en bleu semble être bien classifiée.

```python
data_image = img imp to numpy(). transpose(1, 2, 0). reshape(img imp shape[1] * img imp shape[2], 4)
y_image = parrallepipedPredict(clf, data_image)
y_image = y_image . reshape(img imp shape[1], img imp shape[2])

fig, ax = plt . subplots(nrows=1, ncols=1, figsize=(8, 6))
ax . show()  # This line is important to display the image
```
6 Classifications d'images supervisées

<figure><img src="image_88.png" /></figure>

On peut calculer quelques mesures de performance sur l'ensemble d'entrainement :

1 y_pred= parrallepipedPredict(clf, X)

2 nom_classes2= [nom_classes[c] for c in np.unique(y).list()]]

3 print(classification_report(y_new, y_pred, target_names=nom_classes2, zero_division(np.nan)))

<table><thead><tr><th></th><th>precision</th><th>recall</th><th>f1-score</th><th>support</th></tr></thead><tbody><tr><td>Commercial</td><td>1.00</td><td>0.06</td><td>0.11</td><td>100</td></tr><tr><td>Foret</td><td>1.00</td><td>0.09</td><td>0.17</td><td>100</td></tr><tr><td>Faible_végétation</td><td>1.00</td><td>0.02</td><td>0.04</td><td>100</td></tr><tr><td>Sol_nu</td><td>nan</td><td>0.00</td><td>0.00</td><td>100</td></tr><tr><td>Roche</td><td>0.00</td><td>0.00</td><td>0.00</td><td>100</td></tr><tr><td>Route</td><td>0.00</td><td>0.00</td><td>0.00</td><td>100</td></tr><tr><td>Urbain</td><td>0.08</td><td>0.08</td><td>0.08</td><td>100</td></tr><tr><td>Eau</td><td>0.83</td><td>0.88</td><td>0.85</td><td>100</td></tr><tr><td>Végétation éparse</td><td>1.00</td><td>0.01</td><td>0.02</td><td>100</td></tr><tr><td>Roche avec végétation</td><td>0.13</td><td>1.00</td><td>0.23</td><td>100</td></tr><tr><td>accuracy</td><td></td><td></td><td>0.21</td><td>1000</td></tr><tr><td>macro avg</td><td>0.56</td><td>0.21</td><td>0.15</td><td>1000</td></tr><tr><td>weighted avg</td><td>0.56</td><td>0.21</td><td>0.15</td><td>1000</td></tr></tbody></table>
### 6.5.1.1 La malédiction de la haute dimension

Augmenter le nombre de dimension ou de caractéristiques des données permet de résoudre des problèmes complexes comme la classification d'image. Cependant, cela amène beaucoup de contraintes sur le volume des données. Supposons que nous avons N points occupant un segment linéaire de taille d. La densité de points est $N/d$. Si nous augmentons le nombre de dimension D, la densité de points va diminuer exponentiellement en $1/d^D$. Par conséquent, pour garder une densité constante et donc une bonne estimation des parallélépipèdes, il nous faudrait augmenter le nombre de points en puissance de D. Ceci porte le nom de la malédiction de la dimensionnalité (*dimensionality curse*). En résumé, l'espace vide augmente plus rapidement que le nombre de données d'entraînement et l'espace des données devient de plus en plus parcimonieux (*sparse*). Pour contrecarrer ce problème, on peut sélectionner les meilleures caractéristiques ou appliquer une réduction de dimension comme une ACP (Analyse en composantes principales).

### 6.5.2 Plus proches voisins

La méthode des plus proches voisins (*K-Nearest-Neighbors*) est certainement la plus simple des méthodes pour classifier des données. Elle consiste à comparer une nouvelle donnée avec ses voisins les plus proches en fonction d'une simple distance Euclidienne. Si une majorité de ces *K* voisins appartiennent à une classe majoritaire alors cette classe est sélectionnée. Afin de permettre un vote majoritaire, on choisira un nombre impair pour la valeur de *K*. Mallgré sa simplicité, cette technique peut devenir assez demandante en temps de calcul pour un nombre important de points et un nombre élevé de dimensions.

Reprenons l'ensemble d'entraînement formé à partir de notre image RGBNIR précédente :

```
df= pd.read_csv('sampling_points.csv')
# Extraire la colonne 'value'.
# 'value' est une chaîne de caractères comme représentation d'une liste de valeurs.
# Nous devons la convertir en données numériques réelles.
X = df['value'].apply(lambda x: np.fromstring(x[1:-1], dtypeanteer, sep=' ')).to_list()
# on obtient une liste de numpy array qu'il faut convertir en un numpy array 2D
X = np.array([row.tolist() for row in X])
idx= X.sum[1]>0 # il se peut qu'il y ait des valeurs erronées
X = X[idx,...]
y = df['class'].to numpy()
y= y[idx]
class labels = np.unique(y).tost ( ) # on cherche à savoir combien de classes uniques
n_classes = len(class labels)
if max(class labels) > n_classes: # il se peut que certaines classes soit absentes
    y_new= []
    for i,l in enumerate(class labels):
        y_new.extend([i]*sum(y==l))
    y_new = np.array(y_new)
nom_classes2= [nom_classes[c] for c in np.unique(y).tost ()]
```
## 6 Classifications d'images supervisées

Il importe de préalablement centrer (moyenne = 0) et de réduire (variance = 1) les données avant d'appliquer la méthode K-NN; avec cette méthode de normalisation, on dit parfois que l'on blanchit les données. Puisque la variance de chaque dimension est égale à 1 (et donc l'inertie totale est égale au nombre de bandes), on s'assure qu'elle ait le même poids ait le même poids dans le calcul des distances entre points. Cette opération porte le nom de **StandardScaler** dans **scikit-learn**. On peut alors former un pipeline de traitement combinant les deux opérations :

```groovy
1 conf = Pipeline(
2     steps=[("scaler", StandardScaler()), ("knn", KNeighborsClassifier(n_neighbors=1))]
3 )
```

Avant d'effectuer un entraînement, on met généralement une portion des données pour valider les performances :

```groovy
1 X_train, X_test, y_train, y_test = train_test_split(X, y_new, test_size=0.2, random_state=0)
```

On peut visualiser les frontières de décision du K-NN pour différentes valeurs de *K* lorsque seulement deux bandes sont utilisées (Rouge et proche infra-rouge ici) :

```
Number of mislabeled points out of a total 200 points : 143
```

```
Number of mislabeled points out of a total 200 points : 141
```

```
Number of mislabeled points out of a total 200 points : 136
```

```
Number of mislabeled points out of a total 200 points : 130
```
6 Classifications d'images supervisées

<figure><img src="image_89.png" /><figcaption>FIGURE 6.3 – Frontières de décision pour le classificateur K-NN</figcaption></figure>

On peut voir comment les différentes frontières de décision se forment dans l'espace des bandes Rouge-NIR.
L'augmentation de K rend ces frontières plus complexes et le calcul plus long.

1. `clf.set_params(knn__weights='distance', knn__n_neighbors = 7).fit(X_train, y_train)`

2. `y_pred =clfPredict(X_test)`
6 Classifications d'images supervisées

3 print("Nombre de points misclassifiés sur %d points : %d")

4 % (X_test.shape[0], (y_test != y_pred).sum())

Nombre de points misclassifiés sur 200 points : 117

Le rapport de performance est le suivant :

1 nom_classes2= [nom_classes[c] for c in np.unique(y).list() ]

2 print(classification_report(y_test, y_pred, target_names=nom_classes2, zero_division= np.nan))

<table>
  <thead>
    <tr>
      <th></th>
      <th>precision</th>
      <th>recall</th>
      <th>f1-score</th>
      <th>support</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td>Commercial</td>
      <td>0.38</td>
      <td>0.40</td>
      <td>0.39</td>
      <td>15</td>
    </tr>
    <tr>
      <td>Foret</td>
      <td>0.45</td>
      <td>0.82</td>
      <td>0.58</td>
      <td>11</td>
    </tr>
    <tr>
      <td>Faible_végétation</td>
      <td>0.29</td>
      <td>0.15</td>
      <td>0.20</td>
      <td>27</td>
    </tr>
    <tr>
      <td>Sol_nu</td>
      <td>0.53</td>
      <td>0.45</td>
      <td>0.49</td>
      <td>22</td>
    </tr>
    <tr>
      <td>Roche</td>
      <td>0.38</td>
      <td>0.26</td>
      <td>0.31</td>
      <td>23</td>
    </tr>
    <tr>
      <td>Route</td>
      <td>0.16</td>
      <td>0.17</td>
      <td>0.16</td>
      <td>18</td>
    </tr>
    <tr>
      <td>Urbain</td>
      <td>0.25</td>
      <td>0.20</td>
      <td>0.22</td>
      <td>20</td>
    </tr>
    <tr>
      <td>Eau</td>
      <td>0.96</td>
      <td>0.96</td>
      <td>0.96</td>
      <td>24</td>
    </tr>
    <tr>
      <td>Végétation éparse</td>
      <td>0.26</td>
      <td>0.53</td>
      <td>0.35</td>
      <td>15</td>
    </tr>
    <tr>
      <td>Roche avec végétation</td>
      <td>0.40</td>
      <td>0.40</td>
      <td>0.40</td>
      <td>25</td>
    </tr>
    <tr>
      <td>accuracy</td>
      <td></td>
      <td></td>
      <td>0.41</td>
      <td>200</td>
    </tr>
    <tr>
      <td>macro avg</td>
      <td>0.40</td>
      <td>0.43</td>
      <td>0.40</td>
      <td>200</td>
    </tr>
    <tr>
      <td>weighted avg</td>
      <td>0.42</td>
      <td>0.41</td>
      <td>0.40</td>
      <td>200</td>
    </tr>
  </tbody>
</table>

La matrice de confusion peut-être affichée de manière graphique :

1 disp= ConfusionMatrixDisplay.from(y_test, y_pred, display labels=nom_classes2,
  → xticks_rotate='vertical')
## 6 Classifications d'images supervisées

<figure><img src="image_90.png" /></figure>

L'application du modèle (la prédiction) peut se faire sur toute l'image en transposant l'image sous forme d'une
matrice avec Largeur x Hauteur lignes et 4 colonnes :

1. `data_image = img_rgbnir.to numpy().transpose(1,2,0).reshape(img_rgbnir.shape[1]*img_rgbnir.shape[2],4)`

2. `y_classe = clf Predict(data_image)`

3. `y_classe = y_classe.reshape(img_rgbnir.shape[1], img_rgbnir.shape[2])`

1. `fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(8, 6))`

2. `plt.imshow(y_classe, cmap=cmap_classes2)`

3. `ax.set_title("Carte d'occupation des sols avec K-NN", fontsize="small")`

4. `plt.show()`
<figure><img src="image_91.png" /></figure>

### 6.5.3 Méthodes par arbre de décision

La méthode par arbre de décision consiste à construire une cascade de règles de décision sur chaque caractéristique du jeu de donnée (Breiman et C. Stone 1984). On pourra trouver plus de détails dans la documentation de scikit-learn (Decision Trees). Les arbres de décision on tendance à surapprendre surtout si le nombre de dimensions est élevé. Il est donc conseillé d'avoir un bon ratio entre le nombre d'échantillons et le nombre de dimensions.

1. `X_train`, `X_test`, `y_train`, `y_test` = `train_test_split(X, y_new, test_size=0.2, random_state=0)`

Number of mislabeled points out of a total 200 points : 167

Number of mislabeled points out of a total 200 points : 154

Number of mislabeled points out of a total 200 points : 143

Number of mislabeled points out of a total 200 points : 128
6 Classifications d'images supervisées

<figure><img src="image_92.png" /><figcaption>FIGURE 6.4 – Frontières de décision pour des arbres de décision de différente profondeur</figcaption></figure>

On peut observer que les frontières de décision sont formées d'un ensemble de plans simple. Chaque plan étant issu d'une règle de décision formé d'un seuil sur chacune des dimensions. On entraine un arbre de décision avec une profondeur maximale de 5:
6 Classifications d'images supervisées

```python
1   clf = tree.DecisionTreeClassifier(max_depth=5)
2   clf.fit(X_train, y_train)
3   y_pred = clf Predict(X_test)
4   print("Nombre de points misclassifiés sur %d points : %d")
5   % (X_test.shape[0], (y_test != y_pred).sum())
```

Nombre de points misclassifiés sur 200 points : 130

Le rapport de performance et la matrice de confusion:

<table><thead><tr><th></th><th>precision</th><th>recall</th><th>f1-score</th><th>support</th></tr></thead><tbody><tr><td>Commercial</td><td>0.37</td><td>0.47</td><td>0.41</td><td>15</td></tr><tr><td>Foret</td><td>0.57</td><td>0.73</td><td>0.64</td><td>11</td></tr><tr><td>Faible_végétation</td><td>0.19</td><td>0.19</td><td>0.19</td><td>27</td></tr><tr><td>Sol_nu</td><td>0.57</td><td>0.18</td><td>0.28</td><td>22</td></tr><tr><td>Roche</td><td>0.40</td><td>0.09</td><td>0.14</td><td>23</td></tr><tr><td>Route</td><td>0.32</td><td>0.44</td><td>0.37</td><td>18</td></tr><tr><td>Urbain</td><td>nan</td><td>0.00</td><td>0.00</td><td>20</td></tr><tr><td>Eau</td><td>0.95</td><td>0.79</td><td>0.86</td><td>24</td></tr><tr><td>Végétation éparse</td><td>nan</td><td>0.00</td><td>0.00</td><td>15</td></tr><tr><td>Roche avec végétation</td><td>0.20</td><td>0.68</td><td>0.31</td><td>25</td></tr><tr><td>accuracy</td><td></td><td></td><td>0.35</td><td>200</td></tr><tr><td>macro avg</td><td>0.45</td><td>0.36</td><td>0.32</td><td>200</td></tr><tr><td>weighted avg</td><td>0.44</td><td>0.35</td><td>0.31</td><td>200</td></tr></tbody></table>

<figure><img src="image_93.png" /></figure>

L'application du modèle (la prédiction) peut se faire sur toute l'image en transposant l'image sous forme d'une
matrice avec Largeur x Hauteur lignes et 4 colonnes:
## 6 Classifications d'images supervisées

<figure><img src="image_94.png" /></figure>

Il est possible de visualiser l'arbre avec l'outil `SuperTree` mais cela contient beaucoup d'information

1

```python
from supertree import SuperTree # <- import supertree :)
```

2

```python
super_tree = SuperTree(clf, X_train, y_train, ['Bleu', 'Vert', 'Rouge', 'NIR'], nom_classes=2)
```

3

4

```python
super_tree.show_tree()
```

```
<IPython.core.display.HTML object>
```

On peut voir que le nœud le plus haut utilise la bande proche-infrarouge pour séparer l'eau de la forêt comme première décision.

## 6.6 Méthodes paramétriques

Les méthodes paramétriques se basent sur des modélisations statistiques des données pour permettre une classification. Contrairement au méthodes non paramétriques, elles ont un nombre fixe de paramètres qui ne dépend pas de la taille du jeu de données. Par contre, des hypothèses sont faites a priori sur le comportement
statistique des données. La classification consiste alors à trouver la classe la plus vraisemblable dont le modèle statistique décrit le mieux les valeurs observées. L'ensemble d'entraînement permettra alors de calculer les paramètres de chaque Gaussienne pour chacune des classes d'intérêt.

### 6.6.1 Méthode Bayésienne naïve

La méthode Bayésienne naïve Gaussienne consiste à poser des hypothèses simplificatrices sur les données, en particulier l'indépendance des données et des dimensions. Ceci permet un calcul plus simple.

```python
1 from sklearn.naive_bayes import GaussianNB
2 gnb = GaussianNB()
3 y_pred = gnb.fit(X_train, y_train).predict(X_test)
4 print("Nombre de points erronés sur %d points : %d")
5 % (X_test.shape[0], (y_test != y_pred).sum())
```

Nombre de points erronés sur 200 points : 131

<figure><img src="image_95.png" /><figcaption>FIGURE 6.5 – Frontières de décision pour un classificateur Bayésien naif</figcaption></figure>

On observe que les frontières de décision sont beaucoup plus régulières que pour K-NN.

```python
1 gnb.fit(X_train, y_train)
2 y_pred = gnb Predict(X_test)
3 print("Nombre de points misclassifiés sur %d points : %d")
4 % (X_test.shape[0], (y_test != y_pred).sum())
```
# 6 Classifications d'images supervisées

Nombre de points misclassifiés sur 200 points : 131

De la même manière, la prédiction peut s'appliquer sur toute l'image:

<figure><img src="image_96.png" /></figure>

## 6.6.2 Analyse discriminante quadratique (ADQ)

L'analyse discriminante quadratique peut-être vue comme une généralisation de l'approche Bayésienne naive qui suppose des modèles Gaussiens indépendants pour chaque dimension et chaque point. Ici, on va considérer un modèle Gaussien multivarié.

1. qda = QuadraticDiscriminantAnalysis(store_covariance=True)

2. qda.fit(X_train, y_train)

3. y_pred = qdaPredict(X_test)

4. print("Nombre de points misclassifiés sur %d points : %d")

5. % (X_test.shape[0], (y_test != y_pred).sum())

Nombre de points misclassifiés sur 200 points : 124
# 6 Classifications d'images supervisées

Les Gaussiennes multivariées peuvent être visualiser sous forme d'ellipses décrivant le domaine des valeurs de chaque classe:

<figure><img src="image_97.png" /></figure>

De la même manière, la prédiction peut s'appliquer sur toute l'image:
6 Classifications d'images supervisées

<figure><img src="image_98.png" /></figure>
# Bibliographie

Achanta, Kevin Smith adhakrishna, Appu Shaji et Sabine Süsstrunk. 2012. « SLIC Superpixels Compared to State-of-the-art Superpixel Methods. » *TPAMI*: 636-643. https://doi.org/10.1109/TPAMI.2012.120.

Breiman, Friedman, L. et C. C. Stone. 1984. *Classification and Regression Trees*. Wadsworth, Belmont, CA.

Cooley, James W. et John W. Tukey. 1965. « An algorithm for the machine calculation of complex Fourier series. » *Math. Comput.*: 297-301. https://web.stanford.edu/class/cme324/classics/cooley-tukey.pdf.

Harris, Millman, C. R. 2020. « Array programming with NumPy. » *Nature*: 357-362. https://doi.org/10.1038/s41586-020-2649-2.

Hoyer, S. et J. Hamman. 2017. « xarray: N-D labeled Arrays and Datasets in Python. » *Journal of Open Research Software* 5 (1): 10. https://doi.org/10.5334/jors.148.

Jahne, Scharr, B et Korkel S. 1999. *Principles of filter design*. Handbook of Computer Vision; Applications; Academic Press.

Jaworek-Korjakowska, P., J.; Kleczek. 2018. « Region Adjacency Graph Approach for Acral Melanocytic Lesion Segmentation. » *Applied Sciences* 8: 1430. 10.3390/app8091430.

Jensen, J. R. 2016. *Introductory digital image processing: A remote sensing perspective*. Pearson Education, Inc.

Kokaly, Clark, R. F. et A. J. Klein. 2017. « USGS Spectral Library Version 7 Data: U.S. Geological Survey data release. » *Applied Sciences*. 10.5066/F7RR1WDJ.

Lee, J. S. 1986. « Speckle suppression and analysis for synthetic aperture radar images. » *Opt. Eng.* 25 (5): 636-643. https://doi.org/10.1117/12.7973877.

OGC. 2019. « OGC GeoTIFF Standard. » https://docs.ogc.org/is/19-008r4/19-008r4.html/.