# BFASTm-L2, an unsupervised LULCC detection based on seasonal change detection - an application to large-scale land acquisitions in Senegal

Y.N. Scarpetta, V. Lebourgeois, Anne-Elisabeth Laques, M. Dieye, J.

Bourgoin, A. Begue

► To cite this version:

Y.N. Scarpetta, V. Lebourgeois, Anne-Elisabeth Laques, M. Dieye, J. Bourgoin, et al., BFASTm-L2, an unsupervised LULCC detection based on seasonal change detection - an application to large-scale land acquisitions in Senegal. International Journal of Applied Earth Observation and Geoinformation, 2023, 121, 103379 [17 p.]. 10.1016/j.jag.2023.103379. hal-04164102

HAL Id: hal-04164102

https://hal.science/hal-04164102

Submitted on 1 Feb 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# BFASTm-L2, an unsupervised LULCC detection based on seasonal change detection - an application to large-scale land acquisitions in Senegal

Y.N. Scarpetta, V. Lebourgeois, Anne-Elisabeth Laques, M. Dieye, J.

Bourgoin, A. Begue

► To cite this version:

Y.N. Scarpetta, V. Lebourgeois, Anne-Elisabeth Laques, M. Dieye, J. Bourgoin, et al., BFASTm-L2, an unsupervised LULCC detection based on seasonal change detection - an application to large-scale land acquisitions in Senegal. International Journal of Applied Earth Observation and Geoinformation, 2023, 121, 103379 [17 p.]. 10.1016/j.jag.2023.103379. hal-04164102

HAL Id: hal-04164102

https://hal.science/hal-04164102v1

Submitted on 1 Feb 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_1.png" /></figure>

BFASTm-L2, an unsupervised LULCC detection based on seasonal change
detection – An application to large-scale land acquisitions in Senegal

Yasmine Ngadi Scarpetta<sup>a,b,d,*</sup>, Valentine Lebourgeois<sup>a,b</sup>, Anne-Elisabeth Laques<sup>b,d</sup>,
Mohamadou Dieye<sup>c</sup>, Jérémy Bourgoin<sup>a,b</sup>, Agnès Bégué<sup>a,b</sup>

ARTICLE INFO

Keywords:

SITS

Land use land cover change

BFAST

MODIS NDVI

LSLA

A B S T R A C T

In the context of Global Change Research, detection, monitoring and characterization of land use/land cover (LULC) changes are of prime importance. The increasing availability of dense satellite image time series (SITS) has led to a shift in the change detection paradigm, with algorithms able to exploit the full temporal information laid down in SITS. So far, most of these algorithms have focused on the detection of abrupt and gradual changes, and thus developed breakpoint detection based on significant deviations from the mean. However, LULC changes may manifest themselves in other patterns, particularly changes in seasonality (amplitude, number and length of the growing seasons) that are harder to detect. In this paper, we propose a simple method to automatically select the breakpoint linked to the biggest seasonal change in long and dense SITS with multiple breakpoints. This approach - BFASTm-L2 - relies on linking a high-speed algorithm (BFAST monitor) with a time series similarity metric (Euclidian distance L2) sensitive to seasonal changes. The capacity of BFASTm-L2 to identify the date of change in different situations was tested on two data sets, and compared to the performances of three other algorithms (BFAST monitor, BFAST lite, and Edyn). The data sets are 1. a published benchmark data set composed of 25 200 simulated SITS with different change types and change magnitudes, and 2. the 2000-2020 MODIS NDVI SITS over a 200x200 pixels area in Senegal including different study sites which have undergone recent LULC changes due to agricultural large-scale land acquisitions (LSLAs) (as reported in the ground field database used in this study). The results show that BFASTm-L2 is efficient in accurately detecting in time most of the changes, and, in contrast with BFAST Lite and BFASTmonitor, to spatially highlight LSLAs-induced changes without the need of any prior knowledge. The automatic proposed approach, faster than BFAST Lite and Edyn, and with very few tuneable parameters, may thus be easily implemented in unsupervised pipelines to map and analyse generic LULC changes at regional scale.

1. Introduction

In the context of Global Change Research, detection, monitoring, and characterization of land use/land cover (LULC) changes are of prime importance. Global satellite-based Earth observation, with its repetitive coverage at short intervals and consistent image quality, has led to major advances in the field by providing insights into the land dynamics of large and remote areas. Recent studies have benefited from the increasing availability of free remote sensing data and the

overwhelming increase in computing power to identify changes that
occur over time. In particular, the availability of dense Satellite Image
Time Series (SITS) has led to a shift in the change detection paradigm,
with ever-increasing approaches and algorithms exploiting the full
temporal information contained in SITS. Changes detected are usually
categorized as “abrupt”, “gradual” and “seasonal”. While “abrupt” refers
to short-term, large magnitude date-to-date changes (e.g., deforestation,
fire, or urbanization), “gradual” (also referred as trends) refers to long-
term (i.e. inter-annual), small magnitude date-to-date changes (e.g. land

<footnote><sup>a</sup> Cirad, UMR TETIS, Maison de la Télédétection, 500 rue Jean-François Breton, 34093 Montpellier, France</footnote>

<footnote><sup>b</sup> Université de Montpellier, 163 rue Auguste Broussonnet, 34 090 Montpellier, France</footnote>

<footnote><sup>c</sup> ISRA-BAME, Route des Hydrocarbures, BP 3120 Dakar, Senegal</footnote>

<footnote><sup>d</sup> IRD UMR ESPACE-DEV, Maison de la Télédétection, 500 rue Jean-François Breton, 34090 Montpellier, France</footnote>

<footnote>* Corresponding author at: Cirad, UMR TETIS, Maison de la Télédétection, 500 rue Jean-François Breton, 34093 Montpellier, France.
*E-mail addresses:* yasmine.ngadi_scarpetta@cirad.fr (Y. Ngadi Scarpetta), valentine.lebourgeois@cirad.fr (V. Lebourgeois), anne-elisabeth.laques@ird.fr (A.-E. Laques), jeremy.bourgoin@cirad.fr (J. Bourgoin), agnes.begue@cirad.fr (A. Bégué).</footnote>
<figure><img src="image_2.png" /></figure>

BFASTm-L2, an unsupervised LULCC detection based on seasonal change
detection – An application to large-scale land acquisitions in Senegal

Yasmine Ngadi Scarpetta<sup>a,b,d,*</sup>, Valentine Lebourgeois<sup>a,b</sup>, Anne-Elisabeth Laques<sup>b,d</sup>,
Mohamadou Dieye<sup>c</sup>, Jérémy Bourgoin<sup>a,b</sup>, Agnès Bégué<sup>a,b</sup>

ARTICLE INFO

Keywords:

SITS

Land use land cover change

BFAST

MODIS NDVI

LSLA

A B S T R A C T

In the context of Global Change Research, detection, monitoring and characterization of land use/land cover (LULC) changes are of prime importance. The increasing availability of dense satellite image time series (SITS) has led to a shift in the change detection paradigm, with algorithms able to exploit the full temporal information laid down in SITS. So far, most of these algorithms have focused on the detection of abrupt and gradual changes, and thus developed breakpoint detection based on significant deviations from the mean. However, LULC changes may manifest themselves in other patterns, particularly changes in seasonality (amplitude, number and length of the growing seasons) that are harder to detect. In this paper, we propose a simple method to automatically select the breakpoint linked to the biggest seasonal change in long and dense SITS with multiple breakpoints. This approach - BFASTm-L2 - relies on linking a high-speed algorithm (BFAST monitor) with a time series similarity metric (Euclidian distance L2) sensitive to seasonal changes. The capacity of BFASTm-L2 to identify the date of change in different situations was tested on two data sets, and compared to the performances of three other algorithms (BFAST monitor, BFAST lite, and Edyn). The data sets are 1. a published benchmark data set composed of 25 200 simulated SITS with different change types and change magnitudes, and 2. the 2000–2020 MODIS NDVI SITS over a 200x200 pixels area in Senegal including different study sites which have undergone recent LULC changes due to agricultural large-scale land acquisitions (LSLAs) (as reported in the ground field database used in this study). The results show that BFASTm-L2 is efficient in accurately detecting in time most of the changes, and, in contrast with BFAST Lite and BFASTmonitor, to spatially highlight LSLAs-induced changes without the need of any prior knowledge. The automatic proposed approach, faster than BFAST Lite and Edyn, and with very few tuneable parameters, may thus be easily implemented in unsupervised pipelines to map and analyse generic LULC changes at regional scale.

1. Introduction

In the context of Global Change Research, detection, monitoring, and characterization of land use/land cover (LULC) changes are of prime importance. Global satellite-based Earth observation, with its repetitive coverage at short intervals and consistent image quality, has led to major advances in the field by providing insights into the land dynamics of large and remote areas. Recent studies have benefited from the increasing availability of free remote sensing data and the

overwhelming increase in computing power to identify changes that
occur over time. In particular, the availability of dense Satellite Image
Time Series (SITS) has led to a shift in the change detection paradigm,
with ever-increasing approaches and algorithms exploiting the full
temporal information contained in SITS. Changes detected are usually
categorized as “abrupt”, “gradual” and “seasonal”. While “abrupt” refers
to short-term, large magnitude date-to-date changes (e.g., deforestation,
fire, or urbanization), “gradual” (also referred as trends) refers to long-
term (i.e. inter-annual), small magnitude date-to-date changes (e.g. land

<footnote><sup>a</sup> Cirad, UMR TETIS, Maison de la Télédétection, 500 rue Jean-François Breton, 34093 Montpellier, France</footnote>

<footnote><sup>b</sup> Université de Montpellier, 163 rue Auguste Broussonnet, 34 090 Montpellier, France</footnote>

<footnote><sup>c</sup> ISRA-BAME, Route des Hydrocarbures, BP 3120 Dakar, Senegal</footnote>

<footnote><sup>d</sup> IRD UMR ESPACE-DEV, Maison de la Télédétection, 500 rue Jean-François Breton, 34090 Montpellier, France</footnote>

<footnote>* Corresponding author at: Cirad, UMR TETIS, Maison de la Télédétection, 500 rue Jean-François Breton, 34093 Montpellier, France.
*E-mail addresses:* yasmine.ngadi_scarpetta@cirad.fr (Y. Ngadi Scarpetta), valentine.lebourgeois@cirad.fr (V. Lebourgeois), anne-elisabeth.laques@ird.fr (A.-E. Laques), jeremy.bourgoin@cirad.fr (J. Bourgoin), agnes.begue@cirad.fr (A. Bégué).</footnote>
degradation, forest recovery) (Zhu, 2017). “Seasonal” changes are those affecting time series seasonality (i.e., vegetation phenology), and refer more explicitly to changes in the number of growing cycles (i.e., number of “peaks” of vegetation activity), in the season amplitude and in the length of season (i.e., growing season baseline length).

So far, change detection algorithms have mostly been used to detect significant deviations from the mean in forest ecosystems (Ochtyra et al., 2020). This is explained by the fact that discrimination between phenological changes driven by climatic variability and disturbances is easier in stable seasonal environments such as evergreen temperate or tropical forests than, for example, in drier ecosystems (Browning et al., 2017; Gao et al., 2021; Zhu and Woodcock, 2014). Less focus has been given to other land surface dynamics such as changes in land use (Verburg et al., 2009). Land use, which refers to the purposes for which human exploit the land cover, is hardly inferred directly from remote sensing images and very often needs ground-knowledge to be accurately assessed. For this task, temporal series of Vegetation Indices (such as the Normalised Difference Vegetation Index or NDVI) are often used. While remote assessment of a single land use may be difficult to achieve, land use changes may however be assessed by detecting persistent seasonal changes (e.g. changes in the amplitude, length of season and/or number of seasons) within VI's time series, as homogeneous land practices are expected to present typical intra-annual patterns (Setiawan and Yoshino, 2014).

This is how, Hentze et al. (2017), focusing on the detection of specific seasonal changes (from unimodal to bimodal distributions and vice-versa) in the NDVI, jointly with a seasonal-trend analysis, were able to identify agricultural land tenure transitions (from large-scale to small-scale and vice-versa) in Zimbabwe using BFAST (Breaks For Additive Season and Trend; (Verbesselt et al., 2010a)). However, their method has the drawback of not being up-scalable because of the algorithm used (BFAST is computationally expensive), and because their need of external data (crop mask).

Speed, automaticity and accuracy of the detection of seasonal changes are thus crucial for the development of *generic* approaches to land use change detection. This is particularly relevant when using long and dense temporal series, in which a mix of different types of changes with different intensities may occur. As a consequence, approaches that perform temporal segmentation (such as DBEST (Jamali et al. (2015)) or LandTrendr (Kennedy et al. (2010))) which aim to find major abrupt and gradual changes without considering the seasonality are out of the scope of this study. The same applies to the deep learning methods, despite the growing interest in the remote sensing community due to their ability to automatically extract spectral-spatial features from satellite imagery (Tuia et al., 2021; Yuan et al., 2020; Zhu et al., 2017). Because 1- the unsupervised/ semi-supervised deep learning approaches for LULC change detection are still at early stage of development (Leenstra et al.; Meshkini et al., 2021), in addition to the 2- poor interpretability of the models, and the 3- low generalization performance linked to the highly contextual-dependent methods, the labeled-data scarcity (hampering the method's scalability at regional scale), and the huge spectral-spatial variability of the targeted object (land use changes), these approaches are not considered in this study. They however should be the object of further research.

To continue with the change detection characteristics sought in this study, to be applicable on large scale, change detection approaches need to be easy to use (i.e number of parameters), and preferably integrated on cloud-computing platforms. Amongst the well-established statistical-based algorithms implemented on big data platforms, BFASTmonitor (Verbesselt et al., 2012) stands out from the rest because of its simple configuration, speed and massively-parallel GPU implementation (Gieseke et al., 2020), and its relative good performance in detecting seasonal changes (Awty-Carroll, 2019). Its major drawback lays in its high commission error (false positives). Minimization of false positives in long and dense time series is often performed using a threshold on the breakpoint magnitude (Gao et al., 2021), or by using some statistical

tests such as the Chow test (Bullock et al., 2020). When more than one true breakpoint is found, selection of a unique breakpoint for mapping purpose may be difficult, unless when looking for specific changes (in magnitude, sign or pattern) in a specific period of time.

Selection of a unique breakpoint may be more difficult when using algorithms that do not perform a season-trend decomposition (e.g. BFAST monitor, EWMACD). These latter, usually faster (and thus suitable for large-scale applications), are often able to detect any type of change, including the seasonal ones, within the same non-decomposed time series. However, because the computed breakpoint's magnitudes are linked to deviations from the model, selection of a unique breakpoint based on its magnitude may disfavour the selection of seasonal changes, particularly those modifying the time series “shape” (i.e. number of growing cycles per year, length of season) without heavily impacting its amplitude.

In this study, we propose a simple, fast, generic and unsupervised approach (hereafter referred to as BFASTm-L2) to select, in long and dense NDVI time series with multiple breakpoints, the optimal breakpoint linked to the most important land use change (i.e., linked to the most important seasonal or “pattern” change within the time series (in amplitude, length of season, NOS, or a mix of them)).

Two sets of research questions aimed to be answered in this paper: 1- The first one is broad and concerns the temporal accuracy (does the breakpoint correspond to the year of change?) and the sensitivity (how likely are different *type* of changes detected?) to different *types* and *intensities* of change of BFASTm-L2 and three change detection algorithms tested for comparison purposes: BFASTmonitor, BFAST Lite, Edyn. 2- The second set of questions concerns an application case, which is the detection (are the different change detection algorithms able to detect in long and dense time series, with likely multiple changes of different types and intensities, the land use changes related to LSLAs?) and mapping (which types of changes are most likely to be highlighted in maps using the breakpoint magnitude as mapping variable? Can LSLAs be pinpointed in such maps?) of Large-Scale Land Acquisitions (LSLAs) in Senegal.

More details on the material and methods used are given in the next section.

## 2. Material and methods

### 2.1. Global approach

In this study BFASTm-L2 is proposed as a method to select, in long and dense NDVI time series with multiple breakpoints, the optimal breakpoint linked to the most important land use change, that is the one related to the most important “pattern” change within time series. Because the breakpoint selection is based on the breakpoint magnitude, we tested a magnitude metric based on the Euclidean distance (L2), a time series similarity metric proposed by Lhermitte et al. (2011), which is prone to be more sensitive to seasonal changes than to abrupt and gradual changes. BFASTm-L2 relies thus on the sequential running of BFASTmonitor for the fast detection of breakpoints, jointly with the Euclidean distance (L2) for the breakpoint selection.

Because long and dense time series most of the time include multiple changes of different types (often combined) and intensities, a benchmark dataset of simulated seasonal time series including a unique change was used to answer the first research set of questions, aiming to assess the temporal accuracy and sensitivity of the tested change detection algorithms to different *types* and *intensities* of change. More specifically, a sub-sample of the benchmark dataset provided by Awty-Carroll et al. (2019) composed of 25 000 simulated SITS was used, including different single change types and multiple noise/change intensity levels (see Table 1). To answer the second set of questions, related to the performance of the change detection algorithms in detecting and mapping LSLAs in Senegal, the MODIS 16-day NDVI time series (MOD13Q1 v.6) acquired over the 2000–2021 period in different
degradation, forest recovery) (Zhu, 2017). “Seasonal” changes are those
affecting time series seasonality (i.e., vegetation phenology), and refer
more explicitly to changes in the number of growing cycles (i.e., number
of “peaks” of vegetation activity), in the season amplitude and in the
length of season (i.e., growing season baseline length).

So far, change detection algorithms have mostly been used to detect
significant deviations from the mean in forest ecosystems (Ochtyra et al.,
2020). This is explained by the fact that discrimination between
phenological changes driven by climatic variability and disturbances is
easier in stable seasonal environments such as evergreen temperate or
tropical forests than, for example, in drier ecosystems (Browning et al.,
2017; Gao et al., 2021; Zhu and Woodcock, 2014). Less focus has been
given to other land surface dynamics such as changes in land use (Ver-
burg et al., 2009). Land use, which refers to the purposes for which
human exploit the land cover, is hardly inferred directly from remote
sensing images and very often needs ground-knowledge to be accurately
assessed. For this task, temporal series of Vegetation Indices (such as the
Normalised Difference Vegetation Index or NDVI) are often used. While
remote assessment of a single land use may be difficult to achieve, land
use changes may however be assessed by detecting persistent seasonal
changes (e.g. changes in the amplitude, length of season and/or number
of seasons) within VI's time series, as homogeneous land practices are
expected to present typical intra-annual patterns (Setiawan and Yosh-
ino, 2014).

This is how, Hentze et al. (2017), focusing on the detection of specific
seasonal changes (from unimodal to bimodal distributions and vice-
versa) in the NDVI, jointly with a seasonal-trend analysis, were able to
identify agricultural land tenure transitions (from large-scale to small-
scale and vice-versa) in Zimbabwe using BFAST (Breaks For Additive
Season and Trend; (Verbesselt et al., 2010a)). However, their method
has the drawback of not being up-scalable because of the algorithm used
(BFAST is computationally expensive), and because their need of
external data (crop mask).

Speed, automaticity and accuracy of the detection of seasonal
changes are thus crucial for the development of *generic* approaches to
land use change detection. This is particularly relevant when using long
and dense temporal series, in which a mix of different types of changes
with different intensities may occur. As a consequence, approaches that
perform temporal segmentation (such as DBEST (Jamali et al. (2015)) or
LandTrendr (Kennedy et al. (2010))) which aim to find major abrupt and
gradual changes without considering the seasonality are out of the scope
of this study. The same applies to the deep learning methods, despite the
growing interest in the remote sensing community due to their ability to
automatically extract spectral-spatial features from satellite imagery
(Tuia et al., 2021; Yuan et al., 2020; Zhu et al., 2017). Because 1- the
unsupervised/ semi-supervised deep learning approaches for LULC
change detection are still at early stage of development (Leenstra et al.;
Meshkini et al., 2021), in addition to the 2- poor interpretability of the
models, and the 3- low generalization performance linked to the highly
contextual-dependent methods, the labeled-data scarcity (hampering
the method's scalability at regional scale), and the huge spectral-spatial
variability of the targeted object (land use changes), these approaches
are not considered in this study. They however should be the object of
further research.

To continue with the change detection characteristics sought in this
study, to be applicable on large scale, change detection approaches need
to be easy to use (i.e number of parameters), and preferably integrated
on cloud-computing platforms. Amongst the well-established statistical-
based algorithms implemented on big data platforms, BFASTmonitor
(Verbesselt et al., 2012) stands out from the rest because of its simple
configuration, speed and massively-parallel GPU implementation (Gie-
seke et al., 2020), and its relative good performance in detecting sea-
sonal changes (Awty-Carroll, 2019). Its major drawback lays in its high
commission error (false positives). Minimization of false positives in
long and dense time series is often performed using a threshold on the
breakpoint magnitude (Gao et al., 2021), or by using some statistical

tests such as the Chow test (Bullock et al., 2020). When more than one
true breakpoint is found, selection of a unique breakpoint for mapping
purpose may be difficult, unless when looking for specific changes (in
magnitude, sign or pattern) in a specific period of time.

Selection of a unique breakpoint may be more difficult when using
algorithms that do not perform a season-trend decomposition (e.g.
BFAST monitor, EWMACD). These latter, usually faster (and thus suit-
able for large-scale applications), are often able to detect any type of
change, including the seasonal ones, within the same non-decomposed
time series. However, because the computed breakpoint's magnitudes
are linked to deviations from the model, selection of a unique breakpoint
based on its magnitude may disfavour the selection of seasonal changes,
particularly those modifying the time series "shape" (i.e. number of
growing cycles per year, length of season) without heavily impacting its
amplitude.

In this study, we propose a simple, fast, generic and unsupervised
approach (hereafter referred to as BFASTm-L2) to select, in long and
dense NDVI time series with multiple breakpoints, the optimal break-
point linked to the most important land use change (i.e., linked to the
most important seasonal or "pattern" change within the time series (in
amplitude, length of season, NOS, or a mix of them)).

Two sets of research questions aimed to be answered in this paper: 1-
The first one is broad and concerns the temporal accuracy (does the
breakpoint correspond to the year of change?) and the sensitivity (how
likely are different *type* of changes detected?) to different *types* and *in-
tensities* of change of BFASTm-L2 and three change detection algorithms
tested for comparison purposes: BFASTmonitor, BFAST Lite, Edyn. 2-
The second set of questions concerns an application case, which is the
detection (are the different change detection algorithms able to detect in
long and dense time series, with likely multiple changes of different
types and intensities, the land use changes related to LSLAs?) and
mapping (which types of changes are most likely to be highlighted in
maps using the breakpoint magnitude as mapping variable? Can LSLAs
be pinpointed in such maps?) of Large-Scale Land Acquisitions (LSLAs)
in Senegal.

More details on the material and methods used are given in the next
section.

2. Material and methods

2.1. *Global approach*

In this study BFASTm-L2 is proposed as a method to select, in long and dense NDVI time series with multiple breakpoints, the optimal breakpoint linked to the most important land use change, that is the one related to the most important “pattern” change within time series. Because the breakpoint selection is based on the breakpoint magnitude, we tested a magnitude metric based on the Euclidean distance (L2), a time series similarity metric proposed by Lhermitte et al. (2011), which is prone to be more sensitive to seasonal changes than to abrupt and gradual changes. BFASTm-L2 relies thus on the sequential running of BFASTmonitor for the fast detection of breakpoints, jointly with the Euclidean distance (L2) for the breakpoint selection.

Because long and dense time series most of the time include multiple
changes of different types (often combined) and intensities, a bench-
mark dataset of simulated seasonal time series including a unique
change was used to answer the first research set of questions, aiming to
assess the temporal accuracy and sensitivity of the tested change
detection algorithms to different *types* and *intensities* of change. More
specifically, a sub-sample of the benchmark dataset provided by Awty-
Carroll et al. (2019) composed of 25 000 simulated SITS was used,
including different single change types and multiple noise/change in-
tensity levels (see Table 1). To answer the second set of questions,
related to the performance of the change detection algorithms in
detecting and mapping LSLAs in Senegal, the MODIS 16-day NDVI time
series (MOD13Q1 v.6) acquired over the 2000–2021 period in different
<caption>**Table 1**</caption>

Types of change, levels of intensities and number of samples present in the simulated data set (Awty-Carroll et al., 2019). For the break/trend set, each abrupt change in NDVI is followed by either no trend or one of the six levels of trend present in the trend only set.

<table><thead><tr><th>Type of changes</th><th>Levels (units)</th><th>Number of simulations</th></tr></thead><tbody><tr><td>- No change</td><td>-</td><td>400</td></tr><tr><td>Trend Trend only</td><td>±[0.001, 0.0015, 0.002] (NDVI / year)</td><td>2 400</td></tr><tr><td>Break (abrupt) / trend</td><td>±[0.1, 0.2, 0.3] (NDVI)</td><td>16 800</td></tr><tr><td>Amplitude</td><td>±[0.1, 0.2, 0.3] (NDVI)</td><td>2 400</td></tr><tr><td>Seasonal LOS (length of season)</td><td>-[13, 22, 30, 37, 43, 49] (days)</td><td>2 400</td></tr><tr><td>NOS (number of season)</td><td>1 to 2, 2 to 1</td><td>800</td></tr><tr><td>Total</td><td>-</td><td>25 200</td></tr></tbody></table>

study areas with known LSLAs were used. Results were validated with ground observations. The flowchart of the data and methods is illustrated here below in Fig. 1.

## 2.2. Data

### 2.2.1. Benchmark dataset and preprocessing

This study made use of the simulated NDVI time series (10-year 2006–2015 at 16-day temporal resolution) data set created by Awty-Carroll (2019) and available at https://osf.io/taf9y/. The 16-day temporal frequency makes this dataset suitable for the purposes of this study which uses the MODIS satellite imagery. In this study, only the gap-free simulations were used consisting of 25 200 seasonal time series generated to represent a large range of ecosystem dynamics based on the previous work of Verbesselt et al. (2010b). This dataset includes a large range of unique trend, abrupt (positive and negative), seasonal (amplitude, length of season (LOS) and number of seasons (NOS)) changes, with eight levels of noise (random value from a normal distribution with a mean of 0 and a standard deviation ∈ [0: 0.01: 0.07]), and 50 replicates (Table 1). The LOS changes are changes that move back the start of the season from 13 to 49 days. All changes were placed in January 2011.

Because of the presence of noise hindering the detection of change, in this study the simulated time series were, as for the real MODIS time series smoothed using the Savitzky-Golay smoothing filter. A moving window length of 9 observations and a polynomial order of 3 were used.

### 2.2.2. Study case dataset

#### 2.2.2.1. Study area and ground dataset.
Senegal is in the west part of the Sudano-Sahelian zone and is characterized by an overall low average annual rainfall, but with high inter- and intra-annual variability that constrains the vegetation growth. The precipitation shows an increasing gradient along the North-South direction. A distinct seasonality is present, with a long dry season and a short rainy season spanning from late June to early October (Abel et al., 2019).

As many African countries concerned by Large Scale Land Acquisitions (LSLAs) (https://landmatrix.org), Senegal had in 2016 around 3% of its total arable land (270 908 ha) declared under contract by foreign investors from 12 countries (Harding et al., 2016). However, as >50% of those deals have an area under contract smaller than 5 000 ha, the country is less affected by the so called “megadeals” (>50 000 ha, only 2 out of 19 in 2016). Because of those climatic and LSLA characteristics, Senegal represents an interesting and difficult study case for the detection of land use system changes such as those induced by LSLAs.

The Senegalese Institute of Agricultural Research (ISRA) conducted in 2019 an extensive field campaign on LSLA. >700 records and corresponding attributes were initially recorded in a database (M. Dieye, personal communication, 2022). Attributes consist in all kinds of information related to the identified deals, such as: name, deal type, coordinates, negotiation status, implementation status, year of transaction, size, previous land use, previous land tenure etc. From this database, and using observations from Google Earth/Sentinel Hub images, a selection was done to keep only active agricultural deals with growing crops (some deals are abandoned, others are still in negotiation), and with an implementation date after 2003 (year as from which changes may be detected from MODIS SITS, considering the training period length needed by the detection algorithms).

From this ground-field database subsample, four individual (points) LSLA study cases with different land processes under different ecorigions (see Fig. 2) were selected for algorithm testing purposes. Point 1 is an example of a conversion from small agriculture to LSLA, points 2 to 4 are examples of conversions from natural vegetation to LSLA.

More specifically, study case 1 (16.1118°N, 15.9954°W) and 2 (16.1159°N, 16.0230°W) are in the Senegal river valley (rainfall: 150–600 mm; (Tappan et al., 2004)) and are within pivot irrigation areas belonging to a concession growing vegetables. Changes occurred from the end of 2011 for point 1, and the start of 2016 for point 2. Point 3 (14.6444°N, 17.0271°W), in the “Agricultural Expansion Region” as defined in Tappan et al. (2004) (rainfall: 600–700 mm), is within one of the production blocks belonging to a concession specialized in vegetable production, and changes are observed as from 2008. Finally, the southern point 4 (13.0104°N, 14.1766°W), in the Anambé basin (Casamance; rainfall: 800–1400 mm), is in the area of a concession specialized in irrigated rice. Field preparation is observed as from the start of 2007.

In addition to the LSLA study cases, an area of 200x200 MODIS pixels close to Dakar (red box in Fig. 2), including six production blocks of different concessions (see red polygons in Fig. 7), was selected for mapping purposes (see section 2.4 for more details). The area is spatially contrasted, with urban/rural areas and natural vegetation (see Fig. 7).

#### 2.2.2.2. MODIS NDVI data and pre-processing.
With its global coverage, moderate spatial resolution (250 m) and high temporal resolution (1 to 2 days), the Moderate Resolution Imaging Spectroradiometer (MODIS) sensor allows for the detection of subtle changes in land cover. A set of MODIS NDVI 16-day composites at 250 m resolution (MOD13Q1, collection 6) acquired over Senegal for the period 2000–2021 was pre-processed in Google Earth Engine. The 16-day composite NDVI product was chosen to reduce NDVI variability due to meteorological distortions like clouds. Series were gap-filled with linear interpolations and smoothed using an optimized weighted Savitzky-Golay filter (Chen et al., 2004) in order to reduce the noise. Weights were computed following the approach developed in Piou et al. (2013), and are function of the pixel's reliability (i.e. quality flag, view zenith angle), and position in the moving window (exponentially decreasing weights with the distance to the window's centre). After some testing, a moving window length of 13 observations, and a polynomial order of 3 (in order to keep the ratio w/p close to 3–4) were used.

## 2.3. Methods

### 2.3.1. Statistical-based change detection algorithms

Up to now, the LULC change community has benefited from an ever-increasing emergence of change detection algorithms exploiting the full temporal information contained in SITS (Molinier et al., 2021). Some of them are now well-established and implemented on big data platforms enabling fast processing of large volume data.

Change detection algorithms mainly differ in their approach to process time series and detect changes in them. In their approach to
<caption>**Table 1**</caption>

Types of change, levels of intensities and number of samples present in the simulated data set (Awty-Carroll et al., 2019). For the break/trend set, each abrupt change in NDVI is followed by either no trend or one of the six levels of trend present in the trend only set.

<table><thead><tr><th>Type of changes</th><th>Levels (units)</th><th>Number of simulations</th></tr></thead><tbody><tr><td>- No change</td><td>-</td><td>400</td></tr><tr><td>Trend Trend only</td><td>±[0.001, 0.0015, 0.002] (NDVI / year)</td><td>2 400</td></tr><tr><td>Break (abrupt) / trend</td><td>±[0.1, 0.2, 0.3] (NDVI)</td><td>16 800</td></tr><tr><td>Amplitude</td><td>±[0.1, 0.2, 0.3] (NDVI)</td><td>2 400</td></tr><tr><td>Seasonal LOS (length of season)</td><td>-[13, 22, 30, 37, 43, 49] (days)</td><td>2 400</td></tr><tr><td>NOS (number of season)</td><td>1 to 2, 2 to 1</td><td>800</td></tr><tr><td>Total</td><td>-</td><td>25 200</td></tr></tbody></table>

study areas with known LSLAs were used. Results were validated with ground observations. The flowchart of the data and methods is illustrated here below in Fig. 1.

## 2.2. Data

### 2.2.1. Benchmark dataset and preprocessing

This study made use of the simulated NDVI time series (10-year 2006–2015 at 16-day temporal resolution) data set created by Awty-Carroll (2019) and available at https://osf.io/taf9y/. The 16-day temporal frequency makes this dataset suitable for the purposes of this study which uses the MODIS satellite imagery. In this study, only the gap-free simulations were used consisting of 25 200 seasonal time series generated to represent a large range of ecosystem dynamics based on the previous work of Verbesselt et al. (2010b). This dataset includes a large range of unique trend, abrupt (positive and negative), seasonal (amplitude, length of season (LOS) and number of seasons (NOS)) changes, with eight levels of noise (random value from a normal distribution with a mean of 0 and a standard deviation ∈ [0: 0.01: 0.07]), and 50 replicates (Table 1). The LOS changes are changes that move back the start of the season from 13 to 49 days. All changes were placed in January 2011.

Because of the presence of noise hindering the detection of change, in this study the simulated time series were, as for the real MODIS time series smoothed using the Savitzky-Golay smoothing filter. A moving window length of 9 observations and a polynomial order of 3 were used.

### 2.2.2. Study case dataset

#### 2.2.2.1. Study area and ground dataset.
Senegal is in the west part of the Sudano-Sahelian zone and is characterized by an overall low average annual rainfall, but with high inter- and intra-annual variability that constrains the vegetation growth. The precipitation shows an increasing gradient along the North-South direction. A distinct seasonality is present, with a long dry season and a short rainy season spanning from late June to early October (Abel et al., 2019).

As many African countries concerned by Large Scale Land Acquisitions (LSLAs) (https://landmatrix.org), Senegal had in 2016 around 3% of its total arable land (270 908 ha) declared under contract by foreign investors from 12 countries (Harding et al., 2016). However, as >50% of those deals have an area under contract smaller than 5 000 ha, the country is less affected by the so called “megadeals” (>50 000 ha, only 2 out of 19 in 2016). Because of those climatic and LSLA characteristics, Senegal represents an interesting and difficult study case for the detection of land use system changes such as those induced by LSLAs.

The Senegalese Institute of Agricultural Research (ISRA) conducted in 2019 an extensive field campaign on LSLA. >700 records and corresponding attributes were initially recorded in a database (M. Dieye, personal communication, 2022). Attributes consist in all kinds of information related to the identified deals, such as: name, deal type, coordinates, negotiation status, implementation status, year of transaction, size, previous land use, previous land tenure etc. From this database, and using observations from Google Earth/Sentinel Hub images, a selection was done to keep only active agricultural deals with growing crops (some deals are abandoned, others are still in negotiation), and with an implementation date after 2003 (year as from which changes may be detected from MODIS SITS, considering the training period length needed by the detection algorithms).

From this ground-field database subsample, four individual (points) LSLA study cases with different land processes under different ecorigions (see Fig. 2) were selected for algorithm testing purposes. Point 1 is an example of a conversion from small agriculture to LSLA, points 2 to 4 are examples of conversions from natural vegetation to LSLA.

More specifically, study case 1 (16.1118°N, 15.9954°W) and 2 (16.1159°N, 16.0230°W) are in the Senegal river valley (rainfall: 150–600 mm; (Tappan et al., 2004)) and are within pivot irrigation areas belonging to a concession growing vegetables. Changes occurred from the end of 2011 for point 1, and the start of 2016 for point 2. Point 3 (14.6444°N, 17.0271°W), in the “Agricultural Expansion Region” as defined in Tappan et al. (2004) (rainfall: 600–700 mm), is within one of the production blocks belonging to a concession specialized in vegetable production, and changes are observed as from 2008. Finally, the southern point 4 (13.0104°N, 14.1766°W), in the Anambé basin (Casamance; rainfall: 800–1400 mm), is in the area of a concession specialized in irrigated rice. Field preparation is observed as from the start of 2007.

In addition to the LSLA study cases, an area of 200x200 MODIS pixels close to Dakar (red box in Fig. 2), including six production blocks of different concessions (see red polygons in Fig. 7), was selected for mapping purposes (see section 2.4 for more details). The area is spatially contrasted, with urban/rural areas and natural vegetation (see Fig. 7).

#### 2.2.2.2. MODIS NDVI data and pre-processing.
With its global coverage, moderate spatial resolution (250 m) and high temporal resolution (1 to 2 days), the Moderate Resolution Imaging Spectroradiometer (MODIS) sensor allows for the detection of subtle changes in land cover. A set of MODIS NDVI 16-day composites at 250 m resolution (MOD13Q1, collection 6) acquired over Senegal for the period 2000–2021 was pre-processed in Google Earth Engine. The 16-day composite NDVI product was chosen to reduce NDVI variability due to meteorological distortions like clouds. Series were gap-filled with linear interpolations and smoothed using an optimized weighted Savitzky-Golay filter (Chen et al., 2004) in order to reduce the noise. Weights were computed following the approach developed in Piou et al. (2013), and are function of the pixel’s reliability (i.e. quality flag, view zenith angle), and position in the moving window (exponentially decreasing weights with the distance to the window’s centre). After some testing, a moving window length of 13 observations, and a polynomial order of 3 (in order to keep the ratio w/p close to 3–4) were used.

## 2.3. Methods

### 2.3.1. Statistical-based change detection algorithms

Up to now, the LULC change community has benefited from an ever-increasing emergence of change detection algorithms exploiting the full temporal information contained in SITS (Molinier et al., 2021). Some of them are now well-established and implemented on big data platforms enabling fast processing of large volume data.

Change detection algorithms mainly differ in their approach to process time series and detect changes in them. In their approach to
<figure><img src="image_2.png" /><figcaption>Fig. 1. Flowchart of the data and methods.</figcaption></figure>

process time-series, two main groups of algorithms exist: While offline algorithms operate retrospectively on the complete time series, online (or real-time) algorithms aim to detect changes as soon as they occur (Bullock et al., 2020). Because offline methods make use of the entire time series, there are often more robust than online methods. In turn, online methods are often faster than offline methods as they only use some of the data preceding the real-time observation (i.e., training period). They are however prone to false positives and require a stable training period (limiting its use in places frequently disturbed (e.g. agriculture (Zhu et al., 2020)). As the length of the training period may have an impact on the quality of the fit - risk of overfitting when training periods are too short vs risk to include breakpoints in the training period when those are too long (Brooks et al., 2017)), algorithms that are able to automatically select the optimal training period length (such as Edyn (Brooks et al., 2017), BFAST monitor (Verbesselt et al., 2012) or CCDC (Zhu and Woodcock, 2014)) are preferred.

In their approach to detect changes, while some algorithms provide

the option of selecting an optimal model with one unique breakpoint (e.g. BFAST Lite (Masiliunas et al., 2021)), almost all of them detect multiple breakpoints. Changes may be detected in the trend and seasonal component separately (e.g. BFAST (Verbesselt et al., 2010a), BEAST (Zhao et al., 2019)) or in the undecomposed time series (e.g. BFAST Lite, BFAST monitor, EWMACD (Brooks et al., 2014)), using temporal segmentation approaches based on residual-errors and angle criterion (e.g. DBEST (Jamali et al. (2015)), LandTrendr (Kennedy et al. (2010)), or model-deviation seeking approaches (e.g. BFAST monitor, CCDC (Zhu and Woodcock (2014)), EWMACD). Decomposition may be interesting for end-users with an a priori knowledge of the type of change foreseen. As an example, Mardian et al. (2021) assuming that pasture/rangeland conversions to cropland mostly impact VIs' seasonal component, applied a modified version of BFAST that constrained the model to detect only seasonal changes and obtained higher change detection accuracy compared to BFAST and BEAST. However, it is worth noting that algorithms that perform decomposition are usually slower than those that do
<figure><img src="image_3.png" /><figcaption>Fig. 1. Flowchart of the data and methods.</figcaption></figure>

process time-series, two main groups of algorithms exist: While offline algorithms operate retrospectively on the complete time series, online (or real-time) algorithms aim to detect changes as soon as they occur (Bullock et al., 2020). Because offline methods make use of the entire time series, there are often more robust than online methods. In turn, online methods are often faster than offline methods as they only use some of the data preceding the real-time observation (i.e., training period). They are however prone to false positives and require a stable training period (limiting its use in places frequently disturbed (e.g. agriculture (Zhu et al., 2020)). As the length of the training period may have an impact on the quality of the fit - risk of overfitting when training periods are too short vs risk to include breakpoints in the training period when those are too long (Brooks et al., 2017)), algorithms that are able to automatically select the optimal training period length (such as Edyn (Brooks et al., 2017), BFAST monitor (Verbesselt et al., 2012) or CCDC (Zhu and Woodcock, 2014)) are preferred.

In their approach to detect changes, while some algorithms provide

the option of selecting an optimal model with one unique breakpoint (e.g. BFAST Lite (Masiliunas et al., 2021)), almost all of them detect multiple breakpoints. Changes may be detected in the trend and seasonal component separately (e.g. BFAST (Verbesselt et al., 2010a), BEAST (Zhao et al., 2019)) or in the undecomposed time series (e.g. BFAST Lite, BFAST monitor, EWMACD (Brooks et al., 2014)), using temporal segmentation approaches based on residual-errors and angle criterion (e.g. DBEST (Jamali et al. (2015)), LandTrendr (Kennedy et al. (2010)), or model-deviation seeking approaches (e.g. BFAST monitor, CCDC (Zhu and Woodcock (2014)), EWMACD). Decomposition may be interesting for end-users with an a priori knowledge of the type of change foreseen. As an example, Mardian et al. (2021) assuming that pasture/rangeland conversions to cropland mostly impact VIs' seasonal component, applied a modified version of BFAST that constrained the model to detect only seasonal changes and obtained higher change detection accuracy compared to BFAST and BEAST. However, it is worth noting that algorithms that perform decomposition are usually slower than those that do
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Study area (red box) and study cases (red dots) in Senegal located within agro-industrial concessions (red polygons) (source: ISRA-BAME field data base). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

not. In addition, errors in the decomposition (because of using an inappropriate model) may be translated in errors in the accuracy of the changes detected in both SITS components (Mardian et al., 2021; Zhao et al., 2019). To avoid selection of a single-model algorithm, focus is recently given to ensemble learning algorithms such as BEAST (Zhao et al., 2019), that are however computational cost expensive. As one can see, selection of a single-model algorithm highly depends on the application scale, targeted type of change, environment, SITS source and frequency.

The approach proposed in this study for the selection of a single breakpoint linked to the biggest seasonal change in long and dense time series, is based on existing algorithms. For the selection of the change detection algorithm, focus was first given on speed, followed by the algorithm's sensitivity to seasonal changes and its ease of use (number of tuneable parameters). On these criteria, BFASTmonitor was chosen as the base algorithm for the change detection approach. For comparison purposes, in addition to BFASTmonitor, two other algorithms were selected: BFAST Lite, because of its robustness and speed (Masiliūnas et al., 2021), and Edyn, because of its speed and ability to capture seasonal changes (Awty-Carroll et al., 2019).

More details on the different change detection algorithms used and

developed in this study are given in the next subsections and summarized in the table below (Table 2).

2.3.1.1. **BFASTmonitor (BFASTm)**. The BFASTmonitor (Breaks For Additive Season and Trend Monitor; (Verbesselt et al., 2012)) is an online (near real-time) unsupervised change detection algorithm that flags abnormal observations within a monitoring period, based on a stable history period. More specifically, once the start of the monitoring period has been defined, a stable history period is automatically selected using the reversed-ordered-cumulative sum (CUSUM) of residuals (default approach). Then, a regression model (here a linear harmonic regression model) is fitted based on the history period. Finally, the moving sums (MOSUM) of residuals are used (bandwidth defined by the h parameter) in the monitoring period to determine whether the model remains stable for new observations. A break is detected when the absolute value of the moving sums exceeds a probability threshold. The magnitude of change recorded represents the median of the difference between the data and the model prediction in the monitoring period. Because BFASTmonitor only needs a single observation to exceed the threshold, the algorithm is prone to false positives (Awty-Carroll, 2019; Ghaderpour and Vujadinović, 2020) and magnitude thresholds are often applied to minimize
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Study area (red box) and study cases (red dots) in Senegal located within agro-industrial concessions (red polygons) (source: ISRA-BAME field data base). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

not. In addition, errors in the decomposition (because of using an inappropriate model) may be translated in errors in the accuracy of the changes detected in both SITS components (Mardian et al., 2021; Zhao et al., 2019). To avoid selection of a single-model algorithm, focus is recently given to ensemble learning algorithms such as BEAST (Zhao et al., 2019), that are however computational cost expensive. As one can see, selection of a single-model algorithm highly depends on the application scale, targeted type of change, environment, SITS source and frequency.

The approach proposed in this study for the selection of a single breakpoint linked to the biggest seasonal change in long and dense time series, is based on existing algorithms. For the selection of the change detection algorithm, focus was first given on speed, followed by the algorithm's sensitivity to seasonal changes and its ease of use (number of tuneable parameters). On these criteria, BFASTmonitor was chosen as the base algorithm for the change detection approach. For comparison purposes, in addition to BFASTmonitor, two other algorithms were selected: BFAST Lite, because of its robustness and speed (Masiliūnas et al., 2021), and Edyn, because of its speed and ability to capture seasonal changes (Awty-Carroll et al., 2019).

More details on the different change detection algorithms used and

developed in this study are given in the next subsections and summarized in the table below (Table 2).

2.3.1.1. **BFASTmonitor (BFASTm)**. The BFASTmonitor (Breaks For Adaptive Season and Trend Monitor; (Verbesselt et al., 2012)) is an online (near real-time) unsupervised change detection algorithm that flags abnormal observations within a monitoring period, based on a stable history period. More specifically, once the start of the monitoring period has been defined, a stable history period is automatically selected using the reversed-ordered-cumulative sum (CUSUM) of residuals (default approach). Then, a regression model (here a linear harmonic regression model) is fitted based on the history period. Finally, the moving sums (MOSUM) of residuals are used (bandwidth defined by the h parameter) in the monitoring period to determine whether the model remains stable for new observations. A break is detected when the absolute value of the moving sums exceeds a probability threshold. The magnitude of change recorded represents the median of the difference between the data and the model prediction in the monitoring period. Because BFASTmonitor only needs a single observation to exceed the threshold, the algorithm is prone to false positives (Awty-Carroll, 2019; Ghaderpour and Vujadinović, 2020) and magnitude thresholds are often applied to minimize
<caption>**Table 2**</caption>

Main characteristics of the change-detection algorithms used in this study, their strengths and limitations.

<table><thead><tr><th></th><th>Algorithm</th><th>Speed</th><th>Parameter number *</th><th>Model</th><th>Decomposition</th><th>Strengths</th><th>Limitations (others than speed and number of parameters)</th></tr></thead><tbody><tr><td>Offline</td><td>BFAST Lite (Masiliūnas et al., 2021)</td><td>++</td><td>+</td><td>Model-deviation</td><td>None</td><td>Detection of abrupt and trend changes. Possibility to select a model with a unique (the highest-magnitude) breakpoint.</td><td>Low performance in capturing seasonal changes</td></tr><tr><td>Online</td><td>BFAST monitor (Verbesselt et al., 2012)<br>EDYN (Brooks et al., 2017)</td><td>+++</td><td>+</td><td>Model-deviation</td><td>None</td><td>Sensitivity to seasonal changes (Awty-Carroll, 2019). Detection of subtle (sub-pixel) changes in long time series (model dynamically retrained); Able to capture seasonal changes</td><td>High false positive rate (Masiliūnas et al., 2021)<br>High sensitivity to algorithm parameters (Saxena et al., 2018)</td></tr></tbody></table>

* Number of key parameters (others than (if applied) the order of the harmonic term, the probability threshold/ statistical significance level and the training period length):

$$ \leq 2: +, \geq 3: - $$

them (Gao et al., 2021; Hamunyela, 2017). However, its massively-parallel GPU implementation (Gieseke et al., 2020), that makes it 14.5 times faster than the newly launched BFAST Lite (Masiliūnas et al., 2021), its implementation on Google Earth Engine (Hamunyela et al., 2020) and its good performance in detecting seasonal changes (Awty-Carroll, 2019), make of BFASTmonitor a potential change detection algorithm for the purposes of this study. This study made use of the Python package bfast0.7, available at https://pypi.org/project/bfast/, that provides a parallel implementation of the BFASTmonitor algorithm. The parameters used were: h = 0.25 (the MOSUM bandwidth), k = 3 (default and minimal number of harmonic terms possible), and threshold level = 0.05. To detect multiple breakpoints, the algorithm was iteratively run on the MODIS NDVI 2000–2021 image stack, each 3 months, using a 3-year training period (after testing of two training period lengths hl = 2 vs. hl = 3 years) and a monitoring period of one year given that the BFAST monitor change magnitude is relative to the monitoring period length. No “penalty” period after each detected breakpoint (as opposed to Awty-Carroll (2019)) was applied to avoid missing significant changes. It is worth to note that this implementation of BFASTmonitor does not include the automatic determination of a stable training period yet.

### 2.3.1.2. BFAST Lite.
The newly BFAST Lite unsupervised algorithm is built upon the BFAST algorithm (Verbesselt et al., 2010a) with the aim to improve its speed and flexibility (Masiliūnas et al., 2021). Compared to BFAST, this offline approach avoids the seasonal-trend decomposition and performs the model fitting in a single step, using a multivariate piecewise linear harmonic regression. In addition, it provides more robust statistics for breakpoint magnitude calculation, such as the Root Mean Squared Deviation (RMSD) and the Mean Absolute Deviation (MAD), which are computed between the predicted values of the adjacent segments over the time span of one year before and after the detected break, as opposed to BFAST that computes the difference of the fitted value immediately before and after the break (Masiliūnas et al., 2021).

This study used the BFAST 1.6.1 R package available at https://cran.r-project.org/web/packages/bfast/index.html. Default parameters were used, with the number of harmonic terms equals to 3. The biggest breakpoint in magnitude was selected using the root mean squared deviation (RMSD).

### 2.3.1.3. Edyn: The dynamic version of EWMACD.
EWMACD (Exponentially Weighted Moving Average Change Detection) is an online monitoring algorithm that aims to detect persistent subtle changes, such as forest degradation or thinning (Brooks et al., 2014). As BFASTmonitor, EWMACD uses a statistical control chart (here, the EWMA) on the residuals to detect deviations from the mean. In the original version, the algorithm does not retrain after having flagged a breakpoint. In Edyn,

the dynamic version of EWMACD (Brooks et al., 2017), the harmonic model coefficients are dynamically updated, and the optimal training period length is automatically found based on the quality of the model fit.

Global parameters defined by the end-user are: 1- the parameter $\lambda$ which defines the algorithm's robustness to low signal-to-noise ratios, and 2- the persistence that refers to the number of times flagged deviations must be successively detected for a change to be considered (in Edyn it is given as a proportion of a year). Finally, and with respect to the breakpoint magnitude, while originally given as a standardized value (i.e. the residual value divided by the chart's control limits, and then rounded down to the nearest integer value), the magnitudes in this study correspond to the residual values at the breakpoints' location. The reason for this choice is that, when ran continuously, more than one breakpoint may have the same standardized value impeding the selection of a unique breakpoint.

In this study the same values as in Awty-Carroll (2019) were applied ($\lambda$ = 0.3, persistence = 6). Their R version of Edyn (available at: https://github.com/klh5/season-trend-comparison/tree/master/ewmacd)

was used, adapted from the original R script (Brooks et al., 2014) and Edyn (Brooks et al., 2017). This version enables selection of stable (without breakpoints) sliding windows of fixed length (2 years) that are used as training periods, thus allowing continuous monitoring. Compared to the original Edyn this adaptation does not find the optimal training period length based on the model fit quality.

### 2.3.2. BFASTm-L2: A new approach

The procedure to select the BFASTmonitor breakpoint linked to the biggest pattern change in long-term NDVI SITS is detailed in this section. The approach (hereafter BFASTm-L2) is based on the Euclidean distance (L2 distance) and is schematized in Fig. 3.

In step 1 of Fig. 3, BFASTmonitor is successively run each 3 months over the entire time series. Then, for each detected breakpoint, the time series segments L2-w (L2-w = 3) before and after the breakpoint are extracted (step 2), and monthly averaged in annual subsamples (step 3). In step 4, the Euclidean distance between the 2 annual subsamples is computed using the Python numpy.linalg.norm function. Finally, the breakpoint with the highest L2 distance is selected (step 5).

This selection procedure is somehow similar to the method in Setiawan and Yoshino (2012), where a mean-based distance measure is computed between each two successive annual segments. It however differs in many aspects: 1- The distance metric is computed only where BFASTmonitor breakpoints are detected; 2- The segment length used here is 3 years (instead of 1 year), enabling to skip non-persistent changes mainly due to climate variability; 3- The (L2) distance used is more representative of an overall pattern change than the mean-based distance used in Setiawan and Yoshino (2012).
<caption>**Table 2**</caption>

Main characteristics of the change-detection algorithms used in this study, their strengths and limitations.

<table><thead><tr><th></th><th>Algorithm</th><th>Speed</th><th>Parameter number *</th><th>Model</th><th>Decomposition</th><th>Strengths</th><th>Limitations (others than speed and number of parameters)</th></tr></thead><tbody><tr><td>Offline</td><td>BFAST Lite (Masiliūnas et al., 2021)</td><td>++</td><td>+</td><td>Model-deviation</td><td>None</td><td>Detection of abrupt and trend changes. Possibility to select a model with a unique (the highest-magnitude) breakpoint.</td><td>Low performance in capturing seasonal changes</td></tr><tr><td>Online</td><td>BFAST monitor (Verbesselt et al., 2012)<br>EDYN (Brooks et al., 2017)</td><td>+++</td><td>+</td><td>Model-deviation</td><td>None</td><td>Sensitivity to seasonal changes (Awty-Carroll, 2019). Detection of subtle (sub-pixel) changes in long time series (model dynamically retrained); Able to capture seasonal changes</td><td>High false positive rate (Masiliūnas et al., 2021)<br>High sensitivity to algorithm parameters (Saxena et al., 2018)</td></tr></tbody></table>

* Number of key parameters (others than if applied) the order of the harmonic term, the probability threshold/ statistical significance level and the training period length):

$$ \leq 2: +, \geq 3: - $$

them (Gao et al., 2021; Hamunyela, 2017). However, its massively-parallel GPU implementation (Gieseke et al., 2020), that makes it 14.5 times faster than the newly launched BFAST Lite (Masiliūnas et al., 2021), its implementation on Google Earth Engine (Hamunyela et al., 2020) and its good performance in detecting seasonal changes (Awty-Carroll, 2019), make of BFASTmonitor a potential change detection algorithm for the purposes of this study. This study made use of the Python package bfast0.7, available at https://pypi.org/project/bfast/, that provides a parallel implementation of the BFASTmonitor algorithm. The parameters used were: h = 0.25 (the MOSUM bandwidth), k = 3 (default and minimal number of harmonic terms possible), and threshold level = 0.05. To detect multiple breakpoints, the algorithm was iteratively run on the MODIS NDVI 2000–2021 image stack, each 3 months, using a 3-year training period (after testing of two training period lengths hl = 2 vs. hl = 3 years) and a monitoring period of one year given that the BFAST monitor change magnitude is relative to the monitoring period length. No “penalty” period after each detected breakpoint (as opposed to Awty-Carroll (2019)) was applied to avoid missing significant changes. It is worth to note that this implementation of BFASTmonitor does not include the automatic determination of a stable training period yet.

### 2.3.1.2. BFAST Lite.
The newly BFAST Lite unsupervised algorithm is built upon the BFAST algorithm (Verbesselt et al., 2010a) with the aim to improve its speed and flexibility (Masiliūnas et al., 2021). Compared to BFAST, this offline approach avoids the seasonal-trend decomposition and performs the model fitting in a single step, using a multivariate piecewise linear harmonic regression. In addition, it provides more robust statistics for breakpoint magnitude calculation, such as the Root Mean Squared Deviation (RMSD) and the Mean Absolute Deviation (MAD), which are computed between the predicted values of the adjacent segments over the time span of one year before and after the detected break, as opposed to BFAST that computes the difference of the fitted value immediately before and after the break (Masiliūnas et al., 2021).

This study used the BFAST 1.6.1 R package available at https://cran.r-project.org/web/packages/bfast/index.html. Default parameters were used, with the number of harmonic terms equals to 3. The biggest breakpoint in magnitude was selected using the root mean squared deviation (RMSD).

### 2.3.1.3. Edyn: The dynamic version of EWMACD.
EWMACD (Exponentially Weighted Moving Average Change Detection) is an online monitoring algorithm that aims to detect persistent subtle changes, such as forest degradation or thinning (Brooks et al., 2014). As BFASTmonitor, EWMACD uses a statistical control chart (here, the EWMA) on the residuals to detect deviations from the mean. In the original version, the algorithm does not retrain after having flagged a breakpoint. In Edyn,

the dynamic version of EWMACD (Brooks et al., 2017), the harmonic model coefficients are dynamically updated, and the optimal training period length is automatically found based on the quality of the model fit.

Global parameters defined by the end-user are: 1- the parameter $\lambda$ which defines the algorithm's robustness to low signal-to-noise ratios, and 2- the persistence that refers to the number of times flagged deviations must be successively detected for a change to be considered (in Edyn it is given as a proportion of a year). Finally, and with respect to the breakpoint magnitude, while originally given as a standardized value (i.e. the residual value divided by the chart's control limits, and then rounded down to the nearest integer value), the magnitudes in this study correspond to the residual values at the breakpoints' location. The reason for this choice is that, when ran continuously, more than one breakpoint may have the same standardized value impeding the selection of a unique breakpoint.

In this study the same values as in Awty-Carroll (2019) were applied ($\lambda$ = 0.3, persistence = 6). Their R version of Edyn (available at: https://github.com/klh5/season-trend-comparison/tree/master/ewmacd)

was used, adapted from the original R script (Brooks et al., 2014) and Edyn (Brooks et al., 2017). This version enables selection of stable (without breakpoints) sliding windows of fixed length (2 years) that are used as training periods, thus allowing continuous monitoring. Compared to the original Edyn this adaptation does not find the optimal training period length based on the model fit quality.

### 2.3.2. BFASTm-L2: A new approach

The procedure to select the BFASTmonitor breakpoint linked to the biggest pattern change in long-term NDVI SITS is detailed in this section. The approach (hereafter BFASTm-L2) is based on the Euclidean distance (L2 distance) and is schematized in Fig. 3.

In step 1 of Fig. 3, BFASTmonitor is successively run each 3 months over the entire time series. Then, for each detected breakpoint, the time series segments L2-w (L2-w = 3) before and after the breakpoint are extracted (step 2), and monthly averaged in annual subsamples (step 3). In step 4, the Euclidean distance between the 2 annual subsamples is computed using the Python numpy.linalg.norm function. Finally, the breakpoint with the highest L2 distance is selected (step 5).

This selection procedure is somehow similar to the method in Setiawan and Yoshino (2012), where a mean-based distance measure is computed between each two successive annual segments. It however differs in many aspects: 1- The distance metric is computed only where BFASTmonitor breakpoints are detected; 2- The segment length used here is 3 years (instead of 1 year), enabling to skip non-persistent changes mainly due to climate variability; 3- The (L2) distance used is more representative of an overall pattern change than the mean-based distance used in Setiawan and Yoshino (2012).
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Flowchart of the BFASTm-L2 approach (in this example with L2-w = 3 years). Blue lines represent the NDVI time series subsamples (dark blue before breakpoint, and light blue after breakpoint). Detected BFASTmonitor breakpoints are presented with dashed red vertical lines. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

### 2.3.3. Evaluation of the algorithms' performance

#### 2.3.3.1. Temporal accuracy to different types and intensities of change.

The temporal accuracy of the single breakpoint detected by the four algorithms (BFASTmonitor, BFAST Lite, Edyn, and BFASTm-L2) was performed as following: 1- Each algorithm was run on each simulated time series, 2- the highest-magnitude breakpoint was selected (no selection performed for BFAST Lite), and 3- the absolute difference between the breakpoint date and the date of change (January of 2011) was recorded. Finally, bar charts of the selected breakpoints' date per type of change relative to the true date of change were plotted to assess the performance of each algorithm to accurately detect in time the different types of change.

#### 2.3.3.2. Breakpoint-magnitude sensitivity to different types and intensities of change.

The breakpoint magnitude is often used to spatially detect significant hotspots of change. To unravel the type of change most likely detected by each algorithm (and thus spatially highlighted), the

sensitivity of each algorithm to different types and intensities of change was assessed. For that, the distributions of the magnitude of the detected breakpoint were calculated per type of change, and represented using violin plots. Absolute values were used and normalized for each algorithm to allow comparisons between them. For all the change types (except the change in trend only, and the no-change category), only the breakpoints comprised between 2010 and 2012 (included) were used.

#### 2.3.3.3. Running times of the algorithms.

The 5-run average times of BFASTmonitor, BFAST Lite, Edyn and BFASTm-L2, but also of the L2 distance metric alone computed sequentially each 3 months, were reported for areas of different sizes: 50x50, 100x100, 150x150, 200x200 pixels (the whole red box Fig. 1), using the MODIS NDVI 2000–2021 data set. Six CPU cores (parallel processing) were used on a 64 GB RAM computer.

#### 2.3.3.4. Detection and mapping of LSLAs induced LULC changes: Individual study cases.

Evaluation of the change detection accuracy on the four
<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Flowchart of the BFASTm-L2 approach (in this example with L2-w = 3 years). Blue lines represent the NDVI time series subsamples (dark blue before breakpoint, and light blue after breakpoint). Detected BFASTmonitor breakpoints are presented with dashed red vertical lines. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

### 2.3.3. Evaluation of the algorithms' performance

#### 2.3.3.1. Temporal accuracy to different types and intensities of change.

The temporal accuracy of the single breakpoint detected by the four algorithms (BFASTmonitor, BFAST Lite, Edyn, and BFASTm-L2) was performed as following: 1- Each algorithm was run on each simulated time series, 2- the highest-magnitude breakpoint was selected (no selection performed for BFAST Lite), and 3- the absolute difference between the breakpoint date and the date of change (January of 2011) was recorded. Finally, bar charts of the selected breakpoints' date per type of change relative to the true date of change were plotted to assess the performance of each algorithm to accurately detect in time the different types of change.

#### 2.3.3.2. Breakpoint-magnitude sensitivity to different types and intensities of change.

The breakpoint magnitude is often used to spatially detect significant hotspots of change. To unravel the type of change most likely detected by each algorithm (and thus spatially highlighted), the

sensitivity of each algorithm to different types and intensities of change was assessed. For that, the distributions of the magnitude of the detected breakpoint were calculated per type of change, and represented using violin plots. Absolute values were used and normalized for each algorithm to allow comparisons between them. For all the change types (except the change in trend only, and the no-change category), only the breakpoints comprised between 2010 and 2012 (included) were used.

#### 2.3.3.3. Running times of the algorithms.

The 5-run average times of BFASTmonitor, BFAST Lite, Edyn and BFASTm-L2, but also of the L2 distance metric alone computed sequentially each 3 months, were reported for areas of different sizes: 50x50, 100x100, 150x150, 200x200 pixels (the whole red box Fig. 1), using the MODIS NDVI 2000–2021 data set. Six CPU cores (parallel processing) were used on a 64 GB RAM computer.

#### 2.3.3.4. Detection and mapping of LSLAs induced LULC changes: Individual study cases.

Evaluation of the change detection accuracy on the four
real study cases presented in 2.1 was performed for: BFASTmonitor (training period of 3 years), BFAST Lite, Edyn (with lambda = 0.3) and BFASTm-L2 (using a training period length and L2-w size of 3 years), which were run continuously every 3 months. All the detected breakpoints, and selected highest-magnitude breakpoint were recorded (date and magnitudes).

Change maps based on the highest-magnitude breakpoint found in the 2000–2020 MODIS NDVI time series by BFASTm-L2, BFAST Lite, BFASTmonitor and Edyn were produced over the 200 × 200 MODIS pixels study area. The magnitudes for each map were then normalized between 0 and 1 in order to allow comparison. To quantitatively assess the performance of each method in detecting LSLAs induced changes, the average breakpoint magnitude within and outside the LSLAs (represented by red polygons in the study area) were computed and are presented in Table 3.

# 3. Results

## 3.1. Temporal accuracy and magnitude sensitivity of the single breakpoint

### 3.1.1. Temporal accuracy of the breakpoint

BFASTmonitor, BFAST Lite, Edyn and BFASTm-L2 were tested on the benchmark dataset (cf. section 2.1) to assess the performance of each method to accurately detect in time the different types of change (occurring in January of 2011). The breakpoints with the highest magnitude were selected, and their temporal distributions relative to the true change are presented in Fig. 4 (cf. section 3.3.1).

As a first observation one can see that most of the highest-magnitude breakpoints (of almost any change intensity), and with the exception of BFASTmonitor (4.a), are close in time to the true real change (less than 6 months). The algorithms have however slightly different sensitivity to the different types of changes:

- Regarding the seasonal changes (amplitude, LOS and NOS), BFASTm-L2 breakpoint magnitudes (4.d) are the most responsive to this type of change (>76% of the breakpoints in each subcategory of change are located ± 6 months of the true change). When looking more in detail, and compared to the other algorithms, BFASTm-L2 breakpoint magnitude is particularly sensitive to LOS changes (84.6% vs 70.1% for BFAST Lite and 41.6% for Edyn, at ± 6 months), and is as good as BFAST Lite in detecting NOS changes (98% vs. 100% at 6 ± months) and amplitude changes (76.8% for both at ± 6 months).

- In what concerns gradual changes (Trend), all the algorithms have at least one breakpoint induced by this type of change, mostly occurring one year before/after the true date of change.

- BFAST Lite (4.a) breakpoint magnitude is very responsive to the break/trend changes (99% of the changes detected within 6 months vs. 83.4% for both BFASTm-L2 and Edyn).

- Finally, regarding the no-change type, Edyn (4.c) was the only algorithm to correctly detect no change in almost 80% of the cases (as seen with the NA class).

<caption>**Table 3**</caption>

Average breakpoint magnitude (mean ± standard deviation) in LSLAs (red polygons in the study area) and the whole study area.

<table><thead><tr><th></th><th>BFAST Lite</th><th>BFASTmonitor</th><th>BFASTm-L2</th><th>Edyn</th></tr></thead><tbody><tr><td>Breakpoint magnitude average in LSLAs (mean ± std)</td><td>0.22 ± 0.14 (n = 258)</td><td>0.60 ± 0.15 (n = 262)</td><td>0.50 ± 0.13 (n = 262)</td><td>0.46 ± 0.16 (n = 258)</td></tr><tr><td>Breakpoint magnitude average outside LSLAs (mean ± std)</td><td>0.22 ± 0.10 (n = 37,653)</td><td>0.38 ± 0.13 (n = 37,343)</td><td>0.23 ± 0.07 (n = 37,343)</td><td>0.28 ± 0.10 (n = 37,684)</td></tr><tr><td>Difference (in-out) (mean ± std)</td><td>0.0 ± 0.17</td><td>0.22 ± 0.20</td><td><strong>0.27 ± 0.15</strong></td><td>0.18 ± 0.19</td></tr></tbody></table>

Those first results show that selecting the BFASTm-L2 highest-magnitude breakpoint to correctly detect in time seasonal changes is effective. However, as one can expect multiple breakpoints detected in long and dense time series triggered by different types of change, the next question that comes up is: how comparable are the breakpoint magnitudes induced by the different change types?

As a small parenthesis, it is worth to note that the performances of Edyn, BFAST Lite and BFASTmonitor are different than in Awty-Carroll et al. (2019). This is explained by the fact that: 1- A subsample of the original dataset was used (dataset without missing data), and smoothed; 2- BFASTmonitor was run continuously without considering any “penalty” period after each detected breakpoint; 3- Only the distributions of the highest-magnitude breakpoints were evaluated; 4- The “correct” detection period considered here for the abrupt changes is larger than in Awty-Carroll et al. (2019), who considered a maximum period of 3 months for the abrupt changes and 1 year for the seasonal ones.

### 3.1.2. Breakpoint-magnitude sensitivity to the type and intensity of change

Fig. 5 shows the violin plots of the breakpoint's magnitude for each algorithm and type of change, as explained in 2.4.2. This type of representation helps in the identification of the type of change that will most likely be highlighted on a map when using the breakpoint magnitude as the mapping variable.

Some general findings can be drawn from Fig. 5. First, one can observe that all the algorithms, with the exception of Edyn, have their highest breakpoint magnitudes (>0.4) induced by medium to large abrupt changes (i.e. the trend/break data subset without any associated trend change) (>0.2 NDVI units; cf. Table 1). Detection of seasonal changes will therefore very likely be hindered in time series with large abrupt changes. Second, and with the exception of the no-change category, the lowest magnitudes are globally associated to the LOS changes. As such, one can expect LOS changes to be hardly detected in breakpoint-magnitude change maps.

Going more in depth with the algorithms' analysis, one can see that the BFASTm-L2 distributions means are the highest for two out of three (LOS and NOS) of the seasonal changes, with a remarkably high average for the NOS changes. In addition, the LOS and NOS distributions are well separated from the “No-change” category. This is important as it ensures the ability of the algorithm to highlight seasonal changes in breakpoint magnitude-based change maps. As a reminder, this remains possible as long as there are no large abrupt changes, or moderate to larger gradual changes (trends). Indeed BFASTm-L2, contrarily to BFASTmonitor or BFAST Lite, is sensitive to trends. Regarding the seasonal changes, BFASTmonitor, BFAST Lite and Edyn respond much better (higher magnitudes) to the amplitude changes than to NOS/LOS changes. This is particularly true for Edyn and even more for BFASTmonitor, which distribution base (for the smallest change intensity) is at almost 0.3.

Finally, it is important to remember that the values represented in Fig. 5 were obtained from short and smoothed simulated time series, containing a unique change. As a result, the “No change” values for the different algorithms could in fact be higher in real time series with multiple small types of change, thus impacting the performance of each algorithm. To assess this, the different algorithms were tested on real time series and over a small area. Results are presented in the next section.

## 3.2. Detection of LSLAs LULC driven changes

### 3.2.1. Assessment on individual (pixel) study cases

In this section, the capacity of the four algorithms to detect real seasonal changes on the four study cases (distributed from arid to humid conditions, see Fig. 1) is evaluated. The overview of Fig. 6, on which are overlaid the 2000–2021 MODIS NDVI time series and the period of change (grey boxes) for each study case, shows the diversity of type and intensity of changes found in the dataset: seasonal changes (all study cases present changes in amplitude, and study cases 1 to 3 show
real study cases presented in 2.1 was performed for: BFASTmonitor (training period of 3 years), BFAST Lite, Edyn (with lambda = 0.3) and BFASTm-L2 (using a training period length and L2-w size of 3 years), which were run continuously every 3 months. All the detected breakpoints, and selected highest-magnitude breakpoint were recorded (date and magnitudes).

Change maps based on the highest-magnitude breakpoint found in the 2000–2020 MODIS NDVI time series by BFASTm-L2, BFAST Lite, BFASTmonitor and Edyn were produced over the 200 × 200 MODIS pixels study area. The magnitudes for each map were then normalized between 0 and 1 in order to allow comparison. To quantitatively assess the performance of each method in detecting LSLAs induced changes, the average breakpoint magnitude within and outside the LSLAs (represented by red polygons in the study area) were computed and are presented in Table 3.

# 3. Results

## 3.1. Temporal accuracy and magnitude sensitivity of the single breakpoint

### 3.1.1. Temporal accuracy of the breakpoint

BFASTmonitor, BFAST Lite, Edyn and BFASTm-L2 were tested on the benchmark dataset (cf. section 2.1) to assess the performance of each method to accurately detect in time the different types of change (occurring in January of 2011). The breakpoints with the highest magnitude were selected, and their temporal distributions relative to the true change are presented in Fig. 4 (cf. section 3.3.1).

As a first observation one can see that most of the highest-magnitude breakpoints (of almost any change intensity), and with the exception of BFASTmonitor (4.a), are close in time to the true real change (less than 6 months). The algorithms have however slightly different sensitivity to the different types of changes:

- Regarding the seasonal changes (amplitude, LOS and NOS), BFASTm-L2 breakpoint magnitudes (4.d) are the most responsive to this type of change (>76% of the breakpoints in each subcategory of change are located ± 6 months of the true change). When looking more in detail, and compared to the other algorithms, BFASTm-L2 breakpoint magnitude is particularly sensitive to LOS changes (84.6% vs 70.1% for BFAST Lite and 41.6% for Edyn, at ± 6 months), and is as good as BFAST Lite in detecting NOS changes (98% vs. 100% at 6 ± months) and amplitude changes (76.8% for both at ± 6 months).

- In what concerns gradual changes (Trend), all the algorithms have at least one breakpoint induced by this type of change, mostly occurring one year before/after the true date of change.

- BFAST Lite (4.a) breakpoint magnitude is very responsive to the break/trend changes (99% of the changes detected within 6 months vs. 83.4% for both BFASTm-L2 and Edyn).

- Finally, regarding the no-change type, Edyn (4.c) was the only algorithm to correctly detect no change in almost 80% of the cases (as seen with the NA class).

<caption>**Table 3**</caption>

Average breakpoint magnitude (mean ± standard deviation) in LSLAs (red polygons in the study area) and the whole study area.

<table><thead><tr><th></th><th>BFAST Lite</th><th>BFASTmonitor</th><th>BFASTm-L2</th><th>Edyn</th></tr></thead><tbody><tr><td>Breakpoint magnitude average in LSLAs (mean ± std)</td><td>0.22 ± 0.14 (n = 258)</td><td>0.60 ± 0.15 (n = 262)</td><td>0.50 ± 0.13 (n = 262)</td><td>0.46 ± 0.16 (n = 258)</td></tr><tr><td>Breakpoint magnitude average outside LSLAs (mean ± std)</td><td>0.22 ± 0.10 (n = 37,653)</td><td>0.38 ± 0.13 (n = 37,343)</td><td>0.23 ± 0.07 (n = 37,343)</td><td>0.28 ± 0.10 (n = 37,684)</td></tr><tr><td>Difference (in-out) (mean ± std)</td><td>0.0 ± 0.17</td><td>0.22 ± 0.20</td><td>0.27 ± 0.15</td><td>0.18 ± 0.19</td></tr></tbody></table>

Those first results show that selecting the BFASTm-L2 highest-magnitude breakpoint to correctly detect in time seasonal changes is effective. However, as one can expect multiple breakpoints detected in long and dense time series triggered by different types of change, the next question that comes up is: how comparable are the breakpoint magnitudes induced by the different change types?

As a small parenthesis, it is worth to note that the performances of Edyn, BFAST Lite and BFASTmonitor are different than in Awty-Carroll et al. (2019). This is explained by the fact that: 1- A subsample of the original dataset was used (dataset without missing data), and smoothed; 2- BFASTmonitor was run continuously without considering any “penalty” period after each detected breakpoint; 3- Only the distributions of the highest-magnitude breakpoints were evaluated; 4- The “correct” detection period considered here for the abrupt changes is larger than in Awty-Carroll et al. (2019), who considered a maximum period of 3 months for the abrupt changes and 1 year for the seasonal ones.

### 3.1.2. Breakpoint-magnitude sensitivity to the type and intensity of change

Fig. 5 shows the violin plots of the breakpoint's magnitude for each algorithm and type of change, as explained in 2.4.2. This type of representation helps in the identification of the type of change that will most likely be highlighted on a map when using the breakpoint magnitude as the mapping variable.

Some general findings can be drawn from Fig. 5. First, one can observe that all the algorithms, with the exception of Edyn, have their highest breakpoint magnitudes (>0.4) induced by medium to large abrupt changes (i.e. the trend/break data subset without any associated trend change) (>0.2 NDVI units; cf. Table 1). Detection of seasonal changes will therefore very likely be hindered in time series with large abrupt changes. Second, and with the exception of the no-change category, the lowest magnitudes are globally associated to the LOS changes. As such, one can expect LOS changes to be hardly detected in breakpoint-magnitude change maps.

Going more in depth with the algorithms' analysis, one can see that the BFASTm-L2 distributions means are the highest for two out of three (LOS and NOS) of the seasonal changes, with a remarkably high average for the NOS changes. In addition, the LOS and NOS distributions are well separated from the “No-change” category. This is important as it ensures the ability of the algorithm to highlight seasonal changes in breakpoint magnitude-based change maps. As a reminder, this remains possible as long as there are no large abrupt changes, or moderate to larger gradual changes (trends). Indeed BFASTm-L2, contrarily to BFASTmonitor or BFAST Lite, is sensitive to trends. Regarding the seasonal changes, BFASTmonitor, BFAST Lite and Edyn respond much better (higher magnitudes) to the amplitude changes than to NOS/LOS changes. This is particularly true for Edyn and even more for BFASTmonitor, which distribution base (for the smallest change intensity) is at almost 0.3.

Finally, it is important to remember that the values represented in Fig. 5 were obtained from short and smoothed simulated time series, containing a unique change. As a result, the “No change” values for the different algorithms could in fact be higher in real time series with multiple small types of change, thus impacting the performance of each algorithm. To assess this, the different algorithms were tested on real time series and over a small area. Results are presented in the next section.

## 3.2. Detection of LSLAs LULC driven changes

### 3.2.1. Assessment on individual (pixel) study cases

In this section, the capacity of the four algorithms to detect real seasonal changes on the four study cases (distributed from arid to humid conditions, see Fig. 1) is evaluated. The overview of Fig. 6, on which are overlaid the 2000–2021 MODIS NDVI time series and the period of change (grey boxes) for each study case, shows the diversity of type and intensity of changes found in the dataset: seasonal changes (all study cases present changes in amplitude, and study cases 1 to 3 show
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

**Fig. 4.** Temporal distributions of the highest-absolute magnitude breakpoint detected by BFASTmonitor (4.a), BFAST Lite (4.b), Edyn (4.c) and BFASTm-L2 (4.d), for six types of change (amplitude, length of season (LOS), number of season (NOS), monotonic trend (Trend), abrupt change (Break/trend) and no-change (None)). The x-axis represents the absolute difference in months between the breakpoint's date of change and the real date of change (January of 2011). The y-axis represents the proportion of samples detected at each x-axis unit. The stacked bars colours indicate the different intensities of change (see Table 1). The second y-axis (in orange) represents the cumulative percentage of the number of samples within each change type. A bar at 6 months was arbitrarily added to allow comparisons (% of samples with date of change ≤ 6 months from real date of change). The NA class represents the cases for which no breakpoint was detected. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

**Fig. 4.** Temporal distributions of the highest-absolute magnitude breakpoint detected by BFASTmonitor (4.a), BFAST Lite (4.b), Edyn (4.c) and BFASTm-L2 (4.d), for six types of change (amplitude, length of season (LOS), number of season (NOS), monotonic trend (Trend), abrupt change (Break/trend) and no-change (None)). The x-axis represents the absolute difference in months between the breakpoint's date of change and the real date of change (January of 2011). The y-axis represents the proportion of samples detected at each x-axis unit. The stacked bars colours indicate the different intensities of change (see Table 1). The second y-axis (in orange) represents the cumulative percentage of the number of samples within each change type. A bar at 6 months was arbitrarily added to allow comparisons (% of samples with date of change ≤ 6 months from real date of change). The NA class represents the cases for which no breakpoint was detected. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
<figure><img src="image_9.png" /><figcaption>**Fig. 5.** Violin plots of the normalized highest-magnitude breakpoint, per type of change for: BFASTm-L2 (green), BFAST Lite (orange), BFASTmonitor (white) and Edyn (blue). Dots represent the mean of each distribution. “Abrupt” refers to the Trend/break data subset without any trend change (refer to section 2.2.1). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 6.** Application of the change detection algorithms to four Senegalese study cases showing agricultural LSLA implementations. For each study case, the corresponding pre-processed MODIS 2000–2021 NDVI time series with the different breakpoints detected (BFAST Lite)/selected (all others) are presented in the a subparts of the figures at the left. Figures at the right correspond to Google Earth snapshots closest in time with the LSLA implementation: the top/bottom snapshot for the closest date available before/after the change. In the left figures: dashed lines correspond to the highest-magnitude breakpoints: green for BFASTm-L2, orange for BFAST Lite (1 break), black for BFASTmonitor and blue for Edyn (break dates are also given in the same colors to facilitate identification when breakpoints are superimposed). Grey shaded areas correspond to the period of change (1 year before and after the observed date of change). Grey subplots (subparts b) at the bottom present all the breakpoints detected by BFASTmonitor (black), BFASTm-L2 (green) and Edyn (blue), along with their magnitude (height of vertical lines). Red dots pinpoint the highest-magnitude breakpoint (in absolute value). Hatched areas correspond to periods without monitoring (because of the initial period of time needed for training, and for computation of L2). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

important NOS changes), an extreme abrupt change (~+0.25 NDVI units) for study case 1 in 2018, and a strong positive trend (study case 2, start in 2016).

On Fig. 6 (subpart b of each subplot), the multiple breakpoints detected by BFASTmonitor, BFASTm-L2 and Edyn are represented with vertical lines. The selected breakpoints are the ones having the highest magnitude, and are highlighted with a red dot at the extremity. BFAST

Lite, BFASTm-L2 and Edyn performances (position of the unique breakpoint) are overall similar, with the following exceptions: i) a gap of 3/4 years is observed between BFASTm-L2 and BFAST Lite/Edyn in study case 2 and, ii) Edyn selected breakpoint is far away (>10 years) BFASTm-L2/BFAST Lite breakpoints in study case 4. Compared to BFASTmonitor, Edyn found significant fewer breakpoints, mainly due to the used-fixed period length needed after each detected breakpoint to
<figure><img src="image_10.png" /><figcaption>**Fig. 5.** Violin plots of the normalized highest-magnitude breakpoint, per type of change for: BFASTm-L2 (green), BFAST Lite (orange), BFASTmonitor (white) and Edyn (blue). Dots represent the mean of each distribution. “Abrupt” refers to the Trend/break data subset without any trend change (refer to section 2.2.1). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 6.** Application of the change detection algorithms to four Senegalese study cases showing agricultural LSLA implementations. For each study case, the corresponding pre-processed MODIS 2000–2021 NDVI time series with the different breakpoints detected (BFAST Lite)/selected (all others) are presented in the a subparts of the figures at the left. Figures at the right correspond to Google Earth snapshots closest in time with the LSLA implementation: the top/bottom snapshot for the closest date available before/after the change. In the left figures: dashed lines correspond to the highest-magnitude breakpoints: green for BFASTm-L2, orange for BFAST Lite (1 break), black for BFASTmonitor and blue for Edyn (break dates are also given in the same colors to facilitate identification when breakpoints are superimposed). Grey shaded areas correspond to the period of change (1 year before and after the observed date of change). Grey subplots (subparts b) at the bottom present all the breakpoints detected by BFASTmonitor (black), BFASTm-L2 (green) and Edyn (blue), along with their magnitude (height of vertical lines). Red dots pinpoint the highest-magnitude breakpoint (in absolute value). Hatched areas correspond to periods without monitoring (because of the initial period of time needed for training, and for computation of L2). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

important NOS changes), an extreme abrupt change (~+0.25 NDVI units) for study case 1 in 2018, and a strong positive trend (study case 2, start in 2016).

On Fig. 6 (subpart b of each subplot), the multiple breakpoints detected by BFASTmonitor, BFASTm-L2 and Edyn are represented with vertical lines. The selected breakpoints are the ones having the highest magnitude, and are highlighted with a red dot at the extremity. BFAST

Lite, BFASTm-L2 and Edyn performances (position of the unique breakpoint) are overall similar, with the following exceptions: i) a gap of 3/4 years is observed between BFASTm-L2 and BFAST Lite/Edyn in study case 2 and, ii) Edyn selected breakpoint is far away (>10 years) BFASTm-L2/BFAST Lite breakpoints in study case 4. Compared to BFASTmonitor, Edyn found significant fewer breakpoints, mainly due to the used-fixed period length needed after each detected breakpoint to
PART II

<figure><img src="image_11.png" /><figcaption>Fig. 7. PART I) From left to right, top to bottom: 1) Snapshot view of the study area (Map data ©2015 Google), (2) MODIS NDVI average over 2000–2020, 3–4.5–6) Breakpoints (normalized) magnitude maps using: (3) BFASTm-L2, (4) BFAST Lite, (5) BFASTmonitor, (6) Edyn. Red polygons represent the active agricultural LSLAs reported in the ISRA ground-field database. PART II) Google Earth zoom-in of points A, B, C, D plotted in all the maps of PART I. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

ensure stability for model retraining.

In study case 1, two major events occur: a central-pivot irrigation system was installed at the end of 2011 (conversion from small to large-scale agriculture), and as from 2018 a huge increase in productivity is observed. The first change mostly translates as seasonal changes (NOS, amplitude), while the second induces a sharp abrupt change in 2018. In this case, BFAST Lite, BFASTm-L2, but also Edyn, output a higher-magnitude breakpoint for the abrupt change than the seasonal one, which is in line with the findings made in section 4.1.2 (Fig. 5).

Study case 2 (conversion from natural vegetation), is also concerned by the installation of a central-pivot irrigation system at the beginning of 2016. It translates in the signal as seasonal changes (NOS, amplitude), with a positive trend. In this case BFAST Lite and Edyn selected breakpoints are closely located in time (~1 year of difference), and both differs (>3 years) from BFASTm-L2 selected breakpoint. While BFAST Lite and Edyn breakpoints are closely located from big changes in amplitude, BFASTm-L2 highest magnitude breakpoint is located at the beginning of 2017, when seasonal changes are most marked, and most probably related to the concession implementation and start of activities.

Study case 3 is another example of transition from natural vegetation to intensive agriculture, but in a less arid environment. The transition is visually observed in 2008 and mostly related to seasonal changes (NOS). Despite the high variability exhibited in the time series (particularly after 2008), BFAST Lite and BFASTm-L2 detects the same breakpoint, while Edyn differs from 9 months away in the middle of the second growing cycle of the first year of production.

Lastly, study case 4 located in the most humid environment, shows more stability than the other study cases. Indeed, the land use change does not induce a significant change in the time series shape. The change being sought translates in a significant change in amplitude and a negative trend. Here again, BFAST Lite and BFASTm-L2 agreed in the breakpoint detected, inside the sought period of change. Edyn diverges from many years away, and outputs its highest-magnitude breakpoint where an amplitude change is observed. As a reminder, the Edyn breakpoint magnitude used here is the residual of the detected breakpoint. While the standardized magnitude should be more adapted for the comparison of breakpoint magnitudes, it was not used here as more than one breakpoint may share the same maximum value thus hindering the selection of a unique breakpoint. Edyn sensitivity to the parameter lambda was tested, with the parameter set to 0.6 (Appendix A). All the highest-magnitude breakpoints found differed from those presented in

Fig. 6, including for study case 1 which contains a huge abrupt change.

3.2.2. Detection of LSLAs induced LULC changes in breakpoint magnitude maps

So far, it has been shown that the performances of BFAST Lite and BFASTm-L2 in detecting in time seasonal changes are similar. However, because of the different sensitivities of the breakpoint magnitudes to the different types of change (Fig. 5), differences in the change maps based on the breakpoint magnitudes are expected. The change maps obtained with the different algorithms were evaluated over a study area close to Dakar, in which several LSLAs have been implemented since 2003 (Fig. 7).

The change maps based on the breakpoint magnitude display different spatial patterns. While LSLAs (red polygons) in the BFASTm-L2 (subplot I.3) change map are clearly highlighted and have significantly higher magnitude values than the background (highest difference for BFASTm-L2 as shown in Table 3), this is not the case for BFAST Lite (subplot I.4) which difference between the values inside and outside LSLAs is null.

BFAST Lite (subplot I.4) highest magnitudes are mainly located in two areas of the map: 1- in the centre of the map, in the wooded vegetation of the reserve of Bandia (see Appendix B.1) and, 2- in the top left of the map, forming a stripe pattern of high values related to a highway construction in 2016 (see Appendix B.2). When looking more in detail the time series of these two places, the first change corresponds to an amplitude change (what seem to probably be the end of a recovery phase), while the second correspond to an abrupt change. This is in line with the results shown in Fig. 5, which indicate that BFAST Lite breakpoint magnitude is mostly sensitive to abrupt changes and in a lesser extent to amplitude changes. Other seasonal changes (NOS and LOS) are not expected to be well discriminated from the background, particularly if, as it can be observed in I.4, the real background values are higher than what was expected from the no-change category in the simulated data set (less than 0.1 vs. 0.1–0.3).

At the opposite, and as mentioned above, BFASTm-L2 change map (Figure I.3) based on the highest-breakpoint magnitude is able to highlight reported LSLAs such as those represented by polygons A (corresponding to the study case 3) and C (also see zoom-ins II.A and II. C), but also to detect newer LULC changes located in points B and D, and induced by agricultural activities as observed in the subplots II.B and II. D of Fig. 7. “Pattern” changes induced by activities other than agriculture are also detected, such as the ones related to the highway
PART II

<figure><img src="image_12.png" /><figcaption>Fig. 7. PART I) From left to right, top to bottom: 1) Snapshot view of the study area (Map data ©2015 Google), (2) MODIS NDVI average over 2000–2020, 3–4.5–6) Breakpoints (normalized) magnitude maps using: (3) BFASTm-L2, (4) BFAST Lite, (5) BFASTmonitor, (6) Edyn. Red polygons represent the active agricultural LSLAs reported in the ISRA ground-field database. PART II) Google Earth zoom-in of points A, B, C, D plotted in all the maps of PART I. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

ensure stability for model retraining.

In study case 1, two major events occur: a central-pivot irrigation system was installed at the end of 2011 (conversion from small to large-scale agriculture), and as from 2018 a huge increase in productivity is observed. The first change mostly translates as seasonal changes (NOS, amplitude), while the second induces a sharp abrupt change in 2018. In this case, BFAST Lite, BFASTm-L2, but also Edyn, output a higher-magnitude breakpoint for the abrupt change than the seasonal one, which is in line with the findings made in section 4.1.2 (Fig. 5).

Study case 2 (conversion from natural vegetation), is also concerned by the installation of a central-pivot irrigation system at the beginning of 2016. It translates in the signal as seasonal changes (NOS, amplitude), with a positive trend. In this case BFAST Lite and Edyn selected breakpoints are closely located in time (~1 year of difference), and both differs (>3 years) from BFASTm-L2 selected breakpoint. While BFAST Lite and Edyn breakpoints are closely located from big changes in amplitude, BFASTm-L2 highest magnitude breakpoint is located at the beginning of 2017, when seasonal changes are most marked, and most probably related to the concession implementation and start of activities.

Study case 3 is another example of transition from natural vegetation to intensive agriculture, but in a less arid environment. The transition is visually observed in 2008 and mostly related to seasonal changes (NOS). Despite the high variability exhibited in the time series (particularly after 2008), BFAST Lite and BFASTm-L2 detects the same breakpoint, while Edyn differs from 9 months away in the middle of the second growing cycle of the first year of production.

Lastly, study case 4 located in the most humid environment, shows more stability than the other study cases. Indeed, the land use change does not induce a significant change in the time series shape. The change being sought translates in a significant change in amplitude and a negative trend. Here again, BFAST Lite and BFASTm-L2 agreed in the breakpoint detected, inside the sought period of change. Edyn diverges from many years away, and outputs its highest-magnitude breakpoint where an amplitude change is observed. As a reminder, the Edyn breakpoint magnitude used here is the residual of the detected breakpoint. While the standardized magnitude should be more adapted for the comparison of breakpoint magnitudes, it was not used here as more than one breakpoint may share the same maximum value thus hindering the selection of a unique breakpoint. Edyn sensitivity to the parameter lambda was tested, with the parameter set to 0.6 (Appendix A). All the highest-magnitude breakpoints found differed from those presented in

Fig. 6, including for study case 1 which contains a huge abrupt change.

3.2.2. Detection of LSLAs induced LULC changes in breakpoint magnitude maps

So far, it has been shown that the performances of BFAST Lite and BFASTm-L2 in detecting in time seasonal changes are similar. However, because of the different sensitivities of the breakpoint magnitudes to the different types of change (Fig. 5), differences in the change maps based on the breakpoint magnitudes are expected. The change maps obtained with the different algorithms were evaluated over a study area close to Dakar, in which several LSLAs have been implemented since 2003 (Fig. 7).

The change maps based on the breakpoint magnitude display different spatial patterns. While LSLAs (red polygons) in the BFASTm-L2 (subplot I.3) change map are clearly highlighted and have significantly higher magnitude values than the background (highest difference for BFASTm-L2 as shown in Table 3), this is not the case for BFAST Lite (subplot I.4) which difference between the values inside and outside LSLAs is null.

BFAST Lite (subplot I.4) highest magnitudes are mainly located in two areas of the map: 1- in the centre of the map, in the wooded vegetation of the reserve of Bandia (see Appendix B.1) and, 2- in the top left of the map, forming a stripe pattern of high values related to a highway construction in 2016 (see Appendix B.2). When looking more in detail the time series of these two places, the first change corresponds to an amplitude change (what seem to probably be the end of a recovery phase), while the second correspond to an abrupt change. This is in line with the results shown in Fig. 5, which indicate that BFAST Lite breakpoint magnitude is mostly sensitive to abrupt changes and in a lesser extent to amplitude changes. Other seasonal changes (NOS and LOS) are not expected to be well discriminated from the background, particularly if, as it can be observed in I.4, the real background values are higher than what was expected from the no-change category in the simulated data set (less than 0.1 vs. 0.1–0.3).

At the opposite, and as mentioned above, BFASTm-L2 change map (Figure I.3) based on the highest-breakpoint magnitude is able to highlight reported LSLAs such as those represented by polygons A (corresponding to the study case 3) and C (also see zoom-ins II.A and II. C), but also to detect newer LULC changes located in points B and D, and induced by agricultural activities as observed in the subplots II.B and II. D of Fig. 7. “Pattern” changes induced by activities other than agriculture are also detected, such as the ones related to the highway
construction at the north-west of the map.

Regarding BFASTmonitor change map (I.5), one can see that the highest-breakpoint magnitudes are globally positively correlated to average NDVI (I.2). Even if the average of the magnitudes within LSLAs is high (0.6, the highest of the set; cf. Table 3), the magnitudes are not particularly sensitive to the changes induced by LSLAs. As an example, and in contrast with the BFASTm-L2 change map, magnitudes are overall low in points C and B. From Fig. 5, and considering that the correlation with NDVI is lower for the other algorithms, we can advance that the BFASTmonitor magnitudes are mostly responsive to amplitude changes (instead of abrupt changes which would also be highlighted by BFAST Lite), or as observed in the previous section, to trend changes.

Finally, one can see in subplot I.6 that LSLAs are somehow discriminated in the Edyn residual-based change map (I.6). This is corroborated by the average difference of the magnitudes within and outside LSLAs, which even if lower than the BFASTm-L2 one, is still non-negligible (0.18 vs. 0.27, cf. Table 3). Despite the fact that the breakpoint residuals are not adapted for the breakpoint selection approach applied in this study (as the standardized magnitudes would be), the algorithm proved to respond to seasonal changes (in particular to amplitude changes). The drawback of the algorithm, and in particular for large-scale applications, is however its high sensitivity to the lambda parameter (Saxena et al., 2018) and low speed. This last point will certainly be improved in the future through the implementation of the algorithm on cloud platforms.

### 3.2.3. Running times of the algorithms

To assess the algorithm speed on real time series, the four algorithms and the L2 distance alone were applied on varying-size areas: from 10x10 to 200x200 pixels. Results are presented in Fig. 8.

From Fig. 8, BFASTm-L2 appears to be very fast (16.4 min for a 200x200 pixels area), just after BFASTmonitor (2.0 min for a 200x200 pixels area, with a median of numbers of breakpoints of 37). The “force-brute” method (the continuous computation of L2 at a 3-month step) took 21.1 min for the same area size, followed by BFAST Lite (86.5 min), and finally Edyn (439.8 min). Worth is to note that the R implementation of Edyn was called from python using the rpy2 interface, which may slower the entire process. Faster python implementation of Edyn may be further tested (as the vey recent pyEWMACD available in github https://github.com/lewistrotter/pyEWMACD).

<figure><img src="image_12.png" /><figcaption>**Fig. 8.** Running time mean and standard-deviations (5 runs) over areas with varying-size, using 6-cores parallel processing (64 Go RAM) for the 5 change detection approaches.</figcaption></figure>

# 4. Discussion

## 4.1. BFASTm-L2, an efficient method to detect seasonal changes

The results obtained with the single-change simulated data set showed that BFASTm-L2 was the approach with the overall best performance in accurately detecting seasonal changes, thereby demonstrating the efficiency of combining BFASTmonitor and L2 distance for timely breakpoint selection (Fig. 4). The algorithm was particularly good at detecting, in a benchmark dataset with single changes, LOS changes (84.8%) and NOS changes (98%), but also changes in amplitude (76.8%) and break/trend changes (83.4%), which were however better detected by BFAST Lite (99%).

When tested over simulated time series with a unique change, all the algorithms showed increased sensitivity (i.e. higher breakpoint magnitude) to abrupt changes (Fig. 5). However, when considering the seasonal distributions alone in the violin plots of Fig. 5 BFASTm-L2 presented the highest magnitude means for LOS and NOS. This was also a goal of this study, as it enhances the probability of being able to spatially identify large-scale persistent seasonal changes.

When applied to real MODIS NDVI time series of individual case studies, BFAST Lite and BFASTm-L2 performed particularly well, producing similar performances in detecting the different types of change (Fig. 6). However, on the study area, BFAST Lite in contrast to BFASTm-L2 failed to spatially highlight the changes in LULC caused by the setting up of the different types of agro-industrial concessions (SLAS) (Fig. 7). This is mostly explained by the lower breakpoint magnitudes of BFAST Lite associated with seasonal changes (SLAS mostly seemed to induce changes in NOS, as shown in Fig. 6), but also because the real background values of BFAST Lite were higher than, the simulated data led us to expect (see Table 3). On the contrary, BFASTm-L2 and in a lesser extent Edyn, efficiently spatially capture LSLA driven changes. To better identify the type of change induced by LSLAs, the violin plots of Fig. 5 along with the change maps of Fig. 7 were helpful in this task. Considering the performance similarities in breakpoint detection between BFAST Lite and BFASTm-L2, and the responsiveness of BFAST Lite magnitudes to abrupt and amplitude changes, in case LSLAs induced changes were of any of these types (i.e. amplitude or abrupt), they would certainly be highlighted in the BFAST Lite’s change map. This is not the case. Agro-industrial LSLAs induced changes are therefore probably of seasonal type (NOS and/or LOS). This interpretation is strengthened by the different specific study cases presented in Fig. 6. In addition, because Edyn (which was also able to capture LSLAs changes) outputs breakpoint magnitudes for gradual changes of the same order as the no-change category, the probability that the LSLAs changes captured by BFASTm-L2 but also by Edyn are due to trend changes is low.

In terms of process speed, BFASTm-L2 produced the second-best performance after BFASTmonitor, despite its high rate of false positives, probably enhanced in this study by the absence of a “penalty” period applied after each detected breakpoint, and the non-automatic determination of a stable training period. Indeed, it is worth to remember that the BFASTm-L2 approach relies on the L2-distance, computed whenever a breakpoint is detected by BFASTmonitor. As a result, the speed of this approach depends on the number of breakpoints found in long time series. With a step frequency of three months, and a monitoring period of 14 years, BFASTm-L2 proved to be faster than the force-brute method represented by L2 computed continuously over areas of >100x100 pixels (Fig. 8). Worth is to notice that running times may be improved by changing the step frequency, by removing some of the detected breakpoints if too close one to the other, or by using a more optimized algorithm (that effectively does implement the automatic detection of the training period).

To resume, the main contributions brought to the change detection community through this study are:
construction at the north-west of the map.

Regarding BFASTmonitor change map (I.5), one can see that the highest-breakpoint magnitudes are globally positively correlated to average NDVI (I.2). Even if the average of the magnitudes within LSLAs is high (0.6, the highest of the set; cf. Table 3), the magnitudes are not particularly sensitive to the changes induced by LSLAs. As an example, and in contrast with the BFASTm-L2 change map, magnitudes are overall low in points C and B. From Fig. 5, and considering that the correlation with NDVI is lower for the other algorithms, we can advance that the BFASTmonitor magnitudes are mostly responsive to amplitude changes (instead of abrupt changes which would also be highlighted by BFAST Lite), or as observed in the previous section, to trend changes.

Finally, one can see in subplot I.6 that LSLAs are somehow discriminated in the Edyn residual-based change map (I.6). This is corroborated by the average difference of the magnitudes within and outside LSLAs, which even if lower than the BFASTm-L2 one, is still non-negligible (0.18 vs. 0.27, cf. Table 3). Despite the fact that the breakpoint residuals are not adapted for the breakpoint selection approach applied in this study (as the standardized magnitudes would be), the algorithm proved to respond to seasonal changes (in particular to amplitude changes). The drawback of the algorithm, and in particular for large-scale applications, is however its high sensitivity to the lambda parameter (Saxena et al., 2018) and low speed. This last point will certainly be improved in the future through the implementation of the algorithm on cloud platforms.

### 3.2.3. Running times of the algorithms

To assess the algorithm speed on real time series, the four algorithms and the L2 distance alone were applied on varying-size areas: from 10x10 to 200x200 pixels. Results are presented in Fig. 8.

From Fig. 8, BFASTm-L2 appears to be very fast (16.4 min for a 200x200 pixels area), just after BFASTmonitor (2.0 min for a 200x200 pixels area, with a median of numbers of breakpoints of 37). The “force-brute” method (the continuous computation of L2 at a 3-month step) took 21.1 min for the same area size, followed by BFAST Lite (86.5 min), and finally Edyn (439.8 min). Worth is to note that the R implementation of Edyn was called from python using the rpy2 interface, which may slower the entire process. Faster python implementation of Edyn may be further tested (as the vey recent pyEWMACD available in github https://github.com/lewistrotter/pyEWMACD).

<figure><img src="image_13.png" /><figcaption>**Fig. 8.** Running time mean and standard-deviations (5 runs) over areas with varying-size, using 6-cores parallel processing (64 Go RAM) for the 5 change detection approaches.</figcaption></figure>

# 4. Discussion

## 4.1. BFASTm-L2, an efficient method to detect seasonal changes

The results obtained with the single-change simulated data set showed that BFASTm-L2 was the approach with the overall best performance in accurately detecting seasonal changes, thereby demonstrating the efficiency of combining BFASTmonitor and L2 distance for timely breakpoint selection (Fig. 4). The algorithm was particularly good at detecting, in a benchmark dataset with single changes, LOS changes (84.8%) and NOS changes (98%), but also changes in amplitude (76.8%) and break/trend changes (83.4%), which were however better detected by BFAST Lite (99%).

When tested over simulated time series with a unique change, all the algorithms showed increased sensitivity (i.e. higher breakpoint magnitude) to abrupt changes (Fig. 5). However, when considering the seasonal distributions alone in the violin plots of Fig. 5 BFASTm-L2 presented the highest magnitude means for LOS and NOS. This was also a goal of this study, as it enhances the probability of being able to spatially identify large-scale persistent seasonal changes.

When applied to real MODIS NDVI time series of individual case studies, BFAST Lite and BFASTm-L2 performed particularly well, producing similar performances in detecting the different types of change (Fig. 6). However, on the study area, BFAST Lite in contrast to BFASTm-L2 failed to spatially highlight the changes in LULC caused by the setting up of the different types of agro-industrial concessions (SLAS) (Fig. 7). This is mostly explained by the lower breakpoint magnitudes of BFAST Lite associated with seasonal changes (SLAS mostly seemed to induce changes in NOS, as shown in Fig. 6), but also because the real background values of BFAST Lite were higher than, the simulated data led us to expect (see Table 3). On the contrary, BFASTm-L2 and in a lesser extent Edyn, efficiently spatially capture LSLA driven changes. To better identify the type of change induced by LSLAs, the violin plots of Fig. 5 along with the change maps of Fig. 7 were helpful in this task. Considering the performance similarities in breakpoint detection between BFAST Lite and BFASTm-L2, and the responsiveness of BFAST Lite magnitudes to abrupt and amplitude changes, in case LSLAs induced changes were of any of these types (i.e. amplitude or abrupt), they would certainly be highlighted in the BFAST Lite’s change map. This is not the case. Agro-industrial LSLAs induced changes are therefore probably of seasonal type (NOS and/or LOS). This interpretation is strengthened by the different specific study cases presented in Fig. 6. In addition, because Edyn (which was also able to capture LSLAs changes) outputs breakpoint magnitudes for gradual changes of the same order as the no-change category, the probability that the LSLAs changes captured by BFASTm-L2 but also by Edyn are due to trend changes is low.

In terms of process speed, BFASTm-L2 produced the second-best performance after BFASTmonitor, despite its high rate of false positives, probably enhanced in this study by the absence of a “penalty” period applied after each detected breakpoint, and the non-automatic determination of a stable training period. Indeed, it is worth to remember that the BFASTm-L2 approach relies on the L2-distance, computed whenever a breakpoint is detected by BFASTmonitor. As a result, the speed of this approach depends on the number of breakpoints found in long time series. With a step frequency of three months, and a monitoring period of 14 years, BFASTm-L2 proved to be faster than the force-brute method represented by L2 computed continuously over areas of >100x100 pixels (Fig. 8). Worth is to notice that running times may be improved by changing the step frequency, by removing some of the detected breakpoints if too close one to the other, or by using a more optimized algorithm (that effectively does implement the automatic detection of the training period).

To resume, the main contributions brought to the change detection community through this study are:
- The proposition of BFASTm-L2 as a change detection method faster than BFAST Lite, already known for its speed, making easier its application at larger scales;

- A change detection approach with breakpoint magnitudes more sensitive to seasonal changes. This allows to highlight agricultural-induced LULC changes, supporting the hypothesis that generic LULC changes are very often seasonal, and that more importance should be given to the detection of this type of change.

- More understandable BFASTmonitor, BFAST Lite and Edyn magnitude-based change detection maps. Because of the sensitivity analysis made to the different types of change, insights on the dynamics behind the changes observed in a change map are gained.

## 4.2. Recommendations for using BFASTm-L2 to detect and map LULC changes at regional/ national scale

While it is clear that changes in land use/land cover (LULC) involve changes in the composition of the vegetation which, in turn, result in changes in phenology and seasonality (e.g. a transition from grassland to croplands (Mardian et al., 2021)), the relationship between a change in land use and a change in seasonality is less straightforward. To check this link, we focussed on changes in the system of land use caused by agroindustrial LSLAs installation in the study area. Foreseen specific changes were transitions from natural vegetation or small-scale agriculture to large-scale agricultural systems. These changes were identified without the use of any threshold, mask, or any prior information on the type of land cover or on the direction of the change, meaning the method is suitable for use in unsupervised change detection pipelines. Most of the known LSLAs-related land use changes in the study area were spatially highlighted with BFASTm-L2 and in a lesser extent with Edyn, but not with BFAST Lite, suggesting that these specific LULC changes may occur without necessarily causing any abrupt change.

However, even if BFASTm-L2 produced promising results, some limitations should be mentioned. First, and as pointed out in the violin plots of Fig. 5, the breakpoint highest-magnitude based selection approach hinders the identification of seasonal changes when time series contain abrupt changes. Second, BFASTm-L2 aims to detect seasonal changes based on the use of the time signal canopy greenness (estimated through the NDVI) as a disturbance indicator. In greener environments such as in Casamance (Southern Senegal), where the difference in seasonality may be less marked, the method may be less sensitive to changes in land use. In such cases, other indicators more sensitive to vegetation biomass should be tested. Third, BFASTm-L2 is highly sensitive (i.e. high breakpoint magnitudes) to trends, which are often linked to vegetation recovery/degradation land processes. While trends induced by these land processes are frequently accompanied by changes in amplitude, other seasonal changes such as NOS changes are less expected. Because of this, and because amplitude changes do not profoundly affect the time series shape, shape metrics such as the Procrustes distance could be used in the future to minimize the detection of these type of changes. A last limitation is the inability of BFASTm-L2 to detect recent changes (earlier than the defined L2-w period, in our case 3 years).

When thinking of using this method with other sensors, it is useful to remember that BFASTm-L2 relies on time series that: 1- have a high temporal frequency sufficient to properly represent the phenology, 2- are long enough (8 years as minimum) and 3- are gap-free and smoothed, in order to minimise false detections. As such, applications with Sentinel are currently hampered because of the short temporal depth. On another hand, under tropical conditions and with frequent cloud coverage, highly temporal frequency time series are hardly obtained with Sentinel and Landsat. Finally, because of those sensors' higher spatial resolutions, running times would increase, thus hindering the application of BFASTm-L2 at larger scales. With MODIS data, BFASTm-L2 is able to detect other LULC changes occurring at large scale than those related to LSLAs, such as the changes related to mining and urbanization

processes.

Lastly, it is worth emphasising that the method proposed here is purely pixel-based, and that it needs to be completed by a spatial analysis, in order to better identify the drivers of changes in land use and to better interpret the changes. Because of the positive correlation between the NDVI and the breakpoint magnitudes, spatial analysis should be limited to areas that do not encompass many different ecoregions. Further research on spatial analysis at larger scale will be done.

# 5. Conclusions and perspectives

We developed a simple, automatic and rapid approach to select the breakpoint linked to the largest seasonal change in long and dense NDVI MODIS real time series with multiple breakpoints. The method, named BFASTm-L2, is based on the combined use of BFASTmonitor algorithm and the L2 euclidean distance for breakpoint selection, and was shown to accurately detect most of the single change types included in a sub-sample of the Awty-Carroll et al. (2019) benchmark set. Applied to a study area in Senegal using 20 years of MODIS satellite imagery the algorithm, through the spatialization of its single breakpoint magnitude, proved to be able to spatially identify LULC changes induced by the implementation of agro-industrial concessions in Senegal. This task, performed automatically without the need for any prior knowledge, is fit to be included in unsupervised pipelines to map and analyse generic LULC changes at regional scale. This was also possible because of the absence of any abrupt changes, supporting the hypothesis that generic LULC changes are very often seasonal, and that more importance should be given to the detection of this type of change.

To improve the detection of changes in the LULC at regional and national scales, an operational tool will be developed on a platform such as Google Earth Engine. This will enable supporting land monitoring initiatives such as the Land Matrix in detecting and monitoring anthropogenic changes such as those driven by LSLAs, for which much information remains to be gathered to help ground local teams.

## CRediT authorship contribution statement

**Yasmine Ngadi Scarpetta:** Conceptualization, Methodology, Data curation, Formal analysis, Validation, Visualization, Writing – original draft. **Valentine Lebourgeois:** Writing – review & editing, Supervision, Resources, Funding acquisition. **Anne-Elisabeth Laques:** Resources, Supervision. **Mohamadou Dieye:** Investigation, Data curation. **Jérémy Bourgoin:** Investigation, Data curation. **Agnès Begue:** Writing – review & editing, Supervision.

## Declaration of Competing Interest

The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: 'Yasmine Ngadi Scarpetta reports financial support was provided by French National Research Agency. Valentine Lebourgeois reports financial support was provided by French Space Agency.'

## Data availability

Data will be made available on request.

## Acknowledgments

The authors wish to thank the ISRA-BAME and the Land Matrix Initiative for data gathering and support. This work was supported by the French Space Agency (CNES) through the funding of the TOSCA-VISAGE project. Y. Ngadi received a fellowship from Montpellier University. This work was also supported by the French National Research Agency under the Investments for the Future Program #DigitAg referred as ANR-16-CONV-0004.
- The proposition of BFASTm-L2 as a change detection method faster than BFAST Lite, already known for its speed, making easier its application at larger scales;

- A change detection approach with breakpoint magnitudes more sensitive to seasonal changes. This allows to highlight agricultural-induced LULC changes, supporting the hypothesis that generic LULC changes are very often seasonal, and that more importance should be given to the detection of this type of change.

- More understandable BFASTmonitor, BFAST Lite and Edyn magnitude-based change detection maps. Because of the sensitivity analysis made to the different types of change, insights on the dynamics behind the changes observed in a change map are gained.

## 4.2. Recommendations for using BFASTm-L2 to detect and map LULC changes at regional/ national scale

While it is clear that changes in land use/land cover (LULC) involve changes in the composition of the vegetation which, in turn, result in changes in phenology and seasonality (e.g. a transition from grassland to croplands (Mardian et al., 2021)), the relationship between a change in land use and a change in seasonality is less straightforward. To check this link, we focussed on changes in the system of land use caused by agroindustrial LSLAs installation in the study area. Foreseen specific changes were transitions from natural vegetation or small-scale agriculture to large-scale agricultural systems. These changes were identified without the use of any threshold, mask, or any prior information on the type of land cover or on the direction of the change, meaning the method is suitable for use in unsupervised change detection pipelines. Most of the known LSLAs-related land use changes in the study area were spatially highlighted with BFASTm-L2 and in a lesser extent with Edyn, but not with BFAST Lite, suggesting that these specific LULC changes may occur without necessarily causing any abrupt change.

However, even if BFASTm-L2 produced promising results, some limitations should be mentioned. First, and as pointed out in the violin plots of Fig. 5, the breakpoint highest-magnitude based selection approach hinders the identification of seasonal changes when time series contain abrupt changes. Second, BFASTm-L2 aims to detect seasonal changes based on the use of the time signal canopy greenness (estimated through the NDVI) as a disturbance indicator. In greener environments such as in Casamance (Southern Senegal), where the difference in seasonality may be less marked, the method may be less sensitive to changes in land use. In such cases, other indicators more sensitive to vegetation biomass should be tested. Third, BFASTm-L2 is highly sensitive (i.e. high breakpoint magnitudes) to trends, which are often linked to vegetation recovery/degradation land processes. While trends induced by these land processes are frequently accompanied by changes in amplitude, other seasonal changes such as NOS changes are less expected. Because of this, and because amplitude changes do not profoundly affect the time series shape, shape metrics such as the Procrustes distance could be used in the future to minimize the detection of these type of changes. A last limitation is the inability of BFASTm-L2 to detect recent changes (earlier than the defined L2-w period, in our case 3 years).

When thinking of using this method with other sensors, it is useful to remember that BFASTm-L2 relies on time series that: 1- have a high temporal frequency sufficient to properly represent the phenology, 2- are long enough (8 years as minimum) and 3- are gap-free and smoothed, in order to minimise false detections. As such, applications with Sentinel are currently hampered because of the short temporal depth. On another hand, under tropical conditions and with frequent cloud coverage, highly temporal frequency time series are hardly obtained with Sentinel and Landsat. Finally, because of those sensors' higher spatial resolutions, running times would increase, thus hindering the application of BFASTm-L2 at larger scales. With MODIS data, BFASTm-L2 is able to detect other LULC changes occurring at large scale than those related to LSLAs, such as the changes related to mining and urbanization

processes.

Lastly, it is worth emphasising that the method proposed here is purely pixel-based, and that it needs to be completed by a spatial analysis, in order to better identify the drivers of changes in land use and to better interpret the changes. Because of the positive correlation between the NDVI and the breakpoint magnitudes, spatial analysis should be limited to areas that do not encompass many different ecoregions. Further research on spatial analysis at larger scale will be done.

# 5. Conclusions and perspectives

We developed a simple, automatic and rapid approach to select the breakpoint linked to the largest seasonal change in long and dense NDVI MODIS real time series with multiple breakpoints. The method, named BFASTm-L2, is based on the combined use of BFASTmonitor algorithm and the L2 euclidean distance for breakpoint selection, and was shown to accurately detect most of the single change types included in a sub-sample of the Awty-Carroll et al. (2019) benchmark set. Applied to a study area in Senegal using 20 years of MODIS satellite imagery the algorithm, through the spatialization of its single breakpoint magnitude, proved to be able to spatially identify LULC changes induced by the implementation of agro-industrial concessions in Senegal. This task, performed automatically without the need for any prior knowledge, is fit to be included in unsupervised pipelines to map and analyse generic LULC changes at regional scale. This was also possible because of the absence of any abrupt changes, supporting the hypothesis that generic LULC changes are very often seasonal, and that more importance should be given to the detection of this type of change.

To improve the detection of changes in the LULC at regional and national scales, an operational tool will be developed on a platform such as Google Earth Engine. This will enable supporting land monitoring initiatives such as the Land Matrix in detecting and monitoring anthropogenic changes such as those driven by LSLAs, for which much information remains to be gathered to help ground local teams.

## CRediT authorship contribution statement

**Yasmine Ngadi Scarpetta:** Conceptualization, Methodology, Data curation, Formal analysis, Validation, Visualization, Writing – original draft. **Valentine Lebourgeois:** Writing – review & editing, Supervision, Resources, Funding acquisition. **Anne-Elisabeth Laques:** Resources, Supervision. **Mohamadou Dieye:** Investigation, Data curation. **Jérémy Bourgoin:** Investigation, Data curation. **Agnès Begue:** Writing – review & editing, Supervision.

## Declaration of Competing Interest

The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: 'Yasmine Ngadi Scarpetta reports financial support was provided by French National Research Agency. Valentine Lebourgeois reports financial support was provided by French Space Agency.'

## Data availability

Data will be made available on request.

## Acknowledgments

The authors wish to thank the ISRA-BAME and the Land Matrix Initiative for data gathering and support. This work was supported by the French Space Agency (CNES) through the funding of the TOSCA-VISAGE project. Y. Ngadi received a fellowship from Montpellier University. This work was also supported by the French National Research Agency under the Investments for the Future Program #DigitAg referred as ANR-16-CONV-0004.
# Appendix A. Supplementary material

Supplementary data to this article can be found online at https://doi.org/10.1016/j.jag.2023.103379.

## References

Awty-Carroll, K., Bunting, P., Hardy, A., Bell, G., 2019. An evaluation and comparison of four dense time series change detection methods using simulated data. Remote Sens. (Basel) 11 (23), 2779.

Awty-Carroll, K., 2019. Simulated NDVI time series repository. Available from: <https://osf.io/taf9y/>.

Brooks, E., Wynne, R.H., Thomas, V.A., Blinn, C.E., Coulston, J.W., 2014. On-the-fly massively multitemporal change detection using statistical quality control charts and landsat data. IEEE Trans. Geosci. Remote Sensing 52 (6), 3316-3332.

Brooks, E., Yang, Z., Thomas, V., Wynne, R., 2017. Edyn: dynamic signaling of changes to forests using exponentially weighted moving average charts. Forests 8 (9), 304.

Browning, D.M., Maynard, J.J., Karl, J.W., Peters, D.C., 2017. Breaks in MODIS time series portend vegetation change: verification using long-term data in an arid grassland ecosystem. Ecol. Appl: Publ. Ecol. Soc. Am. 27 (5), 1677-1693.

Bullock, E.L., Woodcock, C.E., Holden, C.E., 2020. Improved change monitoring using an ensemble of time series algorithms. Remote Sens. Environ. 238, 111165.

Gao, Y., Solórzano, J.V., Quevedo, A., Loya-Carrillo, J.O., 2021. How BFAST trend and seasonal model components affect disturbance detection in tropical dry forest and temperate forest. Remote Sens. (Basel) 13 (11), 2033.

Ghaderpour, E., Vujadinovic, T., 2020. Change detection within remotely sensed satellite image time series via spectral analysis. Remote Sens. (Basel) 12 (23), 4001.

Gieseke, F., Rosca, S., Henriksen, T., Verbesselt, J., Oancea, C.E., 2020. Massively-parallel change detection for satellite time series data with missing values. In: 2020 IEEE 36th International Conference on Data Engineering (ICDE). 2020 IEEE 36th International Conference on Data Engineering (ICDE), Dallas, TX, USA. 20/04/2020 - 24/04/2020. IEEE, pp. 385-396.

Hamunyela, E., Rosca, S., Mirt, A., Engle, E., Herold, M., Gieseke, F., Verbesselt, J., 2020. Implementation of BFASTmonitor algorithm on google earth engine to support large-area and sub-annual change monitoring using earth observation data. Remote Sens. (Basel) 12 (18), 2953.

Hamunyela, E., 2017. Space-time monitoring of tropical forest changes using observations from multiple satellites. Ph.D. Thesis, Wageningen University and Research, Laboratory of Geo-information Science and Remote Sensing, Wageningen, The Netherlands, 202 pp.

Harding, A., Chamberlain, W., Anseeuw, W., Manco, G., Niassy, S., 2016. Large-scale land acquisitions profile. Senegal. https://landmatrix.org/documents/60/LM_Country Profile Senegal English.pdf. Accessed 1 June 2023.

Hentze, K., Thonfeld, F., Menz, G., 2017. Beyond trend analysis: How a modified breakpoint analysis enhances knowledge of agricultural production after Zimbabwe's fast track land reform. Int. J. Appl. Earth Obs. Geoinf. 62, 78-87.

Jamali, S., Jónsson, P., Eklundh, L., Ardó, J., Seaquist, J., 2015. Detecting changes in vegetation trends using time series segmentation. Remote Sens. Environ. 156, 182-195.

Kennedy, R.E., Yang, Z., Cohen, W.B., 2010. Detecting trends in forest disturbance and recovery using yearly Landsat time series: 1. LandTrendr — temporal segmentation algorithms. Remote Sens. Environ. 114 (12), 2897-2910.

Leenstra, M., Marcos, D., Bovolo, F., Tuia, D. Self-supervised pre-training enhances change detection in Sentinel-2 imagery. Part of the Lecture Notes in Computer Science book series (LNCS.

Lhermitte, S., Verbesselt, J., Verstraeten, W.W., Coppin, P., 2011. A comparison of time series similarity measures for classification and change detection of ecosystem dynamics. Remote Sens. Environ. 115 (12), 3129-3152.

Mardian, J., Berg, A., Daneshfar, B., 2021. Evaluating the temporal accuracy of grassland to cropland change detection using multitemporal image analysis. Remote Sens. Environ. 255, 112292.

Masiliunas, D., Tsendbazar, N.-E., Herold, M., Verbesselt, J., 2021. BFAST lite: a lightweight break detection method for time series analysis. Remote Sens. (Basel) 13 (16), 3308.

Meshkini, K., Bovolo, F., Bruzzone, L., 2021. An Unsupervised Change Detection Approach for Dense Satellite Image Time Series Using 3D CNN. In: 2021 IEEE International Geoscience and Remote Sensing Symposium IGARSS. IGARSS 2021 - 2021 IEEE International Geoscience and Remote Sensing Symposium, Brussels, Belgium. 11/07/2021 - 16/07/2021. IEEE, pp. 4336-4339.

Molinier, M., Miettinen, J., Ienco, D., Qiu, S., zhu, Z., 2021. Optical satellite image time series analysis for environment applications: from classical methods to deep learning and beyond. In: Change Detection and Image Time Series Analysis 2. John Wiley & Sons, Ltd, pp. 109-154.

Ochtyra, A., Marcinkowska-Ochtyra, A., Raczko, E., 2020. Threshold- and trend-based vegetation change monitoring algorithm based on the inter-annual multi-temporal normalized difference moisture index series: A case study of the Tatra Mountains. Remote Sensing of Environment 249, 112026.

Piou, C., Lebourgeois, V., Benahi, A.S., Bonnal, V., Jaavar, M.e.H., Lecoq, M., Vassal, J.-M., 2013. Coupling historical prospection data and a remotely-sensed vegetation index for the preventative control of Desert locusts. Appendix A. Basic Appl. Ecol. 14 (7), 593-604.

Saxena, R., Watson, L.T., Wynne, R.H., Brooks, E.B., Thomas, V.A., Zhiqiang, Y., Kennedy, R.E., 2018. Towards a polyalgorithm for land use change detection. ISPRS J. Photogramm. Remote Sens. 144, 217-234.

Setiawan, Y., Yoshino, K., 2012. Change detection in land-use and land-cover dynamics at a regional scale from MODIS time-series imagery. ISPRS Ann. Photogramm. Remote Sens. Spatial Inf. Sci. I-7, 243-248.

Setiawan, Y., Yoshino, K., 2014. Detecting land-use change from seasonal vegetation dynamics on regional scale with MODIS EVI 250-m time-series imagery. J. Land Use Sci. 9 (3), 304-330.

Tappan, G., Sall, M., Wood, E., Cushing, M., 2004. Ecoregions and land cover trends in Senegal. J. Arid Environ. 59 (3), 427-462.

Tuia, D., Roscher, R., Wegner, J.D., Jacobs, N., Zhu, X., Camps-Valls, G., 2021. Toward a collective agenda on AI for earth science data analysis. IEEE Geosci. Remote Sens. Mag. 9 (2), 88-104.

Verbesselt, J., Hyndman, R., Newnham, G., Culvenor, D., 2010a. Detecting trend and seasonal changes in satellite image time series. Remote Sens. Environ. 114 (1), 106-115.

Verbesselt, J., Hyndman, R., Zeileis, A., Culvenor, D., 2010b. Phenological change detection while accounting for abrupt and gradual trends in satellite image time series. Remote Sens. Environ. 114 (12), 2970-2980.

Verbesselt, J., Zeileis, A., Herold, M., 2012. Near real-time disturbance detection using satellite image time series. Remote Sens. Environ. 123, 98-108.

Verburg, P.H., van de Steeg, J., Veldkamp, A., Willemen, L., 2009. From land cover change to land function dynamics: a major challenge to improve land characterization. J. Environ. Manage. 90 (3), 1327-1335.

Yuan, Q., Shen, H., Li, T., Li, Z., Li, S., Jiang, Y., Xu, H., Tan, W., Yang, Q., Wang, J., Gao, J., Zhang, L., 2020. Deep learning in environmental remote sensing: Achievements and challenges. Remote Sens. Environ. 241, 111716.

Zhao, K., Wulder, M.A., Hu, T., Bright, R., Wu, Q., Qin, H., Li, Y., Toman, E., Mallick, B., Zhang, X., Brown, M., 2019. Detecting change-point, trend, and seasonality in satellite time series data to track abrupt changes and nonlinear dynamics: a Bayesian ensemble algorithm. Remote Sens. Environ. 232, 111181.

Zhu, Z., 2017. Change detection using landsat time series: A review of frequencies, preprocessing, algorithms, and applications. ISPRS J. Photogramm. Remote Sens. 130, 370-384.

Zhu, X.X., Tuia, D., Mou, L., Xia, G.-S., Zhang, L., Xu, F., Fraundorfer, F., 2017. Deep learning in remote sensing: a comprehensive review and list of resources. IEEE Geosci. Remote Sens. Mag. 5 (4), 8-36.

Zhu, Z., Woodcock, C.E., 2014. Continuous change detection and classification of land cover using all available Landsat data. Remote Sens. Environ. 144, 152-171.

Zhu, Z., Zhang, J., Yang, Z., Aljaddani, A.H., Cohen, W.B., Qiu, S., Zhou, C., 2020. Continuous monitoring of land disturbance based on Landsat time series. Remote Sens. Environ. 238, 111116.
# Appendix A. Supplementary material

Supplementary data to this article can be found online at https://doi.org/10.1016/j.jag.2023.103379.

## References

Awty-Carroll, K., Bunting, P., Hardy, A., Bell, G., 2019. An evaluation and comparison of four dense time series change detection methods using simulated data. Remote Sens. (Basel) 11 (23), 2779.

Awty-Carroll, K., 2019. Simulated NDVI time series repository. Available from: <https://osf.io/taf9y/>.

Brooks, E., Wynne, R.H., Thomas, V.A., Blinn, C.E., Coulston, J.W., 2014. On-the-fly massively multitemporal change detection using statistical quality control charts and landsat data. IEEE Trans. Geosci. Remote Sensing 52 (6), 3316-3332.

Brooks, E., Yang, Z., Thomas, V., Wynne, R., 2017. Edyn: dynamic signaling of changes to forests using exponentially weighted moving average charts. Forests 8 (9), 304.

Browning, D.M., Maynard, J.J., Karl, J.W., Peters, D.C., 2017. Breaks in MODIS time series portend vegetation change: verification using long-term data in an arid grassland ecosystem. Ecol. Appl: Publ. Ecol. Soc. Am. 27 (5), 1677-1693.

Bullock, E.L., Woodcock, C.E., Holden, C.E., 2020. Improved change monitoring using an ensemble of time series algorithms. Remote Sens. Environ. 238, 111165.

Gao, Y., Solórzano, J.V., Quevedo, A., Loya-Carrillo, J.O., 2021. How BFAST trend and seasonal model components affect disturbance detection in tropical dry forest and temperate forest. Remote Sens. (Basel) 13 (11), 2033.

Ghaderpour, E., Vujadinovic, T., 2020. Change detection within remotely sensed satellite image time series via spectral analysis. Remote Sens. (Basel) 12 (23), 4001.

Gieseke, F., Rosca, S., Henriksen, T., Verbesselt, J., Oancea, C.E., 2020. Massively-parallel change detection for satellite time series data with missing values. In: 2020 IEEE 36th International Conference on Data Engineering (ICDE). 2020 IEEE 36th International Conference on Data Engineering (ICDE), Dallas, TX, USA. 20/04/2020 - 24/04/2020. IEEE, pp. 385-396.

Hamunyela, E., Rosca, S., Mirt, A., Engle, E., Herold, M., Gieseke, F., Verbesselt, J., 2020. Implementation of BFASTmonitor algorithm on google earth engine to support large-area and sub-annual change monitoring using earth observation data. Remote Sens. (Basel) 12 (18), 2953.

Hamunyela, E., 2017. Space-time monitoring of tropical forest changes using observations from multiple satellites. Ph.D. Thesis, Wageningen University and Research, Laboratory of Geo-information Science and Remote Sensing, Wageningen, The Netherlands, 202 pp.

Harding, A., Chamberlain, W., Anseeuw, W., Manco, G., Niassy, S., 2016. Large-scale land acquisitions profile. Senegal. https://landmatrix.org/documents/60/LM_Country Profile Senegal English.pdf. Accessed 1 June 2023.

Hentze, K., Thonfeld, F., Menz, G., 2017. Beyond trend analysis: How a modified breakpoint analysis enhances knowledge of agricultural production after Zimbabwe's fast track land reform. Int. J. Appl. Earth Obs. Geoinf. 62, 78-87.

Jamali, S., Jónsson, P., Eklundh, L., Ardó, J., Seaquist, J., 2015. Detecting changes in vegetation trends using time series segmentation. Remote Sens. Environ. 156, 182-195.

Kennedy, R.E., Yang, Z., Cohen, W.B., 2010. Detecting trends in forest disturbance and recovery using yearly Landsat time series: 1. LandTrendr — temporal segmentation algorithms. Remote Sens. Environ. 114 (12), 2897-2910.

Leenstra, M., Marcos, D., Bovolo, F., Tuia, D. Self-supervised pre-training enhances change detection in Sentinel-2 imagery. Part of the Lecture Notes in Computer Science book series (LNCS.

Lhermitte, S., Verbesselt, J., Verstraeten, W.W., Coppin, P., 2011. A comparison of time series similarity measures for classification and change detection of ecosystem dynamics. Remote Sens. Environ. 115 (12), 3129-3152.

Mardian, J., Berg, A., Daneshfar, B., 2021. Evaluating the temporal accuracy of grassland to cropland change detection using multitemporal image analysis. Remote Sens. Environ. 255, 112292.

Masiliunas, D., Tsendbazar, N.-E., Herold, M., Verbesselt, J., 2021. BFAST lite: a lightweight break detection method for time series analysis. Remote Sens. (Basel) 13 (16), 3308.

Meshkini, K., Bovolo, F., Bruzzone, L., 2021. An Unsupervised Change Detection Approach for Dense Satellite Image Time Series Using 3D CNN. In: 2021 IEEE International Geoscience and Remote Sensing Symposium IGARSS. IGARSS 2021 - 2021 IEEE International Geoscience and Remote Sensing Symposium, Brussels, Belgium. 11/07/2021 - 16/07/2021. IEEE, pp. 4336-4339.

Molinier, M., Miettinen, J., Ienco, D., Qiu, S., zhu, Z., 2021. Optical satellite image time series analysis for environment applications: from classical methods to deep learning and beyond. In: Change Detection and Image Time Series Analysis 2. John Wiley & Sons, Ltd, pp. 109-154.

Ochtyra, A., Marcinkowska-Ochtyra, A., Raczko, E., 2020. Threshold- and trend-based vegetation change monitoring algorithm based on the inter-annual multi-temporal normalized difference moisture index series: A case study of the Tatra Mountains. Remote Sensing of Environment 249, 112026.

Piou, C., Lebourgeois, V., Benahi, A.S., Bonnal, V., Jaavar, M.e.H., Lecoq, M., Vassal, J.-M., 2013. Coupling historical prospection data and a remotely-sensed vegetation index for the preventative control of Desert locusts. Appendix A. Basic Appl. Ecol. 14 (7), 593-604.

Saxena, R., Watson, L.T., Wynne, R.H., Brooks, E.B., Thomas, V.A., Zhiqiang, Y., Kennedy, R.E., 2018. Towards a polyalgorithm for land use change detection. ISPRS J. Photogramm. Remote Sens. 144, 217-234.

Setiawan, Y., Yoshino, K., 2012. Change detection in land-use and land-cover dynamics at a regional scale from MODIS time-series imagery. ISPRS Ann. Photogramm. Remote Sens. Spatial Inf. Sci. I-7, 243-248.

Setiawan, Y., Yoshino, K., 2014. Detecting land-use change from seasonal vegetation dynamics on regional scale with MODIS EVI 250-m time-series imagery. J. Land Use Sci. 9 (3), 304-330.

Tappan, G., Sall, M., Wood, E., Cushing, M., 2004. Ecoregions and land cover trends in Senegal. J. Arid Environ. 59 (3), 427-462.

Tuia, D., Roscher, R., Wegner, J.D., Jacobs, N., Zhu, X., Camps-Valls, G., 2021. Toward a collective agenda on AI for earth science data analysis. IEEE Geosci. Remote Sens. Mag. 9 (2), 88-104.

Verbesselt, J., Hyndman, R., Newnham, G., Culvenor, D., 2010a. Detecting trend and seasonal changes in satellite image time series. Remote Sens. Environ. 114 (1), 106-115.

Verbesselt, J., Hyndman, R., Zeileis, A., Culvenor, D., 2010b. Phenological change detection while accounting for abrupt and gradual trends in satellite image time series. Remote Sens. Environ. 114 (12), 2970-2980.

Verbesselt, J., Zeileis, A., Herold, M., 2012. Near real-time disturbance detection using satellite image time series. Remote Sens. Environ. 123, 98-108.

Verburg, P.H., van de Steeg, J., Veldkamp, A., Willemen, L., 2009. From land cover change to land function dynamics: a major challenge to improve land characterization. J. Environ. Manage. 90 (3), 1327-1335.

Yuan, Q., Shen, H., Li, T., Li, Z., Li, S., Jiang, Y., Xu, H., Tan, W., Yang, Q., Wang, J., Gao, J., Zhang, L., 2020. Deep learning in environmental remote sensing: Achievements and challenges. Remote Sens. Environ. 241, 111716.

Zhao, K., Wulder, M.A., Hu, T., Bright, R., Wu, Q., Qin, H., Li, Y., Toman, E., Mallick, B., Zhang, X., Brown, M., 2019. Detecting change-point, trend, and seasonality in satellite time series data to track abrupt changes and nonlinear dynamics: a Bayesian ensemble algorithm. Remote Sens. Environ. 232, 111181.

Zhu, Z., 2017. Change detection using landsat time series: A review of frequencies, preprocessing, algorithms, and applications. ISPRS J. Photogramm. Remote Sens. 130, 370-384.

Zhu, X.X., Tuia, D., Mou, L., Xia, G.-S., Zhang, L., Xu, F., Fraundorfer, F., 2017. Deep learning in remote sensing: a comprehensive review and list of resources. IEEE Geosci. Remote Sens. Mag. 5 (4), 8-36.

Zhu, Z., Woodcock, C.E., 2014. Continuous change detection and classification of land cover using all available Landsat data. Remote Sens. Environ. 144, 152-171.

Zhu, Z., Zhang, J., Yang, Z., Aljaddani, A.H., Cohen, W.B., Qiu, S., Zhou, C., 2020. Continuous monitoring of land disturbance based on Landsat time series. Remote Sens. Environ. 238, 111116.