Generating NDVI animations¶
Products used: ndvi_anomaly
Keywords: data used; ndvi anomaly, agriculture
Background¶
This notebook shows how to generate an animation to visualise the NDVI across space, mean NDVI through time, and NDVI anomalies. It also generates a similar animation for monthly rainfall and rainfall anomalies.
Animations can help visualise changes in patterns through time. The visualisation of both mean values and anomalies means that seasonal and inter-annual aspects of time can be assessed. Complementing this with spatial images enhances interpretability of earth observation data.
Description¶
The workflow for generating an animation progresses through the following steps: 1. Select an area of interest. 2. Load DE Africa Mean NDVI and Anomalies product. 3. Mask to cropland. 4. Generate NDVI animation. 5. Generate rainfall animation. ***
Getting started¶
To run this analysis, run all the cells in the notebook, starting with the “Load packages and apps” cell.
Load packages¶
[1]:
%matplotlib inline
import datacube
import datetime as dt
import xarray as xr
import geopandas as gpd
from datacube.utils.geometry import Geometry
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from IPython.display import HTML
import matplotlib.dates as mdates
from deafrica_tools.plotting import display_map
from deafrica_tools.areaofinterest import define_area
Connect to the datacube.
[2]:
dc = datacube.Datacube(app='NDVI-animation')
Analysis parameters¶
The default area is cropland in Madagascar.
The following cell sets the parameters, which define the area of interest and the length of time to conduct the analysis over. The parameters are
lat
: The central latitude to analyse (e.g.10.338
).lon
: The central longitude to analyse (e.g.-1.055
).buffer
: The number of square degrees to load around the central latitude and longitude. For reasonable loading times, set this as0.1
or lower.
Select location¶
To define the area of interest, there are two methods available:
By specifying the latitude, longitude, and buffer. This method requires you to input the central latitude, central longitude, and the buffer value in square degrees around the center point you want to analyze. For example,
lat = 10.338
,lon = -1.055
, andbuffer = 0.1
will select an area with a radius of 0.1 square degrees around the point with coordinates (10.338, -1.055).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 buttonin 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.
To use one of these methods, you can uncomment the relevant line of code and comment out the other one. To comment out a line, add the "#"
symbol before the code you want to comment out. By default, the first option which defines the location using latitude, longitude, and buffer is being used.
If running the notebook for the first time, keep the default settings below. This will demonstrate how the analysis works and provide meaningful results.
[3]:
# # Define the area of interest for the analysis
year = "2022"
# Method 1: Specify the latitude, longitude, and buffer
aoi = define_area(lat=-17.71, lon=48.21, buffer=0.05)
# 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])
display_map(lon_range, lat_range)
[3]:
We will define the start and end date as the beginning and end of our year of interest.
[4]:
start_date = dt.datetime.strptime(year, "%Y").replace(month=1, day=1)
end_date = dt.datetime.strptime(year, "%Y").replace(month=12, day=31)
Below, the analysis parameters are used to construct a query which loads the DE Africa mean NDVI and anomalies product.
[5]:
# Construct the data cube query
query = {
'x': lon_range,
'y': lat_range,
'time': (start_date, end_date),
'output_crs': 'EPSG:6933',
'resolution': (-30, 30)
}
ds = dc.load(product="ndvi_anomaly", **query)
ds
[5]:
<xarray.Dataset> Dimensions: (time: 12, y: 406, x: 323) Coordinates: * time (time) datetime64[ns] 2022-01-16T11:59:59.999999 ... 20... * y (y) float64 -2.218e+06 -2.218e+06 ... -2.23e+06 -2.231e+06 * x (x) float64 4.647e+06 4.647e+06 ... 4.656e+06 4.656e+06 spatial_ref int32 6933 Data variables: ndvi_mean (time, y, x) float32 0.4654 0.4267 0.3652 ... nan nan nan ndvi_std_anomaly (time, y, x) float32 -0.1621 0.2036 0.3302 ... nan nan nan clear_count (time, y, x) int8 4 4 4 4 4 4 4 4 4 ... 0 0 0 0 0 0 0 0 0 Attributes: crs: epsg:6933 grid_mapping: spatial_ref
Load cropland mask¶
This is an optional step for masking the NDVI data to cropland. In this instance, we are interested in inspecting crop vegetation patterns only, so we use the mask to do this.
[6]:
query = {
'x': lon_range,
'y': lat_range,
'time': '2019',
'output_crs': 'EPSG:6933',
'resolution': (-30, 30),
'measurements': 'mask'
}
cm = dc.load(product='crop_mask', **query).squeeze()
cm.mask.where(cm.mask<255).plot()
[6]:
<matplotlib.collections.QuadMesh at 0x7f9e3304e710>

Mask NDVI to cropland¶
We mask the NDVI data to the cropland extent shown above.
[7]:
ds = ds.where(cm.mask ==1)
Prepare animation data¶
We need to prepare three objects for the animation: 1. Mean NDVI with spatial attributes. 2. Mean NDVI through time, aggregated as an average across space. 3. NDVI anomalies for each month, aggregated as an average across space.
[8]:
ds_mean = ds.ndvi_mean
ds_mean_1D = ds_mean.mean(dim=['x', 'y'], skipna=True)
ndvi_anom = ds.ndvi_std_anomaly.mean(dim=['x', 'y']).to_dataframe().drop(['spatial_ref'], axis=1)
Run animation¶
To create the animation, the cell below follows a few key steps: 1. Define the layout of the three axes and their size. 2. Set axis parameters (titles etc) 3. Start each axis with the first image of the animation (Januray). 4. Define an update function which adds each month sequentially to each of the three charts. 5. Run the animation.
Note the animation renders in its own window. The bottom right toggle enables download of the animation.
[9]:
# create a figure and axes
fig = plt.figure()
fig.set_figheight(5)
fig.set_figwidth(10)
ax1 = plt.subplot2grid(shape=(2, 2), loc=(0, 1), colspan=2)
ax2 = plt.subplot2grid(shape=(2, 2), loc=(0, 0), rowspan=2)
ax3 = plt.subplot2grid(shape=(2, 2), loc=(1, 1), colspan=2)
# setup starting plot
cax = ds.ndvi_mean[0,:,:].plot(
add_colorbar=True,
cmap='RdYlGn',
vmin=-1, vmax=1,
cbar_kwargs={
'extend':'neither'
},
ax=ax2
)
ax1.set_ylim([0, 1])
ax1.xaxis.set_major_formatter(mdates.DateFormatter("%b"))
ax1.set_title("Mean NDVI")
ax1.set_xlabel("Date")
ax1.set_ylabel("NDVI")
ax3.set_ylim([-4, 4])
ax3.xaxis_date()
ax3.xaxis.set_major_formatter(mdates.DateFormatter("%b"))
ax3.set_title("NDVI Anomaly")
ax3.set_xlabel("Date")
ax3.set_ylabel("Std Deviations from monthly mean")
ax3.axhline(0, color='black', linestyle='--')
x = ds_mean_1D.time
y = ds_mean_1D
x2 = ndvi_anom.index
y2 = [0] * ndvi_anom.ndvi_std_anomaly
bars = ax3.bar(x2, y2, align='center', width=20,color=(ndvi_anom.ndvi_std_anomaly > 0).map({True: 'g', False: 'brown'}))
line, = ax1.plot(x, y, marker='*',
color='blue')
def update(num, x, y, x2, y2, bars, line):
y2[num] = ndvi_anom.ndvi_std_anomaly[num]
bars[num].set_height(y2[num])
cax.set_array(ds.ndvi_mean[(num),:,:].values.flatten())
ax2.set_title("Time = " + str(ds.ndvi_mean.coords['time'].values[(num)])[:12])
line.set_data(x[:num+1], y[:num+1])
return line,
plt.tight_layout()
ani = animation.FuncAnimation(fig, update, len(x), fargs=[x, y, x2, y2, bars, line],
interval=400, blit=False, repeat=True)
#ani.save('ndvi_anim.gif')
plt.close()
HTML(ani.to_html5_video())
[9]:
Make an NDVI animation with rainfall¶
The animation above summarises cropland NDVI in space and time, within a period of interest. We can add associated information to this visualisation, such as rainfall. The next few steps demonstrate how to load rainfall data and calculate anomalies for plotting.
Rainfall anomaly parameters¶
We need to set some parameters for calculating the rainfall anomaly. The time_m
period is a baseline period while time_x
is the year of interest.
[10]:
# Set the range of dates for the climatology, this will be the reference period (m) for the anomaly calculation.
# Standard practice is to use a 30 year period, so we've used 1981 to 2011 in this example.
time_m = ('2000', '2020')
# time period for monthly anomaly (x)
time_x = ('2022')
# CHIRPS has a spatial resolution of ~5x5 km
resolution = (-5000, 5000)
#size of dask chunks
dask_chunks = dict(x=500,y=500)
query_m = {
'x': lon_range,
'y': lat_range,
'time': time_m,
'output_crs': 'EPSG:6933',
'resolution': (-5000, 5000)
}
query_x = {
'x': lon_range,
'y': lat_range,
'time': time_x,
'output_crs': 'EPSG:6933',
'resolution': (-5000, 5000)
}
Load rainfall data¶
[11]:
ds_rf_m = dc.load(product='rainfall_chirps_monthly', **query_m)
ds_rf_x = dc.load(product='rainfall_chirps_monthly', **query_x)
Below, we calculate the monthly mean rainfall and anomaly.
[12]:
# monthly means
climatology_mean = ds_rf_m.groupby('time.month').mean('time').compute()
#calculate monthly std dev
climatology_std = ds_rf_m.groupby('time.month').std('time').compute()
stand_anomalies = xr.apply_ufunc(
lambda x, m, s: (x - m) / s,
ds_rf_x.groupby("time.month"),
climatology_mean,
climatology_std,
output_dtypes=[ds_rf_x.rainfall.dtype],
dask="allowed"
).compute()
Prepare animation data¶
As for the NDVI animation above, we need to prepare dataframes for the rainfall animation.
[13]:
spatial_mean_rain = ds_rf_x.rainfall.mean(['x','y'])
spatial_mean_rain.time
spatial_mean_anoms = stand_anomalies.rainfall.mean(['x','y'], skipna=True).sel(time=slice('2021-12-31', end_date)).to_dataframe().drop(
['spatial_ref', 'month'], axis=1).fillna(0)
Run the animation¶
We follow the steps listed for the NDVI animation above to produce the rainfall animation.
[14]:
# create a figure and axes
fig = plt.figure()
fig.set_figheight(5)
fig.set_figwidth(10)
ax1 = plt.subplot2grid(shape=(2, 2), loc=(0, 1), colspan=2)
ax2 = plt.subplot2grid(shape=(2, 2), loc=(0, 0), rowspan=2)
ax3 = plt.subplot2grid(shape=(2, 2), loc=(1, 1), colspan=2)
cax = ds.ndvi_mean[0,:,:].plot(
add_colorbar=True,
cmap='RdYlGn',
vmin=-1, vmax=1,
cbar_kwargs={
'extend':'neither'
},
ax=ax2
)
ax1.set_title("Monthly Total Rainfall")
ax1.set_xlabel("Date")
ax1.xaxis.set_major_formatter(mdates.DateFormatter("%b"))
ax1.set_ylabel("Monthly Rainfall (mm)")
ax3.clear()
ax3.set_ylim([-4, 4])
ax3.xaxis_date()
ax3.xaxis.set_major_formatter(mdates.DateFormatter("%b"))
ax3.set_title("Rainfall Anomaly")
ax3.set_xlabel("Date")
ax3.set_ylabel("Std Deviations from monthly mean")
ax3.axhline(0, color='black', linestyle='--')
def update(num, x, y, line):
y2[num] = spatial_mean_anoms.rainfall[num]
bars[num].set_height(y2[num])
cax.set_array(ds.ndvi_mean[(num),:,:].values.flatten())
ax2.set_title("Time = " + str(ds.ndvi_mean.coords['time'].values[(num)])[:12])
line.set_data(x[:num+1], y[:num+1])
return line,
x = spatial_mean_rain.time
y = spatial_mean_rain
x2 = spatial_mean_anoms.index
y2 = [0] * spatial_mean_anoms.rainfall
bars = ax3.bar(x2, y2, align='center', width=20,color=(spatial_mean_anoms.rainfall > 0).map({True: 'b', False: 'r'}))
line, = ax1.plot(x, y, marker='*',
color='blue')
plt.tight_layout()
ani = animation.FuncAnimation(fig, update, len(x), fargs=[x, y, line],
interval=400, blit=True)
plt.close()
HTML(ani.to_html5_video())
[14]:
Drawing conclusions¶
The animations enable us to inspect how vegetation changes within a year, and by showing anomalies, how the year of interest compares to a historical baseline.
This kind of analysis can assist in assessment of crop production in a given year, and can feed into early warning systems and forecasts on crop production and food security.
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 report an issue with this notebook, you can file one on
Github.
Compatible datacube version:
[15]:
print(datacube.__version__)
1.8.15
Last Tested:
[16]:
from datetime import datetime
datetime.today().strftime('%Y-%m-%d')
[16]:
'2023-08-14'