Phenological analysis for image time series

1 Introduction

This module implements a several algorithms allowing to extract phenological information from time profiles. These time profiles should represent vegetation status as for instance NDVI, LAI, etc.

2 Library

The library provides tools for fitting parametric models to time profiles:

  • double logistics
  • assymetric gaussians

Only the double logistic has been thoroughly tested. A double logistic has this form:

\begin{equation} g(x) = A\left(f_1(x)-f_2(x)\right)+B = A\left(\frac{1}{1+e\frac{x_0-x{x_1}}}-\frac{1}{1+e\frac{x_2-x{x_3}}}\right)+B, \end{equation}

where $A+B$ is the maximum value and $B$ is the minimum.

From the double logistic fitting, some key parametres can be obtained. The 6 parameters of the model can be used to define the following phenological metrics:


\begin{equation} \frac{df(x)}{dx} = \frac{e\frac{x_0-x{x_1}}}{x_1\left(1+e\frac{x_0-x{x_1}} \right)^2} \end{equation}

we have

\begin{equation} g’(x)=\frac{dg(x)}{dx} = A\left( \frac{e\frac{x_0-x{x_1}}}{x_1\left(1+e\frac{x_0-x{x_1}} \right)^2} - \frac{e\frac{x_2-x{x_3}}}{x_3\left(1+e\frac{x_2-x{x_3}} \right)^2}\right) \end{equation}

set grid on
set xrange [1:365]
set yrange [-0.2:1.2]
set xlabel "DoY"
set ylabel "Vegetation Index"
set arrow from x0,0 to x0,1 nohead lw 2 lc rgb 'blue'
set arrow from x2,0 to x2,1 nohead lw 2 lc rgb 'blue'
set label "x_0" at x0-10,-0.05
set label "x_2" at x2-10,-0.05
set arrow from 56,0 to 90,1.05 nohead lw 2 lc rgb 'green'
set label "t_{0}" at 45,-0.05
set label "t_1" at 85,1.1
set arrow from 90,1.15 to 225,1.15 heads lw 2 lc rgb 'green'
set label "L" at 150,1.10
set arrow from 230,1.05 to 275,0 nohead lw 2 lc rgb 'green'
set label "t_2" at 230,1.1
set label "t_3" at 275,-0.05
plot A*(1.0/(1.0+exp((x0-x)/x1))-1.0/(1.0+exp((x2-x)/x3)))+B t ""


2.1 Date of the maximum positive gradient

By definition, this is $x_0$.

2.2 Starting date

The date for which the straight line with the slope of $x_0$ intercepts the horizontal axis. If

\begin{equation} g’(x_0) = m \end{equation}

this line’s equation is

\begin{equation} y = m x + b \end{equation}

and verifies that

\begin{equation} \begin{split} 0 = & m t_0 + b
g(x_0) = & m x_0 + b \end{split} \end{equation}

which gives

\begin{equation} y = g’(x_0) x + (g(x_0) - g’(x_0) x_0) \end{equation}

and therefore,

\begin{equation} t_0 = \frac{m x_0 - g(x_0) }{m} = x_0 - \frac{g(x_0)}{g’(x_0)} \end{equation}

2.3 Length of the plateau

We define $t_1$ as the date for which the previous straight line reaches the maximum value.

\begin{equation} t_1 = \frac{A+B - (g(x_0) - g’(x_0) x_0)}{g’(x_0)} \end{equation}

Similarly, we can use the straight line associated to the descending slope:

\begin{equation} y = g’(x_2) x + (g(x_2)-g’(x_2)x_2) \end{equation}

and define

\begin{equation} t_2 = \frac{A+B - (g(x_2) - g’(x_2) x_2)}{g’(x_2)} \end{equation}

And the length of the plateau is:

\begin{equation} L = t_2 - t_1 \end{equation}

2.4 Senescense slope

By definition this is g’(x_2).

3 Applications

Only one application is provided at this time. It allows fitting double logistics to each pixel of an image time series. The output contains 2 double logistics, one for the main phenological cycle and another one for a secondary cycle. This secondary cycle may not be present in the input data. This should not have any impact in the estimation of the main cycle.

The application can generate a 12 band image, where each band is one of the parameters of the 2 double logistics in the order $A$, $B$, $x_0$, $x_1$, $x_2$, $x_3$.

The application can also instead generate an output image with the same number of bands as the input image, but replacing the value of each date for each pixel by the value taken by the logistic fitting.

Finally, the application can output an image where each band is one of the phenological metrics for the 2 cycles. The order of the metrics is $g’(x_0)$, $t_0$, $t_1$, $t_2$, $t_3$, $g’(x_2)$.