Cartographie des changements à long terme dans l’étendue annuelle des eaux du delta de l’Okavango
Produits utilisés : wofs_ls_summary_annual,
Mots-clés : données utilisées ; WOfS, eau ; étendue, analyse ; séries chronologiques, visualisation ; animation
Description
Le bloc-notes montre comment charger, visualiser et analyser le produit « résumé annuel WOfS <https://docs.digitalearthafrica.org/en/latest/data_specs/Landsat_WOfS_specs.html> » pour recueillir des informations sur l’étendue à long terme du delta de l’Okavango.
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
Importez les packages Python utilisés pour l’analyse.
[1]:
%matplotlib inline
# Force GeoPandas to use Shapely instead of PyGEOS
# In a future release, GeoPandas will switch to using Shapely by default.
import os
os.environ['USE_PYGEOS'] = '0'
import datacube
import numpy as np
import xarray as xr
import seaborn as sns
import geopandas as gpd
import matplotlib.pyplot as plt
from IPython.display import Image
from matplotlib.colors import ListedColormap
from matplotlib.patches import Patch
from deafrica_tools.bandindices import calculate_indices
from deafrica_tools.plotting import map_shapefile, xr_animation
from deafrica_tools.dask import create_local_dask_cluster
from deafrica_tools.spatial import xr_rasterize
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 activer Dask, configurez le cluster de calcul local à l’aide de la cellule ci-dessous.
[ ]:
create_local_dask_cluster()
Se connecter au datacube
Activez la base de données Datacube, qui fournit des fonctionnalités pour le chargement et l’affichage des données d’observation de la Terre stockées.
[3]:
dc = datacube.Datacube(app='water_extent')
Paramètres d’analyse
La cellule suivante définit les paramètres qui définissent la zone d’intérêt et la durée de réalisation de l’analyse.
Les paramètres sont :
vector_file: Le chemin vers le shapefile ou geojson qui définira la zone d’analyse de l’étudestart_yearetend_year: la plage de dates à analyser (par exemple('1990', '2020').
[4]:
# Define the area of interest
vector_file = 'data/okavango_delta_outline.geojson'
# Define the start year and end year
start_year = '1990'
end_year = '2020'
Visualisez les données vectorielles d’intérêt sur une carte interactive
[5]:
#read shapefile
gdf = gpd.read_file(vector_file)
# map_shapefile(gdf, attribute='GRID_CODE')
Télécharger les résumés annuels du WOfS
[6]:
#Create a query object
bbox=list(gdf.total_bounds)
lon_range = (bbox[0], bbox[2])
lat_range = (bbox[1], bbox[3])
query = {
'x': lon_range,
'y': lat_range,
'resolution': (-30, 30),
'output_crs':'EPSG:6933',
'time': (start_year, end_year),
'dask_chunks':dict(x=1000,y=1000)
}
#load wofs
ds = dc.load(product="wofs_ls_summary_annual",
**query)
print(ds)
<xarray.Dataset>
Dimensions: (time: 31, y: 7774, x: 6955)
Coordinates:
* time (time) datetime64[ns] 1990-07-02T11:59:59.999999 ... 2020-07...
* y (y) float64 -2.279e+06 -2.279e+06 ... -2.512e+06 -2.513e+06
* x (x) float64 2.097e+06 2.097e+06 ... 2.305e+06 2.305e+06
spatial_ref int32 6933
Data variables:
count_wet (time, y, x) int16 dask.array<chunksize=(1, 1000, 1000), meta=np.ndarray>
count_clear (time, y, x) int16 dask.array<chunksize=(1, 1000, 1000), meta=np.ndarray>
frequency (time, y, x) float32 dask.array<chunksize=(1, 1000, 1000), meta=np.ndarray>
Attributes:
crs: epsg:6933
grid_mapping: spatial_ref
[7]:
#create mask
mask = xr_rasterize(gdf,ds)
#mask data
ds = ds.where(mask)
Animation de séries chronologiques
Dans la cellule suivante, nous traçons l’ensemble de données que nous avons chargé ci-dessus sous forme de GIF d’animation, en utilisant la fonction `xr_animation <../Frequently_used_code/Animated_timeseries.ipynb>`__. Le output_path sera enregistré dans le répertoire où se trouve le script et vous pouvez modifier les noms pour éviter l’écrasement des fichiers.
[8]:
out_path = 'results/okavango_annual_water_frequency.gif'
xr_animation(ds=ds,
output_path=out_path,
interval=400,
bands=['frequency'],
show_text='WOfS Annual Summary',
show_date = '%Y',
width_pixels=700,
annotation_kwargs={'fontsize': 15},
imshow_kwargs={'cmap': sns.color_palette("mako_r", as_cmap=True), 'vmin': 0.0, 'vmax': 0.9},
colorbar_kwargs={'colors': 'black'},
show_colorbar=False)
# Plot animated gif
plt.close()
Image(filename=out_path)
Exporting animation to results/okavango_annual_water_frequency.gif
[8]:
<IPython.core.display.Image object>
Calculer la superficie annuelle de l’étendue d’eau
Le nombre de pixels peut être utilisé pour la surface du plan d’eau si la surface en pixels est connue. Exécutez la cellule suivante pour générer les constantes nécessaires à la réalisation de cette conversion.
[9]:
pixel_length = query["resolution"][1] # in metres
m_per_km = 1000 # conversion from metres to kilometres
area_per_pixel = pixel_length**2 / m_per_km**2
Seuil de fréquence annuelle WOfS pour classer eau/non-eau
Calcule la surface des pixels classés comme eau (c’est-à-dire si « ds.frequency » est > 0,1, alors le pixel sera considéré comme de l’eau libre pendant l’année)
[10]:
water_threshold = 0.1
[11]:
#threshold
water_extent = (ds.frequency > water_threshold).persist()
#calculate area
ds_valid_water_area = water_extent.sum(dim=['x', 'y']) * area_per_pixel
Tracer la superficie annuelle des eaux libres
[12]:
plt.figure(figsize=(18, 4))
ds_valid_water_area.plot(marker='o', color='#9467bd')
plt.title(f'Observed Annual Area of Water from {start_year} to {end_year}')
plt.xlabel('Dates')
plt.ylabel('Waterbody area (km$^2$)')
plt.tight_layout()
Déterminer l’étendue d’eau minimale et maximale
La cellule suivante extrait l’étendue minimale et maximale de l’eau de l’ensemble de données à l’aide des fonctions « min » et « max », nous ajoutons ensuite les dates à un « xarray.DataArray ».
[13]:
min_water_area_date, max_water_area_date = min(ds_valid_water_area), max(ds_valid_water_area)
time_xr = xr.DataArray([min_water_area_date.time.values, max_water_area_date.time.values], dims=["time"])
print(time_xr)
<xarray.DataArray (time: 2)>
array(['1996-07-01T23:59:59.999999000', '2011-07-02T11:59:59.999999000'],
dtype='datetime64[ns]')
Dimensions without coordinates: time
Tracez les dates auxquelles se produisent les étendues d’eau minimales et maximales
Tracez le pixel classé de l’eau pour les deux dates où nous avons l’étendue minimale et maximale de l’eau de surface.
[14]:
water_extent.sel(time=time_xr).plot.imshow(col="time", col_wrap=2, figsize=(14, 6));
Comparer deux périodes
Les cellules suivantes déterminent l’étendue maximale de l’eau pour deux années différentes.
baseline_year: l’année de référence pour l’analyseanalysis_year: l’année à comparer à l’année de référence
[15]:
baseline_time = '2019'
analysis_time = '2020'
baseline_ds, analysis_ds = ds_valid_water_area.sel(time=baseline_time), ds_valid_water_area.sel(time=analysis_time)
Traçage
Tracez l’étendue de l’eau pour les deux périodes choisies.
[16]:
compare = water_extent.sel(time=[baseline_ds.time.values[0], analysis_ds.time.values[0]])
compare.plot.imshow(col="time",col_wrap=2,figsize=(15, 8), cmap='viridis', add_colorbar=False);
Calcul de la variation pour les deux périodes désignées
Les cellules ci-dessous calculent la quantité d’eau gagnée, perdue et stable pour les deux périodes
[17]:
analyse_total_value = compare.isel(time=1).astype(int)
change = analyse_total_value - compare.isel(time=0).astype(int)
water_appeared = change.where(change == 1)
permanent_water = change.where((change == 0) & (analyse_total_value == 1))
permanent_land = change.where((change == 0) & (analyse_total_value == 0))
water_disappeared = change.where(change == -1)
La cellule ci-dessous calcule la superficie de l’étendue de l’eau pour la perte d’eau, le gain d’eau, l’eau permanente et le sol
[18]:
total_area = analyse_total_value.count().values * area_per_pixel
water_apperaed_area = water_appeared.count().values * area_per_pixel
permanent_water_area = permanent_water.count().values * area_per_pixel
water_disappeared_area = water_disappeared.count().values * area_per_pixel
Traçage
Les variables de l’eau sont tracées pour visualiser le résultat
[19]:
water_appeared_color = "Green"
water_disappeared_color = "Yellow"
stable_color = "Blue"
land_color = "Brown"
fig, ax = plt.subplots(1, 1, figsize=(15, 15))
compare[1].plot.imshow(cmap="Pastel1",
add_colorbar=False,
add_labels=False,
ax=ax)
water_appeared.plot.imshow(
cmap=ListedColormap([water_appeared_color]),
add_colorbar=False,
add_labels=False,
ax=ax,
)
water_disappeared.plot.imshow(
cmap=ListedColormap([water_disappeared_color]),
add_colorbar=False,
add_labels=False,
ax=ax,
)
permanent_water.plot.imshow(cmap=ListedColormap([stable_color]),
add_colorbar=False,
add_labels=False,
ax=ax)
plt.legend(
[
Patch(facecolor=stable_color),
Patch(facecolor=water_disappeared_color),
Patch(facecolor=water_appeared_color),
Patch(facecolor=land_color),
],
[
f"Water to Water {round(permanent_water_area, 2)} km2",
f"Water to No Water {round(water_disappeared_area, 2)} km2",
f"No Water to Water: {round(water_apperaed_area, 2)} km2",
],
loc="lower left",
)
plt.title("Change in water extent: " + baseline_time + " to " + analysis_time);
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 :
[20]:
print(datacube.__version__)
1.8.15
Dernier test :
[21]:
from datetime import datetime
datetime.today().strftime('%Y-%m-%d')
[21]:
'2023-08-21'