Génération de composites GeoMAD

Mots clés données utilisées ; sentinel-2, données utilisées ; landsat 8, index:données utilisées ; MAD, méthodes de données ; géomédiane, analyse ; composites, dask, méthodes de données ; rééchantillonnage, méthodes de données ; groupby

Aperçu

Les images individuelles de télédétection peuvent être affectées par des données bruyantes, notamment des nuages, des ombres de nuages et de la brume. Pour produire des images plus nettes qui peuvent être comparées plus facilement dans le temps, nous pouvons créer des images « récapitulatives » ou « composites » qui combinent plusieurs images en une seule image pour révéler l’apparence médiane ou « typique » du paysage pendant une certaine période.

Une approche consiste à créer un « géomédian » <https://github.com/daleroberts/hdmedians>. Un « géomédian » est basé sur une statistique de haute dimension appelée « médiane géométrique » (Small 1990) <https://www.jstor.org/stable/1403809), qui échange efficacement une pile temporelle d’observations de mauvaise qualité contre un seul composite de pixels de haute qualité avec un bruit spatial réduit (Roberts et al. 2017) <https://ieeexplore.ieee.org/abstract/document/8004469). Contrairement à une médiane standard, un géomédian maintient la relation entre les bandes spectrales. Cela nous permet d’effectuer des analyses plus approfondies sur les images composites comme nous le ferions sur les images satellites originales (par exemple en permettant le calcul d’indices de bande communs, comme le NDVI).

Le produit GEoMAD de DE Africa (Geomedian et Median Absolute Deviations), en plus de contenir les géomédianes des bandes de réflectance de surface, GeoMAD comprend également trois couches de déviation absolue médiane (MAD). Il s’agit de mesures statistiques d’ordre supérieur sur la variation par rapport à la géomédiane. 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.

Toutes les données de la période sélectionnée doivent être chargées pour calculer un composite, de sorte que le calcul des géomédianes peut nécessiter beaucoup de calculs, en particulier sur de grandes zones ou de longues échelles de temps. Pour faciliter ces analyses, DE Africa héberge un certain nombre de produits geoMAD pré-calculés :

Cela réduit le temps et les ressources nécessaires pour calculer la géomédiane si vous effectuez une analyse sur une échelle de temps annuelle. Des instructions sur la façon d’utiliser la géomédiane de Sentinel-2 GeoMAD se trouvent dans le bloc-notes Datasets/GeoMAD.ipynb.

Pour les analyses sur d’autres échelles de temps, comme l’étude des changements au fil des saisons, il n’est pas possible d’utiliser le produit géomédian annuel. Dans ces cas, il peut être utile de calculer GeoMAD pour cette période spécifique.

Description

Dans ce cahier, nous prendrons des séries chronologiques d’images satellites bruyantes et calculerons un composite GeoMAD qui est en grande partie exempt de nuages et d’autres données bruyantes.

Les calculs GeoMAD sont coûteux en termes de mémoire, de bande passante de données et d’utilisation du processeur. L’ODC a une fonction utile, geomedian_with_mads qui permet à dask d’effectuer le calcul en parallèle sur plusieurs threads pour accélérer les choses. Dans ce notebook, un cluster dask local est utilisé, mais la même approche devrait fonctionner en utilisant un cluster dask plus grand et distribué.

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]:

%matplotlib inline

import datacube
import numpy as np
import geopandas as gpd
import matplotlib.pyplot as plt
from odc.algo import geomedian_with_mads
from datacube.utils.geometry import Geometry
from deafrica_tools.plotting import rgb
from deafrica_tools.datahandling import load_ard
from deafrica_tools.dask import create_local_dask_cluster
from deafrica_tools.areaofinterest import define_area

Configurer un cluster Dask

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

[ ]:

create_local_dask_cluster()

Se connecter au datacube

[3]:

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

Charger Sentinel-2 à partir du datacube

Nous chargeons ici une série chronologique d’images satellites Sentinel-2 masquées par des nuages via l’API datacube à l’aide de la fonction load_ard. Cela nous fournira des données avec lesquelles travailler. Pour limiter le calcul et la mémoire, cet exemple n’utilise que trois bandes optiques (rouge, vert, bleu).

Pour définir la zone d’intérêt, deux méthodes sont disponibles :

  1. En spécifiant la latitude, la longitude et la zone tampon. Cette méthode nécessite que vous saisissiez la latitude centrale, la longitude centrale et la valeur de la zone tampon en degrés carrés autour du point central que vous souhaitez analyser. Par exemple, « lat = 10,338 », « lon = -1,055 » et « buffer = 0,1 » sélectionneront une zone avec un rayon de 0,1 degré carré autour du point avec les coordonnées (10,338, -1,055).

  2. By uploading a polygon as a GeoJSON or Esri Shapefile. If you choose this option, you will need to upload the geojson or ESRI shapefile into the Sandbox using Upload Files button a257f4515f9a4ef893b7e994d8962679 in the top left corner of the Jupyter Notebook interface. ESRI shapefiles must be uploaded with all the related files (.cpg, .dbf, .shp, .shx). Once uploaded, you can use the shapefile or geojson to define the area of interest. Remember to update the code to call the file you have uploaded.

Pour utiliser l’une de ces méthodes, vous pouvez décommenter la ligne de code concernée et commenter l’autre. Pour commenter une ligne, ajoutez le symbole "#" avant le code que vous souhaitez commenter. Par défaut, la première option qui définit l’emplacement à l’aide de la latitude, de la longitude et du tampon est utilisée.

[4]:

# Set the area of interest

# Method 1: Specify the latitude, longitude, and buffer
aoi = define_area(lat=-31.535, lon=18.2702, buffer=0.035)

# Method 2: Use a polygon as a GeoJSON or Esri Shapefile.
# aoi = define_area(vector_path='aoi.shp')

#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])

# Create a reusable query
query = {
    'x': lon_range,
    'y': lat_range,
    'time': ('2019-01', '2019-12'),
    'measurements': ['green',
                     'red',
                     'blue'],
    'resolution': (-10, 10),
    'group_by': 'solar_day',
    'output_crs': 'EPSG:6933'
}

Par rapport à l’utilisation typique de « load_ard » qui renvoie par défaut des données avec des nombres à virgule flottante contenant « NaN » (c’est-à-dire « float32 »), dans cet exemple, nous allons définir le « dtype » sur « native ». Cela conservera nos données dans leur type de données entier d’origine (c’est-à-dire « Int16 »), avec des valeurs nodata marquées avec « -999 ». Cela réduira de moitié la quantité de mémoire occupée par nos données, ce qui peut être extrêmement utile lors de la réalisation d’analyses à grande échelle.

[5]:

# Load available data
ds = load_ard(dc=dc,
              products=['s2_l2a'],
              dask_chunks={'x':1000, 'y':1000},
              mask_filters=(['opening',5], ['dilation',5]), #improve cloud mask
              dtype='native',
              **query)

# Print output data
print(ds)

Using pixel quality parameters for Sentinel 2
Finding datasets
    s2_l2a
Applying morphological filters to pq mask (['opening', 5], ['dilation', 5])
Applying pixel quality/cloud mask
Returning 72 time steps as a dask array
<xarray.Dataset>
Dimensions:      (time: 72, y: 765, x: 676)
Coordinates:
  * time         (time) datetime64[ns] 2019-01-04T08:59:00 ... 2019-12-30T08:...
  * y            (y) float64 -3.824e+06 -3.824e+06 ... -3.831e+06 -3.831e+06
  * x            (x) float64 1.759e+06 1.759e+06 ... 1.766e+06 1.766e+06
    spatial_ref  int32 6933
Data variables:
    green        (time, y, x) uint16 dask.array<chunksize=(1, 765, 676), meta=np.ndarray>
    red          (time, y, x) uint16 dask.array<chunksize=(1, 765, 676), meta=np.ndarray>
    blue         (time, y, x) uint16 dask.array<chunksize=(1, 765, 676), meta=np.ndarray>
Attributes:
    crs:           EPSG:6933
    grid_mapping:  spatial_ref

Tracez les pas de temps en vraies couleurs

Pour visualiser les données, utilisez la fonction utilitaire « rgb » préchargée pour tracer une image en vraies couleurs pour une série d’intervalles de temps. Les zones noires indiquent les endroits où les nuages ou autres pixels non valides de l’image ont été définis sur « -999 » pour indiquer l’absence de données.

Le code ci-dessous tracera trois étapes temporelles de la série temporelle que nous venons de charger.

[6]:

# Set the timesteps to visualise
timesteps = [1, 5, 7]

# Generate RGB plots at each timestep
rgb(ds, index=timesteps)

/usr/local/lib/python3.10/dist-packages/rasterio/warp.py:344: NotGeoreferencedWarning: Dataset has no geotransform, gcps, or rpcs. The identity matrix will be returned.
  _reproject(
../../../_images/sandbox_notebooks_Frequently_used_code_Generating_geomedian_composites_16_1.png

Générer une géomédiane

Comme vous pouvez le voir ci-dessus, la plupart des images satellites auront au moins quelques zones masquées en raison de nuages ou d’autres interférences entre le satellite et le sol. L’exécution de la fonction geomedian_with_mads générera un composite géomédian en combinant toutes les observations de notre xarray.Dataset en une seule image complète (ou presque complète) représentant la médiane géométrique de la période.

Remarque : étant donné que nos données ont été chargées paresseusement avec « dask », l’algorithme géomédian lui-même ne sera pas déclenché tant que nous n’aurons pas appelé la méthode « .compute() » à l’étape suivante.

[7]:

# geomedian = int_geomedian(ds)
geomedian = geomedian_with_mads(ds,
                               compute_mads=False,
                               compute_count=False)

Exécuter le calcul

La méthode « .compute() » déclenchera le calcul de l’algorithme géomédian ci-dessus. L’exécution prendra environ quelques minutes ; consultez le « tableau de bord dask » pour vérifier la progression.

[8]:

%%time
geomedian = geomedian.compute()

CPU times: user 5.17 s, sys: 551 ms, total: 5.72 s
Wall time: 1min 32s

Si nous imprimons notre résultat, vous verrez que la dimension « temps » a maintenant été supprimée et qu’il nous reste une seule image qui représente la médiane géométrique de toutes les images satellites de notre série temporelle initiale :

[9]:

geomedian

[9]:

<xarray.Dataset>
Dimensions:  (y: 765, x: 676)
Coordinates:
  * y        (y) float64 -3.824e+06 -3.824e+06 ... -3.831e+06 -3.831e+06
  * x        (x) float64 1.759e+06 1.759e+06 1.759e+06 ... 1.766e+06 1.766e+06
Data variables:
    green    (y, x) uint16 1425 1471 1510 1525 1456 ... 1881 1846 1859 1851 1821
    red      (y, x) uint16 2219 2265 2335 2352 2230 ... 3178 3145 3154 3112 3059
    blue     (y, x) uint16 811 843 864 876 827 786 ... 1125 1109 1113 1109 1100
Attributes:
    units:         1
    nodata:        0
    crs:           EPSG:6933
    grid_mapping:  spatial_ref

Tracer le composite géomédian

En traçant le résultat, nous pouvons voir que l’image géomédiane est beaucoup plus complète que n’importe laquelle des images individuelles. Nous pouvons également utiliser ces données dans l’analyse en aval, car les relations entre les bandes spectrales sont maintenues par la statistique médiane géométrique.

[10]:

# Plot the result
rgb(geomedian, size=8)
../../../_images/sandbox_notebooks_Frequently_used_code_Generating_geomedian_composites_24_0.png

Générer des écarts absolus médians

En plus d’exécuter un composite géomédian, qui renverra les observations médianes pour chaque bande, nous pouvons également calculer la variation que subit chaque pixel au cours de la série temporelle. La fonction geomedian_with_mads renvoie trois mesures de variabilité : l’écart absolu médian spectral (smad), l’écart absolu médian euclidien (emad) et l’écart absolu médian Bray-Curtis (bcmad).

Pour calculer les trois mesures de variabilité, transmettez simplement « compute_mads=True » à la fonction « geomedian_with_mads ». Cette fonction prend également en charge le renvoi de « count », qui renvoie le nombre d’observations claires pour chaque pixel.

[11]:

geomads = geomedian_with_mads(ds,
                              compute_mads=True,
                              compute_count=True)

geomads = geomads.compute()
print(geomads)

<xarray.Dataset>
Dimensions:  (y: 765, x: 676)
Coordinates:
  * y        (y) float64 -3.824e+06 -3.824e+06 ... -3.831e+06 -3.831e+06
  * x        (x) float64 1.759e+06 1.759e+06 1.759e+06 ... 1.766e+06 1.766e+06
Data variables:
    green    (y, x) uint16 1425 1471 1510 1525 1456 ... 1881 1846 1859 1851 1821
    red      (y, x) uint16 2219 2265 2335 2352 2230 ... 3178 3145 3154 3112 3059
    blue     (y, x) uint16 811 843 864 876 827 786 ... 1125 1109 1113 1109 1100
    smad     (y, x) float32 0.0002015 0.0001259 ... 0.0001583 0.0001578
    emad     (y, x) float32 370.7 378.3 371.5 372.2 ... 287.2 224.4 256.9 321.8
    bcmad    (y, x) float32 0.06463 0.06477 0.06227 ... 0.03043 0.03468 0.04292
    count    (y, x) uint16 48 48 48 48 47 47 47 47 ... 44 44 44 45 45 45 45 45
Attributes:
    units:         1
    nodata:        0
    crs:           EPSG:6933
    grid_mapping:  spatial_ref

Tracer les écarts absolus médians et effacer le décompte

[12]:

fig,ax=plt.subplots(1,4, sharey=True, figsize=(20,5))
geomads['smad'].plot(ax=ax[0], cmap='plasma', robust=True)
geomads['emad'].plot(ax=ax[1], cmap='magma', robust=True)
geomads['bcmad'].plot(ax=ax[2], cmap='cividis', robust=True)
geomads['count'].plot(ax=ax[3], cmap='viridis', robust=True)
ax[0].set_title('SMAD')
ax[1].set_title('EMAD')
ax[2].set_title('BCMAD')
ax[3].set_title('Clear Count');
../../../_images/sandbox_notebooks_Frequently_used_code_Generating_geomedian_composites_28_0.png

Exécution de GeoMADs sur des séries chronologiques groupées ou rééchantillonnées

Dans le bloc-notes Generating composites, des fonctions intégrées telles que mean et median sont exécutées sur des données de séries chronologiques qui ont été rééchantillonnées ou groupées.
Nous pouvons utiliser les mêmes techniques avec la fonction géomédiane.

Rééchantillonnage

Nous allons d’abord diviser la série temporelle en groupes souhaités. Le rééchantillonnage peut être utilisé pour créer un nouvel ensemble de temps à intervalles réguliers :

  • groupé = da_scaled.resample(time=1M)

    • « nD » - nombre de jours (par exemple, « 7D » pour sept jours)

    • « nM » - nombre de mois (par exemple, « 6M » pour six mois)

    • « nY » - nombre d’années (par exemple, « 2Y » pour deux ans)

Grouper par

Le regroupement fonctionne en examinant une partie de la date, mais en ignorant les autres parties. Par exemple, « time.month » regrouperait toutes les données de janvier, quelle que soit l’année à laquelle elles appartiennent.

  • groupé = da_scaled.groupby('time.month')

    • ``”time.day””` - regroupe par jour du mois (1-31)

    • 'time.dayofyear' - regroupe par jour de l’année (1-365)

    • 'time.week' - groupes par semaine (1-52)

    • 'time.month' - groupes par mois (1-12)

    • 'time.season' - regroupe en saisons de 3 mois :

      • 'DJF' Décembre, Janvier, Février

      • 'MAM' mars, avril, mai

      • « JJA » juin, juillet, août

      • 'SON' Septembre, octobre, novembre

    • 'time.year' - regroupe par année

Ici, nous allons rééchantillonner en trois groupes de deux mois (« 2M »), le groupe commençant au début du mois (représenté par le « S » à la fin).

[13]:

grouped = ds.resample(time='2MS')
grouped

[13]:

DatasetResample, grouped over '__resample_dim__'
6 groups with labels 2019-01-01, ..., 2019-11-01.

Nous allons maintenant appliquer la fonction « geomedian_with_mads » à chaque groupe rééchantillonné en utilisant la méthode « map ».

Au lieu d’appeler « geomedian_with_mads(ds) » sur l’ensemble du tableau, nous passons la fonction « geomedian_with_mads » à « map » pour l’appliquer séparément à chaque groupe rééchantillonné.

[14]:

geomedian_grouped = grouped.map(geomedian_with_mads)

Nous pouvons maintenant déclencher le calcul et observer la progression à l’aide du tableau de bord Dask.

[15]:

geomedian_grouped = geomedian_grouped.compute()

Nous pouvons tracer les géomédianes de sortie et voir l’évolution du paysage au cours de l’année :

[16]:

rgb(geomedian_grouped, col='time', col_wrap=3)
../../../_images/sandbox_notebooks_Frequently_used_code_Generating_geomedian_composites_36_0.png

Avancé : GeoMADs sur données Landsat

Ci-dessous, nous calculons les GeoMADs à l’aide de Landsat 8. Pour que le GeoMAD corresponde à Sentinel-2, nous devons franchir une étape supplémentaire. Les données de réflectance de surface Landsat sont mises à l’échelle de 0 à 1, tandis que les données Sentinel-2 sont mises à l’échelle de 0 à 10 000. Pour que les deux ensembles de données correspondent, nous devons multiplier les bandes Landsat par 10 000 (cela fera que les geoMADs produits dans ce bloc-notes correspondront à l’ensemble de données « gm_ls8_annual » disponible dans le datacube). Les bandes « bcmad » et « smad » n’ont pas besoin d’être multipliées par 10 000.

[17]:

# Load available data for Landsat 8
ds = load_ard(dc=dc,
              x=lon_range,
              y=lat_range,
              time=('2019-01', '2019-12'),
              measurements=['green','red','blue'],
              products=['ls8_sr'],
              dask_chunks={'x':1000, 'y':1000},
              mask_filters=(['opening',3], ['dilation',3]), #improve cloud mask
              resolution=(-30,30),
              group_by='solar_day',
              output_crs='EPSG:6933'
             )

# Print output data
print(ds)

Using pixel quality parameters for USGS Collection 2
Finding datasets
    ls8_sr
Applying morphological filters to pq mask (['opening', 3], ['dilation', 3])
Applying pixel quality/cloud mask
Re-scaling Landsat C2 data
Returning 23 time steps as a dask array
<xarray.Dataset>
Dimensions:      (time: 23, y: 256, x: 226)
Coordinates:
  * time         (time) datetime64[ns] 2019-01-08T08:40:51.853294 ... 2019-12...
  * y            (y) float64 -3.824e+06 -3.824e+06 ... -3.831e+06 -3.831e+06
  * x            (x) float64 1.759e+06 1.759e+06 ... 1.766e+06 1.766e+06
    spatial_ref  int32 6933
Data variables:
    green        (time, y, x) float32 dask.array<chunksize=(1, 256, 226), meta=np.ndarray>
    red          (time, y, x) float32 dask.array<chunksize=(1, 256, 226), meta=np.ndarray>
    blue         (time, y, x) float32 dask.array<chunksize=(1, 256, 226), meta=np.ndarray>
Attributes:
    crs:           EPSG:6933
    grid_mapping:  spatial_ref

Calculer les GeoMADs sur les données Landsat

[18]:

geomads = geomedian_with_mads(ds, compute_mads=True, compute_count=True).compute()
print(geomads)

<xarray.Dataset>
Dimensions:  (y: 256, x: 226)
Coordinates:
  * y        (y) float64 -3.824e+06 -3.824e+06 ... -3.831e+06 -3.831e+06
  * x        (x) float64 1.759e+06 1.759e+06 1.76e+06 ... 1.766e+06 1.766e+06
Data variables:
    green    (y, x) float32 0.1419 0.1439 0.1352 0.1399 ... 0.1597 0.1509 0.1506
    red      (y, x) float32 0.2158 0.2185 0.2056 0.2149 ... 0.2627 0.2512 0.2384
    blue     (y, x) float32 0.07643 0.07741 0.07278 ... 0.0825 0.07558 0.07652
    smad     (y, x) float32 1.876e-05 3.31e-05 3.873e-05 ... 0.0002149 0.0003207
    emad     (y, x) float32 0.03394 0.02931 0.03076 ... 0.01675 0.02122 0.01753
    bcmad    (y, x) float32 0.05951 0.0552 0.05722 ... 0.02686 0.03528 0.03135
    count    (y, x) uint16 15 15 14 14 14 14 14 14 ... 14 14 14 14 14 14 14 14
Attributes:
    crs:           EPSG:6933
    grid_mapping:  spatial_ref

Multipliez les bandes par 10 000 pour correspondre à Sentinel-2

[19]:

exclude = ['bcmad', 'smad', 'count']

for band in geomads.data_vars:
    if band not in exclude:
        geomads[band] = geomads[band]*10000

Tracer le Landsat GeoMAD

[20]:

fig,ax=plt.subplots(1,5, sharey=True, figsize=(22,5))
rgb(geomads, ax=ax[0])
geomads['smad'].plot(ax=ax[1], cmap='plasma', robust=True)
geomads['emad'].plot(ax=ax[2], cmap='magma', robust=True)
geomads['bcmad'].plot(ax=ax[3], cmap='cividis', robust=True)
geomads['count'].plot(ax=ax[4], cmap='viridis', robust=True)
ax[0].set_title('Geomedian')
ax[1].set_title('SMAD')
ax[2].set_title('EMAD')
ax[3].set_title('BCMAD')
ax[4].set_title('Clear Count');
../../../_images/sandbox_notebooks_Frequently_used_code_Generating_geomedian_composites_44_0.png

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 :

[21]:

print(datacube.__version__)

1.8.15

Dernier test :

[22]:

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

[22]:

'2023-08-14'