Monthly vegetation condition anomalies

• Products used: ndvi_climatology_ls, crop_mask, wofs_ls_summary_alltime, s2_ls2, ls8_sr, ls9_sr

Keywords: data used; ndvi_climatology_ls, vegetation; anomalies, band index; NDVI;

Contexte

Understanding how the vegetated landscape responds to longer-term environmental drivers such as the El Nino Southern Oscillation (ENSO) or climate change, requires the calculation of standardised anomalies. Standardised anomalies subtract the long-term mean from an observation of interest and then divide the result by the long-term standard deviation, thus removing seasonal variability and highlighting change related to longer-term drivers.

Description

This notebook will calculate monthly standardised NDVI anomalies for any given month and year, and relies on DE Africa’s 30 metre resolution NDVI climatology product to define the baseline vegetation condition. Sentinel-2, and Landsat 8/9 are used to calculate the monthly mean NDVI for the month of interest. The standardised anaomaly is then calculated as:

\begin{equation} \text{Standardised anomaly }=\frac{x-m}{s} \end{equation}

where \(x\) is NDVI for the month of interest, \(m\) is the long-term mean, and \(s\) is the long-term standard deviation.

The following steps are taken in the notebook: 1. Load a polygon for the country of interest (a geojson with country boundaries is already provided) 2. Load cloud masked Sentinel-2 and Landsat satellite data 3. Preprocess the satellite data so we have one continuous, noise free time-series of NDVI. 4. Create mask of the country boundary 5. Load DE Africa’s NDVI climatology product 6. Use WOfS to remove waterbodies from the region 7. Calculate standardised NDVI anomalies 8. Plot the result 9. Load DE Africa’s crop-mask for the region 10. Plot NDVI anomalies for regions only within the crop-mask

Getting started

To run this analysis, run all the cells in the notebook, starting with the « Load packages » cell.

After finishing the analysis, return to the « Analysis parameters » cell, modify some values (e.g. choose a different location or time period to analyse) and re-run the analysis.

Load packages

Load key Python packages and supporting functions for the analysis.

%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 xarray as xr
import numpy as np
import geopandas as gpd
import matplotlib.pyplot as plt
from datacube.utils import masking
from datacube.utils.geometry import Geometry, CRS
from deafrica_tools.spatial import xr_rasterize

from deafrica_tools.datahandling import load_ard
from deafrica_tools.bandindices import calculate_indices
from deafrica_tools.plotting import display_map
from deafrica_tools.dask import create_local_dask_cluster

Set up a Dask cluster

Dask can be used to better manage memory use and conduct the analysis in parallel. For an introduction to using Dask with Digital Earth Africa, see the Dask notebook.

Note: We recommend opening the Dask processing window to view the different computations that are being executed; to do this, see the Dask dashboard in DE Africa section of the Dask notebook.

To use Dask, set up the local computing cluster using the cell below.

create_local_dask_cluster()

Connect to the datacube

Activate the datacube database, which provides functionality for loading and displaying stored Earth observation data.

dc = datacube.Datacube(app="Vegetation_anomalies")

Analysis parameters

The following cell sets important parameters for the analysis:

• country: The country over which to load data, e.g. 'Uganda'. The full list of Africa countries can be found in cell [6]

• year: The year for which we will calculate the standardised anomaly, e.g. '2022'

• month: The month for which we will calculate the standardised anomaly, written as a lower-case three letter abbreviation e.g 'jan','feb' etc. This is because the ndvi_climatology_ls product’s measurements are stored as 'mean_jan', 'mean_feb' etc.

• crop_mask: Depending on where you are loading data, you can change the crop mask to load with this parameter. Use the explorer to find which crop mask covers which region of Africa.

• resolution: The x and y cell resolution of the satellite data and ndvi_climatology_ls data, e.g. (-30, 30) will load data at its native 30 x 30m resolution. If loading a large area, then increase the resolution so the data will fit into memory. For example, if trying to find the NDVI anomaly for a large country, consider increasing the resolution to (-210,210).

• dask_chunks: the size of the dask chunks, dask breaks data into manageable chunks that can be easily stored in memory, e.g. dict(x=1000,y=1000)

If running the notebook for the first time, keep the default settings below. This will demonstrate how the analysis works and provide meaningful results.

# Select a country
country = "Uganda"

# year and month for anomaly
year = "2022"
month = "jan"

# the crop mask to use
crop_mask = "crop_mask"

# pixel resolution can got as low as (-30,30) metres,
# but make larger for big areas.
resolution = (-120, 120)

# how to chunk the dataset for use with dask
dask_chunks = dict(x=1000, y=1000)

Load the African Countries shapefile

This shapefile contains polygons for the boundaries of African countries and will allows us to calculate NDVI anomalies within a chosen country

african_countries = gpd.read_file('../Supplementary_data/Rainfall_anomaly_CHIRPS/african_countries.geojson')

List countries

You can change the country in the analysis parameters cell to any African country. A complete list of countries is printed below.

print(np.unique(african_countries.COUNTRY))

['Algeria' 'Angola' 'Benin' 'Botswana' 'Burkina Faso' 'Burundi' 'Cameroon'
 'Cape Verde' 'Central African Republic' 'Chad' 'Comoros'
 'Congo-Brazzaville' 'Cote d`Ivoire' 'Democratic Republic of Congo'
 'Djibouti' 'Egypt' 'Equatorial Guinea' 'Eritrea' 'Ethiopia' 'Gabon'
 'Gambia' 'Ghana' 'Guinea' 'Guinea-Bissau' 'Kenya' 'Lesotho' 'Liberia'
 'Libya' 'Madagascar' 'Malawi' 'Mali' 'Mauritania' 'Morocco' 'Mozambique'
 'Namibia' 'Niger' 'Nigeria' 'Rwanda' 'Sao Tome and Principe' 'Senegal'
 'Sierra Leone' 'Somalia' 'South Africa' 'Sudan' 'Swaziland' 'Tanzania'
 'Togo' 'Tunisia' 'Uganda' 'Western Sahara' 'Zambia' 'Zimbabwe']

Setup polygon

The country selected needs to be transformed into a geometry object to be used in the load_ard function.

idx = african_countries[african_countries['COUNTRY'] == country].index[0]
geom = Geometry(geom=african_countries.iloc[idx].geometry, crs=african_countries.crs)

Load and process cloud-masked satellite data

The first step in this analysis is to load in Landsat 8/9 and Sentinel-2 data for the country we’ve picked, and for the month of the anomaly we want to assess. The code below uses the load_ard function to load in data from Landsat 8/9 and Sentinel-2 satellites for the area and time specified. For more information, see the Using load_ard notebook.

Firstly, we will set up a dictionary that maps the abbreviated month name, month to a string number, allowing us to load data for that month from the datacube. The second dictionary, query maps the parameters we provided above to the datacube.load() parameters.

# map months to string numbers
months = {
    "jan": "01",
    "feb": "02",
    "mar": "03",
    "apr": "04",
    "may": "05",
    "jun": "06",
    "jul": "07",
    "aug": "08",
    "sep": "09",
    "oct": "10",
    "nov": "11",
    "dec": "12",
}

# Create the 'query' dictionary object for loading data form the datacube
query = {
    "geopolygon": geom,
    "time": year + "-" + months[month],
    "resolution": resolution,
    "output_crs": "epsg:6933",
    "dask_chunks": dask_chunks,
}

Now let’s load the satellite data.

# Load available data Landsat 8 and 9
ds_ls = load_ard(dc=dc,
              products=['ls8_sr', 'ls9_sr'],
              group_by='solar_day',
              measurements= ['red','nir'],
              resampling={'pixel_quality':'nearest', '*':'average'},
              mask_filters=[("opening", 3),("dilation", 3)], #improve cloud-mask
              verbose=False,
              **query)

# Load available sentinel2
ds_s2 = load_ard(dc=dc,
              products=['s2_l2a'],
              measurements=['red','nir_2'], #use nir narrow to match with LS8 nir
              group_by='solar_day',
              resampling={'pixel_quality':'nearest', '*':'average'},
              mask_filters=[("opening", 3),("dilation", 3)], #improve cloud-mask
              verbose=False,
              **query)

Preprocess satellite data

Here we perform a number of steps to get the satellite data into the format we need for calculating an NDVI anomaly. 1. Rename the 'nir_2' band in Sentinel-2 to 'nir' so the calculate_indices function will understand to calculate NDVI using 'nir_2'. 2. Calculate NDVI on both Landsat and Sentinel-2 data, then we concatenate the two time-series together. 3. Now that we one continous time-series of NDVI, we can clean the time-series by removing data that is outside the range 0-1, and minimize the impact of noise/missed cloud by running a temporal rolling mean filter with a small window size. 4. Calculate the mean NDVI for the month.

Note that this cell will take several minutes to run (depending on size and resolution of the datasets) as here we trigger all of the calculations we’ve listed up to this point with the .compute() method

#rename to trick 'calculate_indices'
ds_s2 = ds_s2.rename({'nir_2':'nir'})

#calculate NDVI
ndvi_s2 = calculate_indices(ds_s2, 'NDVI', satellite_mission='s2', drop=True)
ndvi_ls = calculate_indices(ds_ls, 'NDVI', satellite_mission='ls', drop=True)

#concatente together
ndvi = xr.concat([ndvi_ls,ndvi_s2], dim='time').sortby('time')

# Remove NDVI's that aren't between 0 and 1
ndvi = ndvi.where((ndvi >= 0) & (ndvi <= 1))

# smooth out noise from missed cloud
ndvi = ndvi.rolling(time=3, min_periods=1).mean()

# summarise into a monthly mean
month_mean = ndvi.NDVI.mean('time').compute()

Dropping bands ['red', 'nir']
Dropping bands ['red', 'nir']

/usr/local/lib/python3.10/dist-packages/rasterio/warp.py:344: NotGeoreferencedWarning: Dataset has no geotransform, gcps, or rpcs. The identity matrix will be returned.
  _reproject(

Create a rasterized mask of the country boundary

african_countries = african_countries.to_crs('epsg:6933')
mask = xr_rasterize(african_countries[african_countries['COUNTRY'] == country], ds_ls)

Load DE Africa’s NDVI climatologies

ndvi_clim = dc.load(
    product="ndvi_climatology_ls",
    like=ds_ls.geobox,
    measurements=['mean_'+month, 'stddev_'+month, 'count_'+month],
    dask_chunks=dask_chunks,
).squeeze()

Make a quality assurance mask where clear observation count is low

Due to the inconsistent data availability of Landsat 5 over equatorial Africa, and due to persistent cloud coverage over these same regions, the quality of the long-term NDVI baseline is poor in in some parts of Africa. We recommend not using/relying on the product in locations where the clear obervation count is less than ~20-30 observations. Below this, the standard deviation layers are prone to data artefacts owing to temporal smoothing operations performing poorly on sparse datasets, and residual cloudy pixels not being sufficiently ‘averaged out’ by a decent volume of input data.

min_num_obs = 20

qa_mask = ndvi_clim['count_'+month] >= min_num_obs

# remove pixels where obs are < min_num_obs
ndvi_clim = ndvi_clim.where(qa_mask).compute()

Mask both datasets with country mask

#mask with country boundary
ndvi_clim = ndvi_clim.where(mask==1)
month_mean = month_mean.where(mask==1)

Remove water bodies

This will constrain all the NDVI anomaly observations to just the land mass.

#load wofs
wofs = dc.load(product='wofs_ls_summary_alltime',
               measurements=['frequency'],
               like=ds_ls.geobox,
               dask_chunks=dask_chunks
              ).frequency.squeeze()

#set masked terrain regions to 0
wofs = xr.where(np.isnan(wofs), 0, wofs)

#threshold to create waterbodies mask
wofs = (wofs < 0.5)

#mask with country boundary
wofs = wofs.where(mask).compute()

#mask with wofs boundary
ndvi_clim = ndvi_clim.where(wofs)
mean_clim = ndvi_clim['mean_'+month].where(wofs)
std_clim = ndvi_clim['stddev_'+month].where(wofs)
month_mean = month_mean.where(wofs)

Calculate standardised anomalies

Now we can calculate the standardised anomalies by subtracting the long-term mean and dividing by the long-term standard deviation.

stand_anomalies = xr.apply_ufunc(
    lambda x, m, s: (x - m) / s,
    month_mean,
    mean_clim,
    std_clim,
    output_dtypes=[month_mean.dtype],
    dask="allowed"
)

Plot NDVI climatology, monthly mean, and standardised anomaly

# plot all layers
plt.rcParams["axes.facecolor"] = "gray"  # makes transparent pixels obvious
fig, ax = plt.subplots(2, 2, sharey=True, sharex=True, figsize=(18, 15))

mean_clim.plot.imshow(ax=ax[0, 0],cmap="RdYlGn",vmin=0,vmax=0.8)
ax[0, 0].set_title("NDVI: " + month + " mean climatology")

std_clim.plot.imshow(ax=ax[0, 1], cmap='magma', vmin=0, vmax=0.25)
ax[0, 1].set_title("NDVI: " + month + " std dev climatology")

month_mean.plot.imshow(ax=ax[1, 0], cmap="RdYlGn",vmin=0,vmax=0.8)
ax[1, 0].set_title("NDVI: " + year + " " + month + " mean")

stand_anomalies.plot.imshow(ax=ax[1, 1], cmap="BrBG", vmin=-3, vmax=3)
ax[1, 1].set_title("NDVI: " + year + " " + month + " standardised anomaly")
plt.tight_layout();

Incorporating DE Africa’s cropland extent map

Load the cropland mask over the region of interest. The default region we’re analysing is in Ethiopia, so we need to load either the crop_mask product which covers the entire African continent, or the crop_mask_eastern product, which cover the countries of Ethiopia, Kenya, Tanzania, Rwanda, and Burundi. If you change the analysis region from the default one, you may need to load a different crop mask - see the docs page to find out more.

cm = dc.load(product=crop_mask,
            time=('2019'),
            measurements='mask',
            resampling='nearest',
            like=ds_ls.geobox).mask.squeeze()

# convert the missing values (255) to NaN
cm = cm.where(cm != 255)

cm.plot.imshow(add_colorbar=False, figsize=(7,7))
plt.title('Cropland Extent');

Plot NDVI climatolgies, monthly mean, and anomalies for cropping regions only

Below we mask out the regions that aren’t cropping, revealing only the condition of the cropped regions.

#mask layers with crop-mask
mean_clim_crop=mean_clim.where(cm==1, np.nan)
std_clim_crop=std_clim.where(cm==1, np.nan)
month_mean_crop=month_mean.where(cm==1, np.nan)
stand_anomalies_crop=stand_anomalies.where(cm==1, np.nan)

#plot al layers
plt.rcParams['axes.facecolor'] = 'gray' # makes transparent pixels obvious
fig,ax = plt.subplots(2,2, sharey=True, sharex=True, figsize=(15,15))

mean_clim_crop.plot.imshow(ax=ax[0,0], cmap='RdYlGn' ,vmin=0, vmax=0.8)
ax[0,0].set_title('NDVI: '+month+' mean climatology')

std_clim_crop.plot.imshow(ax=ax[0,1], cmap='magma', vmin=0, vmax=0.2)
ax[0,1].set_title('NDVI: '+month+' std dev climatology')

month_mean_crop.plot.imshow(ax=ax[1,0], cmap="RdYlGn",vmin=0,vmax=0.8)
ax[1,0].set_title('NDVI: '+year+" "+month+' mean')

stand_anomalies_crop.plot.imshow(ax=ax[1,1], cmap='BrBG',vmin=-3, vmax=3)
ax[1,1].set_title('NDVI: '+year+" "+month+' standardised anomaly')
plt.tight_layout();

---

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:

print(datacube.__version__)

1.8.15

Last Tested:

from datetime import datetime
datetime.today().strftime('%Y-%m-%d')

'2023-08-11'