Phénologie de la végétation

• Products used: s2_l2a_c1, crop_mask

Keywords data used; sentinel-2 c1, data used; crop_mask, band index; NDVI, band index; EVI, phenology, analysis; time series

Aperçu

La phénologie est l’étude des cycles de vie des plantes et des animaux dans le contexte des saisons. Elle peut être utile pour comprendre les tendances du cycle de vie des cultures et la manière dont les saisons de croissance sont affectées par les changements climatiques. Pour plus d’informations, consultez la page de l’USGS sur la dérivation de la phénologie à partir des séries chronologiques NDVI <https://www.usgs.gov/land-resources/eros/phenology/science/deriving-phenological-metrics-ndvi?qt-science_center_objects=0#qt-science_center_objects>`__.

Description

Ce bloc-notes montre comment calculer les statistiques de phénologie de la végétation à l’aide de la fonction DE Africa xr_phenology. Pour détecter les changements dans la vie végétale pour Sentinel-2, le script utilise soit l”indice de végétation par différence normalisée (NDVI), soit l”indice de végétation amélioré (EVI), qui sont des indicateurs courants de la croissance et de la santé de la végétation.

Les résultats de ce carnet peuvent être utilisés pour évaluer les différences spatio-temporelles dans les saisons de croissance des champs agricoles ou de la végétation indigène.

Ce cahier illustre les étapes suivantes :

1. Load cloud-masked Sentinel 2 c1 data for an area of interest.

2. Calculer un indice proxy de végétation (NDVI ou EVI).

3. Générer une série chronologique zonale de la santé de la végétation.

4. Complétez et lissez la série temporelle de la végétation pour supprimer les lacunes et le bruit.

5. Calculer les statistiques phénologiques sur une série chronologique de végétation 1D simple.

6. Calculer les statistiques de phénologie par pixel.

7. Optionnel : calcul de statistiques temporelles génériques à l’aide de la bibliothèque hdstats.

---

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

Chargez les principaux packages Python et les fonctions de support pour l’analyse.

%matplotlib inline

import datetime as dt
from datetime import datetime
import os

import datacube
import geopandas as gpd
import matplotlib as mpl
import matplotlib.pyplot as plt
import mpl_toolkits.axisartist as AA
import numpy as np
import pandas as pd
import xarray as xr
from datacube.utils.aws import configure_s3_access
from odc.geo.geom import Geometry
from mpl_toolkits.axes_grid1 import host_subplot

from deafrica_tools.areaofinterest import define_area
from deafrica_tools.bandindices import calculate_indices
from deafrica_tools.datahandling import load_ard
from deafrica_tools.plotting import display_map, rgb
from deafrica_tools.spatial import xr_rasterize
from deafrica_tools.temporal import temporal_statistics, xr_phenology

configure_s3_access(aws_unsigned=True, cloud_defaults=True)

Se connecter au datacube

Connectez-vous au datacube pour que nous puissions accéder aux données de DE Africa. Le paramètre « app » est un nom unique pour l’analyse qui est basé sur le nom du fichier du notebook.

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

Paramètres d’analyse

La cellule suivante définit des paramètres importants pour l’analyse :

• « veg_proxy » : indice de bande à utiliser comme proxy pour la santé de la végétation, par exemple « NDVI » ou « EVI ».

• « lat » : la latitude centrale à analyser (par exemple « -10,6996 »).

• « lon » : la longitude centrale à analyser (par exemple « 35,2708 »).

• « buffer » : le nombre de degrés carrés à charger autour de la latitude et de la longitude centrales. Pour des temps de chargement raisonnables, définissez cette valeur sur « 0,1 » ou moins.

• time_range : la plage d’années à analyser (par exemple ('2019-01', '2019-06')).

Sélectionnez l’emplacement

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

   Vous pouvez également fournir des valeurs de tampon distinctes pour la latitude et la longitude d’une zone rectangulaire. Par exemple, « lat = 10.338 », « lon = -1.055 » et « lat_buffer = 0.1 » et « lon_buffer = 0.08 » sélectionneront une zone rectangulaire s’étendant sur 0,1 degré au nord et au sud et 0,08 degré à l’est et à l’ouest à partir du point « (10.338, -1.055) ».

   Pour des temps de chargement raisonnables, définissez la mémoire tampon sur « 0,1 » ou moins.

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

# Set the vegetation proxy to use
veg_proxy = "NDVI"

time_range = ("2019-01-01", "2020-12-20")

#To use crop mask to filter result
use_crop_mask = False

# Method 1: Specify the latitude, longitude, and buffer
aoi = define_area(lat=9.775557, lon=8.947609, buffer=0.02)

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

Afficher l’emplacement sélectionné

La cellule suivante affichera la zone sélectionnée sur une carte interactive. N’hésitez pas à zoomer et dézoomer pour mieux comprendre la zone que vous allez analyser. En cliquant sur n’importe quel point de la carte, vous découvrirez les coordonnées de latitude et de longitude de ce point.

display_map(x=lon_range, y=lat_range)

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

Charger les données Sentinel-2 masquées par le cloud

La première étape consiste à charger les données Sentinel-2 pour la zone d’intérêt et la plage de temps spécifiées. La fonction « load_ard » est utilisée ici pour charger les données qui ont été masquées pour les filtres de nuages, d’ombre et de qualité, les rendant ainsi prêtes à être analysées.

# Create a reusable query
query = {
    "y": lat_range,
    "x": lon_range,
    "time": time_range,
    "measurements": ["red", "green", "blue", "nir"],
    "resolution": (-20, 20),
    "output_crs": "epsg:6933",
    "group_by": "solar_day",
}

# Load available data from Sentinel-2
ds = load_ard(
    dc=dc,
    products=["s2_l2a_c1"],
    mask_filters=[("opening", 3), ("dilation", 3)],
    **query,
)

print(ds)

Using pixel quality parameters for Sentinel 2
Finding datasets
    s2_l2a_c1
Applying morphological filters to pq mask [('opening', 3), ('dilation', 3)]
Applying pixel quality/cloud mask
Re-scaling Sentinel-2 C1 data
Loading 126 time steps
<xarray.Dataset> Size: 99MB
Dimensions:      (time: 126, y: 252, x: 194)
Coordinates:
  * time         (time) datetime64[ns] 1kB 2019-03-26T09:58:13.015000 ... 202...
  * y            (y) float64 2kB 1.244e+06 1.244e+06 ... 1.239e+06 1.239e+06
  * x            (x) float64 2kB 8.614e+05 8.614e+05 ... 8.652e+05 8.652e+05
    spatial_ref  int32 4B 6933
Data variables:
    red          (time, y, x) float32 25MB 1.322e+03 1.828e+03 ... 1.794e+03
    green        (time, y, x) float32 25MB 957.0 1.334e+03 ... 1.158e+03
    blue         (time, y, x) float32 25MB 514.0 750.0 524.0 ... 806.0 806.0
    nir          (time, y, x) float32 25MB 2.104e+03 2.564e+03 ... 2.81e+03
Attributes:
    crs:           epsg:6933
    grid_mapping:  spatial_ref

Découpez les ensembles de données selon la forme de la zone d’intérêt

Un géopolygone représente les limites et non la forme réelle, car il est conçu pour représenter l’étendue de l’entité géographique cartographiée, plutôt que sa forme exacte. En d’autres termes, le géopolygone est utilisé pour définir la limite extérieure de la zone d’intérêt, plutôt que les entités et caractéristiques internes.

Il est important de découper les données selon la forme exacte de la zone d’intérêt, car cela permet de garantir que les données utilisées sont pertinentes pour la zone d’étude spécifique. Bien qu’un géopolygone fournisse des informations sur la limite de l’entité géographique représentée, il ne reflète pas nécessairement la forme ou l’étendue exacte de la zone d’intérêt.

# Rasterise the area of interest polygon
aoi_raster = xr_rasterize(gdf=geopolygon_gdf, da=ds, crs=ds.crs)
# Mask the dataset to the rasterised area of interest
ds = ds.where(aoi_raster == 1)

Masquer la région avec la carte de l’étendue des terres cultivées de l’Afrique de l’Est

Chargez le masque de terres cultivées sur la région d’intérêt.

if use_crop_mask:
    cm = dc.load(
    product="crop_mask",
    time=("2019"),
    measurements="filtered",
    resampling="nearest",
    like=ds.geobox,
    ).filtered.squeeze()

    cm.where(cm < 255).plot.imshow(
        add_colorbar=False, figsize=(6, 6)
    )  # we filter to <255 to omit missing data
    plt.title("Cropland Extent");

Nous allons maintenant utiliser la carte des terres cultivées pour masquer les régions des données Sentinel-2 qui ne contiennent que des cultures.

# Filter out the no-data pixels (255) and non-crop pixels (0) from the cropland map and
# mask the Sentinel-2 data.
if use_crop_mask:
    ds = ds.where(cm == 1)

Une fois le chargement terminé, nous pouvons tracer certaines images sous forme d’image en vraies couleurs en utilisant la fonction « rgb ».

rgb(ds, index=[1, 22, 44, 70], col_wrap=4, size=3)

Calculer les indices de bande

Cette étude mesure la présence de végétation soit par l’intermédiaire de « l’indice de végétation à différence normalisée (NDVI) », soit par l’intermédiaire de « l’indice de végétation amélioré (EVI) ». L’indice qui sera utilisé est déterminé par le paramètre « veg_proxy » qui a été défini dans la section « Paramètres d’analyse ».

L’indice de végétation différentiel normalisé (NDVI) nécessite les bandes « rouge » et « nir » (proche infrarouge). La formule est

\[\begin{split}\begin{aligned} \text{NDVI} & = \frac{(\text{NIR} - \text{Rouge})}{(\text{NIR} + \text{Rouge})} \\ \end{aligned}\end{split}\]

L’indice de végétation amélioré nécessite les bandes « rouge », « nir » et « bleue ». La formule est

\[\begin{split}\begin{aligned} \text{EVI} & = \frac{2,5 \fois (\text{NIR} - \text{Rouge})}{(\text{NIR} + 6 \fois \text{Rouge} - 7,5 \fois \text{Bleu} + 1)} \\ \end{aligned}\end{split}\]

Les deux indices sont disponibles via la fonction calculate_indices, importée depuis deafrica_tools.bandindices. Ici, nous utilisons satellite_mission='s2' puisque nous travaillons avec les données Sentinel-2.

# Calculate the chosen vegetation proxy index and add it to the loaded data set
ds = calculate_indices(ds, index=veg_proxy, satellite_mission="s2")

L’indice proxy de végétation doit maintenant apparaître comme une variable de données, avec les mesures chargées, dans l’objet « ds ».

Tracer l’indice de végétation au fil du temps

Pour avoir une idée de l’évolution de la santé de la végétation au cours des années, nous pouvons tracer une série chronologique zonale sur la région d’intérêt. Nous allons d’abord réaliser un tracé simple de la moyenne zonale des données.

ds.NDVI.mean(["x", "y"]).plot.line("b-^", figsize=(11, 4))
plt.title("Zonal mean of vegetation timeseries");

Lissage/interpolation de séries chronologiques de végétation

Ici, nous allons lisser et interpoler les données pour garantir que nous travaillons avec une série chronologique cohérente. Il s’agit d’une étape très importante du flux de travail et il existe de nombreuses façons de lisser, d’interpoler, de combler les lacunes, de supprimer les valeurs aberrantes ou d’ajuster les courbes des données pour garantir une série chronologique utilisable. Si vous n’utilisez pas l’exemple par défaut, vous devrez peut-être définir des méthodes supplémentaires à celles utilisées ici.

Pour ce faire, nous procédons en deux étapes :

1. Rééchantillonner les données par pas de temps bimensuels en utilisant la médiane bimensuelle

2. Calculer une moyenne mobile avec une fenêtre de 4 étapes

resample_period = "2W"
window = 4

veg_smooth = (
    ds[veg_proxy]
    .resample(time=resample_period)
    .median()
    .rolling(time=window, min_periods=1)
    .mean()
)

Alternativement, complétons la série chronologique en utilisant la fonction DE Africa deafrica_temporal_statistics.fast_complete()

veg_smooth_1D = veg_smooth.mean(["x", "y"])
veg_smooth_1D.plot.line("b-^", figsize=(15, 5))
_max = veg_smooth_1D.max()
_min = veg_smooth_1D.min()
plt.vlines(np.datetime64("2019-01-01"), ymin=_min, ymax=_max)
plt.vlines(np.datetime64("2020-01-01"), ymin=_min, ymax=_max)
plt.vlines(np.datetime64("2021-01-01"), ymin=_min, ymax=_max)
plt.title(veg_proxy + " time-series, year start/ends marked with vertical lines")
plt.ylabel(veg_proxy);

Calculer les statistiques phénologiques à l’aide de « xr_phenology »

La fonction DE Africa « xr_phenology » peut calculer un certain nombre de statistiques sur la phénologie de la surface terrestre qui décrivent ensemble les caractéristiques du cycle de vie d’une plante. La fonction peut calculer les statistiques suivantes soit sur une série temporelle zonale (comme celle ci-dessus), soit sur une base par pixel :

SOS = Date of start of season
POS = Date of peak of season
EOS = Date of end of season
vSOS = Value at start of season
vPOS = Value at peak of season
vEOS = Value at end of season
Trough = Minimum value of season
LOS = Length of season (Days)
AOS = Amplitude of season (in value units)
ROG = Rate of greening
ROS = Rate of senescence

Par défaut, la fonction renverra toutes les statistiques sous la forme d’un « xarray.Dataset », pour renvoyer uniquement un sous-ensemble de ces statistiques, transmettez une liste des statistiques souhaitées à la fonction, par exemple « stats=[“SOS”, “EOS”, “ROG”] ».

La fonction « xr_phenology » permet également d’interpoler et/ou de lisser les séries chronologiques de la même manière que nous l’avons fait ci-dessus, l’interpolation/le lissage se produira avant le calcul des statistiques.

Consultez le script deafrica_tools.temporal pour plus d’informations sur chacun des paramètres de xr_phenology.

Statistiques de phénologie zonale

Pour nous aider à comprendre à quoi ces statistiques font référence, passons d’abord la série chronologique moyenne zonale plus simple (moyenne de tous les pixels de l’image) à la fonction et traçons les résultats sur les mêmes courbes que ci-dessus.

Tout d’abord, fournissez une liste de statistiques à calculer avec le paramètre « pheno_stats ».

method_sos : Si “first” alors vSOS est estimé comme la première pente positive du côté verdoyant de la courbe. Si “median”, alors vSOS est estimé comme la valeur médiane des pentes positives du côté verdoyant de la courbe.

method_eos : Si “last” alors vEOS est estimé comme la dernière pente négative du côté sénescent de la courbe. Si “median”, alors vEOS est estimé comme la valeur “médiane” des pentes négatives du côté sénescent de la courbe.

basic_pheno_stats = [
    "SOS",
    "vSOS",
    "POS",
    "vPOS",
    "EOS",
    "vEOS",
    "Trough",
    "LOS",
    "AOS",
    "ROG",
    "ROS",
]
method_sos = "first"
method_eos = "last"

Tracez les résultats avec nos statistiques annotées sur le graphique

# find all the years to assist with plotting
years = veg_smooth_1D.groupby("time.year")

# store results in dict
pheno_results = {}

# loop through years and calculate phenology
for y, year in years:
    # calculate stats
    stats = xr_phenology(
        year,
        method_sos=method_sos,
        method_eos=method_eos,
        stats=basic_pheno_stats,
        verbose=False,
    )
    # add results to dict
    pheno_results[str(y)] = stats

df_dict = {}
for key, value in pheno_results.items():
    df_dict_1 = {}
    for b in value.data_vars:
        if value[b].dtype == np.dtype("<M8[ns]") or value[b].dtype == np.dtype("int16"):
            result = pd.to_datetime(value[b].values)
        else:
            result = round(float(value[b].values), 3)
        df_dict_1[b] = result
    df_dict[key] = df_dict_1

df = (pd.DataFrame(df_dict)).T
df

	SOS	vSOS	POS	vPOS	EOS	vEOS	Trough	LOS	AOS	ROG	ROS
2019	2019-04-14 00:00:00	0.085	2019-09-01 00:00:00	0.456	2019-12-22 00:00:00	0.237	0.085	252.0	0.371	0.003	-0.002
2020	2020-04-12 00:00:00	0.076	2020-09-13 00:00:00	0.381	2020-12-20 00:00:00	0.193	0.075	252.0	0.306	0.002	-0.002

Tracez les résultats avec nos statistiques annotées sur le graphique

# Figure to which the subplot will be added
fig = plt.figure(figsize=(15, 7))
# Create a subplot that can act as a host to parasitic axes
host = host_subplot(111, figure=fig, axes_class=AA.Axes)
# fig, ax = plt.subplots()


# Function to use to edit the axes of the plot
def adjust_axes(ax):
    # Set the location of the major and minor ticks.
    ax.xaxis.set_major_locator(mpl.dates.MonthLocator())
    ax.xaxis.set_minor_locator(mpl.dates.MonthLocator(bymonthday=16))

    # Format the major and minor tick labels.

    ax.xaxis.set_major_formatter(mpl.ticker.NullFormatter())
    ax.xaxis.set_minor_formatter(mpl.dates.DateFormatter("%b"))


#     # Turn off unnecessary ticks.
    ax.axis["bottom"].minor_ticks.set_visible(False)

    ax.axis["top"].major_ticks.set_visible(False)
    ax.axis["top"].minor_ticks.set_visible(False)

    ax.axis["right"].major_ticks.set_visible(False)
    ax.axis["right"].minor_ticks.set_visible(False)


# Find all the years to assist with plotting
years = veg_smooth_1D.groupby("time.year")

# Counter to aid in plotting.
counter = 0

for y, year in years:
    # Grab all the values we need for plotting.
    eos = df.loc[str(y)].EOS
    sos = df.loc[str(y)].SOS
    pos = df.loc[str(y)].POS

    veos = df.loc[str(y)].vEOS
    vsos = df.loc[str(y)].vSOS
    vpos = df.loc[str(y)].vPOS

    if counter == 0:
        ax = host
    else:
        # Create the secondary axis.
        ax = host.twiny()

    # Plot the data
    year.plot(ax=ax, label=y)
    # add start of season

    ax.plot(sos, vsos, "or")
    ax.annotate(
        "SOS",
        xy=(sos, vsos),
        xytext=(-15, 20),
        textcoords="offset points",
        arrowprops=dict(arrowstyle="-|>"),
    )

    # add end of season
    ax.plot(eos, veos, "or")
    ax.annotate(
        "EOS",
        xy=(eos, veos),
        xytext=(0, 20),
        textcoords="offset points",
        arrowprops=dict(arrowstyle="-|>"),
    )

    # add peak of season
    ax.plot(pos, vpos, "or")
    ax.annotate(
        "POS",
        xy=(pos, vpos),
        xytext=(-10, -25),
        textcoords="offset points",
        arrowprops=dict(arrowstyle="-|>"),
    )

    # Set the x-axis limits
    min_x = dt.date(y, 1, 1)
    max_x = dt.date(y, 12, 31)
    ax.set_xlim(min_x, max_x)

    adjust_axes(ax)
    counter += 1


host.legend(labelcolor="linecolor")
host.set_ylim([_min - 0.025, _max.values + 0.05])

plt.ylabel(veg_proxy)
plt.xlabel("Month")

plt.title("Yearly " + veg_proxy);

Statistiques de phénologie par pixel

Nous pouvons maintenant calculer les statistiques pour chaque pixel de notre série chronologique et tracer les résultats.

# find all the years to assist with plotting
years = veg_smooth.groupby("time.year")

# get list of years in ts to help with looping
years_int = [y[0] for y in years]

# store results in dict
pheno_results = {}

# loop through years and calculate phenology
for year in years_int:
    # select year
    da = dict(years)[year]

    # calculate stats
    stats = xr_phenology(
        da,
        method_sos=method_sos,
        method_eos=method_eos,
        stats=basic_pheno_stats,
        verbose=False,
    )
    # add results to dict
    pheno_results[str(year)] = stats

Les statistiques phénologiques ont été calculées séparément pour chaque pixel de l’image. Représentons chacun d’eux pour voir les résultats.

Ci-dessous, choisissez une année parmi les résultats phénologiques à tracer.

# Pick a year to plot
year_to_plot = "2019"

En haut du code de traçage, nous re-masquons les résultats de phénologie avec le masque de recadrage. En effet, xr_phenologypossède des méthodes pour gérer les pixels avec uniquement des NaN (comme les régions en dehors du masque de polygone), de sorte que les résultats peuvent avoir des résultats de phénologie pour les régions en dehors du masque. Nous devrons donc masquer à nouveau les données.

# select the year to plot
phen = pheno_results[year_to_plot]

# Define a few items to aid in plotting.
start_date = dt.date(int(year_to_plot), 1, 1)
end_date = dt.date(int(year_to_plot), 10, 27)

date_list = pd.date_range(start_date, end_date, freq="MS")
bounds = [int(i.strftime("%Y%m%d")) for i in date_list]


@mpl.ticker.FuncFormatter
def float_to_date(x, pos):
    tick_str = str(int(x))
    year = tick_str[:4]
    month = tick_str[4:6]
    day = tick_str[6:]
    return f"{year}-{month}-{day}"


# set up figure
fig, ax = plt.subplots(nrows=2,
                       ncols=5,
                       figsize=(18, 8),
                       sharex=True,
                       sharey=True)

# set colorbar size
cbar_size = 0.7

# set aspect ratios
for a in fig.axes:
    a.set_aspect("equal")

# start of season
# Convert SOS values from np.date64 to float values for plotting
cax = phen.SOS.dt.dayofyear.plot(
    ax=ax[0, 0],
    cmap="magma_r",
#     levels=bounds,
    add_colorbar=True,
    cbar_kwargs=dict(shrink=cbar_size, label=None),
)
ax[0, 0].set_title("Start of Season")

phen.vSOS.plot(ax=ax[0, 1],
               cmap="YlGn",
               vmax=0.8,
               cbar_kwargs=dict(shrink=cbar_size, label=None))
ax[0, 1].set_title(veg_proxy + " at SOS")

# peak of season
# # Convert POS values from np.date64 to float values for plotting
cax = phen.POS.dt.dayofyear.plot(
    ax=ax[0, 2],
    cmap="magma_r",
#     levels=bounds,
    add_colorbar=True,
    cbar_kwargs=dict(shrink=cbar_size, label=None),
)
ax[0, 2].set_title("Peak of Season")

phen.vPOS.plot(ax=ax[0, 3],
               cmap="YlGn",
               vmax=0.8,
               cbar_kwargs=dict(shrink=cbar_size, label=None))
ax[0, 3].set_title(veg_proxy + " at POS")

# end of season
# Convert EOS values from np.date64 to float values for plotting
cax = phen.EOS.dt.dayofyear.plot(
    ax=ax[0, 4],
    cmap="magma_r",
#     levels=bounds,
    add_colorbar=True,
    cbar_kwargs=dict(shrink=cbar_size, label=None),
)
ax[0, 4].set_title("End of Season")

phen.vEOS.plot(ax=ax[1, 0],
               cmap="YlGn",
               vmax=0.8,
               cbar_kwargs=dict(shrink=cbar_size, label=None))
ax[1, 0].set_title(veg_proxy + " at EOS")

# Length of Season
phen.LOS.plot(
    ax=ax[1, 1],
    cmap="magma_r",
    vmax=300,
    vmin=0,
    cbar_kwargs=dict(shrink=cbar_size, label=None),
)
ax[1, 1].set_title("Length of Season (Days)")

# Amplitude
phen.AOS.plot(ax=ax[1, 2],
              cmap="YlGn",
              vmax=0.8,
              cbar_kwargs=dict(shrink=cbar_size, label=None))
ax[1, 2].set_title("Amplitude of Season")

# rate of growth
phen.ROG.plot(
    ax=ax[1, 3],
    cmap="coolwarm_r",
    vmin=-0.02,
    vmax=0.02,
    cbar_kwargs=dict(shrink=cbar_size, label=None),
)
ax[1, 3].set_title("Rate of Growth")

# rate of Sensescence
phen.ROS.plot(
    ax=ax[1, 4],
    cmap="coolwarm_r",
    vmin=-0.02,
    vmax=0.02,
    cbar_kwargs=dict(shrink=cbar_size, label=None),
)
ax[1, 4].set_title("Rate of Senescence")
plt.suptitle("Phenology for " + year_to_plot)
plt.tight_layout();

Conclusions

Dans l’exemple ci-dessus, nous pouvons voir que ces quatre champs suivent le même calendrier de culture et qu’il s’agit donc probablement de la même espèce de culture. Nous pouvons également observer des différences intra-champ dans les taux de croissance et dans les valeurs NDVI à différents moments de la saison, qui peuvent être attribuables à des différences de qualité du sol, d’intensité d’arrosage ou d’autres pratiques agricoles.

Les statistiques phénologiques sont un moyen efficace de résumer le cycle saisonnier de la vie d’une plante. Les graphiques par pixel de la phénologie peuvent nous aider à comprendre le calendrier de croissance et de perception de la végétation sur de grandes surfaces et sur diverses espèces végétales, car chaque pixel est traité comme une série indépendante d’observations. Cela pourrait être important, par exemple, si nous voulions évaluer comment les saisons de croissance évoluent à mesure que le climat se réchauffe.

Avancé : Calcul de statistiques temporelles génériques

En plus de la fonction « xr_phenology », le script DE Africa deafrica_tools.temporal contient une autre fonction pour calculer des statistiques de séries chronologiques génériques, temporal_statistics. Cette fonction est basée sur la bibliothèque hdstats (une bibliothèque d’algorithmes statistiques multivariés et à haute dimension). Cette fonction accepte une série chronologique 2 ou 3D de, par exemple, NDVI, et calcule un certain nombre de statistiques récapitulatives, notamment :

• désaccord

• coefficients de la transformée de Fourier discrète (moyenne, écart-type et médiane)

• changement médian

• changement absolu

• complexité

• différence centrale

• nombre de pics (très lent à courir)

Ci-dessous, nous allons calculer un certain nombre de ces statistiques et les représenter graphiquement.

statistics = [
    "f_mean",
    "median_change",
    "abs_change",
    "complexity",
    "central_diff",
#     "discordance",
]

ts_stats = temporal_statistics(veg_smooth, statistics)
print(ts_stats)

   Statistics:
      f_mean
      median_change
      abs_change
      complexity
      central_diff
<xarray.Dataset> Size: 1MB
Dimensions:        (y: 252, x: 194)
Coordinates:
  * x              (x) float64 2kB 8.614e+05 8.614e+05 ... 8.652e+05 8.652e+05
  * y              (y) float64 2kB 1.244e+06 1.244e+06 ... 1.239e+06 1.239e+06
    spatial_ref    int32 4B 6933
Data variables:
    f_mean_n1      (y, x) float32 196kB nan nan nan nan ... 1.198 1.396 1.565
    f_mean_n2      (y, x) float32 196kB nan nan nan nan ... 0.2225 0.2497 0.1253
    f_mean_n3      (y, x) float32 196kB nan nan nan ... 0.09022 0.1026 0.06695
    median_change  (y, x) float32 196kB nan nan nan ... 0.002481 0.005793 0.0
    abs_change     (y, x) float32 196kB nan nan nan ... 0.03155 0.03198 0.0313
    complexity     (y, x) float32 196kB nan nan nan nan ... 4.759 4.142 3.91
    central_diff   (y, x) float32 196kB nan nan nan ... -0.0004274 -0.0006773

if use_crop_mask:
    ts_stats = ts_stats.where(cm == 1)

# set up figure
fig, axes = plt.subplots(2, 4, figsize=(16, 8), sharey=True, sharex=True )

# set aspect ratios
for a in fig.axes:
    a.set_aspect("equal")

# set colorbar size
cbar_size = 0.7
stats = list(ts_stats.data_vars)

# plot
axx = 0
col = 0
for st in stats:
    if axx > 3:
        col += 1
        axx = 0
    ts_stats[st].plot(ax=axes[col][axx], cmap="plasma",cbar_kwargs=dict(shrink=cbar_size, label=None))
    axes[col][axx].set_title(st)
    axx += 1

if len(stats) < 8: fig.delaxes(axes[1][3])

plt.tight_layout();

Prochaines étapes

Une fois terminé, si vous souhaitez exécuter ce code pour une autre région, revenez à la cellule « Paramètres d’analyse », modifiez certaines valeurs (par exemple, « time_range » ou « lat »/« lon ») et réexécutez l’analyse.

Pour les utilisateurs avancés, « xr_phenology » peut être utilisé pour générer des couches de caractéristiques phénologiques dans un classificateur d’apprentissage automatique (voir « Machine Learning with ODC <../Real_world_examples/Machine_learning_with_ODC.ipynb> » pour un exemple d’exécution de modèles ML avec des données ODC). « xr_phenology » peut être passé à l’intérieur du paramètre « custom_func » dans la fonction « deafrica_classificationtools.collect_training_data() », permettant de calculer des statistiques phénologiques lors de la collecte des données d’apprentissage. Un exemple ressemblerait à ceci :

from deafrica_tools.temporal import xr_phenology
from deafrica_tools.classification import collect_training_data

def phenology_stats(da):
    stats = xr_phenology(da, complete='fast_complete')
    return stats

training = collect_training_data(...,
                                 feature_func=phenology_stats)

---

Informations Complémentaires

Licence Le code de ce notebook 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 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 DE Africa <https://digitalearthafrica.slack.com/> ou sur le GIS Stack Exchange <https://gis.stackexchange.com/questions/ask?tags=open-data-cube> en utilisant la balise open-data-cube (vous pouvez consulter les questions posées précédemment `ici <https://gis.stackexchange.com/questions/tagged/open-data-

Si vous souhaitez signaler un problème avec ce notebook, 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-06-04'