Vous êtes géomaticiens  et vous voulez vous initier au machine learning, mais vous rencontrez des difficultés pour démarrer? Cette article est faite pour vous. Vous allez découvrir étape par étape les bases de l’apprentissage automatique sur un projet de classification d’images. 

Machine learning de quoi on parle?

Le machine learning ou apprentissage automatique en Français, est une technologie d’intelligence artificielle qui permet aux ordinateurs d’apprendre par eux-mêmes, sans avoir été programmés à effectuer des tâches spécifiques. Comment ? En reproduisant un comportement grâce à des algorithmes, eux-mêmes alimentés par un grand nombre de données. Confronté à de nombreuses situations, l’algorithme apprend quelle est la décision à adopter et créé un modèle. La machine peut automatiser les tâches en fonction des situations. 

Alors qu’elle est notre objectif? J’ai longtemps travaillé sur des images satellites mais toujours en utilisant des logiciels payants comme Arcgis, IDRISI, ENVI, Ecognition ou dans le domaine de l’open source comme Grass, SAGA, QGis, GDAL ou Monteverdi. Les résultats sont excellents mais on rencontre toujours des problèmes de classifications de précisions et on passe énormément du temps pour produire des modèles peu prédictifs ou non reproductibles. Aujourd’hui le machine learning est devenu une solution incontournable dans le domaine de la télédétection et offre d’immenses possibilités de modélisation de données dans l’espace et dans le temps.  

Pour les besoins de ce tutoriel,  nous allons utiliser une image du 15 mai 2018 à l’ouest de Bordeaux . Sur la dalle que j’ai téléchargé, j’ai également utilisé un shp pour identifier 3 classes d’occupation du sol:  les zones artificialisés (Batiment, chemin, routes..) , la végétation(cultures, prairie, bois…)  et l’eau. Le jeux de données est accessible ici .Dans article précèdent, nous avons déjà vu comment télécharger les images Sentinel. 

L’objectif ici c’est d’apprendre à la machine à prédire ces trois classes d’occupation du sol et d’évaluer la précision du model. Il faut savoir qu’un projet d’apprentissage automatique peut ne pas être linéaire, mais comporte un certain nombre d’étapes bien connues schématisé ci-dessous.

De la modélisation au déploiement d’un projet en production

Source: Wikipedia

Pour résumer il existe différentes étapes qu’on peux décomposer en 2 grandes catégories : Modélisation et Déploiement. Nous allons les décortiquer étape par étape. Vous n’avez pas besoin de tout comprendre (du moins pas pour le moment) . Votre objectif est de parcourir le didacticiel de bout en bout et d’obtenir des résultats. Vous n’avez pas besoin de tout comprendre lors du premier passage. Faites la liste de vos questions au fur et à mesure.

1-Installation de la plateforme Python et SciPy.

Il existe énormément d’outils open source pour mener à bien votre projet. Vous trouverez sur dynacentrix.com les logiciels couramment utilisés. Pour les besoins de ce tuto, vous aurez besoin de ces 5 logiciels :

• scipy : ensemble d’outils de calculs scientifiques sous python
• gdal :  traitement de données rasteurs
• rasterio: exploitation des données rasteur avec python
• numpy : est destinée à manipuler des matrices ou tableaux multidimensionnels
• matplotlib:  est destinée à tracer et visualiser des données sous formes de graphiques. 
• pandas:  pour la manipulation et l’analyse des données.
• sklearn :  une librairie d’apprentissage automatique 

Il existe de nombreuses façons d’installer ces bibliothèques. Mon meilleur conseil est  d’installer Anaconda il fourni l’ensemble des bibliothèques dédiés à la data science : Python, Jupyter Notebook et pleins d’autres  bibliothèques scientifiques pour l’analyse de données. Si vous voulez en savoir plus sur l’installation d’anaconda,  il existe un bon tuto sur openclassroom.

Apres avoir installer les bibliothèques, vous pouvez lancer jupyter avec cette commande : jupyter notebook. Le tuto est également accessible ici.

2- Préparation des données

Apres avoir charger les données,  j’utilise le module raster du package gdal « ReadAsArray » permet d’assembler les informations de géolocalisation du GeoTIFF et les valeurs numériques (DN) sous la forme d’un tableau NumPy. Ce qui nous permet d’avoir le même nombre de lignes pour nos deux images et dans le bon ordre. Puis les données ont été transformées en un tableau bidimensionnel. Ce qui est attendu par la majorité des algorithmes ML, où chaque ligne représente un pixel.

Schémas de restructuration des données

#Importer les images
img_ds = gdal.Open('bordeaux1.tif', gdal.GA_ReadOnly)
roi_ds = gdal.Open('training_data3classes.tif', gdal.GA_ReadOnly)

img = np.zeros((img_ds.RasterYSize, img_ds.RasterXSize, img_ds.RasterCount),
gdal_array.GDALTypeCodeToNumericTypeCode(img_ds.GetRasterBand(1).DataType))
for b in range(img.shape[2]):
img[:, :, b] = img_ds.GetRasterBand(b + 1).ReadAsArray()

roi = roi_ds.GetRasterBand(1).ReadAsArray().astype(np.uint8)

# Affichez-les
plt.subplot(121)
plt.imshow(img[:, :,3 ], cmap='gray')
plt.title('B4')

plt.subplot(122)
plt.imshow(roi, cmap=plt.cm.Spectral)
plt.title('ROI Training Data')

plt.show()

Restructuration des données en un tableau bidimensionnel

Maintenant, nous allons diviser notre tableau  en deux jeux de données:  apprentissage et validation. Nous utiliserons la validation croisée pour éviter un surapprentissage. Mais aussi nous assurer une bonne répartition des classes dans les deux jeux de données.

# #### Association de Y avec X
# Maintenant que nous avons l'image que nous voulons classer (nos entrées de fonctionnalités X),
#et le retour sur investissement avec les étiquettes de couverture terrestre (nos données étiquetées Y), nous devons les coupler
# dans les tableaux NumPy afin que nous puissions les transmettre aux modèles:


# Trouvez combien d'entrées non nulles nous avons : c'est-à-dire combien d'échantillons de données d'entraînement?
n_samples = (roi > 0).sum()
print('Nous avons {n} pixels'.format(n=n_samples))

# What are our classification labels?
labels = np.unique(roi[roi > 0])
print('The training data include {n} classes: {classes}'.format(n=labels.size, 
classes=labels))
# Nous aurons besoin d'une matrice "X" contenant nos caractéristiques, et d'un tableau "y" contenant nos étiquettes
# Ceux-ci auront n_samples lignes

X = img[roi > 0, :] # inclure toutes les bandes
y = roi[roi > 0]

print('Dimensions du tableau en X: {sz}'.format(sz=X.shape))
print('Dimensions du tableauen Y: {sz}'.format(sz=y.shape))

# Scinder les données en deux: apprentissage et test

X_train, X_test, Y_train, Y_test = train_test_split(X, y, test_size=0.20, shuffle=True,random_state=42,stratify=y)

#outil standardisation
from sklearn.preprocessing import StandardScaler
#centrage-réduction des variables
cr = StandardScaler(with_mean=True,with_std=True)
#calcul des paramètres + centrage réduction du train set
XTrainStd = cr.fit_transform(X_train)
#comparaison des moyennes, avant…
print(np.mean(X_train,axis=0))
#... et après CR (centrage-réduction)
print(np.mean(XTrainStd,axis=0))

Données d’apprentissage et Validation

Le test_size (0.3) dans l’extrait de code ci-dessus signifie que la proportion de  d’entraînement est de 70% et 30% pour la validation. Nous nous assurons d’avoir les mêmes proportions d’observations positives dans les deux sous-échantillons par un tirage stratifié. Les variables prédictives étant définies sur des échelles différentes, il est recommandé de les standardiser avant de procéder à un apprentissage. Nous utilisons la classe StandardScaler de la librairie scikit-learn. fit_transform() calcule tout d’abord les paramètres de la transformation (les moyennes et écartstype pour chaque variable), puis l’applique sur ces mêmes données X_Train.

3- Modélisation et apprentissage

Maintenant que tout est en place, nous allons pouvoir construire notre model de prédiction. Nous allons tester le modèl de prédiction:  forêt aléatoire ou (Random Forest) qui permet d’obtenir des modèles prédictifs pour la classification et la régression proposés par Breiman et al. (1984). Cette  méthode permet de prédire l’appartenance d’observations (observations, individus) à une classe d’une variable qualitative, sur la base de variables explicatives quantitatives et/ou qualitatives. Vous pouvez trouver plus d’infos sur Random Forrest ici . 

# ## Formation du modèle
# Maintenant que nous avons notre matrice X d'entrées de caractéristiques (les bandes spectrales) et notre tableau y (les étiquettes), nous pouvons entraîner notre modèle.
# Initialize our models
rf = RandomForestClassifier( n_estimators = 500, criterion = 'gini', max_depth = 4, 
min_samples_split = 2, min_samples_leaf = 1, max_features = 'auto', 
bootstrap = True, oob_score = True, n_jobs = 4, random_state = None, verbose = True) 

#Apprendre notre modèle sur les  données d'entraînement
rf = rf.fit(X_train, Y_train)
print ("Trained model :: ", rf)

Models Random Forrest et SVM

Enfin, nous exécutons le modèle sur X_train et Y_Train . L’ajustement du modèle prendra un certain temps en fonction de la taille de vos données et de la puissance de calcul. Ici le traitement a été parallélisé.

4- Evaluation des deux modèles

L’approche courante d’évaluation consiste à réaliser la prédiction sur l’échantillon test, puis à la confronter avec les valeurs observées de la variable cible. Nous allons analyser ici la précision, le rappel , le taux d’erreurs et les matrices de confusion des deux modèles

##### Random Forest diagnostics
#Grâce à l'ajustement de notre modèle Random Forest, nous pouvons consulter le score de prédiction "Out-of-Bag" (OOB):

# ÉVALUATION DE LA PRÉCISION
print('Our OOB prediction of accuracy is: {oob}%'.format(oob=rf.oob_score_ * 100))
########## Random Forest ########
#Évaluation du modèle avec jeu de données de test
# Prédiction au niveau de l'ensemble de données de test
Y_pred = rf.predict(X_test)
Y_pred_train=rf.predict(X_train)
# Précision de l'entraînement et des tests
print ("Train Accuracy :: ", accuracy_score(Y_train, Y_pred_train))
print ("Test Accuracy :: ", accuracy_score( Y_test, Y_pred))

## Matrice de Confusion
print("Matrice de Confusion:\n%s" %
metrics.confusion_matrix(Y_test, Y_pred))

# Précision de la classification:
print("Rapport de classification:\n%s" %
metrics.classification_report(Y_test, Y_pred))
print("Précision de la classification: %f" %
metrics.accuracy_score(Y_test, Y_pred))

 

Comme nous pouvons le voir dans les matrices de confusions ci-dessus, il y a des milliers de pixels  bien classés et des erreurs de commission ou d ‘omission. Mais la proportion par rapport à la taille totale des données est moindre. La précision et le rappel obtenus sur les données de test sont supérieurs à 0,8 avec une précision globale de 90%. Il est également possible d’intégrer des indices de végétation pour optimiser le model.

3- Déploiement du model 

Une fois que la précision souhaitée est atteinte, nous pouvons sauvegarder notre modèle et l’utiliser pour prédire les nouvelles données. Ce modèle  avec des modifications mineures peut être appliqué pour des applications similaires..

# call your random forest model
rf = "model_rf.pkl" 
clf = joblib.load(rf)

4- Prédiction

L’objectif attendu ici c’est de prédire nos trois classes d’occupation du sol avec les jeux de données de test que nous avons conservées séparément. Mais aussi d’effectuer divers contrôles de précision. Pour cela nous devons d’abord transformer notre images en 3 démentions pour dans un tableaux avant d’effectuer une prédiction

# Take our full image, ignore the Fmask band, and reshape into long 2d array (nrow * ncol, nband) for classification
new_shape = (img.shape[0] * img.shape[1], img.shape[2] )

img_as_array = img[:, :, :7].reshape(new_shape)
print('Reshaped from {o} to {n}'.format(o=img.shape, n=img_as_array.shape))

# Now predict for each pixel
class_prediction = clf.predict(img_as_array)

# Reshape our classification map
class_prediction = class_prediction.reshape(img[:, :, 0].shape)

Visualiser le s résultats

Exporter les données

#Exporter les données
#Commandes pyrsgis:
#https://ddeevv.com/lib/pypi/pyrsgis.html
from pyrsgis import raster 
# Assign file names 
imagpredict = 'bordeaux1.tif'
# Read the rasters as array
ds3, arrayimagpredict = raster.read(imagpredict, bands='all')
# Print the size of the arrays
print("Image classifiée: ", arrayimagpredict.shape)
outFile = 'rf_sentinel_predicted.tif'
raster.export(class_prediction, ds3, filename=outFile, dtype='int')

 

La précision du modèle a déjà été évaluée avec précision et rappel – vous pouvez également effectuer les vérifications comme  le coefficient kappa, la moyenne harmonique sur le nouveau raster prédit. Outre les défis susmentionnés de la classification des données satellitaires, d’autres limitations intuitives incluent l’incapacité du modèle à prédire sur les données acquises à différentes saisons et sur différentes régions, en raison de la variation des signatures spectrales.

A présent vous avez les clefs en main  pour bien structurer votre projet de Machine Learning.  J’ai volontairement utilisé ici le Random forest mais il serait lus pertinent d’utiliser un SVM, j’y reviendrai dans un prochain tuto. Certains modèles sont encore plus complexes comme les réseaux de neurones à convolution (CNN)  et produisent  de meilleurs résultats. Le principal avantage de ces techniques de classification est l’évolutivité une fois que le modèle est formé. Pour commencer avec CNN pour la classification des données satellitaires, j’en parlerai également dans mon prochain article. 

Veuillez trouver les données utilisées et le script complet ici .