GéoMAD roulant

Mots-clés : données utilisées ; rolling geomad, données utilisées ; sentinel-2

Aperçu

L’imagerie satellite nous permet d’observer la Terre de manière répétitive et détaillée. Cependant, les données manquantes, comme les trous causés par la couverture nuageuse, peuvent rendre difficile la création d’une image complète. Afin de produire une vue unique et complète d’une certaine zone, les données satellite peuvent être consolidées, en empilant les mesures de différents moments dans le temps pour créer une image composite.

Le service Sentinel-2 Rolling GeoMAD (Geomedian et Median Absolute Deviations) de Digital Earth Africa (DE Africa) fournit des données mensuelles sur les géomédianes et les MAD calculés à l’aide d’une fenêtre mobile de 3 mois. Il s’agit d’une série chronologique sans nuage qui peut être utilisée pour surveiller les changements sur une base plus fréquente qu’un produit annuel ou semestriel.

Chaque produit combine des mesures collectées sur une période de 3 mois pour produire une image multispectrale représentative de chaque pixel du continent africain pour chaque mois calendaire. Le résultat final est un ensemble de données complet qui peut être utilisé soit pour générer des images en vraies couleurs pour l’inspection visuelle du paysage, soit pour développer des algorithmes plus complexes.

Détails importants :

  • Noms des produits Datacube : « gm_s2_rolling »

  • Produit de réflectance de surface géomédiane

    • Plage de mise à l’échelle valide : « 1 - 10 000 »

    • « 0 » signifie « aucune donnée »

  • Produit d’écart absolu médian

    • Plage de mise à l’échelle valide : MAD spectral : « 0 - 1 », MAD Bray-Curtis « 0 - 1 », MAD euclidien « 0 - 10 000 »

    • « NaN » signifie « nodata »

  • Statut : provisoire

  • Période : janvier 2019 – aujourd’hui

  • Résolution spatiale : 10 m

Pour plus d’informations sur le service GeoMAD de DE Africa, consultez le service GeoMAD de DE Africa <https://docs.digitalearthafrica.org/en/latest/data_specs/GeoMAD_specs.html>`__.

Description

Dans ce carnet, nous travaillerons avec le Sentinel-2 Rolling GeoMAD et démontrerons comment il peut être utilisé pour surveiller les changements dans les paysages au fil du temps.

Les sujets abordés comprennent :

  1. Inspection des produits et mesures Rolling GeoMAD disponibles dans le datacube

  2. Chargement des données GeoMAD roulantes

  3. Affichage RVB du GeoMAD Roulant

  4. Calculer et tracer le NDVI à partir du GeoMAD Roulant

  5. Inspectez les bandes MAD du GeoMAD ***

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 numpy as np
import geopandas as gpd
import matplotlib.pyplot as plt
import matplotlib.dates as mdates

import datacube
from odc.geo.geom import Geometry
from datetime import datetime

from deafrica_tools.datahandling import load_ard
from deafrica_tools.plotting import rgb, display_map
from deafrica_tools.areaofinterest import define_area
from deafrica_tools.bandindices import calculate_indices

Se connecter au datacube

[2]:
dc = datacube.Datacube(app='rolling_geomad')

Liste des mesures

Inspectez les mesures ou les bandes disponibles pour le Sentinel-2 Rolling GeoMAD à l’aide de la fonctionnalité « list_measurements » de Datacube.

[3]:
product_name = 'gm_s2_rolling'

dc_measurements = dc.list_measurements()
dc_measurements.loc[product_name].drop('flags_definition', axis=1)
[3]:
name dtype units nodata aliases scale_factor
measurement
B02 B02 uint16 1 0.0 [band_02, blue] NaN
B03 B03 uint16 1 0.0 [band_03, green] NaN
B04 B04 uint16 1 0.0 [band_04, red] NaN
B05 B05 uint16 1 0.0 [band_05, red_edge_1] NaN
B06 B06 uint16 1 0.0 [band_06, red_edge_2] NaN
B07 B07 uint16 1 0.0 [band_07, red_edge_3] NaN
B08 B08 uint16 1 0.0 [band_08, nir, nir_1] NaN
B8A B8A uint16 1 0.0 [band_8a, nir_narrow, nir_2] NaN
B11 B11 uint16 1 0.0 [band_11, swir_1, swir_16] NaN
B12 B12 uint16 1 0.0 [band_12, swir_2, swir_22] NaN
SMAD SMAD float32 1 NaN [smad, sdev, SDEV] NaN
EMAD EMAD float32 1 NaN [emad, edev, EDEV] NaN
BCMAD BCMAD float32 1 NaN [bcmad, bcdev, BCDEV] NaN
COUNT COUNT uint16 1 0.0 [count] NaN

Définir la zone d’intérêt

Dans cet exemple, nous allons inspecter une zone de plantations de canne à sucre au Zimbabwe. Le GeoMAD est particulièrement adapté aux applications de classification et de détection de changements, comme dans l’agriculture et la foresterie. Nous verrons plus loin dans le carnet comment le GeoMAD peut être utilisé pour tirer des conclusions sur les pratiques agricoles, comme les dates de récolte.

[4]:
#Specify the latitude, longitude, and buffer
aoi = define_area(lat=-21.07, lon=31.53, buffer=0.03)

#Create a geopolygon and geodataframe of the area of interest
geopolygon = Geometry(aoi["features"][0]["geometry"], crs="epsg:4326")
geopolygon_gdf = gpd.GeoDataFrame(geometry=[geopolygon], crs=geopolygon.crs)

# Get the latitude and longitude range of the geopolygon
lat_range = (geopolygon_gdf.total_bounds[1], geopolygon_gdf.total_bounds[3])
lon_range = (geopolygon_gdf.total_bounds[0], geopolygon_gdf.total_bounds[2])

Afficher la zone d’intérêt avec une carte de base en utilisant display_map()

[5]:
display_map(x=lon_range, y=lat_range)
[5]:
Make this Notebook Trusted to load map: File -> Trust Notebook

Chargez les données Sentinel-2 Rolling GeoMAD à l’aide de dc.load().

Pour plus d’informations sur la manière de charger des données à l’aide du datacube, reportez-vous au bloc-notes Introduction au chargement des données.

Gérer le temps

La nature évolutive de ce produit signifie que le comportement de chargement pour des périodes de temps spécifiées est différent des autres produits Digital Earth Africa. La raison en est qu’il existe un chevauchement temporel entre les images. Ceci est illustré ci-dessous dans le graphique qui conceptualise le comportement de chargement pour la période spécifiée « 2021-01-01 - 2021-03-31 » indiquée par des lignes noires en pointillés. Comme cette période couvre trois mois (janvier, février et mars), les utilisateurs peuvent s’attendre à ce que la requête renvoie trois images. Cependant, dans le cas du GeoMAD évolutif, cela renverra cinq images comme indiqué ci-dessous. C’est parce que la fonction « dc.load » importera toutes les images qui incluent des acquisitions d’images dans la période d’intérêt. Par exemple, le GeoMAD de décembre 2020 comprend des images de novembre et décembre 2020, et de janvier 2021, et sera donc chargé dans le cadre de la requête « 2021-01-01 - 2021-03-31 ». L’étiquette de date dans l’ensemble de données chargé sera « 2020-12-16 ».

[6]:
begin = np.array(["2020-11-01", "2020-12-01", "2021-01-01", "2021-02-01", "2021-03-01"])
end =   np.array(["2021-01-31", "2021-02-28", "2021-03-31", "2021-04-30", "2021-05-31"])

begin = [datetime.strptime(i, "%Y-%m-%d") for i in begin]
end = [datetime.strptime(i, "%Y-%m-%d") for i in end]

event = ["December 2020 GeoMAD", "January 2021 GeoMAD", "February 2021 GeoMAD", "March 2021 GeoMAD", "April 2021 GeoMAD"]

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

plt.barh(range(len(begin)),  [(end[i]-begin[i]) for i in range(len(begin))],
         left=begin, color=['orange', 'green', 'green', 'green', 'orange'])

plt.yticks(range(len(begin)), event)
ax.xaxis.set_major_locator(mdates.MonthLocator(interval=1))
ax.xaxis.set_major_formatter(mdates.DateFormatter("%b %Y"))
plt.setp(ax.get_xticklabels(), rotation=30, ha="right")
ax.axvline(x=datetime.strptime("2021-01-01", "%Y-%m-%d"), linestyle='dashed', color='black')
ax.axvline(x=datetime.strptime("2021-03-31", "%Y-%m-%d"), linestyle='dashed', color='black')
plt.show()
../../../_images/sandbox_notebooks_Datasets_Rolling_GeoMAD_16_0.png

Le comportement de chargement décrit ci-dessus est évident dans le « ds » renvoyé ci-dessous.

[7]:
ds = dc.load(product="gm_s2_rolling",
             measurements=['red','green','blue','nir','emad', 'smad', 'bcmad'],
             x=lon_range,
             y=lat_range,
             resolution=(-20, 20),
             output_crs = 'epsg:6933',
             time=("2021-01","2021-12"),

             )
display(ds)
<xarray.Dataset> Size: 29MB
Dimensions:      (time: 14, y: 359, x: 291)
Coordinates:
  * time         (time) datetime64[ns] 112B 2020-12-16T23:59:59.999999 ... 20...
  * y            (y) float64 3kB -2.626e+06 -2.626e+06 ... -2.633e+06 -2.633e+06
  * x            (x) float64 2kB 3.039e+06 3.039e+06 ... 3.045e+06 3.045e+06
    spatial_ref  int32 4B 6933
Data variables:
    red          (time, y, x) uint16 3MB 595 573 543 851 ... 1027 983 779 648
    green        (time, y, x) uint16 3MB 712 698 656 872 ... 1061 1026 929 886
    blue         (time, y, x) uint16 3MB 458 452 421 538 468 ... 728 704 638 597
    nir          (time, y, x) uint16 3MB 2634 2607 2514 2949 ... 3712 3597 3663
    emad         (time, y, x) float32 6MB 1.858e+03 1.863e+03 ... 1.197e+03
    smad         (time, y, x) float32 6MB 0.02137 0.01984 ... 0.001466 0.002157
    bcmad        (time, y, x) float32 6MB 0.1483 0.1527 ... 0.07304 0.07646
Attributes:
    crs:           epsg:6933
    grid_mapping:  spatial_ref

La cellule ci-dessous découpe les heures de début et de fin pour obtenir la plage horaire spécifiée. Notez que l’objet « ds » ci-dessous inclut uniquement les images étiquetées « 2021 ».

[8]:
ds_time = ds.time[1:-1]
ds = ds.sel(time=ds_time.values)
display(ds)
<xarray.Dataset> Size: 25MB
Dimensions:      (time: 12, y: 359, x: 291)
Coordinates:
  * time         (time) datetime64[ns] 96B 2021-01-14T23:59:59.999999 ... 202...
  * y            (y) float64 3kB -2.626e+06 -2.626e+06 ... -2.633e+06 -2.633e+06
  * x            (x) float64 2kB 3.039e+06 3.039e+06 ... 3.045e+06 3.045e+06
    spatial_ref  int32 4B 6933
Data variables:
    red          (time, y, x) uint16 3MB 332 336 339 678 ... 1053 1007 802 636
    green        (time, y, x) uint16 3MB 583 602 557 781 ... 1027 983 900 838
    blue         (time, y, x) uint16 3MB 317 332 320 451 376 ... 704 678 622 561
    nir          (time, y, x) uint16 3MB 3159 3204 2827 3226 ... 3392 3313 3353
    emad         (time, y, x) float32 5MB 883.5 698.1 756.0 ... 725.1 679.8
    smad         (time, y, x) float32 5MB 0.0004319 0.0004552 ... 0.002771
    bcmad        (time, y, x) float32 5MB 0.0666 0.0551 ... 0.05421 0.04862
Attributes:
    crs:           epsg:6933
    grid_mapping:  spatial_ref

Affichage RVB du Sentinel-2 Rolling GeoMAD

Nous pouvons tracer les données que nous avons chargées à l’aide de la fonction « rgb() ». Par défaut, la fonction tracera les données sous forme d’image en vraies couleurs à l’aide des bandes « rouge », « verte » et « bleue ».

[9]:
rgb(ds, col='time', col_wrap=4, size=3)
../../../_images/sandbox_notebooks_Datasets_Rolling_GeoMAD_22_0.png

Calculer le NDVI en utilisant « calculer les indices »

Reportez-vous au bloc-notes « Calcul des indices de bande <../Frequently_used_code/Calculating_band_indices.ipynb> » pour plus d’informations. >Remarque : ce produit convient au calcul de tous les indices, pas seulement du NDVI.

[10]:
ds_NDVI = calculate_indices(ds, index=['NDVI'], satellite_mission='s2')

Graphique NDVI

Le NDVI est représenté ci-dessous pour chaque mois de 2021. Existe-t-il des tendances ou des changements notables dans la distribution du NDVI dans l’espace ?

[11]:
# NDVI Rolling GeoMAD
ds_NDVI.NDVI.plot(x='x',y='y', col='time', col_wrap = 4, vmin=0, vmax=1,
                                           cmap='RdYlGn', robust=True)

# Title NDVI Rolling GeoMAD
plt.suptitle('NDVI from Rolling GeoMAD', y = 1.05, fontproperties={'style':'oblique', 'weight':'bold'})

plt.show()
../../../_images/sandbox_notebooks_Datasets_Rolling_GeoMAD_26_0.png

Terrain SMAD

Les bandes d’écart absolu médian (MAD) fournissent des informations sur les changements au cours de la période de calcul géomédiane. Cela se fait tout en préservant les relations à haute dimension entre les bandes satellites. Par exemple, des valeurs d’écart absolu médian spectral (SMAD) plus élevées indiquent une plus grande quantité de changement au cours d’une période ; dans ce cas, la fenêtre glissante de trois mois.

Ces informations peuvent aider à identifier les zones les plus dynamiques au cours d’une période donnée. Dans les environnements gérés (par exemple agricoles), cela peut concerner des activités telles que la culture, la récolte ou la plantation de cultures.

En examinant les images ci-dessous, pouvez-vous identifier les champs de la plantation de canne à sucre qui ont subi des changements ? Essayez de comparer ces champs avec les images NDVI et RVB ci-dessus pour voir si vous pouvez tirer des conclusions sur ce qui s’est passé.

[12]:
ds.smad.plot(x='x',y='y', col='time', col_wrap=4, cmap='viridis', robust=True)

plt.suptitle('SMAD from Rolling GeoMAD', y = 1.05, fontproperties={'style':'oblique', 'weight':'bold'})

plt.show()
../../../_images/sandbox_notebooks_Datasets_Rolling_GeoMAD_28_0.png

Conclusion

Ce carnet a démontré le chargement et le traçage du Rolling GeoMAD. Il a également montré comment la nature évolutive de ce produit signifie que ses images se chevauchent dans le temps, ce qui a des implications sur les procédures de chargement.


Informations Complémentaires

Licence : Le code de ce carnet 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 par 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 Open Data Cube <http://slack.opendatacube.org/>`__ ou sur le GIS Stack Exchange en utilisant la balise open-data-cube (vous pouvez consulter les questions posées précédemment ici). Si vous souhaitez signaler un problème avec ce bloc-notes, vous pouvez en déposer un sur Github.

Version de Datacube compatible :

[13]:
print(datacube.__version__)
1.8.20

Dernier test :

[14]:
from datetime import datetime
datetime.today().strftime('%Y-%m-%d')
[14]:
'2025-01-14'