GeoMAD
Produits utilisés : gm_s2_annual, gm_s2_semiannual, gm_ls8_ls9_annual, gm_ls8_annual, gm_ls5_ls7_annual,
Mots-clés : données utilisées ; géomédiane sentinel-2, données utilisées ; géomédiane landsat 8, données utilisées ; géomédiane semi-annuelle sentinel-2, index:données utilisées ; MAD, données utilisées ; géomédiane landsat 5 et 7, données utilisées ; géomédiane landsat 8 et 9
Aperçu
L’imagerie satellite nous permet d’observer la Terre avec une précision et un niveau de détail considérables. 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 doivent être consolidées, en empilant les mesures de différents points dans le temps pour créer une image composite.
Le GeoMAD de Digital Earth Africa (DE Africa) (Geomedian et Median Absolute Deviations) est un composite sans nuages de données satellitaires compilées sur des périodes annuelles et semestrielles (six mois) au cours de chaque année civile.
Les produits GeoMAD suivants sont disponibles sur les plateformes DE Africa :
gm_s2_annual: composite annuel (année civile) GeoMAD utilisant l’imagerie Sentinel-2, disponible pour les années 2017 - présentgm_s2_semiannual: composites biannuels (janvier-juin, juillet-décembre) GeoMAD utilisant l’imagerie Sentinel-2, disponibles pour les années 2017 - présentgm_ls8_ls9_annual: composite annuel (année civile) GeoMAD utilisant des images Landsat-8 et Landsat-9, disponible pour les années 2021 - présentgm_ls8_annual: composite annuel (année civile) GeoMAD utilisant l’imagerie Landsat-8, disponible pour les années 2013 - 2020gm_ls5_ls7_annual: composite annuel GeoMAD (année civile) combinant les images Landsat-5 et Landsat-7, disponible pour les années 1984 - 2012
Chaque produit combine des mesures collectées sur une période définie (annuelle ou semestrielle) pour produire une image multispectrale représentative de chaque pixel du continent africain. 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.
Pour produire un composite GeoMAD, pour chaque pixel, les données non valides sont éliminées et les observations restantes sont résumées mathématiquement à l’aide d’une statistique géomédiane de grande dimension. GeoMAD comprend également trois mesures de l’écart absolu médian (MAD). Il s’agit de mesures statistiques d’ordre supérieur sur la variation par rapport à la géomédiane, précalculées à la même échelle de temps annuelle. Ces couches peuvent être utilisées seules ou avec la géomédiane pour obtenir des informations sur la surface terrestre et comprendre son évolution au fil du temps. Pour une description détaillée de la manière dont le GeoMAD est calculé, consultez les spécifications techniques GeoMAD.
Détails importants :
Noms des produits Datacube : « gm_s2_annual », « gm_s2_semiannual », « gm_ls8_ls9_annual », « gm_ls8_annual », « gm_ls5_ls7_annual »
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 : Opérationnel
Période : 1984 – aujourd’hui
Résolution spatiale : 10 m pour les produits S2, 30 m pour les produits Landsat
Remarque : pour une description détaillée du service GeoMAD de DE Africa, consultez les spécifications techniques GeoMAD de DE Africa <https://docs.digitalearthafrica.org/en/latest/data_specs/GeoMAD_specs.html>`__.
Description
Dans ce bloc-notes, nous allons charger les données GeoMAD à l’aide de dc.load() pour renvoyer une série chronologique d’images satellite. Le xarray.Dataset renvoyé contiendra des images prêtes à être analysées.
Les sujets abordés comprennent :
Inspecting the GeoMAD products and measurements available in the datacube
Using the native
dc.load()function to load in GeoMAD dataGeomedian surface reflectance example
Median Absolute Deviations example
A simple example anlaysis using GeoMAD
An interactive widget for understanding the geomedian statistic
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 numpy as np
import matplotlib.pyplot as plt
from deafrica_tools.plotting import rgb
from deafrica_tools.app import geomedian
Se connecter au datacube
[2]:
dc = datacube.Datacube(app='GeoMAD')
Produits et mesures disponibles
Liste des produits
Nous pouvons utiliser la fonctionnalité « list_products » de Datacube pour inspecter les produits DE Africa GeoMAD disponibles dans le Datacube. Le tableau ci-dessous indique les noms des produits que nous utiliserons pour charger les données, une brève description des données et l’instrument satellite qui a acquis les données.
[3]:
# List GeoMAD products available in DE Africa
dc_products = dc.list_products()
display_columns = ['name', 'description']
dc_products[dc_products.name.str.contains(
'gm').fillna(
False)][display_columns].set_index('name')
[3]:
| description | |
|---|---|
| name | |
| gm_ls5_ls7_annual | Surface Reflectance Annual Geometric Median an... |
| gm_ls5_ls7_annual_lowres | Surface Reflectance Annual Geometric Median an... |
| gm_ls8_annual | Surface Reflectance Annual Geometric Median an... |
| gm_ls8_annual_lowres | Surface Reflectance Annual Geometric Median an... |
| gm_ls8_ls9_annual | Surface Reflectance Annual Geometric Median an... |
| gm_ls8_ls9_annual_lowres | Surface Reflectance Annual Geometric Median an... |
| gm_s2_annual | Surface Reflectance Annual Geometric Median an... |
| gm_s2_annual_lowres | Annual Geometric Median, Sentinel-2 - Low Reso... |
| gm_s2_rolling | Surface Reflectance 3 Monthly Rolling Geometri... |
| gm_s2_semiannual | Surface Reflectance Semiannual Geometric Media... |
| gm_s2_semiannual_lowres | Surface Reflectance Semiannual Geometric Media... |
| gmw | Global Mangrove Watch data sourced from the UN... |
Liste des mesures
Nous pouvons inspecter plus en détail les données disponibles pour GeoMAD en utilisant la fonctionnalité « list_measurements » de Datacube. Le tableau ci-dessous répertorie chacune des mesures disponibles dans les données.
Ci-dessous, saisissez le nom du produit dont vous souhaitez voir les mesures. Dans l’exemple par défaut, nous examinerons « gm_s2_annual », le service Sentinel-2 Annual GeoMAD.
[4]:
product_name = 'gm_s2_annual'
dc_measurements = dc.list_measurements()
dc_measurements.loc[product_name].drop('flags_definition', axis=1)
[4]:
| name | dtype | units | nodata | aliases | |
|---|---|---|---|---|---|
| measurement | |||||
| B02 | B02 | uint16 | 1 | 0.0 | [band_02, blue] |
| B03 | B03 | uint16 | 1 | 0.0 | [band_03, green] |
| B04 | B04 | uint16 | 1 | 0.0 | [band_04, red] |
| B05 | B05 | uint16 | 1 | 0.0 | [band_05, red_edge_1] |
| B06 | B06 | uint16 | 1 | 0.0 | [band_06, red_edge_2] |
| B07 | B07 | uint16 | 1 | 0.0 | [band_07, red_edge_3] |
| B08 | B08 | uint16 | 1 | 0.0 | [band_08, nir, nir_1] |
| B8A | B8A | uint16 | 1 | 0.0 | [band_8a, nir_narrow, nir_2] |
| B11 | B11 | uint16 | 1 | 0.0 | [band_11, swir_1, swir_16] |
| B12 | B12 | uint16 | 1 | 0.0 | [band_12, swir_2, swir_22] |
| SMAD | SMAD | float32 | 1 | NaN | [smad, sdev, SDEV] |
| EMAD | EMAD | float32 | 1 | NaN | [emad, edev, EDEV] |
| BCMAD | BCMAD | float32 | 1 | NaN | [bcmad, bcdev, BCDEV] |
| COUNT | COUNT | uint16 | 1 | 0.0 | [count] |
Regardons également les mesures du Landsat 8-9 GeoMAD, gm_ls8_ls9_annual
[5]:
product_name = 'gm_ls8_ls9_annual'
dc_measurements = dc.list_measurements()
dc_measurements.loc[product_name].drop('flags_definition', axis=1)
[5]:
| name | dtype | units | nodata | aliases | |
|---|---|---|---|---|---|
| measurement | |||||
| SR_B2 | SR_B2 | uint16 | 1 | 0.0 | [band_2, blue] |
| SR_B3 | SR_B3 | uint16 | 1 | 0.0 | [band_3, green] |
| SR_B4 | SR_B4 | uint16 | 1 | 0.0 | [band_4, red] |
| SR_B5 | SR_B5 | uint16 | 1 | 0.0 | [band_5, nir] |
| SR_B6 | SR_B6 | uint16 | 1 | 0.0 | [band_6, swir_1] |
| SR_B7 | SR_B7 | uint16 | 1 | 0.0 | [band_7, swir_2] |
| SMAD | SMAD | float32 | 1 | NaN | [smad, sdev, SDEV] |
| EMAD | EMAD | float32 | 1 | NaN | [emad, edev, EDEV] |
| BCMAD | BCMAD | float32 | 1 | NaN | [bcmad, bcdev, BCDEV] |
| COUNT | COUNT | uint16 | 1 | 0.0 | [count] |
Charger les données GeoMAD à l’aide de dc.load()
Maintenant que nous savons quels produits et mesures sont disponibles pour les produits, nous pouvons charger des données à partir du datacube en utilisant « dc.load ».
Exemple 1 : Réflectivité de surface
Dans le premier exemple ci-dessous, nous allons charger les données GeoMAD de l’Œil de l’Afrique dans le désert du Sahara, en Mauritanie, à partir de 2020. Nous allons charger les données de trois bandes spectrales de satellites : « rouge », « vert » et « bleu ».
Remarque : pour une discussion plus générale sur la manière de charger des données à l’aide du datacube, reportez-vous au bloc-notes « Introduction au chargement des données <../Beginners_guide/03_Loading_data.ipynb> ».
[6]:
# load data
ds = dc.load(product="gm_s2_annual",
measurements=['red', 'green', 'blue'],
x=(-11.50, -11.27),
y=(21.00, 21.20),
time=("2020")
)
print(ds)
<xarray.Dataset>
Dimensions: (time: 1, y: 2385, x: 2220)
Coordinates:
* time (time) datetime64[ns] 2020-07-01T23:59:59.999999
* y (y) float64 2.645e+06 2.645e+06 ... 2.621e+06 2.621e+06
* x (x) float64 -1.11e+06 -1.11e+06 ... -1.087e+06 -1.087e+06
spatial_ref int32 6933
Data variables:
red (time, y, x) uint16 1933 1956 1979 2006 ... 4392 4413 4439 4458
green (time, y, x) uint16 1451 1472 1491 1511 ... 2833 2847 2868 2879
blue (time, y, x) uint16 1039 1054 1064 1073 ... 1695 1709 1719 1724
Attributes:
crs: EPSG:6933
grid_mapping: spatial_ref
Tracé des données géomédianes GeoMAD
Nous pouvons tracer les données que nous avons chargées à l’aide de la fonction « deafrica_tools.plotting.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 ».
[7]:
rgb(ds)
Exemple 2 – Écarts absolus médians
Dans cet exemple, nous chargeons les bandes de données de l’écart absolu médian (MAD) et traçons une carte en fausses couleurs à l’aide de ces trois bandes. La zone sélectionnée est un district agricole autour d’Oum Er-Rbia Rover au Maroc.
[8]:
# load data
lat, lon = 32.4493, -7.4236
buffer = 0.05
ds = dc.load(product="gm_s2_annual",
measurements=['red','green','blue','smad','emad', 'bcmad'],
x= (lon-buffer, lon+buffer),
y=(lat+buffer, lat-buffer),
time=("2020")
)
print(ds)
<xarray.Dataset>
Dimensions: (time: 1, y: 1081, x: 965)
Coordinates:
* time (time) datetime64[ns] 2020-07-01T23:59:59.999999
* y (y) float64 3.932e+06 3.932e+06 ... 3.922e+06 3.922e+06
* x (x) float64 -7.211e+05 -7.211e+05 ... -7.115e+05 -7.115e+05
spatial_ref int32 6933
Data variables:
red (time, y, x) uint16 2940 2948 2945 2934 ... 2940 2903 2879 2862
green (time, y, x) uint16 1798 1802 1803 1795 ... 1916 1893 1876 1858
blue (time, y, x) uint16 1114 1118 1118 1110 ... 1275 1253 1237 1224
smad (time, y, x) float32 0.0002043 0.000197 ... 0.0001852 0.0001878
emad (time, y, x) float32 719.1 734.2 745.5 ... 754.0 740.9 744.1
bcmad (time, y, x) float32 0.03218 0.03213 0.03373 ... 0.0337 0.03397
Attributes:
crs: EPSG:6933
grid_mapping: spatial_ref
Traçage des données GeoMAD MAD
Les données MAD comportent trois bandes (SMAD, EMAD, BCMAD) et sont donc bien adaptées à une visualisation en fausses couleurs. Cela signifie que chacun des MAD est attribué à l’un des canaux de couleur rouge, vert et bleu de l’image.
En inspectant le xarray.Dataset ci-dessus, nous pouvons voir que la mise à l’échelle de « smad (0-1) », « emad (0-10,000) » et « bcmad (0-1) » ont tous des ordres de grandeur très différents. Cela signifie que si nous les traçons comme une image RVB à trois bandes (c’est-à-dire en spécifiant l’argument « rgb(bands=[« emad », « smad », « bcmad »]``) pour les fausses couleurs), les valeurs plus élevées dans « emad » sursatureront l’image (essayez !).
Pour compenser les différentes plages de valeurs dans l’ensemble de données, nous pouvons mettre à l’échelle les données de chacune des trois bandes en fonction de la plage de valeurs présente dans cette bande. Cela amène tous les MAD à peu près à la même plage de valeurs. Le graphique représentera alors chaque MAD de manière plus égale, ce qui permettra d’identifier facilement par inspection visuelle les caractéristiques présentant une forte variabilité dans les trois mesures.
Il existe deux types de mise à l’échelle : fixe et dynamique.
La mise à l’échelle fixe est utile lorsque vous comparez plusieurs zones et que vous souhaitez avoir la même échelle sur chacune d’elles. Elle est utilisée dans la couche WMS GeoMAD et dans le portail Digital Earth Africa Maps pour les données MAD. > Sur le portail Cartes, sélectionnez « Ajouter des données » > « Images satellite » > « Réflectance de surface » > « Annuel » > « GeoMAD annuel (Sentinel-2) » pour afficher les données GeoMAD. Pour obtenir des instructions sur la connexion aux services Web SIG DE Africa, consultez ce didacticiel <https://docs.digitalearthafrica.org/en/latest/platform_tools/index.html>`__.
Dans ce cahier, nous allons démontrer une échelle dynamique. Cette échelle est automatiquement ajustée en fonction de la zone d’intérêt sélectionnée, en fonction de la plage de valeurs de données pour chacun des MAD. Cette méthode est plus adaptée à l’étude d’une certaine zone d’intérêt, car l’échelle est adaptée aux données contenues dans cette zone.
L’échelle dynamique utilise ici une fonction « log » pour transformer chaque point de données MAD. La plage est coupée aux 2e et 98e percentiles, supprimant les valeurs aberrantes extrêmes.
[9]:
smad_scaling = [np.log(ds.smad.quantile(0.02).values), np.log(ds.smad.quantile(0.98).values)]
emad_scaling = [np.log(ds.emad.quantile(0.02).values), np.log(ds.emad.quantile(0.98).values)]
bcmad_scaling = [np.log(ds.bcmad.quantile(0.02).values), np.log(ds.bcmad.quantile(0.98).values)]
Dans ce cas, pour chaque MAD, nous avons utilisé la mise à l’échelle :
\begin{align*} \text{MAD}_{\text{mis à l'échelle}} = \frac{\log{\text{MAD}} - \log \text{MAD}_{\text{2 centile}}}{\log \text{MAD}_{\text{98 centile}} - \log \text{MAD}_{\text{2 centile}}} \end{align*}
Il s’agit d’un exemple de mise à l’échelle qui peut être utilisée pour transformer les données MAD en un ordre de grandeur commun, sans fausser les données. La mise à l’échelle peut être ajustée en fonction de l’objectif ou de l’application des données.
[10]:
ds['smad'] = (np.log(ds.smad)-smad_scaling[0])/(smad_scaling[1]- smad_scaling[0])
ds['emad'] = (np.log(ds.emad)-emad_scaling[0])/(emad_scaling[1]- emad_scaling[0])
ds['bcmad'] = (np.log(ds.bcmad)-bcmad_scaling[0])/(bcmad_scaling[1]- bcmad_scaling[0])
Ci-dessous, nous allons tracer l’image en vraies couleurs de la région à côté de notre image RVB mise à l’échelle des trois mesures MAD. Le graphique ci-dessous indique comment la composition des couleurs RVB peut être interprétée. Les zones où les trois mesures ont des valeurs élevées indiquent les endroits où le pixel subit un degré élevé de variation au cours de l’année. Dans cet exemple, les champs irrigués le long de la rivière sont très variables et apparaissent en blanc sur l’image.
[11]:
fig,ax=plt.subplots(1,2, sharey=True, figsize=(15,7))
rgb(ds,ax=ax[0])
rgb(ds, bands=['smad','emad','bcmad'],ax=ax[1])
ax[0].set_title('True Colour')
ax[1].set_title('Ternary MADs');
Exemple d’application : Identification des zones irriguées à l’aide de SMAD
La section suivante présente un flux de travail d’analyse simple basé sur le produit GeoMAD. Dans cet exemple, nous allons définir un seuil pour la bande SMAD afin de distinguer les champs agricoles irrigués des champs agricoles pluviaux.
Tout d’abord, rechargeons nos données sur la même région que ci-dessus. Nous allons charger les bandes « rouge », « verte » et « bleue », ainsi que « smad »
[12]:
# load data
lat, lon = 32.4493, -7.4236
buffer = 0.05
ds = dc.load(product="gm_s2_annual",
measurements=['red','green','blue','smad'],
x= (lon-buffer, lon+buffer),
y=(lat+buffer, lat-buffer),
time=("2020")
).squeeze()
print(ds)
<xarray.Dataset>
Dimensions: (y: 1081, x: 965)
Coordinates:
time datetime64[ns] 2020-07-01T23:59:59.999999
* y (y) float64 3.932e+06 3.932e+06 ... 3.922e+06 3.922e+06
* x (x) float64 -7.211e+05 -7.211e+05 ... -7.115e+05 -7.115e+05
spatial_ref int32 6933
Data variables:
red (y, x) uint16 2940 2948 2945 2934 2888 ... 2940 2903 2879 2862
green (y, x) uint16 1798 1802 1803 1795 1772 ... 1916 1893 1876 1858
blue (y, x) uint16 1114 1118 1118 1110 1092 ... 1275 1253 1237 1224
smad (y, x) float32 0.0002043 0.000197 ... 0.0001852 0.0001878
Attributes:
crs: EPSG:6933
grid_mapping: spatial_ref
Tracez maintenant la bande « smad » pour voir si elle est utile pour faire la distinction entre les champs irrigués et non irrigués
[13]:
fig,ax=plt.subplots(1,2, sharey=True, figsize=(15,7))
rgb(ds,ax=ax[0])
ds.smad.plot.imshow(vmax=0.025, ax=ax[1])
ax[0].set_title('True Colour')
ax[1].set_title('Spectral MAD');
Il semble que « smad » sera utile pour identifier l’irrigation, définissons donc un seuil à l’aide de la méthode « skimage.filters.threshold_li ». Cette méthode trouvera automatiquement le seuil optimal en minimisant l’entropie croisée entre le premier plan et la moyenne du premier plan, et entre l’arrière-plan et la moyenne de l’arrière-plan
[14]:
from skimage.filters import threshold_li
#find the threshold
threshold = threshold_li(ds.smad.values)
print("Threshold: ", threshold)
#plot a histogram of the result
fig, ax = plt.subplots(figsize=(10, 3))
ds.smad.plot.hist(bins=1000, label="SMAD")
plt.xlim(0,0.003)
ax.axvspan(xmin=0, xmax=threshold, alpha=0.25, color="green", label="Non-irrigated"),
ax.axvspan(xmin=threshold,
xmax=0.005,
alpha=0.25,
color="red",
label="Irrigated")
plt.legend();
Threshold: 0.0016189023
Appliquer le seuil et tracer le résultat
[15]:
#Apply the threshold
irrigated = ds.smad >= threshold
#plot
fig,ax=plt.subplots(1,2, sharey=True, figsize=(15,7))
rgb(ds,ax=ax[0])
irrigated.plot.imshow(ax=ax[1], add_colorbar=False)
ax[0].set_title('True Colour')
ax[1].set_title('Irrigated areas');
Widget interactif : Comprendre la géomédiane
Un exemple du widget interactif géomédien.
Le site Web Digital Earth Africa Docs est en lecture seule. Connectez-vous au Digital Earth Africa Sandbox et accédez au dossier Datasets > GeoMAD pour interagir avec le widget.
Dans cet exemple de widget interactif, vous disposez d’un ensemble de données de satellites d’observation de la Terre. Il contient les bandes rouge, verte et bleue, qui sont les bandes généralement utilisées pour générer des images en couleur. Le widget se concentre sur un seul pixel qui contient des données pour 3 pas de temps différents. Nous pouvons composer (combiner) ces 3 pas de temps en un seul en utilisant une méthode de composition statistique telle que « médiane » ou « géomédiane ».
Ce widget montre que la médiane ne tient pas toujours compte des valeurs de bande très grandes ou très petites et peut ne pas être particulièrement représentative de la variation que nous observons dans les trois pas de temps.
La variation est mieux intégrée dans le résultat « Geomedian RGB - All timesteps », car la formule géomédiane traite chaque pas de temps comme un vecteur multidimensionnel. Nous constatons que cela se traduit par des valeurs différentes entre la médiane et la géomédiane.
[16]:
geomedian.run_app()
[16]:
Par exemple, essayez :
Pas de temps 1 : (0,0,0)
Pas de temps 2 : (0, 255, 0)
Pas de temps 3 : (255, 255, 255)
La géomédiane montre une valeur plus représentative qui intègre une partie de la variation à travers les pas de temps.
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 :
[17]:
print(datacube.__version__)
1.8.15
Dernier test :
[18]:
from datetime import datetime
datetime.today().strftime('%Y-%m-%d')
[18]:
'2023-08-11'