Products used: ndvi_climatology_ls
Keywords: data used; ndvi climatology
Climatology refers to conditions averaged over a long period of time, typically greater than 30 years. The Digital Earth Africa NDVI Climatology product represents the long-term average baseline condition of vegetation for every Landsat pixel over the African continent. Both mean and standard deviation NDVI climatologies are available for each calendar month based on calculation over the period 1984-2020. NDVI climatologies may be used for many applications including identifying extremeties (anomalies) in vegetation condition, identifying both long and short-term changes in vegetation condition, and as an input into machine learning processes for land use classification.
Further details on the calculation of the product are available in the NDVI Climatology technical specifications documentation.
Datacube product name:
mean_<month>: These measurements show the mean NDVI calculated from all available NDVI data from 1984-2020 for the given month.
stddev_<month>: These measurements show the standard deviation of NDVI values of all available NDVI data from 1984-2020 for the given month.
count_<month>: These measaurements show the number of clear observations that go into creating the mean and standard deviation measurements. Importantly, caution should be used when applying this product to regions where the clear observation count is less than approximately 20-30. This can often be the case over equatorial Africa due to frequent cloud cover and inconsistent coverage of Landsat-5.
Date-range: The time dimension represents calendar months aggregated across the period 1984-2020.
Spatial resolution: 30m
In this notebook we will load the NDVI Climatology product using
dc.load() to return mean, standard deviation, and clear observation count for each calendar month. A final section explores an example analysis using the product.
Topics covered include: 1. Inspecting the NDVI Climatology product and measurements available in the datacube. 2. Using the native
dc.load() function to load in NDVI Climatology. 3. Inspect perennial and annual vegtation using NDVI mean and standard deviation. 4. Exaple analysis: analyse and plot the climatological (long-term average) phenology of croplands
To run this analysis, run all the cells in the notebook, starting with the “Load packages” cell.
import datacube import numpy as np import xarray as xr import pandas as pd import matplotlib.pyplot as plt from odc.ui import with_ui_cbk from deafrica_tools.plotting import display_map import warnings warnings.simplefilter(action='ignore', category=FutureWarning)
/usr/local/lib/python3.8/dist-packages/geopandas/_compat.py:112: UserWarning: The Shapely GEOS version (3.8.0-CAPI-1.13.1 ) is incompatible with the GEOS version PyGEOS was compiled with (3.10.1-CAPI-1.16.0). Conversions between both will be slow. warnings.warn(
Connect to the datacube¶
dc = datacube.Datacube(app='NDVI_clim')
The table printed below shows the measurements available in the NDVI Climatology product. The mean NDVI, standard deviation of NDVI, and clear obervations count can be loaded for each calendar month.
product_name = 'ndvi_climatology_ls' dc_measurements = dc.list_measurements() dc_measurements.loc[product_name].drop('flags_definition', axis=1)
This section defines the analysis parameters, including:
lat, lon, buffer: center lat/lon and analysis window size for the area of interest
resolution: the pixel resolution to use for loading the
ndvi_climatology_ls. The native resolution of the product is 30 metres i.e.
(-30,30)as the product is Landsat derived.
The default location is an cropping region in Western Cape, South Africa where an irrigation scheme along a river is surrounded by rain-fed cropping.
lat, lon = -34.0602, 20.2800 # capetown buffer = 0.08 resolution=(-30, 30) #join lat, lon, buffer to get bounding box lon_range = (lon - buffer, lon + buffer) lat_range = (lat + buffer, lat - buffer)
View the selected location¶
The next cell will display the selected area on an interactive map. Feel free to zoom in and out to get a better understanding of the area you’ll be analysing. Clicking on any point of the map will reveal the latitude and longitude coordinates of that point.
Below we use the
dc.load function to load all the measurements over the region specified above.
# load data ndvi_clim = dc.load( product="ndvi_climatology_ls", resolution=resolution, x=lon_range, y=lat_range, progress_cbk=with_ui_cbk(), ) print(ndvi_clim)
<xarray.Dataset> Dimensions: (time: 1, y: 567, x: 515) Coordinates: * time (time) datetime64[ns] 2002-07-02T11:59:59.999999 * y (y) float64 -4.091e+06 -4.091e+06 ... -4.108e+06 -4.108e+06 * x (x) float64 1.949e+06 1.949e+06 ... 1.964e+06 1.964e+06 spatial_ref int32 6933 Data variables: (12/36) mean_jan (time, y, x) float32 0.4916 0.4292 0.4075 ... 0.1989 0.2001 mean_feb (time, y, x) float32 0.5281 0.4745 0.4531 ... 0.222 0.2238 mean_mar (time, y, x) float32 0.4731 0.416 0.3949 ... 0.2116 0.2121 mean_apr (time, y, x) float32 0.5536 0.5103 0.5034 ... 0.2504 0.2526 mean_may (time, y, x) float32 0.6418 0.6385 0.6137 ... 0.4022 0.3982 mean_jun (time, y, x) float32 0.6343 0.6148 0.6111 ... 0.5 0.5278 0.5339 ... ... count_jul (time, y, x) int16 30 30 30 30 30 31 33 ... 31 31 31 31 31 31 count_aug (time, y, x) int16 29 29 29 29 29 29 30 ... 23 23 23 23 23 23 count_sep (time, y, x) int16 33 33 33 33 33 33 33 ... 32 32 32 32 31 31 count_oct (time, y, x) int16 32 32 32 33 33 33 32 ... 30 30 30 31 31 31 count_nov (time, y, x) int16 36 36 36 36 37 37 37 ... 30 30 31 31 31 31 count_dec (time, y, x) int16 37 38 38 38 38 38 37 ... 35 35 35 34 35 35 Attributes: crs: EPSG:6933 grid_mapping: spatial_ref
Plotting NDVI Climatology¶
As this dataset has a lot of measurements, it is easier to store them in a list and access them as we need when plotting etc., rather than typing out each measurement
measurements = [i for i in ndvi_clim.data_vars] mean_bands = measurements[0:12] std_bands = measurements[12:24] count_bands = measurements[24:36]
Add a time dimension to the climatologies¶
Adding a time dimension to the dataset by concatenating the months together allows us to quickly make plots or calculate ket statistics without looping through separate datasets.
Important note: The date ‘2002’ is meaningless, its simply the middle point between 1984 and 2020 - the range over which the climatologies were calculated.
months = ['01','02','03','04','05','06','07','08','09','10','11','12'] arrs_mean =  arrs_std =  arrs_count =  for m,s,c,i in zip(mean_bands, std_bands, count_bands, months): #index out the different types of measurements xx_mean=ndvi_clim[m] xx_std=ndvi_clim[s] xx_count=ndvi_clim[c] #add a time dimension corresponding to the right month time = pd.date_range(np.datetime64('2002-'+i+'-15'), periods=1) xx_mean['time'] = time xx_std['time'] = time xx_count['time'] = time #append to lists arrs_count.append(xx_count) arrs_mean.append(xx_mean) arrs_std.append(xx_std) #concatenate into 3D data array along time dimension ndvi_mean = xr.concat(arrs_mean, dim='time').rename('NDVI_Climatology_Mean') ndvi_std = xr.concat(arrs_std, dim='time').rename('NDVI_Climatology_Std') ndvi_count = xr.concat(arrs_count, dim='time').rename('NDVI_Climatology_Count')
Facet plot the mean NDVI climatology¶
The plots below show that, on average, the irrigated fields around the river, and the riparian vegetetation are more persistently green throughout the year, while the rainfed cropping regions are green through the winter months of June, July, August, and September.
ndvi_mean.plot.imshow(col='time', col_wrap=6, cmap="RdYlGn");
Facet plot the standard deviation NDVI climatology¶
Note the irrigated regions tend to show greater variability than the other pixels. By enabling distinction between classes of agricultural land and other vegetation, this might be a good input to land use classification methods such as thresholding or more complex machine learning approaches.
ndvi_std.plot.imshow(col='time', col_wrap=6, cmap="magma", vmax=0.25);
Facet plot the clear count¶
We can use these plots to understand how many observations are going into the climatology calculations. Below (in the default example) we can that in the drier months there are ~25-35 observations (over the 37 year period of the climatology), while in the winter months there are ~20. In this case, because the observation count is relatvely low for the winter months, we should be careful when using the climatology in those months. In this specific example the mean and standard deviation measurements plotted above look okay, probably because the region is at a higher latitude and therefore not so affected by cloud. If this region was in the tropics then the dataset may be of a poorer quality.
ndvi_count.plot.imshow(col='time', col_wrap=6, cmap="viridis");
Example Analysis: Cropland Climatological Phenology¶
We can use the NDVI climatology product to generate long-term average phenological patterns for specific regions, landscapes, or vegetation types. Limiting our analysis of the area above to a crop mask enables us to investigate the long-term average phenology of cropland.
#import deafrica's phenology function from deafrica_tools.temporal import xr_phenology
Load the cropmask dataset for the region
cm = dc.load(product="crop_mask", measurements='filtered', resampling='nearest', like=ndvi_clim.geobox).filtered.squeeze() cm.plot.imshow();
Mask the datasets with the crop mask
ndvi_mean_crop = ndvi_mean.where(cm) ndvi_std_crop = ndvi_std.where(cm)
Plot phenology curve¶
Below we summarise the datasets spatially by taking the mean across the
y dimensions, this will leave us with the average trend through
time for the region we’ve loaded.
When we plot the phenology curve below, we will add +- 1 standard deviation around the mean curve to indicate, on average, how much the trends vary around the long-term mean.
From the phenology plot, we can conclude that crop growth generally commences around May and continues until harvest around October.
ndvi_mean_crop_line = ndvi_mean_crop.mean(['x','y']) ndvi_std_crop_line = ndvi_std_crop.mean(['x','y']) #upper and lower boundaries ndvi_std_upper = ndvi_mean_crop_line+ndvi_std_crop_line ndvi_std_lower = ndvi_mean_crop_line-ndvi_std_crop_line
Plot the mean phenology curve, with +-1 standard deviation envelope
ndvi_mean_crop_line.plot(figsize=(9,5), marker='o') plt.fill_between( ndvi_mean_crop_line.time.values, ndvi_std_upper.values, ndvi_std_lower.values, facecolor="blue", alpha=0.3, ) plt.ylim(0,0.9) plt.title('Climatological phenology averaged over the region') plt.xlabel('Month') plt.ylabel('NDVI') plt.tight_layout()
Per pixel climatological phenology¶
Having plotted in the spatially averaged climatological phenology, we can now use the DE Africa function
xr_phenology to calculate the typical phenology of every pixel in the dataset.
This function computes key phenological statistics expressed as the day-of-the-year (DOY) at: * Start of Season (SOS), * Peak of Season (POS), * End of Season (EOS).
It also gives NDVI values (v) for the times above (e.g. vSOS for the NDVI value at the start of season), and other variables including: * Trough: minimum NDVI value of season, * LOS: Length of season, * AOS: Amplitude of season (difference between the base and maximum NDVI values), * ROG: Rate of greening, * ROS: Rate of senescence.
Important note: Remember that this product is monthly (not daily) so the DOY is simply the middle of the month when the NDVI peaks, for example (we set the time dimension to be the 15th of each month in the “adding time dimension” section above).
#calculate per-pixel phenology phen = xr_phenology(ndvi_mean_crop, method_eos='median', method_sos='median', verbose=False) #mask again with crop-mask phen = phen.where(cm) print(phen)
<xarray.Dataset> Dimensions: (y: 567, x: 515) Coordinates: * y (y) float64 -4.091e+06 -4.091e+06 ... -4.108e+06 -4.108e+06 * x (x) float64 1.949e+06 1.949e+06 ... 1.964e+06 1.964e+06 spatial_ref int32 6933 time datetime64[ns] 2019-07-02T11:59:59.999999 Data variables: SOS (y, x) float64 nan nan nan nan 46.0 ... 46.0 46.0 46.0 46.0 POS (y, x) float64 nan nan nan nan ... 227.0 227.0 227.0 227.0 EOS (y, x) float64 nan nan nan nan ... 288.0 288.0 288.0 319.0 Trough (y, x) float32 nan nan nan nan ... 0.1941 0.1952 0.1989 0.2001 vSOS (y, x) float32 nan nan nan nan ... 0.2217 0.221 0.222 0.2238 vPOS (y, x) float32 nan nan nan nan ... 0.7151 0.7131 0.7291 0.7375 vEOS (y, x) float32 nan nan nan nan ... 0.3536 0.3585 0.3682 0.2471 LOS (y, x) float64 nan nan nan nan ... 242.0 242.0 242.0 273.0 AOS (y, x) float32 nan nan nan nan ... 0.5209 0.518 0.5302 0.5374 ROG (y, x) float32 nan nan nan nan ... 0.002719 0.002801 0.002838 ROS (y, x) float32 nan nan nan ... -0.005813 -0.005916 -0.005331 Attributes: grid_mapping: spatial_ref
Plot per pixel phenology¶
# set up figure fig, ax = plt.subplots(nrows=2, ncols=5, figsize=(20, 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 phen.SOS.plot.imshow(ax=ax[0, 0], cmap='magma_r', vmax=365, vmin=0, cbar_kwargs=dict(shrink=cbar_size, label=None)) ax[0, 0].set_title('Start of Season (DOY)') phen.vSOS.plot.imshow(ax=ax[0, 1], cmap='YlGn', vmax=0.8, cbar_kwargs=dict(shrink=cbar_size, label=None)) ax[0, 1].set_title('NDVI at SOS') # peak of season phen.POS.plot.imshow(ax=ax[0, 2], cmap='magma_r', vmax=365, vmin=0, cbar_kwargs=dict(shrink=cbar_size, label=None)) ax[0, 2].set_title('Peak of Season (DOY)') phen.vPOS.plot.imshow(ax=ax[0, 3], cmap='YlGn', vmax=0.8, cbar_kwargs=dict(shrink=cbar_size, label=None)) ax[0, 3].set_title('NDVI at POS') # end of season phen.EOS.plot.imshow(ax=ax[0, 4], cmap='magma_r', vmax=365, vmin=0, cbar_kwargs=dict(shrink=cbar_size, label=None)) ax[0, 4].set_title('End of Season (DOY)') phen.vEOS.plot.imshow(ax=ax[1, 0], cmap='YlGn', vmax=0.8, cbar_kwargs=dict(shrink=cbar_size, label=None)) ax[1, 0].set_title('NDVI at EOS') # Length of Season phen.LOS.plot.imshow(ax=ax[1, 1], cmap='magma_r', vmax=365, vmin=0, cbar_kwargs=dict(shrink=cbar_size, label=None)) ax[1, 1].set_title('Length of Season (DOY)') # Amplitude phen.AOS.plot.imshow(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.imshow(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.imshow(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('Climatological Phenology') plt.tight_layout();
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
Compatible datacube version:
from datetime import datetime datetime.today().strftime('%Y-%m-%d')