GéoMAD roulant
Produits utilisés : gm_s2_rolling, ndvi_anomaly
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 :
Inspection des produits et mesures Rolling GeoMAD disponibles dans le datacube
Chargement des données GeoMAD roulantes
Affichage RVB du GeoMAD Roulant
Calculer et tracer le NDVI à partir du GeoMAD Roulant
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]:
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()

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)

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()

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()

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'