Sentinel-1 Monthly Mosaic

• Produits utilisés : s1_monthly_mosaic

Mots-clés données utilisées ; sentinel-1,:index:ensembles de données ; mosaïque mensuelle sentinel-1, index:SAR

Aperçu

Les capteurs radar à synthèse d’ouverture (SAR) présentent l’avantage de fonctionner à des longueurs d’onde non gênées par la couverture nuageuse et peuvent acquérir des données sur un site de jour comme de nuit. La mission « Sentinel-1 », qui fait partie de l’initiative conjointe Copernicus de la Commission européenne (CE) et de l’Agence spatiale européenne (ESA), assure une surveillance fiable et répétée de vastes zones grâce à son instrument SAR.

Les mosaïques mensuelles Sentinel-1 sont générées à partir de données de rétrodiffusion corrigées du terrain radiométrique (RTC), les variations dues aux variations géométriques d’observation étant atténuées. Dans la mosaïque mensuelle, les images RTC acquises au cours d’un mois calendaire sont combinées à l’aide d’un algorithme de composition multitemporelle. Cet algorithme calcule une moyenne pondérée des pixels valides, attribuant une pondération plus élevée aux pixels présentant une résolution locale plus élevée (par exemple, les pentes opposées au capteur). Cette approche de pondération en résolution locale minimise le bruit et améliore l’homogénéité spatiale des composites.

Vous trouverez plus d’informations sur la manière dont le produit est généré dans la « Documentation de l’écosystème de l’espace de données Copernicus <https://documentation.dataspace.copernicus.eu/Data/SentinelMissions/Sentinel1.html#sentinel-1-level-3-monthly-mosaics> ».

DE Africa donne accès aux mosaïques mensuelles de Sentinel-1 générées en mode d’acquisition interférométrique large (IW). Ces mosaïques sont organisées selon une grille UTM de 100 x 100 km et présentent une résolution spatiale de 20 mètres.

Description

Dans ce notebook, nous allons charger les mosaïques mensuelles Sentinel-1.

Les sujets abordés comprennent :

1. Inspection du produit mensuel de mosaïques Sentinel-1 et des mesures disponibles dans le datacube

2. Utilisation de la fonction « dc.load() » pour charger les mosaïques mensuelles Sentinel-1

3. Visualisation des mosaïques mensuelles

4. Comparaison des mosaïques mensuelles avec le NDVI mensuel et les scènes individuelles Sentinel-1 ***

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

%matplotlib inline
import datacube
import xarray as xr
import matplotlib.pyplot as plt

from deafrica_tools.plotting import rgb
from deafrica_tools.plotting import display_map
from deafrica_tools.datahandling import mostcommon_crs
from deafrica_tools.bandindices import dualpol_indices

Se connecter au datacube

dc = datacube.Datacube(app="Sentinel_1_mosaic")

Produits et mesures disponibles

Liste des produits

Nous pouvons utiliser la fonctionnalité « list_products » de Datacube pour inspecter les produits Sentinel-1 de DE Africa 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 de celles-ci et l’instrument satellite qui les a acquises.

dc.list_products().loc[dc.list_products()['name'].str.contains('s1')]

	name	description	license	default_crs	default_resolution
name
s1_monthly_mosaic	s1_monthly_mosaic	Monthly Sentinel-1 mosaics of weighted terrain...	None	None	None
s1_rtc	s1_rtc	Sentinel 1 Gamma0 normalised radar backscatter	CC-BY-4.0	EPSG:4326	(-0.0002, 0.0002)

product = "s1_monthly_mosaic"

Liste des mesures

Nous pouvons examiner plus en détail les données disponibles pour chaque produit SAR à l’aide de la fonctionnalité « list_measurements » de Datacube. Le tableau ci-dessous répertorie chacune des mesures disponibles dans les données.

measurements = dc.list_measurements()
measurements.loc[product]

	name	dtype	units	nodata	aliases	flags_definition	add_offset	scale_factor
measurement
VV	VV	float32	1	NaN	[vv]	NaN	NaN	NaN
VH	VH	float32	1	NaN	[vh]	NaN	NaN	NaN

Le produit Mosaïque mensuelle Sentinel-1 comporte deux mesures :

• Rétrodiffusion en polarisation « VV », représentant la moyenne pondérée des valeurs de rétrodiffusion « VV » normalisées valides mesurées au cours du mois.

• Rétrodiffusion en polarisation « VH », représentant la moyenne pondérée des valeurs de rétrodiffusion « VH » normalisées valides mesurées au cours du mois.

Les noms de bande à deux lettres indiquent la polarisation de la lumière transmise et reçue par le satellite. VV fait référence au satellite émettant de la lumière polarisée verticalement et recevant en retour de la lumière polarisée verticalement, tandis que VH fait référence au satellite émettant de la lumière polarisée verticalement et recevant en retour de la lumière polarisée horizontalement.

Charger les mosaïques mensuelles Sentinel-1 à l’aide de « dc.load() »

Maintenant que nous savons quels produits et mesures sont disponibles pour le produit, nous pouvons charger des données à partir du datacube en utilisant « dc.load ».

Dans l’exemple ci-dessous, nous allons charger la mosaïque Sentinel-1 pour une zone autour de la ville de Boma sur le fleuve Congo, d’avril 2023 à mars 2024.

Nous allons charger les données de deux bandes de polarisation : « VV » et « VH ». Les données étant fournies sous forme de grille UTM, nous utiliserons la projection la plus courante pour la zone à charger.

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> ».

# Define the area of interest
latitude = -5.87
longitude = 13.06

buffer = 0.05
time = ('2023-04', '2024-03')

bands = ['vv','vh']

#add spatio-temporal extent to the query
query = {
    'x': (longitude-buffer, longitude+buffer),
    'y': (latitude+buffer, latitude-buffer),
    'time':time,
}

Visualiser la zone sélectionnée

display_map(x=(longitude-buffer, longitude+buffer), y=(latitude+buffer, latitude-buffer))

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

output_crs = mostcommon_crs(dc=dc, product=product, query=query)
print(output_crs)

epsg:32733

# loading the data
ds_monthly = dc.load(
    product=product,
    measurements=bands,
    output_crs=output_crs,
    resolution=(-20, 20),
    **query
)

print(ds_monthly)

<xarray.Dataset> Size: 22MB
Dimensions:      (time: 9, y: 556, x: 557)
Coordinates:
  * time         (time) datetime64[ns] 72B 2023-05-16T11:59:59.500000 ... 202...
  * y            (y) float64 4kB 9.356e+06 9.356e+06 ... 9.345e+06 9.345e+06
  * x            (x) float64 4kB 2.796e+05 2.797e+05 ... 2.908e+05 2.908e+05
    spatial_ref  int32 4B 32733
Data variables:
    vv           (time, y, x) float32 11MB 0.07513 0.06036 ... 0.1843 0.1746
    vh           (time, y, x) float32 11MB 0.02182 0.02375 ... 0.03499 0.03216
Attributes:
    crs:           epsg:32733
    grid_mapping:  spatial_ref

Visualiser la mosaïque Sentinel-1

Pour visualiser les données du capteur SAR Sentinel-1, nous calculons le rapport entre les bandes de polarisation verticale et horizontale.

# adding a ratio band
ds_monthly['vh_over_vv'] = ds_monthly.vh/ds_monthly.vv

Les données sont mises à l’échelle en fonction des 2e et 98e percentiles des observations.

data_min = ds_monthly.quantile(0.02)
data_max = ds_monthly.quantile(0.98)

scaled = (ds_monthly - data_min)/(data_max - data_min)

scaled

<xarray.Dataset> Size: 67MB
Dimensions:      (time: 9, y: 556, x: 557)
Coordinates:
  * time         (time) datetime64[ns] 72B 2023-05-16T11:59:59.500000 ... 202...
  * y            (y) float64 4kB 9.356e+06 9.356e+06 ... 9.345e+06 9.345e+06
  * x            (x) float64 4kB 2.796e+05 2.797e+05 ... 2.908e+05 2.908e+05
    spatial_ref  int32 4B 32733
    quantile     float64 8B 0.02
Data variables:
    vv           (time, y, x) float64 22MB 0.1371 0.1089 ... 0.3455 0.3271
    vh           (time, y, x) float64 22MB 0.2926 0.319 0.2788 ... 0.4726 0.4339
    vh_over_vv   (time, y, x) float64 22MB 0.5152 0.7303 ... 0.3057 0.2938

xarray.Dataset

• Dimensions:
  • time: 9
  • y: 556
  • x: 557
• Coordinates: (5)
  • time
    (time)
    datetime64[ns]
    2023-05-16T11:59:59.500000 ... 2...

    units :

    seconds since 1970-01-01 00:00:00

    array(['2023-05-16T11:59:59.500000000', '2023-06-15T23:59:59.500000000',
           '2023-07-16T11:59:59.500000000', '2023-08-16T11:59:59.500000000',
           '2023-09-15T23:59:59.500000000', '2023-11-15T23:59:59.500000000',
           '2023-12-16T11:59:59.500000000', '2024-02-15T11:59:59.500000000',
           '2024-03-16T11:59:59.500000000'], dtype='datetime64[ns]')

  • y
    (y)
    float64
    9.356e+06 9.356e+06 ... 9.345e+06

    units :

    metre

    resolution :

    -20.0

    crs :

    epsg:32733

    array([9356350., 9356330., 9356310., ..., 9345290., 9345270., 9345250.])

  • x
    (x)
    float64
    2.796e+05 2.797e+05 ... 2.908e+05

    units :

    metre

    resolution :

    20.0

    crs :

    epsg:32733

    array([279650., 279670., 279690., ..., 290730., 290750., 290770.])

  • spatial_ref
    ()
    int32
    32733

    spatial_ref :

    PROJCS["WGS 84 / UTM zone 33S",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",15],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32733"]]

    grid_mapping_name :

    transverse_mercator

    array(32733, dtype=int32)

  • quantile
    ()
    float64
    0.02

    array(0.02)

• Data variables: (3)
  • vv
    (time, y, x)
    float64
    0.1371 0.1089 ... 0.3455 0.3271

    array([[[0.13713501, 0.10893147, 0.09377798, ..., 0.16790584,
             0.15433854, 0.15291156],
            [0.13192096, 0.10956809, 0.1490998 , ..., 0.18366764,
             0.19367417, 0.1661963 ],
            [0.15673033, 0.13886682, 0.15002915, ..., 0.15935644,
             0.21536437, 0.1656329 ],
            ...,
            [0.19324526, 0.27276159, 0.19517083, ..., 0.27447007,
             0.13817161, 0.09046013],
            [0.17106056, 0.13999942, 0.17308267, ..., 0.26345037,
             0.15447046, 0.1138465 ],
            [0.1616256 , 0.13647902, 0.19294172, ..., 0.26208826,
             0.25819347, 0.18127564]],
    
           [[0.11750281, 0.10674454, 0.12690276, ..., 0.20834656,
             0.19346831, 0.1590152 ],
            [0.10129372, 0.11069256, 0.10171512, ..., 0.16414703,
             0.19898853, 0.17662186],
            [0.08837233, 0.09909947, 0.07135853, ..., 0.1405912 ,
             0.17678664, 0.14700703],
    ...
            [0.1804262 , 0.2169798 , 0.18348252, ..., 0.15538716,
             0.12977302, 0.08535825],
            [0.14576377, 0.18308927, 0.18340064, ..., 0.2155028 ,
             0.20291469, 0.1085696 ],
            [0.21676859, 0.20489945, 0.24652384, ..., 0.21794521,
             0.41591443, 0.2828878 ]],
    
           [[0.2694684 , 0.23026168, 0.14300277, ..., 0.20031684,
             0.23932512, 0.22795958],
            [0.21429232, 0.2583397 , 0.16097628, ..., 0.20214331,
             0.27657961, 0.29296577],
            [0.2059946 , 0.30739673, 0.21693559, ..., 0.26092077,
             0.30685448, 0.46742875],
            ...,
            [0.24882351, 0.23280466, 0.22342341, ..., 0.21859729,
             0.19733029, 0.15076844],
            [0.2067028 , 0.26285216, 0.24692865, ..., 0.27696271,
             0.25530749, 0.20696173],
            [0.16817103, 0.24322567, 0.24898484, ..., 0.24747912,
             0.3455302 , 0.32712603]]])

  • vh
    (time, y, x)
    float64
    0.2926 0.319 ... 0.4726 0.4339

    array([[[0.29257128, 0.31904277, 0.27884865, ..., 0.44510139,
             0.33013671, 0.28318317],
            [0.22996209, 0.28052898, 0.25367483, ..., 0.39675328,
             0.27884   , 0.31072764],
            [0.21718775, 0.31972101, 0.29981941, ..., 0.49327833,
             0.33566335, 0.35640815],
            ...,
            [0.37637834, 0.32768326, 0.443328  , ..., 0.19911338,
             0.16167978, 0.16627931],
            [0.35054301, 0.31817919, 0.4608546 , ..., 0.32126297,
             0.37346936, 0.28723695],
            [0.32914253, 0.35638997, 0.56976016, ..., 0.3944613 ,
             0.78012585, 0.54864231]],
    
           [[0.31694968, 0.19421092, 0.20321491, ..., 0.65726963,
             0.30784399, 0.34726003],
            [0.24592618, 0.15378199, 0.19399504, ..., 0.64658464,
             0.34593059, 0.37989922],
            [0.16318548, 0.12553514, 0.14834196, ..., 0.49361827,
             0.45260185, 0.35792698],
    ...
            [0.49949008, 0.52130996, 0.4413568 , ..., 0.43702814,
             0.32299955, 0.2548374 ],
            [0.34517144, 0.43084484, 0.43231707, ..., 0.5314072 ,
             0.40971961, 0.2935609 ],
            [0.32592953, 0.39771393, 0.46955556, ..., 0.42294538,
             0.64554536, 0.63628276]],
    
           [[0.4077815 , 0.38151207, 0.387001  , ..., 0.43465443,
             0.33515069, 0.28876927],
            [0.3384429 , 0.27295872, 0.38913311, ..., 0.47581246,
             0.37207186, 0.35730367],
            [0.35884401, 0.25658635, 0.47320658, ..., 0.46106881,
             0.40190282, 0.53545223],
            ...,
            [0.6691938 , 0.5348687 , 0.53400769, ..., 0.46212438,
             0.4029507 , 0.35669021],
            [0.48594473, 0.38034367, 0.57634953, ..., 0.61478119,
             0.45310622, 0.34486767],
            [0.54194253, 0.54088228, 0.67519623, ..., 0.56876084,
             0.47259322, 0.43393432]]])

  • vh_over_vv
    (time, y, x)
    float64
    0.5152 0.7303 ... 0.3057 0.2938

    array([[[0.51524427, 0.73030194, 0.73754442, ..., 0.66333536,
             0.51852669, 0.43808802],
            [0.40615362, 0.62898736, 0.39573896, ..., 0.52669817,
             0.32409921, 0.44391851],
            [0.30781678, 0.56249345, 0.47882069, ..., 0.78692082,
             0.35823242, 0.52305044],
            ...,
            [0.46735529, 0.257691  , 0.55881253, ..., 0.12220207,
             0.24702972, 0.42709936],
            [0.4946557 , 0.55437216, 0.6671701 , ..., 0.26276393,
             0.59652832, 0.61964296],
            [0.49039221, 0.648206  , 0.75089238, ..., 0.34399035,
             0.77507245, 0.77040501]],
    
           [[0.66858203, 0.42456136, 0.36636589, ..., 0.80927948,
             0.36681265, 0.53147733],
            [0.59061892, 0.30655011, 0.4479415 , ..., 1.02423446,
             0.40856202, 0.52361978],
            [0.42898675, 0.27212399, 0.48703386, ..., 0.8996278 ,
             0.63872212, 0.6004268 ],
    ...
            [0.69765241, 0.59715395, 0.5957514 , ..., 0.7071253 ,
             0.6131216 , 0.73727112],
            [0.58164145, 0.58100954, 0.58216842, ..., 0.6149893 ,
             0.48800698, 0.66827897],
            [0.34270097, 0.46601785, 0.45722359, ..., 0.46643678,
             0.35904422, 0.55632187]],
    
           [[0.34644156, 0.38644848, 0.67566117, ..., 0.53042451,
             0.31388998, 0.27582395],
            [0.36408118, 0.21641561, 0.59709729, ..., 0.58249683,
             0.2987002 , 0.26303841],
            [0.40982201, 0.1532794 , 0.53452957, ..., 0.41848586,
             0.28885706, 0.24251408],
            ...,
            [0.68020356, 0.56818589, 0.59397621, ..., 0.51556236,
             0.49417911, 0.58159777],
            [0.58195433, 0.32748122, 0.57918036, ..., 0.54773456,
             0.42053339, 0.38846725],
            [0.82385111, 0.54771344, 0.68656494, ..., 0.56901928,
             0.30571946, 0.29379778]]])

• Indexes: (3)
  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2023-05-16 11:59:59.500000', '2023-06-15 23:59:59.500000',
                   '2023-07-16 11:59:59.500000', '2023-08-16 11:59:59.500000',
                   '2023-09-15 23:59:59.500000', '2023-11-15 23:59:59.500000',
                   '2023-12-16 11:59:59.500000', '2024-02-15 11:59:59.500000',
                   '2024-03-16 11:59:59.500000'],
                  dtype='datetime64[ns]', name='time', freq=None))

  • y
    PandasIndex

    PandasIndex(Index([9356350.0, 9356330.0, 9356310.0, 9356290.0, 9356270.0, 9356250.0,
           9356230.0, 9356210.0, 9356190.0, 9356170.0,
           ...
           9345430.0, 9345410.0, 9345390.0, 9345370.0, 9345350.0, 9345330.0,
           9345310.0, 9345290.0, 9345270.0, 9345250.0],
          dtype='float64', name='y', length=556))

  • x
    PandasIndex

    PandasIndex(Index([279650.0, 279670.0, 279690.0, 279710.0, 279730.0, 279750.0, 279770.0,
           279790.0, 279810.0, 279830.0,
           ...
           290590.0, 290610.0, 290630.0, 290650.0, 290670.0, 290690.0, 290710.0,
           290730.0, 290750.0, 290770.0],
          dtype='float64', name='x', length=557))

• Attributes: (0)

L’attribution des bandes verticales et horizontales respectivement rouges et vertes, et du rapport à la bande bleue, est un moyen courant de visualiser les données SAR.

rgb(scaled, bands=['vv','vh','vh_over_vv'], col='time', col_wrap=6)

Générer une visualisation en fausses couleurs

La cellule ci-dessous manipule les bandes de double polarisation pour produire une visualisation en fausses couleurs basée sur les classes de couverture terrestre attendues, compte tenu des valeurs de bande. Le code est adapté de Markuse, P., présenté dans Sentinel Hub <https://custom-scripts.sentinel-hub.com/custom-scripts/sentinel-1/sar_false_color_visualization-2/>`__. Le script inclut des valeurs de paramètres d’étirement pour divers objectifs de visualisation, notamment « terres plus rouges », « marées noires » et « eaux claires ». Dans ce cas, nous utilisons la spécification « land_green » car notre zone d’intérêt est caractérisée par une végétation dense et verte.

Tout d’abord, des fonctions sont définies pour produire la visualisation.

def sat_enh(rgb_arr, saturation=1.0):
    avg = sum(rgb_arr) / len(rgb_arr)
    return [avg * (1 - saturation) + band * saturation for band in rgb_arr]

def classify_pixel_xr(VV, VH,
                      saturation=1.0,
                      brightness=1.0,
                      water_threshold=25,
                      manual_correction_water=[0.0, 0.0, 0.0],
                      manual_correction_land=[0.0, 0.0, 0.0]):

    def safe_stretch(val, min_val, max_val):
        return ((val - min_val) / (max_val - min_val)).clip(0.0, 1.0)

    water_normal = [
        safe_stretch(VV * 1, 0.0, 0.99),
        safe_stretch(VV * 8, 0.0, 0.99),
        0.5 + VV * 3 + safe_stretch(VH * 2000, 0.0, 0.99)
    ]
    # adjust stretch values to make land greener
    land_green = [
        safe_stretch(VV * 3.0, 0.0, 0.99),
        safe_stretch(VV * 1.4, 0.0, 0.99) + safe_stretch(VH * 8.75, 0.0, 0.99),
        safe_stretch(VH * 1.75, 0.0, 0.99)
    ]

    # Apply manual corrections and brightness
    watervis = [band + corr for band, corr in zip(water_normal, manual_correction_water)]
    landvis = [(band + corr) * brightness for band, corr in zip(land_green, manual_correction_land)]

    # Saturation enhancement
    watervis = sat_enh(watervis, saturation)
    landvis = sat_enh(landvis, saturation)

    # Water mask
    water_mask = (VV / VH) > water_threshold
    result = [xr.where(water_mask, w, l) for w, l in zip(watervis, landvis)]

    # Stack to RGB array
    rgb = xr.concat(result, dim='band')
    rgb = rgb.assign_coords(band=['R', 'G', 'B'])

    return rgb

La fonction est ensuite appliquée à notre ensemble de données mensuelles Sentinel-1 et le résultat est représenté graphiquement. Nous pouvons voir les zones urbaines, la végétation et l’eau en « fausses couleurs » améliorées.

false_col = classify_pixel_xr(ds_monthly.vv, ds_monthly.vh).isel(time=0)
false_col.plot.imshow()

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers). Got range [0.000113667695..298.3589].

<matplotlib.image.AxesImage at 0x7f4ff403a4b0>

Avantages de l’utilisation des mosaïques mensuelles

L’espacement mensuel régulier permet une analyse facile des séries chronologiques, seules ou combinées avec d’autres ensembles de données.

Dans les régions nuageuses comme celle du sud-ouest du Nigéria, les mosaïques radar constituent un moyen fiable de suivre les changements tout au long de l’année, en particulier pendant la saison des pluies, lorsque les cultures poussent activement.

Ci-dessous, nous comparons les mosaïques radar aux mesures mensuelles du NDVI. Dans une grande partie de cette région, les données NDVI sont soit indisponibles certains mois (par exemple, juillet, septembre, novembre) en raison d’une couverture nuageuse persistante, soit affectées par des nuages ​​résiduelles mal filtrées (par exemple, juin, août). Cela conduit à une représentation inexacte des tendances de croissance, comme le montre le graphique ci-dessous.

En revanche, les données radar ne sont pas affectées par la couverture nuageuse et suivent avec succès le développement des cultures tout au long de l’année, fournissant ainsi un ensemble de données plus cohérent et plus fiable pour surveiller les changements agricoles.

Dans les zones moins touchées par les nuages, le radar peut fournir des informations complémentaires aux données optiques, améliorant ainsi les capacités globales de surveillance.

# Define the area of interest
latitude = 6.86
longitude = 3.33

buffer = 0.02

#add spatio-temporal extent to the query
query = {
    'x': (longitude-buffer, longitude+buffer),
    'y': (latitude+buffer, latitude-buffer),
    'time':time,
}

La zone d’intérêt suivante se situe près de Lagos, au Nigéria. Il s’agit principalement de terres agricoles, traversées par une route principale d’est en ouest, à l’extrême nord de la région.

display_map(x=(longitude-buffer, longitude+buffer), y=(latitude+buffer, latitude-buffer))

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

# loading the monthly NDVI data
ds_monthly = dc.load(
    product=product,
    measurements=bands,
    output_crs=output_crs,
    resolution=(-20, 20),
    **query
)

print(ds_monthly)

<xarray.Dataset> Size: 4MB
Dimensions:      (time: 10, y: 232, x: 232)
Coordinates:
  * time         (time) datetime64[ns] 80B 2023-05-16T11:59:59.500000 ... 202...
  * y            (y) float64 2kB 1.078e+07 1.078e+07 ... 1.077e+07 1.077e+07
  * x            (x) float64 2kB -8.005e+05 -8.004e+05 ... -7.959e+05 -7.958e+05
    spatial_ref  int32 4B 32733
Data variables:
    vv           (time, y, x) float32 2MB 0.159 0.1662 0.1793 ... 0.2324 0.2076
    vh           (time, y, x) float32 2MB 0.04548 0.04205 ... 0.04563 0.05365
Attributes:
    crs:           epsg:32733
    grid_mapping:  spatial_ref

# load monthly mean NDVI from the NDVI anomaly product

ndvi_anom = dc.load(
    product="ndvi_anomaly",
    measurements=["ndvi_mean"],
    **query,
)

# visualise monthly NDVI

ndvi_anom.ndvi_mean.plot.imshow(col='time', col_wrap=6, cmap="RdYlGn");

Dans ce cas, nous calculerons l’indice de végétation radar (RVI) qui est l’un des indices de double polarisation disponibles dans le package d’outils DE Africa <https://docs.digitalearthafrica.org/en/latest/sandbox/notebooks/Tools/gen/deafrica_tools.bandindices.html>`__.

# Calculate the chosen radar vegetation index and add it to the loaded data set

ds_monthly = dualpol_indices(ds_monthly, index='RVI')

# Compare monthly NDVI trend and radar measured vegetation trend

fig, ax1 = plt.subplots(figsize=(11,4))
ndvi_anom.ndvi_mean.mean(['x','y']).isel(time=ndvi_anom.ndvi_mean.isnull().mean(['x','y'])<0.5).plot.line('b:^', ax=ax1, label='NDVI');
ax1.set_title('');
plt.legend(loc='lower left');

# Create a second y-axis
ax2 = ax1.twinx()
ds_monthly.RVI.mean(['x', 'y']).plot.line('r:^', ax=ax2, label='RVI');
ax2.set_title('');

plt.legend(loc='upper right');

La tendance du NDVI ci-dessus est peu fiable en raison de données manquantes et de la présence de nuages ​​résiduels, notamment pendant la saison des pluies. On constate également que certains mois manquent dans la série chronologique du NDVI.

Comparaison avec des scènes individuelles de Sentinel-1

Les mosaïques mensuelles offrent un bruit réduit (comme indiqué ci-dessous) et une homogénéité spatiale améliorée par rapport aux scènes individuelles Sentinel-1.

# reduce the time window to one month for this demonstration
query['time'] = ('2023-04')

# loading individual Sentinel-1 scenes
ds_S1 = dc.load(product='s1_rtc',
                measurements=bands + ['mask', 'area'],
                group_by="solar_day",
                output_crs = ds_monthly.attrs['crs'], resolution = (-20,20),
                **query)

print(ds_S1)

<xarray.Dataset> Size: 4MB
Dimensions:      (time: 5, y: 232, x: 232)
Coordinates:
  * time         (time) datetime64[ns] 40B 2023-04-04T05:30:28.075105 ... 202...
  * y            (y) float64 2kB 1.078e+07 1.078e+07 ... 1.077e+07 1.077e+07
  * x            (x) float64 2kB -8.005e+05 -8.004e+05 ... -7.959e+05 -7.958e+05
    spatial_ref  int32 4B 32733
Data variables:
    vv           (time, y, x) float32 1MB 0.1334 0.1753 0.1383 ... 0.1928 0.1374
    vh           (time, y, x) float32 1MB 0.03488 0.03196 ... 0.05957 0.08005
    mask         (time, y, x) uint8 269kB 1 1 1 1 1 1 1 1 1 ... 1 1 1 1 1 1 1 1
    area         (time, y, x) float32 1MB 1.015 1.031 1.019 ... 1.016 1.041
Attributes:
    crs:           epsg:32733
    grid_mapping:  spatial_ref

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8,4))

# Plot individual Sentinel-1 observation
ds_S1.vh.isel(time=0).plot.imshow(robust=True, ax=ax1)

# monthly mosaic
ds_monthly.isel(time=0).vh.plot.imshow(robust=True, ax=ax2)

# title
ax1.set_title('S1 RTC scene')
ax2.set_title('S1 Monthly Mosaic')

Text(0.5, 1.0, 'S1 Monthly Mosaic')

Avec un bruit réduit, les petites caractéristiques sont plus visibles dans la mosaïque mensuelle, notamment la route/autoroute principale dans la partie nord de la zone d’intérêt.

Autres considérations

La combinaison d’acquisitions provenant de différentes directions orbitales et angles de visée permet de réduire les effets d’ombre radar dans les mosaïques. Cependant, en terrain escarpé, certaines zones peuvent être constamment dans l’ombre radar, rendant impossible toute mesure fiable de la rétrodiffusion. Ces zones sont indiquées par NaN dans la mosaïque.

Lorsque la résolution temporelle élevée est cruciale, les scènes individuelles Sentinel-1 offrent une plus grande flexibilité pour l’analyse.

De plus, une approche de pondération de résolution locale peut être appliquée pour combiner des scènes Sentinel-1 de différentes périodes. Pour les paysages à évolution lente, une fenêtre temporelle supérieure à un mois peut être utilisée pour réduire davantage le bruit.

Génération d’une mosaïque pondérée

# apply local resolution weighting to generate a mosaic

weight = 1 / ds_S1.area.where(ds_S1.mask==1)
ds_S1_weighted = (ds_S1[['vv', 'vh']].where(ds_S1.mask==1)*weight/weight.sum(dim='time')).sum(dim='time')

ds_S1_weighted.vh.plot.imshow(robust=True);

La mosaïque générée ci-dessus diffère légèrement du produit mosaïque mensuel en raison des variations dans les méthodes de projection et de rééchantillonnage utilisées lors du traitement.

---

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 :

print(datacube.__version__)

1.8.20

Dernier test :

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

'2025-04-07'