Radar vegetation phenology using Sentinel-1

Keywords data used; sentinel-1, band index; RVI, phenology, analysis; time series

Background

Phenology is the study of plant and animal life cycles in the context of the seasons. It can be useful in understanding the life cycle trends of crops and how the growing seasons are affected by changes in climate. However, in cloudy regions optical satellites may not view the vegetation cycle during key moments of growth (e.g. during the peak of the growing sesaon) leading to inaccurate estimates of vegetation phenology variables. The Sentinel-1 radar satellite ‘sees’ through clouds and can, in theory, view the entire life cycle of plants. This notebook will use Sentinel-1 to derive key phenology statistics. You can compare these results with the results from the Sentinel-2 vegetation phenology notebook which relies on optical data to calculate phenology over the same region as this notebook

Description

This notebook demonstrates how to calculate vegetation phenology statistics over cropping regions using the DE Africa function xr_phenology. To detect changes in plant life for Sentinel-1, the script uses the Radar Vegetation Index

The outputs of this notebook can be used to assess spatio-temporal differences in the growing seasons of agriculture fields or native vegetation.

This notebook demonstrates the following steps:

  1. Load Sentinel 1 data for an area of interest.

  2. Mask region with DE Africa’s cropland extent map

  3. Calculate the radar vegetation index.

  4. Plot false colour images of the region.

  5. Compute the radar vegetation index

  6. Smooth the vegetation timeseries to minimize noise.

  7. Calculate phenology statistics on a simple 1D vegetation time series

  8. Calculate per-pixel phenology statistics


Getting started

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

Load packages

Load key Python packages and supporting functions for the analysis.

[1]:
%matplotlib inline

import os
import datacube
import numpy as np
import datetime as dt
import matplotlib.pyplot as plt

from deafrica_tools.temporal import xr_phenology
from deafrica_tools.datahandling import load_ard
from deafrica_tools.plotting import display_map, rgb
from deafrica_tools.bandindices import dualpol_indices
from deafrica_tools.classification import HiddenPrints

from datacube.utils.aws import configure_s3_access
configure_s3_access(aws_unsigned=True, cloud_defaults=True)
/env/lib/python3.8/site-packages/datacube/storage/masking.py:7: DeprecationWarning: datacube.storage.masking has moved to datacube.utils.masking
  warnings.warn("datacube.storage.masking has moved to datacube.utils.masking",
/env/lib/python3.8/site-packages/geopandas/_compat.py:106: UserWarning: The Shapely GEOS version (3.8.0-CAPI-1.13.1 ) is incompatible with the GEOS version PyGEOS was compiled with (3.9.1-CAPI-1.14.2). Conversions between both will be slow.
  warnings.warn(

Connect to the datacube

Connect to the datacube so we can access DE Africa data. The app parameter is a unique name for the analysis which is based on the notebook file name.

[2]:
dc = datacube.Datacube(app='Radar_vegetation_phenology')

Analysis parameters

The following cell sets important parameters for the analysis:

  • lat: The central latitude to analyse (e.g. -10.6996).

  • lon: The central longitude to analyse (e.g. 35.2708).

  • buffer: The number of square degrees to load around the central latitude and longitude. For reasonable loading times, set this as 0.1 or lower.

  • time_range: The year range to analyse (e.g. ('2019-01', '2020-12')).

[3]:
# Define area of interest
lat = 8.7186
lon = 40.8646

lon_buffer = 0.02
lat_buffer = 0.02

# Set the range of dates for the analysis
time_range = ('2019-01-01', '2020-12-20')

# Combine central lat,lon with buffer to get area of interest
lat_range = (lat-lat_buffer, lat+lat_buffer)
lon_range = (lon-lon_buffer, lon+lon_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.

[4]:
display_map(x=lon_range, y=lat_range)
[4]:

Load Sentinel-1 data

The first step is to load Sentinel-1 data for the specified area of interest and time range. The load_ard function is used here to load data that has been masked using quality filters, making it ready for analysis.

[5]:
# Create a reusable query
query = {
    'y': lat_range,
    'x': lon_range,
    'time': time_range,
    'measurements': ['vv', 'vh'],
    'resolution': (-20, 20),
    'output_crs': 'epsg:6933',
    'group_by':'solar_day'
}

# Load available data from Sentinel-1
ds = load_ard(
    dc=dc,
    products=['s1_rtc'],
    **query,
)

print(ds)
Using pixel quality parameters for Sentinel 1
Finding datasets
    s1_rtc
Applying pixel quality/cloud mask
Loading 111 time steps
<xarray.Dataset>
Dimensions:      (time: 111, x: 194, y: 253)
Coordinates:
  * time         (time) datetime64[ns] 2019-01-01T03:01:22.069586 ... 2020-12...
  * y            (y) float64 1.111e+06 1.111e+06 ... 1.106e+06 1.106e+06
  * x            (x) float64 3.941e+06 3.941e+06 ... 3.945e+06 3.945e+06
    spatial_ref  int32 6933
Data variables:
    vv           (time, y, x) float32 0.105032705 0.16441293 ... 0.172195
    vh           (time, y, x) float32 0.02470629 0.024684427 ... 0.02133612
Attributes:
    crs:           epsg:6933
    grid_mapping:  spatial_ref

Mask region with 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 the crop_mask_eastern product, which cover the countries of Ethiopia, Kenya, Tanzania, Rwanda, and Burundi

[6]:
cm = dc.load(product='crop_mask_eastern',
            time=('2019'),
            measurements='filtered',
            resampling='nearest',
            like=ds.geobox).filtered.squeeze()

cm.plot.imshow(add_colorbar=False, figsize=(6,6))
plt.title('Cropland Extent');
../../../_images/sandbox_notebooks_Real_world_examples_Phenology_radar_17_0.png

Now we will use the cropland map to mask the regions in the Sentinel-1 data that only have cropping

[7]:
ds = ds.where(cm)

False colour RGB plots

Now can plot some of the Sentinel-1 images as false colour images using the rgb function. Notice that the backscatter intensity of the cropping regions is changing throughout the year as the plants in the fields go through their growth and senescence cycle.

[8]:
# VH/VV will help create an RGB image
ds['vh/vv'] = ds.vh/ds.vv

# median values are used to scale the measurements so they have a similar range for visualization
med_s1 = ds[['vv','vh','vh/vv']].median()
[9]:
# plotting an RGB image for selected timesteps
time_steps = [1,10,22,45]

rgb(ds[['vv', 'vh', 'vh/vv']]/med_s1, #normalize image ranges
    bands=['vv','vh', 'vh/vv'],
    index=time_steps,
    col_wrap=4);
../../../_images/sandbox_notebooks_Real_world_examples_Phenology_radar_22_0.png

Compute the radar vegetation index

This study measures the presence of vegetation through either the Radar vegetation index (RVI).

The formula is

RVI = 4*VH/(VV+VH)

RVI is available through the dualpol_indices function, imported from deafrica_tools.bandindices.

[10]:
# Calculate the chosen vegetation proxy index and add it to the loaded data set
ds = dualpol_indices(ds, index='RVI')

Plot the vegetation index over time

To get an idea of how the vegetation health changes throughout the year(s), we can plot a zonal time series over the region of interest. First we will do a simple plot of the zonal mean of the data.

[11]:
ds.RVI.mean(['x', 'y']).plot.line('b-^', figsize=(11,4))
plt.title('Zonal mean of radar vegetation index timeseries');
../../../_images/sandbox_notebooks_Real_world_examples_Phenology_radar_27_0.png

Smoothing/Interpolating vegetation time-series

Here, we will smooth and interpolate the data to ensure we working with a consistent time-series. This is a very important step in the workflow and there are many ways to smooth, interpolate, gap-fill, remove outliers, or curve-fit the data to ensure a useable time-series. If not using the default example, you may have to define additional methods to those used here.

To do this we take two steps:

  1. Resample the data to fortnightly time-steps using the fortnightly median

  2. Calculate a rolling mean with a window of 4 steps

[12]:
resample_period='2W'
window=4

rvi_smooth=ds.RVI.resample(time=resample_period).median().rolling(time=window, min_periods=1).mean()

Plot the smoothed and interpolated time-series

[13]:
rvi_smooth_1D = rvi_smooth.mean(['x', 'y'])
rvi_smooth_1D.plot.line('b-^', figsize=(15,5))
_max=rvi_smooth_1D.max()
_min=rvi_smooth_1D.min()
plt.vlines(np.datetime64('2019-01-01'), ymin=_min, ymax=_max)
plt.vlines(np.datetime64('2020-01-01'), ymin=_min, ymax=_max)
plt.vlines(np.datetime64('2021-01-01'), ymin=_min, ymax=_max)
plt.title('RVI time-series, year start/ends marked with vertical lines')
plt.ylabel('RVI');
../../../_images/sandbox_notebooks_Real_world_examples_Phenology_radar_31_0.png

Calculate phenology statistics using xr_phenology

The DE Africa function xr_phenology can calculate a number of land-surface phenology statistics that together describe the characteristics of a plant’s lifecycle. The function can calculate the following statistics on either a zonal timeseries (like the one above), or on a per-pixel basis:

SOS = DOY of start of season
POS = DOY of peak of season
EOS = DOY of end of season
vSOS = Value at start of season
vPOS = Value at peak of season
vEOS = Value at end of season
Trough = Minimum value of season
LOS = Length of season (DOY)
AOS = Amplitude of season (in value units)
ROG = Rate of greening
ROS = Rate of senescence

where DOY = day-of-year (Jan 1st = 0, Dec 31st = 365). By default the function will return all the statistics as an xarray.Dataset, to return only a subset of these statistics pass a list of the desired statistics to the function e.g. stats=['SOS', 'EOS', 'ROG'].

See the deafrica_tools.temporal script for more information on each of the parameters in xr_phenology.

Zonal phenology statistics

To help us understand what these statistics refer too, lets first pass the simpler zonal mean (mean of all pixels in the image) time-series to the function and plot the results on the same curves as above.

First, provide a list of statistics to calculate with the parameter, pheno_stats.

method_sos : If ‘first’ then vSOS is estimated as the first positive slope on the greening side of the curve. If ‘median’, then vSOS is estimated as the median value of the postive slopes on the greening side of the curve.

method_eos : If ‘last’ then vEOS is estimated as the last negative slope on the senescing side of the curve. If ‘median’, then vEOS is estimated as the ‘median’ value of the negative slopes on the senescing side of the curve.

[14]:
pheno_stats = ['SOS','vSOS','POS','vPOS','EOS','vEOS','Trough','LOS','AOS','ROG','ROS']
method_sos = 'first'
method_eos = 'last'
[15]:
# find all the years to assist with plotting
years=rvi_smooth_1D.groupby('time.year')

# get list of years in ts to help with looping
years_int=[y[0] for y in years]

#store results in dict
pheno_results = {}

#loop through years and calculate phenology
for year in years_int:
    #select year
    da = dict(years)[year]

    #calculate stats
    stats=xr_phenology(
            da,
            method_sos=method_sos,
            method_eos=method_eos,
            stats=pheno_stats,
            verbose=False
                )
    #add results to dict
    pheno_results[str(year)] = stats

for key,value in pheno_results.items():
    print('Year: ' +key)
    for b in value.data_vars:
        print("   "+b+": ", round(float(value[b].values),3))
Year: 2019
   SOS:  160.0
   vSOS:  0.421
   POS:  286.0
   vPOS:  0.835
   EOS:  356.0
   vEOS:  0.785
   Trough:  0.421
   LOS:  196.0
   AOS:  0.414
   ROG:  0.003
   ROS:  -0.001
Year: 2020
   SOS:  159.0
   vSOS:  0.552
   POS:  299.0
   vPOS:  0.798
   EOS:  355.0
   vEOS:  0.718
   Trough:  0.548
   LOS:  196.0
   AOS:  0.25
   ROG:  0.002
   ROS:  -0.001

Plot the results with our statistcs annotated on the plot

[16]:
# find all the years to assist with plotting
years=rvi_smooth_1D.groupby('time.year')

fig, ax = plt.subplots()
fig.set_size_inches(15,7)
for year, y in zip(years, years_int):
    #grab tall the values we need for plotting
    eos = pheno_results[str(y)].EOS.values
    sos = pheno_results[str(y)].SOS.values
    pos = pheno_results[str(y)].POS.values

    veos = pheno_results[str(y)].vEOS.values
    vsos = pheno_results[str(y)].vSOS.values
    vpos = pheno_results[str(y)].vPOS.values

    #create plot
    #the groupby and mean below doesn't actually do anything, except allow
    #the plots to be on the same x-axis. i.e. veg-data isn't modified, just time units
    year[1].groupby('time.dayofyear').mean().plot(ax=ax, label=year[0])

    #add start of season
    ax.plot(sos, vsos, 'or')
    ax.annotate('SOS',
                xy=(sos, vsos),
                xytext=(-15, 20),
                textcoords='offset points',
                arrowprops=dict(arrowstyle='-|>'))

    #add end of season
    ax.plot(eos, veos, 'or')
    ax.annotate('EOS',
                xy=(eos, veos),
                xytext=(0, 20),
                textcoords='offset points',
                arrowprops=dict(arrowstyle='-|>'))

    #add peak of season
    ax.plot(pos, vpos, 'or')
    ax.annotate('POS',
                xy=(pos, vpos),
                xytext=(-10, -25),
                textcoords='offset points',
                arrowprops=dict(arrowstyle='-|>'))

ax.legend()
plt.ylim(_min-0.025,_max.values+0.05)
month_abbr=['Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec', '']
new_m=[]
for m in month_abbr:
    new_m.append('       %s'%m) #Add spaces before the month name
plt.xticks(np.linspace(0,365,13), new_m, horizontalalignment='left')
plt.xlabel('Month')
plt.ylabel('RVI')
plt.title('Yearly RVI');
../../../_images/sandbox_notebooks_Real_world_examples_Phenology_radar_37_0.png

Per-pixel phenology statistics

We can now calculate the statistics for every pixel in our time-series and plot the results.

[17]:
# find all the years to assist with plotting
years=rvi_smooth.groupby('time.year')

# get list of years in ts to help with looping
years_int=[y[0] for y in years]

#store results in dict
pheno_results = {}

#loop through years and calculate phenology
for year in years_int:
    #select year
    da = dict(years)[year]

    #calculate stats
    stats=xr_phenology(
            da,
            method_sos=method_sos,
            method_eos=method_eos,
            stats=pheno_stats,
            verbose=False
                )
    #add results to dict
    pheno_results[str(year)] = stats

The phenology statistics have been calculated seperately for every pixel in the image. Let’s plot each of them to see the results.

Below, pick a year from the phenology results to plot.

[18]:
#Pick a year to plot
year_to_plot = '2019'

At the top if the plotting code we re-mask the phenology results with the crop-mask. This is because xr_phenologyhas methods for handling pixels with only NaNs (such as those regions outside of the polygon mask), so the results can have phenology results for regions outside the mask. We will therefore have to mask the data again.

[19]:
#select the year to plot
phen = pheno_results[year_to_plot]

#mask again with crop-mask
phen = phen.where(cm)

# set up figure
fig, ax = plt.subplots(nrows=2,
                       ncols=5,
                       figsize=(18, 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(ax=ax[0, 0],
              cmap='magma_r',
              vmax=300,
              vmin=0,
              cbar_kwargs=dict(shrink=cbar_size, label=None))
ax[0, 0].set_title('Start of Season (DOY)')

phen.vSOS.plot(ax=ax[0, 1],
               cmap='YlGn',
               vmax=0.8,
               cbar_kwargs=dict(shrink=cbar_size, label=None))
ax[0, 1].set_title('RVI at SOS')

# peak of season
phen.POS.plot(ax=ax[0, 2],
              cmap='magma_r',
              vmax=300,
              vmin=0,
              cbar_kwargs=dict(shrink=cbar_size, label=None))
ax[0, 2].set_title('Peak of Season (DOY)')
phen.vPOS.plot(ax=ax[0, 3],
               cmap='YlGn',
               vmax=0.8,
               cbar_kwargs=dict(shrink=cbar_size, label=None))
ax[0, 3].set_title('RVI at POS')

# end of season
phen.EOS.plot(ax=ax[0, 4],
              cmap='magma_r',
              vmax=300,
              vmin=0,
              cbar_kwargs=dict(shrink=cbar_size, label=None))
ax[0, 4].set_title('End of Season (DOY)')
phen.vEOS.plot(ax=ax[1, 0],
               cmap='YlGn',
               vmax=0.8,
               cbar_kwargs=dict(shrink=cbar_size, label=None))
ax[1, 0].set_title('RVI at EOS')

# Length of Season
phen.LOS.plot(ax=ax[1, 1],
              cmap='magma_r',
              vmax=300,
              vmin=0,
              cbar_kwargs=dict(shrink=cbar_size, label=None))
ax[1, 1].set_title('Length of Season (DOY)')

# Amplitude
phen.AOS.plot(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(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(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('Radar derived phenology for '+year_to_plot)
plt.tight_layout();
../../../_images/sandbox_notebooks_Real_world_examples_Phenology_radar_43_0.png

Conclusions

In the example above, we can see most of the fields are following a similar cropping schedule and are therefore likely the same species of crop. We can also observe that some fields have not followed this schedule (e.g. the EOY plot shows some fields didn’t follow this schedule, probably because they either remained fallow during the year, or because there were two peaks during the year which may have confused the phenology stats). Differences in the rates of growth, and in the RVI values at different times of the season, may be attributable to differences in soil quality, watering intensity, or other farming practices.

Phenology statistics are a powerful way to summarise the seasonal cycle of a plant’s life. Per-pixel plots of phenology can help us understand the timing of vegetation growth and sensecence across large areas and across diverse plant species as every pixel is treated as an independent series of observations. This could be important, for example, if we wanted to assess how the growing seasons are shifting as the climate warms.


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:

[20]:
print(datacube.__version__)
1.8.5

Last Tested:

[21]:
from datetime import datetime
datetime.today().strftime('%Y-%m-%d')
[21]:
'2021-09-15'