Phénologie de la végétation par pixel dans la réserve de Ruko
Produits utilisés : s2_l2a
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 <https://www.usgs.gov/land-resources/eros/phenology/science/deriving-phenological-metrics-ndvi?qt-science_center_objects=0#qt-science_center_objects>`__.
Description
Ce carnet produira un tracé bidimensionnel (phénologie par pixel) pour une année donnée.
Plusieurs étapes sont nécessaires pour produire les résultats souhaités :
Charger des données satellite pour une région spécifiée par un fichier vectoriel (shapefile ou geojson)
Mettez en mémoire tampon la couche de masquage des nuages pour mieux masquer les nuages dans les données (le masque de nuages Sentinel-2 est assez médiocre)
Préparez ensuite les données pour l’analyse en supprimant les valeurs erronées (infs), en masquant les eaux de surface et en supprimant les valeurs aberrantes dans l’indice de végétation.
Interpolez et lissez les séries chronologiques pour garantir un ensemble de données cohérent avec toutes les lacunes et le bruit supprimés.
Calculer les statistiques de phénologie par pixel
Tracez les résultats et exportez le tracé sur le disque au format .png
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.
[1]:
%matplotlib inline
import datacube
import numpy as np
import pandas as pd
import xarray as xr
import datetime as dt
import geopandas as gpd
import matplotlib.pyplot as plt
import matplotlib as mpl
from odc.geo.geom import Geometry
from deafrica_tools.datahandling import load_ard
from deafrica_tools.bandindices import calculate_indices
from deafrica_tools.plotting import map_shapefile
import deafrica_tools.temporal as ts
from deafrica_tools.dask import create_local_dask_cluster
from deafrica_tools.spatial import xr_rasterize
from deafrica_tools.classification import HiddenPrints
from datacube.utils.aws import configure_s3_access
configure_s3_access(aws_unsigned=True, cloud_defaults=True)
Configurer un cluster Dask
Dask peut être utilisé pour mieux gérer l’utilisation de la mémoire et effectuer l’analyse en parallèle. Pour une introduction à l’utilisation de Dask avec Digital Earth Africa, consultez le Dask notebook.
Remarque : nous vous recommandons d’ouvrir la fenêtre de traitement Dask pour afficher les différents calculs en cours d’exécution ; pour ce faire, consultez la section Tableau de bord Dask en Afrique de l’Ouest du Dask notebook.
Pour utiliser Dask, configurez le cluster de calcul local à l’aide de la cellule ci-dessous.
[2]:
create_local_dask_cluster(spare_mem='2Gb')
/opt/venv/lib/python3.12/site-packages/distributed/node.py:188: UserWarning: Port 8787 is already in use.
Perhaps you already have a cluster running?
Hosting the HTTP server on port 34803 instead
warnings.warn(
Client
Client-b49dd2a1-43c6-11f1-9e1f-ae180e9931c9
| Connection method: Cluster object | Cluster type: distributed.LocalCluster |
| Dashboard: /user/mpho.sadiki@digitalearthafrica.org/proxy/34803/status |
Cluster Info
LocalCluster
35e785f3
| Dashboard: /user/mpho.sadiki@digitalearthafrica.org/proxy/34803/status | Workers: 1 |
| Total threads: 4 | Total memory: 27.14 GiB |
| Status: running | Using processes: True |
Scheduler Info
Scheduler
Scheduler-483a2aa2-4063-4170-98a4-17201d6bb557
| Comm: tcp://127.0.0.1:40127 | Workers: 0 |
| Dashboard: /user/mpho.sadiki@digitalearthafrica.org/proxy/34803/status | Total threads: 0 |
| Started: Just now | Total memory: 0 B |
Workers
Worker: 0
| Comm: tcp://127.0.0.1:42097 | Total threads: 4 |
| Dashboard: /user/mpho.sadiki@digitalearthafrica.org/proxy/46871/status | Memory: 27.14 GiB |
| Nanny: tcp://127.0.0.1:34787 | |
| Local directory: /tmp/dask-scratch-space/worker-motyduvx | |
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 »
product: Le produit satellite à charger. Soit Sentinel-2 :'s2_l2a', soit Landsat-8 :'ls8_cl2'shapefile: Le chemin d’accès au fichier vectoriel délimitant la région d’analyse. Peut être un shapefile ou un geojson« année » : l’année entière à analyser, par exemple « 2017 »
« min_gooddata » : la fraction de bonnes données (non troubles) qu’une scène doit avoir avant d’être renvoyée sous forme d’ensemble de données
« résolution » : la résolution en pixels, en mètres, de l’ensemble de données renvoyé
« dask_chunks » : la taille, en nombre de pixels, des morceaux dask sur chaque dimension.
[3]:
veg_proxy = 'NDVI'
product = 's2_l2a'
shapefile='data/Ruko_conservancy.geojson'
year = '2018'
resolution = (-20,20)
dask_chunks = {'x':750, 'y':750, 'time':1}
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.
[4]:
dc = datacube.Datacube(app='Vegetation_phenology')
Voir la région d’intérêt
La cellule suivante affichera la zone sélectionnée sur une carte Web.
[5]:
#First open the shapefile using geopandas
gdf = gpd.read_file(shapefile)
[6]:
map_shapefile(gdf, attribute='ConsrvName')
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.
La cellule directement ci-dessous créera un objet de requête en utilisant la première géométrie du fichier de formes, ainsi que les paramètres que nous avons définis dans la section Paramètres d’analyse ci-dessus.
[7]:
# Create a reusable query
geom = Geometry(geom=gdf.iloc[0].geometry, crs=gdf.crs)
query = {
"geopolygon": geom,
'time': year,
'measurements': ['red','nir','green','swir_1'],
'resolution': resolution,
'output_crs': 'epsg:6933',
'group_by':'solar_day'
}
Charger les données disponibles depuis S2 :
[8]:
filters=[("opening", 3),("dilation", 2)]
ds = load_ard(
dc=dc,
products=['s2_l2a'],
dask_chunks=dask_chunks,
mask_filters=filters,
**query,
)
print(ds)
Using pixel quality parameters for Sentinel 2
Finding datasets
s2_l2a
Applying morphological filters to pq mask [('opening', 3), ('dilation', 2)]
/opt/venv/lib/python3.12/site-packages/odc/algo/_masking.py:425: FutureWarning: `binary_opening` is deprecated since version 0.26 and will be removed in version 0.28. Use `skimage.morphology.opening` instead.
mask = op(mask, _disk(radius, mask.ndim))
/opt/venv/lib/python3.12/site-packages/odc/algo/_masking.py:425: FutureWarning: `binary_dilation` is deprecated since version 0.26 and will be removed in version 0.28. Use `skimage.morphology.dilation` instead. Note the lack of mirroring for non-symmetric footprints (see docstring notes).
mask = op(mask, _disk(radius, mask.ndim))
Applying pixel quality/cloud mask
/opt/venv/lib/python3.12/site-packages/odc/algo/_masking.py:425: FutureWarning: `binary_opening` is deprecated since version 0.26 and will be removed in version 0.28. Use `skimage.morphology.opening` instead.
mask = op(mask, _disk(radius, mask.ndim))
/opt/venv/lib/python3.12/site-packages/odc/algo/_masking.py:425: FutureWarning: `binary_dilation` is deprecated since version 0.26 and will be removed in version 0.28. Use `skimage.morphology.dilation` instead. Note the lack of mirroring for non-symmetric footprints (see docstring notes).
mask = op(mask, _disk(radius, mask.ndim))
Returning 72 time steps as a dask array
/opt/venv/lib/python3.12/site-packages/deafrica_tools/datahandling.py:565: FutureWarning: In a future version of xarray the default value for compat will change from compat='no_conflicts' to compat='override'. This is likely to lead to different results when combining overlapping variables with the same name. To opt in to new defaults and get rid of these warnings now use `set_options(use_new_combine_kwarg_defaults=True) or set compat explicitly.
ds = xr.merge([ds_data, ds_masks])
<xarray.Dataset> Size: 835MB
Dimensions: (time: 72, y: 1194, x: 607)
Coordinates:
* time (time) datetime64[ns] 576B 2018-01-02T08:07:58 ... 2018-12-2...
* y (y) float64 10kB 9.461e+04 9.459e+04 ... 7.077e+04 7.075e+04
* x (x) float64 5kB 3.479e+06 3.479e+06 ... 3.491e+06 3.491e+06
spatial_ref int32 4B 6933
Data variables:
red (time, y, x) float32 209MB dask.array<chunksize=(1, 750, 607), meta=np.ndarray>
nir (time, y, x) float32 209MB dask.array<chunksize=(1, 750, 607), meta=np.ndarray>
green (time, y, x) float32 209MB dask.array<chunksize=(1, 750, 607), meta=np.ndarray>
swir_1 (time, y, x) float32 209MB dask.array<chunksize=(1, 750, 607), meta=np.ndarray>
Attributes:
crs: EPSG:6933
grid_mapping: spatial_ref
Masquer les données satellite avec la forme
[9]:
#create mask
mask = xr_rasterize(gdf,ds)
#mask data
ds = ds.where(mask)
#convert to float 32 to conserve memory
ds=ds.astype(np.float32)
Calculer les indices de végétation et d’eau
[10]:
# Calculate the chosen vegetation proxy index and add it to the loaded data set
ds = calculate_indices(ds, index=[veg_proxy, 'MNDWI'], satellite_mission='s2', drop=True)
Dropping bands ['red', 'nir', 'green', 'swir_1']
Préparer les données pour l’analyse
Supprimez toutes les valeurs NaN ou infinies, masquez l’eau, supprimez toutes les valeurs aberrantes dans l’indice de végétation. Nous réduisons ensuite les données à une série temporelle 1D en calculant la moyenne sur les dimensions x et y.
Cette cellule prendra quelques minutes à calculer puisque le « water_mask » doit être mis en mémoire (nous l’utiliserons plus tard pour masquer les résultats de la phénologie).
[11]:
# remove any infinite values
ds = ds.where(~np.isinf(ds))
# mask water
ds = ds.where(ds.MNDWI < 0)
#create a all-time water mask for use later
water_mask = ds.MNDWI < 0
water_mask = water_mask.max('time')
water_mask = water_mask.compute()
#remove outliers (if EVI greater than 1.0, set to NaN)
ds[veg_proxy] = xr.where(ds[veg_proxy]>1.0, np.nan, ds[veg_proxy])
ds[veg_proxy] = xr.where(ds[veg_proxy]<0, np.nan, ds[veg_proxy])
# create 1D line plots
veg = ds[veg_proxy]
Lisser et interpoler des séries temporelles
En raison de nombreux facteurs (par exemple, des nuages obscurcissant la région, une couverture nuageuse manquante dans la couche SCL), les données seront incomplètes et bruyantes. Ici, nous allons lisser et interpoler les données pour garantir que nous travaillons avec une série chronologique cohérente.
Pour ce faire, nous procédons en deux étapes :
Rééchantillonner les données à des intervalles de temps bimensuels en utilisant la médiane
Calculer une moyenne mobile avec une fenêtre de 4 étapes
Ces calculs prendront plusieurs minutes à terminer car nous exécuterons .compute(), déclenchant toutes les tâches que nous avons planifiées ci-dessus et mettant les tableaux en mémoire.
[12]:
resample_period='2W'
window=4
# calculate median first, then bring into memory with .compute()
veg_smooth = veg.resample(time=resample_period, label='left').median('time')
# Update the time coordinates of the resampled dataset.
veg_smooth = veg_smooth.assign_coords(time=(veg_smooth.time + np.timedelta64(1, 'W')))
veg_smooth = veg_smooth.compute()
#calculate rolling mean to smooth time-series
veg_smooth=veg_smooth.rolling(time=window, min_periods=1).mean()
Calculer les statistiques de phénologie 2D
Ci-dessous, nous spécifions les statistiques à calculer et la méthode que nous utiliserons pour déterminer les statistiques. Les options sont « first » et « median » pour « method_sos », et « last » et « median » pour « method_eos ».
method_sos : str
If 'first' then vSOS is estimated as the first positive
slope on the greening side of the curve. If 'median',
then vSOS is estimated as the median value of the postive
slopes on the greening side of the curve.
method_eos : str
If 'last' then vEOS is estimated as the last negative slope
on the senescing side of the curve. If 'median', then vEOS is
estimated as the 'median' value of the negative slopes on the
senescing side of the curve.
[13]:
pheno_stats = ['SOS','vSOS','POS','vPOS','EOS','vEOS','Trough','LOS','AOS','ROG','ROS']
method_sos = 'first'
method_eos = 'last'
[14]:
with HiddenPrints():
phen=ts.xr_phenology(
veg_smooth,
method_sos=method_sos,
method_eos=method_eos,
stats=pheno_stats,
)
print(phen)
<xarray.Dataset> Size: 41MB
Dimensions: (y: 1194, x: 607)
Coordinates:
* y (y) float64 10kB 9.461e+04 9.459e+04 ... 7.077e+04 7.075e+04
* x (x) float64 5kB 3.479e+06 3.479e+06 ... 3.491e+06 3.491e+06
spatial_ref int32 4B 0
Data variables:
SOS (y, x) datetime64[ns] 6MB NaT NaT NaT NaT ... NaT NaT NaT NaT
vSOS (y, x) float32 3MB nan nan nan nan nan ... nan nan nan nan nan
POS (y, x) datetime64[ns] 6MB NaT NaT NaT NaT ... NaT NaT NaT NaT
vPOS (y, x) float32 3MB nan nan nan nan nan ... nan nan nan nan nan
EOS (y, x) datetime64[ns] 6MB NaT NaT NaT NaT ... NaT NaT NaT NaT
vEOS (y, x) float32 3MB nan nan nan nan nan ... nan nan nan nan nan
Trough (y, x) float32 3MB nan nan nan nan nan ... nan nan nan nan nan
LOS (y, x) float32 3MB nan nan nan nan nan ... nan nan nan nan nan
AOS (y, x) float32 3MB nan nan nan nan nan ... nan nan nan nan nan
ROG (y, x) float32 3MB nan nan nan nan nan ... nan nan nan nan nan
ROS (y, x) float32 3MB nan nan nan nan nan ... nan nan nan nan nan
Attributes:
grid_mapping: spatial_ref
Remasquer les données
Le code de phénologie dispose de méthodes permettant de gérer les pixels contenant uniquement des NaN (comme les régions situées à l’extérieur du masque de polygone), de sorte que les résultats peuvent avoir des résultats de phénologie pour les régions situées au-dessus de l’eau et à l’extérieur du masque. Nous devrons donc masquer à nouveau les données.
[15]:
#mask with polygon
phen = phen.where(mask)
#mask with water-mask
phen = phen.where(water_mask)
Tracer les statistiques de phénologie 2D
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.
[16]:
year_to_plot = year
# 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();
plt.savefig('results/phenology_2D_'+year+".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 :
[17]:
print(datacube.__version__)
1.9.13
Dernier test :
[18]:
from datetime import datetime
datetime.today().strftime('%Y-%m-%d')
[18]:
'2026-04-29'