Monitoring Deforestation at Sub-Annual Scales as Extreme Events in Landsat Data Cubes

Current methods for monitoring deforestation from satellite data at sub-annual scales require pixel time series to have many historical observations in the reference period to model normal forest dynamics before detecting deforestation. However, in some areas, pixel time series often do not have many historical observations. Detecting deforestation at a pixel with scarce historical observations can be improved by complementing the pixel time series with spatial context information. In this work, we propose a data-driven space-time change detection method that detects deforestation events at sub-annual scales in data cubes of satellite image time series. First we spatially normalised observations in the local space-time data cube to reduce seasonality. Subsequently, we detected deforestation by assessing whether a newly acquired observation in the monitoring period is an extreme when compared against spatially normalised values in a local space-time data cube defined over reference period. We demonstrated our method at two sites, a dry tropical Bolivian forest and a humid tropical Brazilian forest, by varying the spatial and temporal extent of data cube. We emulated a “near real-time” monitoring scenario, implying that observations in the monitoring period were sequentially rather than simultaneously assessed for deforestation. Using Landsat normalised difference vegetation index (NDVI) time series, we achieved a median temporal detection delay of less than three observations, a producer’s accuracy above 70%, a user’s accuracy above 65%, and an overall accuracy above 80% at both sites, even when the reference period of the data cube only contained one year of data. Our results also show that large percentile thresholds (e.g., 5th percentile) achieve higher producer’s accuracy and shorter temporal detection delay, whereas smaller percentiles (e.g., 0.1 percentile) achieve higher user’s accuracy, but longer temporal detection delay. The method is data-driven, not based on statistical assumption on the data distribution, and can be applied on different forest types. However, it may face challenges in mixed forests where, for example, deciduous and evergreen forests coexist within short distances. A pixel to be assessed for deforestation should have a minimum of three temporal observations, the first of which must be known to represent forest. Such short time series allow rapid deployment of newly launched sensors (e.g., Sentinel-2) for detecting deforestation events at sub-annual scales.


Introduction
Monitoring deforestation at sub-annual scales (e.g., weekly or monthly) using satellite data is increasingly becoming an important part of the initiatives that aim to reduce deforestation across the globe.This is because monitoring deforestation at the sub-annual scale, unlike annual monitoring, allows for timely detection of deforestation events, thus providing an opportunity for early interventions to stop illegal deforestation activities [1].For example, the Brazilian Institute for Space Research (INPE) monitors deforestation events at sub-annual scales in the Amazon using a system based on Moderate Resolution Imaging Spectroradiometer and bi-temporal change detection which has played an important role in reducing deforestation in Brazil.However, forest monitoring systems which detect deforestation events at sub-annual scales based on a bi-temporal change detection approach may face challenges in areas where forest has strong seasonality.To address this challenge, methods that detect deforestation at sub-annual scales from satellite image time series while accounting for seasonal variations have been developed in recent years [2][3][4][5][6][7][8].These methods detect deforestation events by testing if a newly acquired observation at a particular pixel is abnormally low when compared to historical temporal dynamics of forest at such pixel [2,3,6,9].However, for the test to be robust, a pixel time series is required to have many historical observations.In some areas, however, pixel time series often do not have enough historical observations, mainly because of persistent cloud cover coupled with a relatively long revisit time in the past [10], especially for satellite sensors which have high and medium spatial resolutions.To remedy the problem of cloud cover, new methods that detect deforestation events at sub-annual scales by combining optical and synthetic aperture radar (SAR) data have been proposed [9,11].However, these methods also require a pixel to have many historical observations.In the past, SAR sensors, which can penetrate the clouds, had limited temporal coverage, also resulting in sparse time series.Such sparse historical observations make it difficult to properly model the normal temporal dynamics of the forest, and may lead to many false detections during deforestation detection.In the near future, however, temporally dense SAR time series from Sentinel sensors will be available, and such dense time series, when combined with data from optical sensors, will address some of the current challenges associated with monitoring deforestation at sub-annual scales in the tropics [12].In particular, it will be possible to detect deforestation events within a few days of occurrence.However, with such dense time series, another challenge related to pre-processing huge amounts of historical data will arise: pre-processing huge amounts of historical data for large areas is likely to take a relatively long time, thus affecting rapid detection of deforestation events.
The challenge associated with monitoring deforestation at sub-annual scales in areas with scarce historical observations can also be addressed by exploiting both temporal and spatial information in satellite image time series.An individual pixel may not have enough observations, but the spatial context derived from neighbouring pixels can provide sufficient information to determine whether a pixel with scarce historical observations is deforested or not.Such an approach may reduce the amount of historical satellite data that needs to be processed when monitoring deforestation at sub-annual scales.
Recently, Huang and colleagues [13] exploited spatiotemporal Landsat data to identify deforested areas, but that work focused on detecting deforestation at the annual scale.Similarly, Hamunyela and co-workers [5] used the spatial context of pixels to reduce seasonal variations in Landsat time series before detecting deforestation at sub-annual scales using a method [3] which relies on individual pixel time series.These studies demonstrate that spatiotemporal information of image time series is useful for deforestation monitoring.However, methods for detecting deforestation at sub-annual scales based on integrated analysis of spatiotemporal data have not been published, to the best of our knowledge.
The approach proposed in this paper is to identify deforested pixels by exploiting spatiotemporal information available in the space-time data cube of satellite image time series.In this way, a pixel with at least three temporal observations in its time series is expected to still allow assessment of deforestation at sub-annual scales.This is because observations in the space-time data cube are likely to be sufficient for deciding whether a newly acquired observation is abnormally low.However, shifting to space-time change analysis requires two major challenges to be overcome.The first challenge concerns dealing with seasonality in the data cube.Seasonal variations may disguise deforestation in the data if not removed, especially in dry forests where seasonality is strong.The second challenge is how to consistently identify anomalies in a space-time data cube.Existing methods for identifying anomalous observations in space-time data cubes of climatic data [14] and global gross primary production [15] nevertheless focus on the temporal perspective, since thresholds are defined without considering spatial context.
This paper describes a data-driven space-time change detection method for monitoring deforestation at sub-annual scales.The method detects deforestation at pixel level as an extreme event in vegetation index values within local space-time cubes of satellite image time series.With this method, spatial-temporal information in satellite image time series is exploited to detect deforestation at sub-annual scales.Our method builds upon the spatial context approach developed in a previous study [5] for reducing seasonal variations in satellite image time series.The method we propose here demonstrates how data cubes of satellite image time series can be used for sub-annual deforestation monitoring.We demonstrated our method at two sites, a dry tropical forest and a humid tropical forest (Section 3), using normalised difference vegetation index (NDVI, [16,17]) time series derived from Landsat-5/TM and Landsat-7/ETM+ data.

Space-Time Approach for Deforestation Detection-The Concept
Methods for detecting deforestation events at sub-annual scales (e.g., monthly) from satellite image time series first model the normal forest cover dynamics before assessing deforestation.Modelling normal forest cover dynamics is challenging for pixels which do not have many historical observations.We can remedy this challenge of insufficient historical observations in individual pixel time series by adding the spatial context of the focal pixel in a so-called neighbouring pixels space-time data cube.To do this, we implement the following steps (Figure 1).First, like in other studies including [4,5,18] we mask non-forest areas in the historical period from the image time series.In this study, we masked non-forest areas using Landsat tree cover continuous fields [19].Second, we define a local space-time data cube (e.g., Figure 2) around each pixel in the image time series.Essentially, we have a spatially moving window with an additional time dimension that moves from one pixel to the other.To avoid the smoothing out effect when deseasonalising the data [5], the spatial extent of the data cube should be larger than the size of deforestation events the user aims to detect.Third, we deseasonalise the observations in the local data cube (e.g., Figure 3) to increase the sensitivity for detecting non-seasonal changes (see Section 2.2).Fourth, we split each local data cube into a reference cube (RC) and monitoring cube (MC).The RC contains historical observations-where non-forest pixels have been masked, whereas the MC contains newly acquired observations-not yet assessed for deforestation.Here, we build upon the idea of specifying the history (reference) and monitoring period as employed in a structural change monitoring framework [3,[20][21][22].Fifth, based on training data and observations in the RC, we compute a threshold percentile for defining an observation as abnormally low (see Section 2.3).Sixth, we assess whether the focal pixel in the MC is below the threshold percentile.If so, the pixel is flagged as deforestation or potential deforestation.
global gross primary production [15] nevertheless focus on the temporal perspective, since thresholds are defined without considering spatial context.
This paper describes a data-driven space-time change detection method for monitoring deforestation at sub-annual scales.The method detects deforestation at pixel level as an extreme event in vegetation index values within local space-time cubes of satellite image time series.With this method, spatial-temporal information in satellite image time series is exploited to detect deforestation at sub-annual scales.Our method builds upon the spatial context approach developed in a previous study [5] for reducing seasonal variations in satellite image time series.The method we propose here demonstrates how data cubes of satellite image time series can be used for sub-annual deforestation monitoring.We demonstrated our method at two sites, a dry tropical forest and a humid tropical forest (Section 3), using normalised difference vegetation index (NDVI, [16,17]) time series derived from Landsat-5/TM and Landsat-7/ETM+ data.

Space-Time Approach for Deforestation Detection-The Concept
Methods for detecting deforestation events at sub-annual scales (e.g., monthly) from satellite image time series first model the normal forest cover dynamics before assessing deforestation.Modelling normal forest cover dynamics is challenging for pixels which do not have many historical observations.We can remedy this challenge of insufficient historical observations in individual pixel time series by adding the spatial context of the focal pixel in a so-called neighbouring pixels space-time data cube.To do this, we implement the following steps (Figure 1).First, like in other studies including [4,5,18] we mask non-forest areas in the historical period from the image time series.In this study, we masked non-forest areas using Landsat tree cover continuous fields [19].Second, we define a local space-time data cube (e.g., Figure 2) around each pixel in the image time series.Essentially, we have a spatially moving window with an additional time dimension that moves from one pixel to the other.To avoid the smoothing out effect when deseasonalising the data [5], the spatial extent of the data cube should be larger than the size of deforestation events the user aims to detect.Third, we deseasonalise the observations in the local data cube (e.g., Figure 3) to increase the sensitivity for detecting non-seasonal changes (see Section 2.2).Fourth, we split each local data cube into a reference cube (RC) and monitoring cube (MC).The RC contains historical observations-where non-forest pixels have been masked, whereas the MC contains newly acquired observations-not yet assessed for deforestation.Here, we build upon the idea of specifying the history (reference) and monitoring period as employed in a structural change monitoring framework [3,[20][21][22].Fifth, based on training data and observations in the RC, we compute a threshold percentile for defining an observation as abnormally low (see Section 2.3).Sixth, we assess whether the focal pixel in the MC is below the threshold percentile.If so, the pixel is flagged as deforestation or potential deforestation.

Deseasonalising Observations in the Space-Time Data Cube
Building on the method of [5], NDVI time series are deseasonalised by dividing the NDVI value for each pixel with the 95th percentile (P95) computed from the NDVI values of the pixels within a defined spatial neighbourhood of such pixel.The use of the P95 is based on the assumption that, for each spatial neighbourhood, forested pixels would occupy the upper tail of the distribution, which is deemed to be properly represented by P95 [5].Note that this assumption may not hold if the forest canopy cover is low, or if there are pixels with high vegetation index (VI) values owing to land cover types other than forests (e.g., agriculture).A good forest mask is needed to make sure seasonality for forest pixels can be captured by P95.In contrast to [5], where P95 was computed for each pixel, here we compute P95 at each time step in the local data cube, and use it to deseasonalise all pixels per time slice in the data cube.The subscript t in P95t denotes a particular time instance t.Within the data cube,

Deseasonalising Observations in the Space-Time Data Cube
Building on the method of [5], NDVI time series are deseasonalised by dividing the NDVI value for each pixel with the 95th percentile (P95) computed from the NDVI values of the pixels within a defined spatial neighbourhood of such pixel.The use of the P95 is based on the assumption that, for each spatial neighbourhood, forested pixels would occupy the upper tail of the distribution, which is deemed to be properly represented by P95 [5].Note that this assumption may not hold if the forest canopy cover is low, or if there are pixels with high vegetation index (VI) values owing to land cover types other than forests (e.g., agriculture).A good forest mask is needed to make sure seasonality for forest pixels can be captured by P95.In contrast to [5], where P95 was computed for each pixel, here we compute P95 at each time step in the local data cube, and use it to deseasonalise all pixels per time slice in the data cube.The subscript t in P95t denotes a particular time instance t.Within the data cube,

Deseasonalising Observations in the Space-Time Data Cube
Building on the method of [5], NDVI time series are deseasonalised by dividing the NDVI value for each pixel with the 95th percentile (P 95 ) computed from the NDVI values of the pixels within a defined spatial neighbourhood of such pixel.The use of the P 95 is based on the assumption that, for each spatial neighbourhood, forested pixels would occupy the upper tail of the distribution, which is deemed to be properly represented by P 95 [5].Note that this assumption may not hold if the forest canopy cover is low, or if there are pixels with high vegetation index (VI) values owing to land cover types other than forests (e.g., agriculture).A good forest mask is needed to make sure seasonality for forest pixels can be captured by P 95 .In contrast to [5], where P 95 was computed for each pixel, here we compute P 95 at each time step in the local data cube, and use it to deseasonalise all pixels per time slice in the data cube.The subscript t in P 95t denotes a particular time instance t.Within the data cube, we calculate P 95t for each time instance, and deseasonalise all pixel values at such time instances by dividing them by P 95t .An example of how observations in the reference cube are distributed (a) before and (b) after reducing the seasonality in a space-time data cube using spatial context is shown in Figure 3.After deseasonalising, the next step (Section 2.3) is to determine whether a newly acquired observation at the focal pixel in the monitoring period is abnormally low.

Detecting Deforestation as Extreme Event in Space-Time Data Cube
We rely on the existing extreme value approach [15] to compute the threshold from RC for determining if an observation in the monitoring period is abnormally low.With the extreme value approach, observations below a specified percentile (e.g., 5th percentile) are regarded as abnormally low [15].Here, a suitable percentile is established by using a training dataset (see Section 3.2).
For each local cube, an observation is abnormally low if it is below the specified percentile, and such an observation is initially flagged as potential deforestation.Definite labelling is postponed since satellite image time series are often noisy, mainly because of remnants of clouds and cloud shadows.Furthermore, some pixels may have naturally low values if they cover forest areas that have lower photosynthetic activity.Values for such pixels might be identified as extreme observations within a local space-time data cube, although no deforestation has occurred.If not accounted for, these pixels could lead to the detection of false deforestation events.
Deforestation is confirmed if the next observation is also below the specified percentile regardless of the size of the time step between the observations.Given this monitoring approach, deforestation events are detected with a delay of one observation.This delay is referred to as temporal detection delay, which is the number of observations between when the deforestation event occurred and when the deforestation is detected [5].We only assess a pixel for deforestation if such a pixel contains a minimum of three observations, one of which must be in the RC.

Monitoring Deforestation in Landsat NDVI Data Cubes-A Case Study
We assessed our space-time change detection method at two study sites, shown in Figure 4.One study site is a dry tropical forest located southeast of Santa Cruz de la Sierra in Bolivia, (centred at 18.388 ˝S, 62.361 ˝W), and the other study site is a humid tropical forest located west of Ariquemes, Rondonia State, Brazil (centred at 10.2952 ˝S, 64.0478 ˝W).The forest at the Bolivian site is characterised by strong seasonality in its photosynthetic activity, whereas the seasonality is less pronounced at the Brazilian site.Each of the study sites covers an area of about 10,000 km 2 .Deforestation at the Bolivian site is dominated mainly by industrial agricultural expansion that resulted in large blocks of deforestation events, whereas deforestation events at the Brazilian site are heterogeneous in size, corresponding mostly to a process of colonisation.With varying degrees of seasonality and different deforestation processes, these study sites are particularly suitable for testing a new method for sub-annual deforestation monitoring.We used the NDVI image time series derived from atmospherically [23] and geometrically corrected Landsat-5/TM and Landsat-7/ETM+ images.Landsat images were obtained from the United State of America's Geological Survey (USGS) Landsat Surface Reflectance (SR) Climate Records (CDR).We used all available (2000-2014) terrain corrected images (L1T).We assumed that co-registration of Landsat-5/TM and Landsat-7/ETM+ images was satisfactory.Clouds and cloud shadows were masked using the Fmask procedure [24].Note that Fmask outputs are distributed with Landsat SR CDR Products.Landsat tree cover continuous fields for 2005 [19] were used to mask non-forest areas and areas with less than 10% tree cover prior to the year 2005.To understand how the spatial extent of the local data cube influences spatial and temporal accuracies, we tested six varying spatial extents of the local data cube (Table 1).We also varied the temporal extent of the RC, herein referred to as temporal extent of the data cube, to understand how it affects spatial and temporal accuracy for deforestation detection.The temporal extent was varied from one to five years of data, at an interval of one year.The RC contained images from 2004 for the one-year data scenario, images from 2003-2004 for the two-year scenario, and images from 2002-2004 for the three-year data scenario.For the four-and five-year data scenarios, the RC contained images from 2001-2004 and 2000-2004, respectively.Figure 5 shows the number of images available in the RC for each temporal extent at each study site.
For each spatial and temporal extent, we trained our method to determine the optimal percentile for detecting deforestation at a sub-annual scale (Section 3.2).Next, we used the optimal percentiles to validate our method (Section 3.3).To understand how the spatial extent of the local data cube influences spatial and temporal accuracies, we tested six varying spatial extents of the local data cube (Table 1).We also varied the temporal extent of the RC, herein referred to as temporal extent of the data cube, to understand how it affects spatial and temporal accuracy for deforestation detection.The temporal extent was varied from one to five years of data, at an interval of one year.The RC contained images from 2004 for the one-year data scenario, images from 2003-2004 for the two-year scenario, and images from 2002-2004 for the three-year data scenario.For the four-and five-year data scenarios, the RC contained images from 2001-2004 and 2000-2004, respectively.Figure 5 shows the number of images available in the RC for each temporal extent at each study site.
For each spatial and temporal extent, we trained our method to determine the optimal percentile for detecting deforestation at a sub-annual scale (Section 3.2).Next, we used the optimal percentiles to validate our method (Section 3.3).

Reference Data
We used 966 sample pixels to validate the final change map for the Bolivian site, and 400 sample pixels were used for validating the Brazilian change map (Table 2).Training was based on 170 and 70 sample pixels for the sites in Bolivia and Brazil, respectively.Similar to [4,18,25], reference data were acquired by visual interpretation of Landsat data along with high spatial resolution imagery available in Google Earth and Bing Maps.High spatial resolution imagery available in Google Earth and Bing Maps were used to determine whether an area is indeed deforested or not.At each study site, we sampled deforested and forested areas by stratified probability sampling [26,27] after manually digitising corresponding areas on Landsat images.The number of sample pixels was proportional to the area of the stratum.The deforested stratum contained areas that had been deforested during the period of 2005-2014, whereas the forested stratum covered areas which were still forested at the end of 2014.For each sample pixel in the deforested area, we estimated the date of deforestation by visually determining the Landsat image in which the deforestation event is first visible.The date of deforestation was used to assess the temporal accuracy.

Training the Space-Time Change Detection Method
To train our method, we generated a series of percentiles (n = 50), ranging from 0.1 to the 5th percentile, at an interval of 0.1 percent.Next, for each spatial and temporal extent, we used each of these percentiles as a threshold for identifying deforestation at a sub-annual scale.A training data set (Table 2) was then used to calculate the overall accuracy, bias (calculated by subtracting the omission error from the commission error [11]), and the median temporal detection delay.Since our monitoring goal was to detect deforestation events as early as possible but with high overall accuracy,

Reference Data
We used 966 sample pixels to validate the final change map for the Bolivian site, and 400 sample pixels were used for validating the Brazilian change map (Table 2).Training was based on 170 and 70 sample pixels for the sites in Bolivia and Brazil, respectively.Similar to [4,18,25], reference data were acquired by visual interpretation of Landsat data along with high spatial resolution imagery available in Google Earth and Bing Maps.High spatial resolution imagery available in Google Earth and Bing Maps were used to determine whether an area is indeed deforested or not.At each study site, we sampled deforested and forested areas by stratified probability sampling [26,27] after manually digitising corresponding areas on Landsat images.The number of sample pixels was proportional to the area of the stratum.The deforested stratum contained areas that had been deforested during the period of 2005-2014, whereas the forested stratum covered areas which were still forested at the end of 2014.For each sample pixel in the deforested area, we estimated the date of deforestation by visually determining the Landsat image in which the deforestation event is first visible.The date of deforestation was used to assess the temporal accuracy.

Training the Space-Time Change Detection Method
To train our method, we generated a series of percentiles (n = 50), ranging from 0.1 to the 5th percentile, at an interval of 0.1 percent.Next, for each spatial and temporal extent, we used each of these percentiles as a threshold for identifying deforestation at a sub-annual scale.A training data set (Table 2) was then used to calculate the overall accuracy, bias (calculated by subtracting the omission error from the commission error [11]), and the median temporal detection delay.Since our monitoring goal was to detect deforestation events as early as possible but with high overall accuracy, for each spatial and temporal extent we selected the percentile, i.e., the optimal percentile, that achieved the shortest median temporal detection delay with the highest overall accuracy.

Validating the Space-Time Change Detection Method
At each study site, we applied our space-time change detection method using the optimal percentiles, determined from the training data, as thresholds for deforestation.We emulated a "near real-time" monitoring scenario, implying that observations in the monitoring period were sequentially rather than simultaneously assessed for deforestation.Although some areas may have experienced multiple deforestation and regrowth regimes between 2005 and 2014, we only considered the first deforestation event per pixel, and once labelled as deforested, we stopped monitoring such a pixel at subsequent time steps.The spatial and temporal accuracy for change detected between 2005 and 2014 were calculated using the test data set (Table 2).More specifically, we calculated the overall accuracy, producer's accuracy, user's accuracy, and the median temporal detection delay.Like [5], the temporal detection delay was calculated at each sample point by counting the number of valid observations available between the image in which deforestation was visually identified, and the image in which deforestation was detected by our method.Confidence intervals for the overall accuracy, as well as producer's and user's accuracies, were calculated using binomial probability of success based on Wilson's method [28].

Results
In this section, first we show the results from the training phase of our method (Section 4.1).In particular, we show how the percentile threshold affects the spatial and temporal accuracy (Section 4.1.1),and how spatial and temporal extents of the data cube influence the optimal percentile for detecting deforestation at a sub-annual scale.Next, we show the validation results (Section 4.2).

Effect of Percentile Threshold on Spatial and Temporal Accuracy
Figure 6 shows an example of how overall accuracy, bias, and median temporal detection delay were changing in relation to the percentile threshold.Here, we only show examples for the smallest (C5) and largest (C25) cubes for each study site.Generally, larger percentiles produced lower overall accuracy (Figure 6a,d) and shorter median temporal detection delay (Figure 6c,f), except for the C25 at Brazilian site.Deforestation events were increasingly being omitted (negative bias) when using smaller percentiles (Figure 6b,e).for each spatial and temporal extent we selected the percentile, i.e., the optimal percentile, that achieved the shortest median temporal detection delay with the highest overall accuracy.

Validating the Space-Time Change Detection Method
At each study site, we applied our space-time change detection method using the optimal percentiles, determined from the training data, as thresholds for deforestation.We emulated a "near real-time" monitoring scenario, implying that observations in the monitoring period were sequentially rather than simultaneously assessed for deforestation.Although some areas may have experienced multiple deforestation and regrowth regimes between 2005 and 2014, we only considered the first deforestation event per pixel, and once labelled as deforested, we stopped monitoring such a pixel at subsequent time steps.The spatial and temporal accuracy for change detected between 2005 and 2014 were calculated using the test data set (Table 2).More specifically, we calculated the overall accuracy, producer's accuracy, user's accuracy, and the median temporal detection delay.Like [5], the temporal detection delay was calculated at each sample point by counting the number of valid observations available between the image in which deforestation was visually identified, and the image in which deforestation was detected by our method.Confidence intervals for the overall accuracy, as well as producer's and user's accuracies, were calculated using binomial probability of success based on Wilson's method [28].

Results
In this section, first we show the results from the training phase of our method (Section 4.1).In particular, we show how the percentile threshold affects the spatial and temporal accuracy (Section 4.1.1),and how spatial and temporal extents of the data cube influence the optimal percentile for detecting deforestation at a sub-annual scale.Next, we show the validation results (Section 4.2).

Effect of Percentile Threshold on Spatial and Temporal Accuracy
Figure 6 shows an example of how overall accuracy, bias, and median temporal detection delay were changing in relation to the percentile threshold.Here, we only show examples for the smallest (C5) and largest (C25) cubes for each study site.Generally, larger percentiles produced lower overall accuracy (Figure 6a,d) and shorter median temporal detection delay (Figure 6c,f), except for the C25 at Brazilian site.Deforestation events were increasingly being omitted (negative bias) when using smaller percentiles (Figure 6b,e).

Effect of Spatial and Temporal Extent of the Data Cube on Optimal Percentile
Figure 7 shows how the optimal percentile changes for different combinations of spatial and temporal extents.For longer temporal extents (TE = 4 or 5 years), the optimal percentile is higher for larger spatial extents (Figure 7).Generally, the optimal percentiles at the Brazilian site were larger than those for the Bolivian site.At the Bolivian site, the optimal percentile was smaller than 3% for each spatial and temporal extent combination.

Effect of Spatial and Temporal Extent of the Data Cube on Optimal Percentile
Figure 7 shows how the optimal percentile changes for different combinations of spatial and temporal extents.For longer temporal extents (TE = 4 or 5 years), the optimal percentile is higher for larger spatial extents (Figure 7).Generally, the optimal percentiles at the Brazilian site were larger than those for the Bolivian site.At the Bolivian site, the optimal percentile was smaller than 3% for each spatial and temporal extent combination.

Effect of Spatial and Temporal Extent of the Data Cube on Optimal Percentile
Figure 7 shows how the optimal percentile changes for different combinations of spatial and temporal extents.For longer temporal extents (TE = 4 or 5 years), the optimal percentile is higher for larger spatial extents (Figure 7).Generally, the optimal percentiles at the Brazilian site were larger than those for the Bolivian site.At the Bolivian site, the optimal percentile was smaller than 3% for each spatial and temporal extent combination.

Validation
Figure 8 shows how the spatial and temporal accuracies varied in relation to spatial and temporal extents.Note that the effect of the spatial extent on spatial and temporal accuracy (Figure 8a,c,e,g), is shown for the results obtained from the longest temporal extent (TE = 5 years).The effect of the temporal extent is shown for the results obtained from the largest spatial extent (56.25 ha, Figure 8b,d,f,h).Varying the spatial extent or temporal extent had a marginal influence on the overall accuracy at both the Bolivian and the Brazilian site (Figure 8a,b).However, the change in either the spatial extent or temporal extent of the data cube had a pronounced effect on the producer's (Figure 8c,d) and user's accuracies (Figure 8e,f).Increasing the spatial or temporal extent of the data cube resulted in higher producer's accuracy and lower user's accuracy at both sites.Data cubes with a smaller spatial extent (e.g., smaller than 20 ha) produced maps with a higher user's accuracy at both study sites than larger data cubes.The reference cube with the shortest temporal extent (one year of data) produced maps with 74.4% user's accuracy at the Bolivian site, and 65.8% at the Brazilian site.Deforestation events detected at Bolivian and Brazilian sites when the reference cube contained one year of data are shown by month (Figure 9) and year (Figure 10) of detection.The variation in producer's accuracy was generally larger when varying the spatial extent of the cube than when varying the temporal extent (Figure 8c,d).
Increasing the spatial or temporal extent of the data cube resulted in a shorter median temporal detection delay at both study sites (Figure 8g,h).However, the change in temporal extent did not have an influence on the median temporal detection delay at the Brazilian site.At this site, reference cubes with a temporal extent of four or five years achieved the shortest median temporal detection delay (one observation).Similar to the producer's accuracy, the variation in median temporal detection delay was generally greater when varying the spatial extent of the cube than when varying the temporal extent (Figure 8g,h).

Validation
Figure 8 shows how the spatial and temporal accuracies varied in relation to spatial and temporal extents.Note that the effect of the spatial extent on spatial and temporal accuracy (Figure 8a,c,e,g), is shown for the results obtained from the longest temporal extent (TE = 5 years).The effect of the temporal extent is shown for the results obtained from the largest spatial extent (56.25 ha, Figure 8b,d,f,h).Varying the spatial extent or temporal extent had a marginal influence on the overall accuracy at both the Bolivian and the Brazilian site (Figure 8a,b).However, the change in either the spatial extent or temporal extent of the data cube had a pronounced effect on the producer's (Figure 8c,d) and user's accuracies (Figure 8e,f).Increasing the spatial or temporal extent of the data cube resulted in higher producer's accuracy and lower user's accuracy at both sites.Data cubes with a smaller spatial extent (e.g., smaller than 20 ha) produced maps with a higher user's accuracy at both study sites than larger data cubes.The reference cube with the shortest temporal extent (one year of data) produced maps with 74.4% user's accuracy at the Bolivian site, and 65.8% at the Brazilian site.Deforestation events detected at Bolivian and Brazilian sites when the reference cube contained one year of data are shown by month (Figure 9) and year (Figure 10) of detection.The variation in producer's accuracy was generally larger when varying the spatial extent of the cube than when varying the temporal extent (Figure 8c,d).
Increasing the spatial or temporal extent of the data cube resulted in a shorter median temporal detection delay at both study sites (Figure 8g,h).However, the change in temporal extent did not have an influence on the median temporal detection delay at the Brazilian site.At this site, reference cubes with a temporal extent of four or five years achieved the shortest median temporal detection delay (one observation).Similar to the producer's accuracy, the variation in median temporal detection delay was generally greater when varying the spatial extent of the cube than when varying the temporal extent (Figure 8g,h).

Discussion
In this paper, we proposed a data-driven space-time change detection method that exploits spatiotemporal information in satellite image time series to detect deforestation at sub-annual scales as extreme events in local data cubes.We demonstrated the space-time change detection method on Landsat NDVI image time series at a dry (Bolivian site) and humid (Brazilian site) tropical forest sites.Our results show that the method is suitable for accurate detection of deforestation events at subannual scales in both dry and humid forests even when image time series contains only one year of historical observations.We were able to achieve a median temporal detection delay of less than three observations, and producer's accuracy above 70%, user's accuracy above 65%, and an overall accuracy above 80% at both dry and humid tropical forest areas when using a data cube with one year of historical observations and a window size of 56.25 ha (Figure 8).A previous study at the same study sites [5], which used a method that only analysed individual pixel time series [3] and used all available Landsat images , also achieved at overall accuracy above 80% and a median temporal detection delay of less than four observations.Other studies [2,4,6], which used different methods to detect deforestation at sub-annual scales at different study sites, also achieved overall accuracies above 80%.These studies expressed the temporal delay in time, whereas here we express the temporal detection delay as number of observations, thus making it difficult to directly compare their temporal delay to ours.
Deforestation events were mapped more accurately at the Bolivian site when using data cubes with a large spatial extent.This is mainly because deforestation events at the Bolivian site were generally large.In areas with large deforestation events, cubes with a larger spatial extent lead to accurate deforestation mapping because the cube is less likely to be entirely within the footprint of a deforestation event.If the spatial extent of a data cube is smaller than the footprint of the deforestation event, the impact of deforestation is likely to be smoothed out when spatially normalising data to reduce seasonal variations.The spatial normalisation approach assumes that there are at least 5% forest pixels within the spatial window, since pixels' values are spatially normalised against the upper 5% tail [5].Using data cubes with small spatial extent can lead to accurate mapping of deforestation events in areas with relatively small deforestation events.With small deforestation events, spatial normalisation is less likely to smooth out the impact of deforestation in the data.This is why deforestation events were mapped accurately at the Brazilian site even when using data cubes with small spatial extent.
The incidents of false detection (low user accuracy) were particularly high for data cubes with large spatial extents because the thresholds calculated from data cubes with large spatial extents or

Discussion
In this paper, we proposed a data-driven space-time change detection method that exploits spatiotemporal information in satellite image time series to detect deforestation at sub-annual scales as extreme events in local data cubes.We demonstrated the space-time change detection method on Landsat NDVI image time series at a dry (Bolivian site) and humid (Brazilian site) tropical forest sites.Our results show that the method is suitable for accurate detection of deforestation events at sub-annual scales in both dry and humid forests even when image time series contains only one year of historical observations.We were able to achieve a median temporal detection delay of less than three observations, and producer's accuracy above 70%, user's accuracy above 65%, and an overall accuracy above 80% at both dry and humid tropical forest areas when using a data cube with one year of historical observations and a window size of 56.25 ha (Figure 8).A previous study at the same study sites [5], which used a method that only analysed individual pixel time series [3] and used all available Landsat images , also achieved at overall accuracy above 80% and a median temporal detection delay of less than four observations.Other studies [2,4,6], which used different methods to detect deforestation at sub-annual scales at different study sites, also achieved overall accuracies above 80%.These studies expressed the temporal delay in time, whereas here we express the temporal detection delay as number of observations, thus making it difficult to directly compare their temporal delay to ours.
Deforestation events were mapped more accurately at the Bolivian site when using data cubes with a large spatial extent.This is mainly because deforestation events at the Bolivian site were generally large.In areas with large deforestation events, cubes with a larger spatial extent lead to accurate deforestation mapping because the cube is less likely to be entirely within the footprint of a deforestation event.If the spatial extent of a data cube is smaller than the footprint of the deforestation event, the impact of deforestation is likely to be smoothed out when spatially normalising data to reduce seasonal variations.The spatial normalisation approach assumes that there are at least 5% forest pixels within the spatial window, since pixels' values are spatially normalised against the upper 5% tail [5].Using data cubes with small spatial extent can lead to accurate mapping of deforestation events in areas with relatively small deforestation events.With small deforestation events, spatial normalisation is less likely to smooth out the impact of deforestation in the data.This is why deforestation events were mapped accurately at the Brazilian site even when using data cubes with small spatial extent.
The incidents of false detection (low user accuracy) were particularly high for data cubes with large spatial extents because the thresholds calculated from data cubes with large spatial extents or longer temporal extents were relatively large.Such large thresholds can lead to many false detections because they might be too sensitive.Such sensitivity could also explain why deforestation events were typically detected with shorter delay when using data cubes with either large spatial extent or longer temporal extent (Figure 8).Sensitive thresholds can also explain why increasing the temporal extent of the data cube at the Bolivian site led to accurate mapping of deforestation events.Increasing the temporal extent of the data cube did not affect the temporal detection delay at the Brazilian site (Figure 8h).This is mainly because, for each temporal extent, the optimal percentile at the Brazilian site was large (Figure 7).
Increasing the spatial extent had more of a major influence on the optimal percentile at the Brazilian site than at the Bolivian site (Figure 7).This inter-site difference can be explained by the number of observations in the reference cube.Temporally, there were fewer images at the Brazilian site than at the Bolivian site (Figure 5).So, optimal percentiles at the Brazilian site were less likely to reach stability at smaller spatial extents because observations in the reference cube were few.In contrast, at the Bolivian site, the optimal percentile was most likely to reach stability at a smaller spatial extents because the images were many, and additional information from spatial context was less likely to have a major influence on the optimal percent.
Our method offers new opportunities to tackle challenges associated with existing methods for monitoring deforestation at sub-annual scales [3][4][5][6]29].In particular, our method exploits both spatial and temporal information in satellite image time series to detect deforestation at a sub-annual scale, thus allowing us to analyse pixels which do not have many historical observations.Results from the two case studies indicate that our method is robust in detecting deforestation events at a sub-annual scale, even when the image time series only contains one year of historical observations (Figure 8).One year of historical observations is often too short to properly differentiate deforestation from normal forest dynamics, especially in forests that exhibit strong seasonality.By combining spatial and temporal information, we can use image time series of high spatial resolution satellite sensors (e.g., RapidEye), whose time series are short, to track small-scale forest disturbances (e.g., selective logging).Similarly, by exploiting spatiotemporal information in image time series, there is no need to wait for image time series from newly launched sensors (e.g., Sentinel-2) to lengthen before exploiting such data to detect deforestation events at sub-annual scales.Since our method remains robust in detecting deforestation at sub-annual scales even when the reference period only contains one year of data, users may not need to pre-process huge amounts of historical data when monitoring deforestation at sub-annual scales.
The method presented in this paper can be applied in different forest areas because it is a data-driven approach using thresholds computed from the data.However, it may face challenges in mixed forests, where deciduous and evergreen forests coexist at short distances.This is mainly because of the way we reduce seasonal variations in the data cube.Normalisation against P 95t (Section 2.2) is not likely to reduce seasonal variations because P 95t would represent evergreen trees in mixed stands.This limitation may be addressed by calculating P 95t for each forest type separately and by deseasonalising pixels from each forest type using the corresponding P 95t .Another limitation is related to how we treat the pixels whose historical observations also qualify as extremes although no deforestation has occurred.
We determined optimal percentiles for identifying deforestation at sub-annual scale, but these percentiles might not be optimal for other areas with different forest types and processes causing deforestation.In areas with gradual changes [18], for example, smaller percentiles (5th percentile) might be preferable.Therefore, users should calibrate the method to identify optimal percentiles for their respective study areas before monitoring for deforestation at sub-annual scales.To do this, users should select sample pixels from both deforested and forested areas in their study areas, and apply space-time change detection method while varying the percentile for deforestation detection (e.g., in Section 3.2).Depending on the monitoring goal, for example shorter temporal detection delay, the user can then select the percentile that achieves the shortest median temporal delay as the optimal percentile.
We tested several window sizes (spatial extents) of the data cube, but identifying a window size appropriate for different parts of the globe is still challenging.This is because a window size which is optimal in one area might not be optimal in another area.Prior knowledge on the size of the deforestation events that typically occur in a particular area can be used to decide on the spatial extent of the data cube.If such prior knowledge is lacking, the user should choose a spatial extent which is larger than the size of deforestation events the user aims to detect.
With the advent of open and free access to data from Sentinel sensors, especially Sentinel-1 and -2, detecting deforestation at small spatial scales within few days of occurrence will be possible.Combining Sentinel and Landsat data will boost monitoring of deforestation at sub-annual scales, allowing agencies responsible for forest protection to timely intervene in areas where illegal deforestation events are occurring.However, such multi-source data would need harmonisation to produce multi-sensor time series which is temporally consistent.

Conclusions
In this paper, we demonstrated how spatial and temporal information can be combined and exploited to detect deforestation from satellite image time series at a sub-annual scale.We proposed a data-driven space-time change detection method that detects deforestation as an extreme event within a space-time data cube of satellite image time series.We detected sub-annual deforestation from Landsat NDVI time series at a dry tropical forest site, where the forest exhibits strong seasonality, and at humid tropical forest site.The method remained robust in detecting deforestation events at a sub-annual scale even when the image time series only contained one year of historical observations.The space-time method we presented in this paper is a novel and robust approach for timely detection of deforestation events in areas where forests exhibit strong seasonality.It provides an opportunity to detect deforestation events using image time series with scarce historical observations.The method can be used in different types of forest, both evergreen and deciduous, but, it may face a challenge in mixed forests, where deciduous and evergreen forests coexist at short distances.Although we used NDVI, the method is expected to be applicable for image time series of any satellite-derived metric that is used for deforestation monitoring.To further improve deforestation monitoring at a sub-annual scale, future research should investigate how data from different satellite sensors (i.e., Landsat 7 and 8, Sentinel-2, RapidEye, and SPOT) can be combined in a space-time change detection framework to facilitate near real-time deforestation detection.

Figure 1 .
Figure 1.The workflow for detecting deforestation at sub-annual as extreme event in space-time data cube.Figure 1.The workflow for detecting deforestation at sub-annual as extreme event in space-time data cube.

Figure 1 .
Figure 1.The workflow for detecting deforestation at sub-annual as extreme event in space-time data cube.Figure 1.The workflow for detecting deforestation at sub-annual as extreme event in space-time data cube.

Figure 2 .
Figure 2.An example of a local space-time data cube for normalised difference vegetation index (NDVI) derived from a stack of Landsat images for the period of 2004-2014.The data cube moves across the image stack (moving window), and it is split into a reference and monitoring cube.The NDVI values in the cube are not yet deseasonalised.The horizontal banding in the cube shows the effect of seasonal variations, with yellowish tones representing the dry season.

Figure 3 .
Figure 3.An example of distributions of observations in the reference cube before (a) and after (b) reducing the seasonality in a space-time data cube using spatial context.sNDVI denotes spatially normalised NDVI-after applying the spatial context method.Note that some values for spatially normalised NDVI are above 1 because we used 95 percentile to normalise the pixel values.Such values represent NDVI values which were above 95 percentile.

Figure 2 .
Figure 2.An example of a local space-time data cube for normalised difference vegetation index (NDVI) derived from a stack of Landsat images for the period of 2004-2014.The data cube moves across the image stack (moving window), and it is split into a reference and monitoring cube.The NDVI values in the cube are not yet deseasonalised.The horizontal banding in the cube shows the effect of seasonal variations, with yellowish tones representing the dry season.

Figure 2 .
Figure 2.An example of a local space-time data cube for normalised difference vegetation index (NDVI) derived from a stack of Landsat images for the period of 2004-2014.The data cube moves across the image stack (moving window), and it is split into a reference and monitoring cube.The NDVI values in the cube are not yet deseasonalised.The horizontal banding in the cube shows the effect of seasonal variations, with yellowish tones representing the dry season.

Figure 3 .
Figure 3.An example of distributions of observations in the reference cube before (a) and after (b) reducing the seasonality in a space-time data cube using spatial context.sNDVI denotes spatially normalised NDVI-after applying the spatial context method.Note that some values for spatially normalised NDVI are above 1 because we used 95 percentile to normalise the pixel values.Such values represent NDVI values which were above 95 percentile.

Figure 3 .
Figure 3.An example of distributions of observations in the reference cube before (a) and after (b) reducing the seasonality in a space-time data cube using spatial context.sNDVI denotes spatially normalised NDVI-after applying the spatial context method.Note that some values for spatially normalised NDVI are above 1 because we used 95 percentile to normalise the pixel values.Such values represent NDVI values which were above 95 percentile.

Figure 4 .
Figure 4.An overview of the location of the study sites south east of Santa Cruz de la Sierra, Bolivia, and west of Ariquemes, Rondonia State, Brazil.The base images are band 1-2-3 composites of Landsat ETM+ images from 15 April 2013 (Bolivia) and 21 July 2014 (Brazil).

Figure 4 .
Figure 4.An overview of the location of the study sites south east of Santa Cruz de la Sierra, Bolivia, and west of Ariquemes, Rondonia State, Brazil.The base images are band 1-2-3 composites of Landsat ETM+ images from 15 April 2013 (Bolivia) and 21 July 2014 (Brazil).

Figure 5 .
Figure 5. Cumulative number of images in the reference cube versus the temporal extent of the reference cube at Bolivian and Brazilian study sites.

Figure 5 .
Figure 5. Cumulative number of images in the reference cube versus the temporal extent of the reference cube at Bolivian and Brazilian study sites.

Figure 6 .
Figure 6.Change in the overall accuracy, bias, and the median temporal detection delay (tD) as a function of the percentile for defining deforestation as extreme event in Landsat space-time data cubes (C5 and C25) at the Bolivian (a-c) and Brazilian (d-f) study sites.The blue lines denote zero bias.Bias was calculated by subtracting the omission error from the commission error ([11]).

Figure 7 .
Figure 7. Change in the optimal percentile for detecting deforestation events at sub-annual scales when varying the spatial extent and temporal extent (TE) of the data cube at the Bolivian and Brazilian site.The optimal percentile was determined by choosing the percentile that achieves the shortest median temporal detection delay and highest overall accuracy based on the training data set.

Figure 6 .
Figure 6.Change in the overall accuracy, bias, and the median temporal detection delay (tD) as a function of the percentile for defining deforestation as extreme event in Landsat space-time data cubes (C5 and C25) at the Bolivian (a-c) and Brazilian (d-f) study sites.The blue lines denote zero bias.Bias was calculated by subtracting the omission error from the commission error ([11]).

Figure 6 .
Figure 6.Change in the overall accuracy, bias, and the median temporal detection delay (tD) as a function of the percentile for defining deforestation as extreme event in Landsat space-time data cubes (C5 and C25) at the Bolivian (a-c) and Brazilian (d-f) study sites.The blue lines denote zero bias.Bias was calculated by subtracting the omission error from the commission error ([11]).

Figure 7 .
Figure 7. Change in the optimal percentile for detecting deforestation events at sub-annual scales when varying the spatial extent and temporal extent (TE) of the data cube at the Bolivian and Brazilian site.The optimal percentile was determined by choosing the percentile that achieves the shortest median temporal detection delay and highest overall accuracy based on the training data set.

Figure 7 .
Figure 7. Change in the optimal percentile for detecting deforestation events at sub-annual scales when varying the spatial extent and temporal extent (TE) of the data cube at the Bolivian and Brazilian site.The optimal percentile was determined by choosing the percentile that achieves the shortest median temporal detection delay and highest overall accuracy based on the training data set.

Figure 8 .
Figure 8. Overall accuracy (a,b); producer's accuracy (c,d); user's accuracy (e,f); and median temporal detection delay (tD) (g,h) of deforestation events at the Bolivian and Brazilian study sites as a function of the spatial (left) and temporal extent (right) of the local data cube.Vertical bars (a-f) indicate 95% confidence intervals.The effect of the spatial extent on spatial and temporal accuracy (Figure 8a,c,e,g), is shown for the results obtained from the longest temporal extent (TE = 5 years).The effect of the temporal extent is shown for the results obtained from the largest spatial extent (56.25 ha, Figure 8b,d,f,h).

Figure 9 .
Figure 9.The month in which deforestation events detected from Landsat time series at Bolivian site (a) and Brazilian site (b) over the period of 2005-2014 occurred.Deforestation was detected as an extreme event in Landsat data cube with a reference period containing one year of data.

Figure 8 . 15 Figure 8 .
Figure 8. Overall accuracy (a,b); producer's accuracy (c,d); user's accuracy (e,f); and median temporal detection delay (tD) (g,h) of deforestation events at the Bolivian and Brazilian study sites as a function of the spatial (left) and temporal extent (right) of the local data cube.Vertical bars (a-f) indicate 95% confidence intervals.The effect of the spatial extent on spatial and temporal accuracy (Figure 8a,c,e,g), is shown for the results obtained from the longest temporal extent (TE = 5 years).The effect of the temporal extent is shown for the results obtained from the largest spatial extent (56.25 ha, Figure 8b,d,f,h).

Figure 9 .
Figure 9.The month in which deforestation events detected from Landsat time series at Bolivian site (a) and Brazilian site (b) over the period of 2005-2014 occurred.Deforestation was detected as an extreme event in Landsat data cube with a reference period containing one year of data.

Figure 9 .
Figure 9.The month in which deforestation events detected from Landsat time series at Bolivian site (a) and Brazilian site (b) over the period of 2005-2014 occurred.Deforestation was detected as an extreme event in Landsat data cube with a reference period containing one year of data.

Figure 10 .
Figure 10.The year in which deforestation events detected from Landsat time series at Bolivian site (a) and Brazilian site (b) over the period of 2005-2014 occurred.Deforestation was detected as an extreme event in Landsat data cube with a reference period containing one-year of data.

Figure 10 .
Figure 10.The year in which deforestation events detected from Landsat time series at Bolivian site (a) and Brazilian site (b) over the period of 2005-2014 occurred.Deforestation was detected as an extreme event in Landsat data cube with a reference period containing one-year of data.

Table 1 .
Data cube names and their spatial extents and dimensions for the Bolivian and Brazilian sites.

Table 1 .
Data cube names and their spatial extents and dimensions for the Bolivian and Brazilian sites.

Table 2 .
Number of sample pixels used for training purposes and final validation of deforestation maps produced at the Bolivian and Brazilian sites.

Table 2 .
Number of sample pixels used for training purposes and final validation of deforestation maps produced at the Bolivian and Brazilian sites.