Waterbodies Long-Term Animations
Products used: DE Africa Waterbodies; DE Africa WOfS; GeoMADs; CHIRPS
Keywords: data used; waterbodies, data used; wofs_ls_summary_annual, data used; CHIRPS, data used; sentinel-2 geomedian, data used; landsat 5 & 7 geomedian, data used; landsat 8 & 9 geomedian
Background
The Digital Earth Africa Continental Waterbodies Monitoring Service identifies more than 700,000 water bodies from over three decades of satellite observations. This service maps persistent and seasonal water bodies and the change in their water surface area over time. Mapped water bodies may include, but are not limited to, lakes, ponds, man-made reservoirs, wetlands, and segments of some river systems. For more information on the waterbodies monitoring service, refer to the Datasets notebook.
Often, it makes sense to visualise changes in waterbodies through time using animations. Animations are a form of exploratory analysis that can help show whether there is an overall trend in water extent or ilustrate repeating seasonal patterns. This notebook demonstrates how to generate animations for the time series component of waterbodies alongside summary data from Water Observations from Space, rainfall data, and true-colour imagery.
Firstly, an animation showing annual changes in water extent is produced for a lake in Kenya with an increasing trend in water extent. Then, an animation is generated for Lake Chad in Nigeria which is known to have a decreasing trend.
Disclaimer: DE Africa Waterbodies Surface Area Change measures the wet surface area of waterbodies as estimated from satellites. This product does not measure depth, volume, purpose of the waterbody, nor the source of the water.
Description
This notebook demonstrates the generation of two animations for inspecting changes in water extent over several years, drawing on the DE Africa Waterbodies service.
Steps taken are:
Loading and preparing annual time series, WOfS, and true-colour data.
Creating an animation for annual changes in water extent.
Getting started
To run this analysis, run all the cells in the notebook, starting with the “Load packages” cell.
Load packages
Import Python packages that are used for the analysis.
[1]:
import matplotlib.pyplot as plt
import datacube
import matplotlib.dates as mdates
import numpy as np
import matplotlib.animation as animation
from deafrica_tools.plotting import display_map
from IPython.display import HTML
from deafrica_tools.spatial import xr_rasterize
from deafrica_tools.waterbodies import (
get_geohashes,
get_waterbodies,
get_waterbody,
get_time_series,
display_time_series,
)
[2]:
dc = datacube.Datacube(app="Waterbody-anim")
Analysis parameters
This section defines the analysis parameters, including:
lat,lon,buffer: center lat/lon and analysis window size for the area of interest (in degrees)
The default area is Lake Baringo in Kenya which has experienced an increase in extent over recent years and has previously been studied with Digital Earth Africa.
[3]:
# Set the central latitude and longitude
lat = 0.62
lon = 36.08
# Set the buffer to load around the central coordinates
buffer = 0.3
# Compute the bounding box coordinates
xlim = (lon - buffer, lon + buffer)
ylim = (lat + buffer, lat - buffer)
# Preview area on a map
display_map(xlim, ylim)
[3]:
Getting data
The deafrica_waterbodies module allows you to query water bodies by location or geohash. For more background, see the waterbodies documentation and the waterbodies Datasets notebook.
List waterbody polygons in an area
We can get a list of waterbody polygons inside a bounding box of coordinates using get_waterbodies.
[4]:
# Create a bounding box from study area coordinates
bbox = (xlim[0], ylim[1], xlim[1], ylim[0])
# Select all water bodies located within the bounding box
polygons = get_waterbodies(bbox, crs="EPSG:4326")
The returned geodataframe includes all the water bodies which are located within the bounding box. This dataset contains metadata for each water body in the dataset, including the ID, UID, WB_UID, area, perimeter and time series. See the Waterbodies Historical Extent documentation for descriptions of each attribute.
Explore and select the polygon
Once the water body polygons are in memory, you can plot them directly, or explore them in an interactive window.
[5]:
# Explore the waterbody polygons located within the bounding box
polygons.explore()
[5]:
Select geohash
Exploring the polygons above, we can see that our waterbody of interest, Lake Baringo, is represented by a single large polygon. There are numerous smaller polygons in the vicinity or immediately surrounding the lake. The geohash for the primary polygon (shown as uid when hovering the mouse) is sb1ekyxw88.
Getting time series data
[6]:
selected_waterbody_geohash = "sb1ekyxw88"
selected_waterbody = get_waterbody(selected_waterbody_geohash)
Get the wet surface area time series for the selected waterbody
For any given geohash or a polygon, we can also use the get_time_series() function to get various measures of the water body surface over time. See the Waterbodies Historical Extent documentation for descriptions of the different surface measures.
The function also calculates a rolling median of the water body surface wet percentage. This is used to visualise the overall trend in the surface wet percentage. The rolling median uses the last three observations to determine the median at a given date.
[7]:
# Get time series for the selected water body
selected_waterbody_timeseries = get_time_series(waterbody=selected_waterbody)
selected_waterbody_timeseries.head()
[7]:
| area_wet_m2 | percent_wet | area_dry_m2 | percent_dry | area_invalid_m2 | percent_invalid | area_observed_m2 | percent_observed | percent_wet_rolling_median | |
|---|---|---|---|---|---|---|---|---|---|
| date | |||||||||
| 1984-05-30 | 142439400.0 | 62.88 | 79190100.0 | 34.96 | 4910400.0 | 2.17 | 226539900.0 | 100.0 | NaN |
| 1984-06-15 | 147272400.0 | 65.01 | 79223400.0 | 34.97 | 44100.0 | 0.02 | 226539900.0 | 100.0 | NaN |
| 1984-07-01 | 147078900.0 | 64.92 | 79350300.0 | 35.03 | 110700.0 | 0.05 | 226539900.0 | 100.0 | 64.92 |
| 1984-07-17 | 147500100.0 | 65.11 | 78994800.0 | 34.87 | 45000.0 | 0.02 | 226539900.0 | 100.0 | 65.01 |
| 1984-12-08 | 132145200.0 | 58.33 | 83509200.0 | 36.86 | 10885500.0 | 4.81 | 226539900.0 | 100.0 | 64.92 |
Display the wet surface area time series for the selected waterbody
After loading the water body time series, we can use the display_time_series() function to create an interactive visualisation of the time series.
The visualisation shows the invalid percentage and the wet percentage.
Limited data availability prior to 2009 is evident, due to the combination of sensor presence (Landsat-8 was launced in 2013 and Landsat-9 in 2021) and data quality.
The increase in water extent from 2010 to present is clearly observable in the time series.
[8]:
display_time_series(selected_waterbody_timeseries)
Load true colour GeoMADs
The cloud-free geomedian annual composite images are a great source of true-colour data to present alongside annual changes in water extent. The cell below loads true colour bands (red, green, and blue) from three annual GeoMAD products: Landsat-5 & Landsat-7 combined, Landsat-8, and Landsat-8 & Landsat-9 combined. This produces a neat time series of annual images for our period of interest, beginning in 2009 due to data availability.
[9]:
lat_range = (selected_waterbody.total_bounds[1], selected_waterbody.total_bounds[3])
lon_range = (selected_waterbody.total_bounds[0], selected_waterbody.total_bounds[2])
ds = dc.load(
product=["gm_ls5_ls7_annual", "gm_ls8_annual", "gm_ls8_ls9_annual"],
measurements=["red", "green", "blue"],
x=lon_range,
y=lat_range,
resolution=(-30, 30),
dask_chunks={"time": 1, "x": 3000, "y": 3000},
time=("2009-01-01", "2022-12-31"),
).compute()
Load Water Observations from Space (WOfS) annual summary
Our first animation will overlay WOfS onto the true colour imagery. As we are producing an annual animation, the WOfS annual summary is a logical product to use.
[10]:
query = {
"x": lon_range,
"y": lat_range,
"time": ("2009-01-01", "2022-12-31"),
"resolution": (-30, 30),
}
wofs = dc.load(product="wofs_ls_summary_annual", output_crs="epsg:6933", **query)
Mask WOfS to waterbody
The WOfS data needs to be masked to the maximum extent of the waterbody to produce a neat overlay. The waterbody polygon is first rasterised then applied as a mask. The resulting WOfS layer is then plotted to check that masking has worked as expected.
[11]:
fig, ax = plt.subplots()
ax.set_aspect('equal')
mask = xr_rasterize(selected_waterbody, wofs)
clipped_wofs = wofs.where(mask)
clipped_wofs.frequency.isel(time=0).plot()
[11]:
<matplotlib.collections.QuadMesh at 0x7f3249232ce0>
Clip and resample waterbody time series
The timeseries needs to be the same length as the raster datasets to produce an intuitive animation. This means we need to resample the time series to annual time steps and limit it to between 2009 and 2022.
[12]:
tsw = (
selected_waterbody_timeseries.percent_wet.loc["2009-01-01":"2022-12-31"]
.resample("1Y")
.mean()
.interpolate()
)
tsw
[12]:
date
2009-12-31 58.107500
2010-12-31 58.945000
2011-12-31 61.466250
2012-12-31 70.880000
2013-12-31 77.166000
2014-12-31 82.610000
2015-12-31 82.553125
2016-12-31 82.785882
2017-12-31 80.502353
2018-12-31 83.196000
2019-12-31 83.286923
2020-12-31 90.783636
2021-12-31 95.685000
2022-12-31 96.086667
Freq: A-DEC, Name: percent_wet, dtype: float64
Printing the length of each time series to be included in the animation ensures they have equal lengths.
[13]:
print(len(ds.time))
print(len(clipped_wofs.time))
print(len(tsw))
14
14
14
Generate annual animation
The code in the cell below produces an animation at annual timesteps. The first sections of the cell set up the plotting structure, then the data for each axis is set. An update function then determines the rules for updating each time series through annual time steps.
The output can be saved as a .gif file by uncommenting (removing the #) the ani.save() line.
The animation shows a marked increase in water extent between 2010 and 2015, then a steadying until 2020 before a steep increase. Both steep increases in water extent are associated with greening of the background landscape.
[ ]:
# create a figure and axes
fig = plt.figure(figsize=(10, 6))
ax1 = plt.subplot(122)
ax2 = plt.subplot(121)
ax1.set_title("Waterbodies Timeseries")
ax1.set_xlabel("Date")
ax1.set_ylabel("Surface water extent (%)")
years = mdates.YearLocator(2)
ax1.xaxis.set_major_locator(years)
ax1.tick_params(axis="x", labelrotation=45)
bands = ["red", "green", "blue"]
cax = (
ds[bands].isel(time=0).to_array().transpose("y", "x", "variable").squeeze().clip(0, 3000)
/ np.max(
ds[bands].isel(time=0).to_array().transpose("y", "x", "variable").squeeze().clip(0, 3000)
)
).plot.imshow(rgb="variable", animated=True, robust=True, ax=ax2)
ax2.set_aspect('equal')
dax = clipped_wofs.frequency[0, :, :].plot(cmap="Blues", ax=ax2)
selected_waterbody.plot(ax=ax2, edgecolor="black", color="none")
def update(num, x, y, line):
dax.set_array(clipped_wofs.frequency[num, :, :])
cax.set_array(
(ds[bands].isel(time=num).to_array().transpose("y", "x", "variable"))
.squeeze()
.clip(0, 3000)
/ np.max(ds[bands].isel(time=num).to_array().transpose("y", "x", "variable"))
.squeeze()
.clip(0, 3000)
)
ax2.set_title("Time = " + str(clipped_wofs.frequency.coords["time"].values[(num)])[:12])
line.set_data(x[:num], y[:num])
return (line,)
x = tsw.index
y = tsw.values
(line,) = ax1.plot(x, y)
plt.tight_layout()
ani = animation.FuncAnimation(
fig, update, len(ds.time), fargs=[x, y, line], interval=800, blit=True
)
# ani.save('LakeBaringo_Kenya_ts.gif')
plt.close()
HTML(ani.to_html5_video())
Obtain data for Lake Chad
We’ll look at another example to inspect a drying trend. Lake Chad in Nigeria has been described as drying by 90% since the 1960s. Of course, the Water Observations from Space service (derived from Landsat) does not go back that far, but we can inspect trends in recent times.
The cells below follow the same process shown for Baringo, with a few less outputs along the way for brevity.
Firstly, we select the polygon. The nature of the area means there are numerous waterbody polygons that comprise the larger area of Lake Chad. We select the largest polygon which covers the largest and most permanent waterbody, which is s66615ywe3.
[ ]:
# Set the central latitude and longitude
lat = 13.63
lon = 13.77
# Set the buffer to load around the central coordinates
buffer = 0.9
# Compute the bounding box coordinates
xlim = (lon - buffer, lon + buffer)
ylim = (lat + buffer, lat - buffer)
# Create a bounding box from study area coordinates
bbox = (xlim[0], ylim[1], xlim[1], ylim[0])
# Select all water bodies located within the bounding box
polygons = get_waterbodies(bbox, crs="EPSG:4326")
# Explore the waterbody polygons located within the bounding box
polygons.explore()
Inspect time series for Lake Chad
We can see that consistent data is available from late 2013. We will look at the period from 2013 - 2023 for trend visualisation.
[ ]:
selected_waterbody_geohash = "s66615ywe3"
selected_waterbody = get_waterbody(selected_waterbody_geohash)
# Get time series for the selected water body
selected_waterbody_timeseries = get_time_series(waterbody=selected_waterbody)
display_time_series(selected_waterbody_timeseries)
Now that we have the period for visualisation, we will load and process WOfS and the waterbodies time series. The printout from this cell ensure each object has the same length in years.
[ ]:
lat_range = (selected_waterbody.total_bounds[1], selected_waterbody.total_bounds[3])
lon_range = (selected_waterbody.total_bounds[0], selected_waterbody.total_bounds[2])
ds = dc.load(
product=["gm_ls5_ls7_annual", "gm_ls8_annual", "gm_ls8_ls9_annual"],
measurements=["red", "green", "blue"],
x=lon_range,
y=lat_range,
resolution=(-30, 30),
dask_chunks={"time": 1, "x": 3000, "y": 3000},
time=("2013-01-01", "2023-12-31"),
).compute()
query = {
"x": lon_range,
"y": lat_range,
"time": ("2013-01-01", "2023-12-31"),
"resolution": (-30, 30),
}
wofs = dc.load(product="wofs_ls_summary_annual", output_crs="epsg:6933", **query)
plt.figure(figsize=(6,6))
mask = xr_rasterize(selected_waterbody, wofs)
clipped_wofs = wofs.where(mask)
tsw = (
selected_waterbody_timeseries.percent_wet.loc["2009-01-01":"2023-12-31"]
.resample("1Y")
.mean()
.interpolate()
)
print(len(ds.time))
print(len(clipped_wofs.time))
print(len(tsw))
Generate Lake Chad animation
Finally, we can run the animation code for Lake Chad. The output shows a declining trend, though it should be noted that the axis ranges between 81% and 88%.
The output can be saved by removing the # from the ani.save('LakeChad_Nigeria_ts.gif') command.
[ ]:
# create a figure and axes
fig = plt.figure(figsize=(10, 6))
ax1 = plt.subplot(122)
ax2 = plt.subplot(121)
ax1.set_title("Waterbodies Timeseries")
ax1.set_xlabel("Date")
ax1.set_ylabel("Surface water extent (%)")
years = mdates.YearLocator(2)
ax1.xaxis.set_major_locator(years)
ax1.tick_params(axis="x", labelrotation=45)
bands = ["red", "green", "blue"]
cax = (
ds[bands].isel(time=0).to_array().transpose("y", "x", "variable").squeeze().clip(0, 3000)
/ np.max(
ds[bands].isel(time=0).to_array().transpose("y", "x", "variable").squeeze().clip(0, 3000)
)
).plot.imshow(rgb="variable", animated=True, robust=True, ax=ax2)
ax2.set_aspect('equal')
dax = clipped_wofs.frequency[0, :, :].plot(cmap="Blues", ax=ax2)
selected_waterbody.plot(ax=ax2, edgecolor="black", color="none")
def update(num, x, y, line):
dax.set_array(clipped_wofs.frequency[num, :, :])
cax.set_array(
(ds[bands].isel(time=num).to_array().transpose("y", "x", "variable"))
.squeeze()
.clip(0, 3000)
/ np.max(ds[bands].isel(time=num).to_array().transpose("y", "x", "variable"))
.squeeze()
.clip(0, 3000)
)
ax2.set_title("Time = " + str(clipped_wofs.frequency.coords["time"].values[(num)])[:12])
line.set_data(x[:num], y[:num])
return (line,)
x = tsw.index
y = tsw.values
(line,) = ax1.plot(x, y)
plt.tight_layout()
ani = animation.FuncAnimation(
fig, update, len(ds.time), fargs=[x, y, line], interval=800, blit=True
)
# ani.save('LakeChad_Nigeria_ts.gif')
plt.close()
HTML(ani.to_html5_video())
Conclusion
This notebook has demonstrated two animations for visualising waterbodies over several years to show medium to long-term trends. The techniques and code shown here can be adapted to user defined purposes. These visualisations can be useful for exploratory analysis of trends in water extent over several years. More formal trend and pattern analysis would need to be undertaken to better understand these phenomena, especially to understand the effects of factors such as rainfall and land cover change.
Additional information
License: The code in this notebook is licensed under the Apache License, Version 2.0. Digital Earth Africa data is licensed under the Creative Commons by Attribution 4.0 license.
Contact: If you need assistance, please post a question on the Open Data Cube Slack channel or on the GIS Stack Exchange using the open-data-cube tag (you can view previously asked questions here). If you would like to repoart an issue with this notebook, you can file one on
Github.
Compatible datacube version:
[ ]:
import datacube
print(datacube.__version__)
Last Tested:
[ ]:
from datetime import datetime
datetime.today().strftime("%Y-%m-%d")