Classification des données satellitaires#

Mots clés données utilisées ; géomédiane sentinel-2, méthodes de données ; classification, méthodes de données ; apprentissage automatique

Description#

Après avoir exécuté avec succès le bloc-notes « 3_Evaluate_optimize_fit_classifier », nous pouvons maintenant utiliser notre modèle de classification pour prédire les valeurs des nouvelles données satellite. Ce bloc-notes vous guidera dans le chargement des données satellite à partir de l’ODC, le calcul des mêmes couches de fonctionnalités que nous avons fait dans le premier bloc-notes lorsque nous avons extrait les données d’entraînement de l’ODC et l’utilisation de notre modèle pour classer les données satellite. Dans un premier temps, nous classons quelques petites régions pour visualiser les performances de notre modèle. Le bloc-notes tentera ensuite de classer une région beaucoup plus grande et d’enregistrer les résultats sur le disque au format Cloud-Optimized GeoTIFF (COG).

Les étapes sont les suivantes :

  1. Ouvrez le modèle que nous avons généré dans le bloc-notes précédent, « 3_Evaluate_optimize_fit_classifier »

  2. Redéfinir la fonction de couche d’entités que nous avons utilisée pour extraire les données d’entraînement de l’ODC dans le premier notebook, « 1_Extract_training_data »

  3. Parcourez un ensemble de petits emplacements de test en extrayant des données satellite de l’ODC, puis calculez les couches de fonctionnalités et classez les données à l’aide de notre modèle.

  4. Représentez graphiquement les résultats de la classification de nos petits emplacements de test

  5. Définir un nouvel emplacement beaucoup plus grand pour charger à partir de l’ODC et classer en utilisant le même modèle

  6. Sauvegarder nos résultats sur le disque en tant que COG

Commencer#

Pour exécuter cette analyse, exécutez toutes les cellules du bloc-notes, en commençant par la cellule « Charger les packages ».

Charger des paquets#

[1]:

import datacube
import xarray as xr
from joblib import load
import matplotlib.pyplot as plt
from odc.geo.xr import write_cog

from deafrica_tools.datahandling import load_ard
from deafrica_tools.bandindices import calculate_indices
from deafrica_tools.dask import create_local_dask_cluster
from deafrica_tools.plotting import rgb, display_map
from deafrica_tools.classification import predict_xr

import warnings
warnings.filterwarnings("ignore")

Configurer un cluster Dask#

Cela nous aidera à réduire l’utilisation de notre mémoire et à effectuer l’analyse en parallèle. Si vous souhaitez afficher le tableau de bord, cliquez sur le lien hypertexte qui s’imprime sous la cellule. Vous pouvez utiliser le tableau de bord pour surveiller la progression des calculs.

[2]:

create_local_dask_cluster()

Client

Client-84ebbb0d-d4a0-11ef-8922-1ede66a139be

Connection method: Cluster object Cluster type: distributed.LocalCluster
Dashboard: /user/victoria@kartoza.com/proxy/37507/status

Cluster Info

Paramètres d’analyse#

  • model_path : le chemin d’accès à l’emplacement où le modèle exporté à partir du bloc-notes précédent est stocké

  • testing_locations : un dictionnaire avec des valeurs contenant des points de latitude et de longitude, et des clés représentant un identifiant unique pour identifier les emplacements. Les points lat et lon définissent le centre des images satellites que nous allons charger pour exécuter de petites classifications de test

  • buffer : la taille, en degrés décimaux, à charger autour des points centraux de latitude et de longitude dans locations. Ce nombre dépendra ici de la capacité de calcul/RAM disponible sur l’instance Sandbox, ainsi que du type et du nombre de couches d’entités calculées. Une valeur de 0,1 (qui donne une étendue d’analyse de 0,2 x 0,2 degré) fonctionne généralement bien sur l’instance Sandbox par défaut.

  • dask_chunks : Dask fonctionne en divisant de grands ensembles de données en morceaux, qui peuvent être lus individuellement. Ce paramètre spécifie la taille des morceaux en nombre de pixels, par exemple {'x':1000,'y':1000}

  • « résultats » : un emplacement de dossier pour stocker les GeoTIFF classifiés

[3]:

model_path = 'results/ml_model.joblib'

testing_locations = {
    '1': (24.4780, 33.0457),
    '2': (22.7712, 28.5223),
    '3': (30.3971, 30.7996),
    '4': (31.1751, 31.2221),
}

buffer = 0.1

dask_chunks = {'x':1000,'y':1000}

results = 'results/'

Se connecter au datacube#

[4]:

dc = datacube.Datacube(app='prediction')

Importer le modèle#

Le code ci-dessous rouvrira également les données de formation que nous avons exportées depuis « 3_Evaluate_optimize_fit_classifier.ipynb » et récupérera les noms de colonnes (fonctionnalités que nous avons sélectionnées).

[5]:

model = load(model_path).set_params(n_jobs=1)

Faire une prédiction#

Redéfinition de la fonction de couche d’entités#

Parce que nous avons choisi d’utiliser toutes les fonctionnalités extraites dans « 1_Extract_training_data.ipynb », nous pouvons simplement copier-coller la fonction « feature_layers » du premier notebook « Extracting_training_data <1_Extract_training_data.ipynb> »__ dans la cellule ci-dessous (cela a déjà été fait pour vous).

Si vous utilisez ce bloc-notes pour exécuter vos propres classifications (c’est-à-dire sans exécuter l’exemple par défaut), vous devrez redéfinir votre fonction de couche d’entités ci-dessous, en prenant soin de faire correspondre les entités du modèle formé. Par exemple, si vous avez effectué une sélection d’entités et supprimé des entités du modèle, vous devrez imiter ce processus ici en supprimant des entités dans les données de prédiction. En bref, les entités du modèle doivent correspondre précisément à celles des données que vous classez.

Remarque : Étant donné que nous utilisons dask pour nous aider à mettre à l’échelle les opérations, nous devons ajouter un paramètre « dask_chunks » à la couche « slope » afin que cette couche interagisse avec les tableaux dask dans les données satellite. Comme « slope » est une couche 2D, nous n’avons pas besoin de « time » dans l’argument chunks, seulement de « x » et « y ».

[6]:

from datacube.testutils.io import rio_slurp_xarray

def feature_layers(query):
    #connect to the datacube
    dc = datacube.Datacube(app='feature_layers')

    #load ls8 geomedian
    ds = dc.load(product='gm_s2_annual',
                 **query)

    #calculate some band indices
    da = calculate_indices(ds,
                           index=['NDVI','LAI','MNDWI'],
                           drop=False,
                           satellite_mission='s2')

    #add slope dataset
    url_slope = "https://deafrica-input-datasets.s3.af-south-1.amazonaws.com/srtm_dem/srtm_africa_slope.tif"
    slope = rio_slurp_xarray(url_slope, gbox=ds.geobox)
    slope = slope.to_dataset(name='slope').chunk(dask_chunks)

    #merge results into single dataset
    result = xr.merge([da, slope],compat='override')

    return result.squeeze()

Configurer une requête Datacube#

Ces options de requête doivent correspondre aux paramètres de requête dans « 1_Extract_training_data.ipynb », à moins qu’il n’y ait des mesures qui n’ont plus besoin d’être chargées car elles ont été supprimées lors d’un processus de sélection de fonctionnalités (ce qui n’a pas été fait dans l’exemple par défaut).

[7]:

# Set up the inputs for the ODC query
time = ('2019')
measurements =  ['blue','green','red','nir','swir_1','swir_2','red_edge_1',
                 'red_edge_2', 'red_edge_3', 'BCMAD', 'EMAD', 'SMAD']
resolution = (-20,20)
output_crs='epsg:6933'

Parcourez les emplacements de test et prédisez#

Pour chaque emplacement que nous avons répertorié dans le dictionnaire « test_locations », nous calculons les couches d’entités, puis utilisons la fonction DE Africa « predict_xr <digitalearthafrica/deafrica-sandbox-notebooks> »__ pour classer les données.

La fonction predict_xr est un wrapper xarray autour des méthodes d’estimation sklearn .predict() et .predict_proba(), et s’appuie sur dask-ml ParallelPostfit pour exécuter les prédictions avec dask. Predict_xr peut calculer des prédictions, des probabilités de prédiction et renvoyer les couches d’entités en entrée. Lisez la documentation pour plus d’informations sur les capacités de cette fonction.

[8]:

predictions = []

for key, value in testing_locations.items():

    print('working on: ' + key)

     # generate a datacube query object
    query = {
        'x': (value[1]-buffer, value[1]+buffer),
        'y': (value[0]+buffer, value[0]-buffer),
        'time': time,
        'measurements': measurements,
        'resolution': resolution,
        'output_crs': output_crs,
        'dask_chunks' : dask_chunks,
    }

    #calculate features
    data = feature_layers(query)

    #predict using the imported model
    predicted = predict_xr(model,
                           data,
                           proba=True,
                           persist=True,
                           clean=True,
                           return_input=True
                          ).compute()

    predictions.append(predicted)

working on: 1
predicting...
   probabilities...
   input features...
working on: 2
predicting...
   probabilities...
   input features...
working on: 3
predicting...
   probabilities...
   input features...
working on: 4
predicting...
   probabilities...
   input features...

Tracé des résultats#

Dans les graphiques ci-dessous, vous verrez à gauche l’image classée (vert = culture, blanc = non-culture), au centre une image en vraies couleurs de la région et à droite une image des probabilités de prédiction. Comme nous utilisons un classificateur de forêt aléatoire, les probabilités de prédiction font référence au pourcentage d’arbres qui ont voté pour la classification résultante. Par exemple, si le modèle avait 200 arbres de décision dans la forêt aléatoire et que 150 de ces arbres ont voté « culture », la probabilité de prédiction est de 150 / 200 x 100 = 75 %

[9]:

for i in range(0, len(predictions)):
    fig, axes = plt.subplots(1, 3, figsize=(24, 8))


    # Plot classified image
    predictions[i].Predictions.plot(ax=axes[0],
                   cmap='Greens',
                   add_labels=False,
                   add_colorbar=False)

    # Plot true colour image
    rgb(predictions[i], bands=['red','green','blue'],
        ax=axes[1], percentile_stretch=(0.01, 0.99))

    predictions[i].Probabilities.plot(ax=axes[2],
                   cmap='magma',
                   vmin=0,
                   vmax=100,
                   add_labels=False,
                   add_colorbar=True)

    # Remove axis on right plot
    axes[2].get_yaxis().set_visible(False)

    # Add plot titles
    axes[0].set_title('Classified Image')
    axes[1].set_title('True Colour Geomedian')
    axes[2].set_title('Probabilities');
../../../../_images/sandbox_notebooks_Real_world_examples_Scalable_machine_learning_4_Classify_satellite_data_22_0.png
../../../../_images/sandbox_notebooks_Real_world_examples_Scalable_machine_learning_4_Classify_satellite_data_22_1.png
../../../../_images/sandbox_notebooks_Real_world_examples_Scalable_machine_learning_4_Classify_satellite_data_22_2.png
../../../../_images/sandbox_notebooks_Real_world_examples_Scalable_machine_learning_4_Classify_satellite_data_22_3.png

Classification à grande échelle#

Si vous êtes satisfait des résultats des emplacements de test, essayez de classer une grande région en saisissant à nouveau une nouvelle latitude, une nouvelle longitude et une taille de tampon plus grande. Vous devrez peut-être ajuster la taille de « dask_chunks » pour optimiser la région plus grande. Si vous modifiez la taille des morceaux, n’oubliez pas d’ajuster les morceaux dans la fonction de couche d’entités ci-dessus (c’est-à-dire dans l’exemple par défaut « feature_layers »)

La cellule directement en dessous effacera d’abord les résultats de l’emplacement du test de la mémoire, afin que nous ayons suffisamment de RAM pour exécuter une prédiction beaucoup plus grande.

[10]:

#clear objects from memory
del predictions
del predicted

Entrez un nouvel ensemble de coordonnées et une taille de tampon plus grande ci-dessous. Vous pouvez utiliser la cellule display_map() pour voir un aperçu de la zone que vous allez classer. L’exemple par défaut est centré sur la ville de Faiyum, en Égypte. Essayez de maintenir la taille de la mémoire tampon en dessous de 0,5, au-delà de cette valeur, le bac à sable par défaut commencera à manquer de RAM, ce qui interrompt le calcul.

[11]:

new_lat, new_lon =  29.3037, 30.825 #centred on Faiyum, Egypt
buf_lat, buf_lon = 0.3, 0.5
dask_chunks={'x':1000,'y':1000}

[12]:

display_map((new_lon-buf_lon, new_lon+buf_lon), (new_lat+buf_lat, new_lat-buf_lat))

[12]:
Make this Notebook Trusted to load map: File -> Trust Notebook

Nous allons maintenant classer la région spécifiée ci-dessus :

[13]:

# generate a datacube query object
query = {
    'x': (new_lon-buf_lon, new_lon+buf_lon),
    'y': (new_lat+buf_lat, new_lat-buf_lat),
    'time': time,
    'measurements': measurements,
    'resolution': resolution,
    'output_crs': output_crs,
    'dask_chunks' : dask_chunks,
}

# calculate features lazily
features = feature_layers(query)


#predict using the imported model
predicted = predict_xr(model,
                       features,
                       proba=True,
                       persist=True,
                       clean=True,
                       return_input=True).compute()

predicting...

/opt/venv/lib/python3.12/site-packages/dask/core.py:133: RuntimeWarning: invalid value encountered in divide
  return func(*(_execute_task(a, cache) for a in args))

   probabilities...
   input features...

Écrire les résultats dans GeoTIFF#

Nos prédictions et probabilités de prédiction sont écrites sur le disque sous forme de fichiers GeoTIFF optimisés pour le cloud. En plus d’exporter les prédictions, nous allons également exporter l’une des couches d’entités, NDVI. Dans le prochain bloc-notes, « 5_Object-based_filtering », nous verrons comment utiliser la segmentation d’image pour nettoyer les résultats basés sur les pixels. La couche NDVI fournira une entrée à l’algorithme de segmentation d’image.

[14]:

write_cog(predicted.Predictions, results+'prediction.tif', overwrite=True)
write_cog(predicted.Probabilities, results+'probabilities.tif', overwrite=True)
write_cog(predicted.NDVI, results+'NDVI.tif', overwrite=True)

[14]:

PosixPath('results/NDVI.tif')

Tracer le résultat#

Ci-dessous, nous allons tracer notre masque de terres cultivées basé sur des pixels pour la région autour de la ville de Fayoum, en Égypte.

Remarque : cela pourrait faire planter le noyau si la région est très grande, mais cela devrait fonctionner correctement si vous utilisez l’exemple par défaut.

[15]:

fig, axes = plt.subplots(1, 3, figsize=(30, 10))

# Plot classified image
predicted.Predictions.plot(ax=axes[0],
               cmap='Greens',
               add_labels=False,
               add_colorbar=False)

# Plot true colour image
rgb(predicted,
        ax=axes[1], percentile_stretch=(0.01, 0.99))

predicted.Probabilities.plot(ax=axes[2],
               cmap='magma',
               vmin=0,
               vmax=100,
               add_labels=False,
               add_colorbar=True)

#Remove axis on right plot
axes[1].get_yaxis().set_visible(False)

#Add plot titles
axes[0].set_title('Classified Image')
axes[1].set_title('True Colour Geomedian')
axes[2].set_title('Probabilities');
../../../../_images/sandbox_notebooks_Real_world_examples_Scalable_machine_learning_4_Classify_satellite_data_33_0.png

Conclusions#

Félicitations, vous avez réussi à créer un modèle d’étendue des terres cultivées pour l’Égypte ! Si vous êtes parfaitement satisfait des résultats, cette classification basée sur les pixels peut être le point final de votre flux de travail. Cependant, en réalité, les flux de travail ML comme celui que vous venez de parcourir sont un processus itératif. Si nous ne sommes pas satisfaits des classifications, nous avons quelques options pour améliorer le modèle :

  1. Effectuez une sélection de fonctionnalités pour supprimer celles qui pourraient perturber notre modèle.

  2. Pensez à ajouter de nouvelles fonctionnalités au modèle. Cela nécessiterait de modifier et de réexécuter la fonction « collect_training_data » dans le bloc-notes « Extracting_training_data <1_Extract_training_data.ipynb> » pour ajouter de nouvelles fonctionnalités à notre ensemble de données d’entraînement.

  3. Essayez d’utiliser un modèle différent (par exemple, au lieu d’utiliser un classificateur de forêt aléatoire, nous pourrions utiliser une machine à vecteurs de support - cela nécessitera de modifier et de réexécuter le bloc-notes Evaluate_optimize_fit_classifier).

  4. Collectez davantage de données d’entraînement dans les régions où notre classificateur obtient de mauvais résultats. Cela peut être fait via les plateformes suggérées dans Extracting_training_data.

Nous pouvons également améliorer nos classifications en passant au prochain bloc-notes de cette série. Le prochain bloc-notes explore la conversion de la classification basée sur les pixels en une classification basée sur les objets à l’aide d’un algorithme de segmentation d’image.

Étapes suivantes recommandées#

Pour continuer à travailler sur les blocs-notes de ce flux de travail « Scalable Machine Learning sur l’ODC », accédez au bloc-notes suivant « 5_Object-based_filtering.ipynb ».

  1. Extraction_des_données_de_formation

  2. Inspecter_les_données_de_formation

  3. Évaluer_optimiser_l’ajustement_classificateur

  4. Classifier_les_données_satellites (ce carnet)

  5. Filtrage basé sur les objets

Informations Complémentaires#

Licence Le code de ce notebook est sous licence Apache, version 2.0 <https://www.apache.org/licenses/LICENSE-2.0>`__.

Les données de Digital Earth Africa sont sous licence Creative Commons Attribution 4.0 <https://creativecommons.org/licenses/by/4.0/>.

Contact Si vous avez besoin d’aide, veuillez poster une question sur le canal Slack DE Africa <https://digitalearthafrica.slack.com/> ou sur le GIS Stack Exchange <https://gis.stackexchange.com/questions/ask?tags=open-data-cube> en utilisant la balise open-data-cube (vous pouvez consulter les questions posées précédemment `ici <https://gis.stackexchange.com/questions/tagged/open-data-

Si vous souhaitez signaler un problème avec ce notebook, vous pouvez en déposer un sur Github.

Version de Datacube compatible

[16]:

print(datacube.__version__)

1.8.20

Dernier test :

[17]:

from datetime import datetime
datetime.today().strftime('%Y-%m-%d')

[17]:

'2025-01-17'