R econstruction of satellite time series with a dynamic smoother

Time series reconstruction methods—used to generate gap-free time series of satellite observations—were historically designed for sensors with frequent image acquisitions. Since 2008, interest in leveraging time series methods has shifted from sensors such as AVHRR and MODIS to Landsat because of free, higher-resolution data availability and improved access to high-performance compute systems. Existing methods are typically designed for speciﬁc applications such as land cover classiﬁcation or for estimating the timing of phenology events. Moreover, approaches developed for speciﬁc ecological systems, such as tropical forests or temperate agriculture, often do not generalize well across land cover, vegetation, and climate types. In this study, we introduce a dynamic temporal smoothing (DTS) method to reconstruct sparse, noisy signals into dense time series at regular intervals. The DTS is a weighted smoother with dynamic parameters that is applied over a signal. The smoother is intended to have wide applicability, with particular focus on applications in vegetation remote sensing. In this paper we present and illustrate the DTS over short-and long-term Landsat (TM, ETM + , and OLI) time series and demonstrate the e ﬀ ectiveness of robust gap-ﬁlling over a range of landscapes in the South American Southern Cone region.


Introduction
Smoothing and curve-fitting techniques for reconstructing time series of remotely sensed observations have been used for decades (Chen et al., 2004;Sakamoto et al., 2005Sakamoto et al., , 2010Ma et al., 2006). The majority of these techniques have been developed for time series derived from instruments such as Advanced Very High Resolution Radiometer and Moderate Resolution Imaging Spectroradiometer. Because these instruments provide daily (or near-daily) repeat frequency, relatively complete time series can be created for many parts of the world using 8-or 16-day composites. Hence, for most applications, widely used reconstruction techniques such as those cited above provide time series that are accurate and gap-free. Gap-filling moderate-resolution data has become easier since the launch of the Sentinel-2 satellites (Drusch et al., 2012). Nevertheless, over the last decade the need for robust and accurate methods to reconstruct historical imagery has become more pronounced as interest in leveraging time series of Landsat imagery to map current and historical changes in Earth surface properties has increased (Hansen et al., 2013).
Techniques designed to overcome the high frequency and duration of data gaps arising from cloud cover in Landsat imagery (and in the case of Landsat 7, from the scan-line corrector failure) are increasingly a focus of research. More specifically, gap-filling (the estimation of missing data at observed dates) and reconstruction (re-gridding of sparse, irregular time series to regular intervals) are used to increase time series density and to standardize metrics for land cover mapping. Until recently, the majority of methods developed for gap-filling or reconstruction of Landsat imagery had been created primarily in support of applications focused on land cover (Kennedy et al., 2010;Griffiths et al., 2013Griffiths et al., , 2019Zhu et al., 2014;DeVries et al., 2015;Roy et al., 2018), phenology (Jönsson et al., 2002(Jönsson et al., , 2018Bolton et al., 2020), or for methodological purposes (Chen et al., 2004;Gao et al., 2006Gao et al., , 2017Roy et al., 2018;Yan et al., 2018Yan et al., , 2020. Many of these approaches generate high-quality datasets that can be used for a variety of applications. However, existing methods do not always generalize well outside their designed purpose. Sinusoidal harmonics, for example, provide a robust basis for classifying forests over long time periods (Arévalo et al., 2019) but have not been proven in more dynamic ecological systems (e.g., agriculture, dry savannas). Similarly, functional curve fitting approaches (Jönsson et al., 2018), which provide high-precision estimates in single-season temperate ecosystems, can be intractable over more dynamic time series. And, while gap-filling by seasonal composites (Flood, 2013;Griffiths et al., 2013;Schmidt et al., 2016) can generate clean, consistent data, these methods may not be sufficiently accurate to estimate land surface phenology.
The goal of this paper is to describe a new method for reconstructing sparse satellite time series that is applicationindependent and, therefore, suitable for a variety of purposes. By providing smooth and gap-free time series of observations at each pixel, the technique we describe provides a general, robust, and accurate basis for creating high-quality, dense time series of satellite observations that are well-suited for analyses of longterm trends and changes in land cover, land use, and ecological processes such as phenology.

Data
Our analysis used Landsat 5 (TM), 7 (ETM+), and 8 (OLI) Collection 1 (C1) surface reflectance data. Our test sites were located in the South American Southern Cone region, which includes a wide range of land cover types and ecological conditions across Argentina, Chile, Paraguay, and Uruguay. Using the on-demand Science Processing Architecture (ESPA), we selected all available imagery from May 1999 to November 2019 that had less than 90% estimated cloud cover over the study region. We used the C1 Quality Assessment (QA) layer to identify all pixels flagged as cloud, cloud shadow, snow, ice, or fill-value in each image (USGS EROS Center, 2021). Scene-wise reflectance variations (bi-directional reflectance distribution function (BRDF) effects) arise from atmospheric conditions and differences across the Landsat track (Toivonen et al., 2006;Hansen et al., 2012;Potapov et al., 2012;Flood et al., 2013). To support corrections for BRDF effects, we created pixel-level sensor and solar azimuth and zenith angle layers using the USGS-provided C1 Angle Coefficient file and the Landsat Angles Creation tool (USGS EROS Center, 2019). The resulting surface reflectance and solar geometry layers were gridded and stored using the South America Albers Equal Area Conic projection with 30 meter spatial resolution.
Using the solar and sensor geometry at each pixel, we corrected C1 surface reflectance values for BRDF effects using the cfactor technique (Roy et al., 2016;Claverie et al., 2018), where the solar zenith angle parameter was derived from a 6 th -order polynomial function of latitude (Zhang et al., 2016). To simplify our presentation, here we focus on reconstructing time series of the two-band enhanced vegetation index (EVI2) (Jiang et al., 2008), which we computed using the BRDF-adjusted surface reflectance data at each pixel as Hereafter, we focus on EVI2 series exclusively, thus any reference to time series refers to times series of EVI2 values at each pixel. Note, however, that our approach is equally valid for application to time series of surface reflectance values or other spectral indices.

Time series pre-processing
The following Sections (2.2.1 -2.2.3) describe the preprocessing steps applied prior to smoothing a time series; the end-to-end processing flow to reconstruct time series is presented in Fig. 1. The core element of our algorithm is a dynamic temporal smoother that can be used with or without the preprocessing steps described in the following sections.

Spatial gap-filling
Data gaps arising from ETM+ scan errors that persist over long time periods (e.g., multiple months), particularly during vegetation greening, can cause irregularities that our DTS-based smoother cannot overcome. Therefore, in the first pre-processing step, we applied spatial interpolation to fill data gaps. The step was designed to address ETM+ scan error gaps, but also conveniently improved instances where a pixel was flagged as nonclear by the C1 cloud mask. To apply this gap-filling procedure, any pixel flagged as either missing or non-clear was filled using local least squares polynomial regression following the method described by Graesser et al. (2018). This method performs spatial gap-filling over variable window sizes that range from 7 to 25 pixels (for details, we refer the reader to Graesser et al. (2018)). We also enhanced the method by including windows from additional years (in this study, ± 1 year). Missing values were not filled if more than 80% of the 25 pixels located around a center pixel were missing or if the r 2 of the local fit was less than 0.5. Hence, in addition to retaining noise (i.e., non-gap pixels were not altered), the resulting time series at each pixel included gaps where the spatial gap-filling was not applied.

Outlier detection
In the second step of the pre-processing phase, we replaced outliers in the time series at each pixel, which are generally caused by undetected clouds or cloud shadows in the C1 product QA layer (Bolton et al., 2020). To identify such outliers, we used a two-step moving window regression procedure. In the first step, an outlier o at time series value y was detected by weighted residuals from linear estimates of three consecutive time series observations distributed over a maximum time period of o d days. For each set of three observations, we estimated the center value by linear interpolation between the two end values, and then flagged the center value as an outlier if its normalized weighted residual exceeded o thresh . The outlier threshold here was determined by simulating a range of parameters over useridentified outliers. More formally, the outlier test was defined as whereŷ d2 is the interpolated center value and o w = w1 × w2. w1 and w2 were designed to apply more weight to values below the time series curve because we assumed most outliers were caused by clouds or cloud shadows. We applied a weight to the time series value using w1, computed as if (y d2 >= 0.45) and min(y d1 , y d3 ) <= 0.15 or 0.03 < y d2 <= 0.15 and 0.03 <= y d1 <= 0.15 or 0.03 <= y d2 <= 0.15 and 0.03 <= y d3 <= 0.15 1 ifŷ d2 >= y d2 otherwise 0.33 otherwise (3) and used w2 to weight the distance between the observation acquisition dates (d) by and L is a logistic function computed by L(y, y 0 , r) = 1 1 + e −r(y−y 0 ) .
Only local outlier spikes were removed by ensuring that y d1 < y d2 > y d3 or y d1 > y d2 < y d3 . A final check was applied to avoid Figure 1: Key steps in the time series smoothing process. Two steps, the input data and dynamic temporal smoothing (DTS), are required. The other steps-spatial gap-filling, outlier detection, and post-DTS filtering-are optional enhancements. The DTS is designed to use as input any time series scaled between 0 and 1. Thus, reflectance bands or normalized indices are viable options. The spatial gap-filling, which is applied as the first pre-processing step, along with the post-DTS filtering are spatial-temporal methods, whereas outlier detection and the DTS are applied on one-dimensional time series at each pixel.
removing values along inflection points: outliers from Eq. 2 were not removed if the proportional difference between either y d1 and y d2 , y d1 and y d3 , y d3 and y d2 , or y d3 and y d1 was below a predetermined threshold (-0.6 in this study). Equation 2 was applied to the moving windows across time series at each pixel and detected outliers were replaced via linear interpolation (i.e., byŷ d2 ).
Outliers located near inflection points can be difficult to detect (i.e., prone to errors of commission) using the linear three-point procedure. Therefore, we applied a second, polynomial regression outlier test to o c consecutive values occurring within o d days of each other. For this second step, w1 was computed using a simpler approach, where w1 = 1 ifŷ o c /2 >= y o c /2 ; otherwise w1 = 0.25, and w2 = 1 − y o c /2 . If the center value of the o c values was an outlier, it was replaced using a polynomial fit (i.e.,ŷ o c /2 = β 0 o c /2 + β 1 (o c /2) 2 + α). The same proportional difference check that was applied with the three-point test in the previous paragraph (using the -0.6 threshold) was also applied here. In this study, we used o d = 120 and o thresh = 0.2 for the three-point linear test, and o d = 120, o c = 7, and o thresh = 0.2 for the polynomial test. The outlier process is illustrated in Fig. 2.

Interpolation to daily time series
In the next step, we interpolated the time series at each pixel to create daily time series. To match the timing of the growing season in the Southern Hemisphere, we defined annual periods for the study region as twelve-month periods beginning and ending on the 1st of July. To reconstruct time series for a single pixel, we extracted 20-month time periods (12 months + 4 months on either end of the series), which ensured continuity in the time series over the 20-year time period we considered (see Bolton et al. (2020)). For example, to process time series for the period 2000-July-1 to 2001-July-1, we extracted all available data for the period from 2000-March-1 to 2001-November-1. We then linearly interpolated the sparse, gap-filled data (i.e., after spa- If the end dates of three consecutive points were within 120 days of each other, the center value was assessed using Eq. 2. If the equation was satisfied, the outlier was removed and replaced with the linearly interpolated value. In the second stage of the outlier test, a seven-point window was used (see Section 2.2.2 for details).
tial gap-filling and outlier detection and adjustment) to a daily series.

Dynamic temporal smoother
The dynamic temporal smoother (DTS) has two key properties. First, it retains sensitivity to high-frequency peaks (e.g., multiple crop cycles over a field) and is, therefore, flexible to different signals. Second, it is robust to sporadic noise in the time series signal, such as outliers missed in pre-processing.

Method overview
The DTS uses a one-dimensional moving window that dynamically adjusts between k min and k max to smooth daily values of a time series (x) on each day (d). Formally, the DTS is defined as Logistic function midpoint to adjust the window size k r Logistic function growth rate to adjust the window size g 0 Logistic function midpoint to adjust the gaussian sigma g r Logistic function growth rate to adjust the gaussian sigma t 0 Logistic function midpoint to adjust the time weight t r Logistic function growth rate to adjust the time weight υ σ Gaussian sigma to adjust the radiometric weight s t σ Gaussian sigma for the spatial-temporal time weight s b σ Gaussian sigma for the spatial-temporal brightness weight s r σ Gaussian sigma for the spatial-temporal radiometric weight s k σ Gaussian sigma for the spatial-temporal distance weight where x is the smoothed time series, DT S is the dynamic temporal smoother function with parameters listed in Tab. 1 that compute a set of dynamically adjusted weights, A is a function (described in Section 2.3.3) that adjusts a smoothed value (i.e., x d +k /2 ) to the upper envelope of the time series, and n is the total number of days in the time series. In this framework, the smoothed value of x is computed as a weighted average of x values by where w t j , w b j , w r j are time, brightness, and radiometric weights, and d and k are the adjusted window starting position and length, respectively. These parameters are adjusted as the moving window passes over the time series. Full details of this process are given in Sections 2.3.2 -2.4.

Window adjustment
With DTS, the size of the window is adjusted dynamically to range between a minimum (k min ) and maximum (k max ) length using a logistic scaler (i.e., Eq. 6). The adjustment is optimized for periods when changes in x are frequent by reducing the window length during periods when values of x are increasing or decreasing rapidly (e.g., vegetation greenup or greendown). To adjust the window size to k , we iterated over x starting at day d with a moving window of length k max , recording the window's center value and standard deviation (x σ d..d+k ), which is normalized from a range of [0, 0.05] to [0, 1] by linear scaling. This window adjustment procedure was designed to adjust values of k to be larger (i.e., closer to k max ) if the time series of x had low variance and was near a low baseline value. To do this, we calculated the proportional difference between the window's center value and 0.1 (x p∆ d+k max /2 ), which was clipped to have a range of [0, 1]. Then, we calculated the window's center weight (x c d+k max /2 ) as 0.5(x σ d..d+k max × x p∆ d+k max /2 ), and then calculated the new window length (k ) using the window's center weight and a scaled logistic function by Thus, the dynamic window extends as x d ..d +k , where d is the adjusted window starting position centered on d + k max /2. Following these window adjustments, window weights were assigned using the new window size k starting at adjusted window starting position d , described in the following section.

Weighting
To accommodate the diversity of temporal dynamics present in time series of remotely sensed imagery, the moving window should account for frequent changes in time series cycles as well as variations in change rates between signal peaks and valleys. Thus, we weighted values within the adjusted window size using a Gaussian or a logistic function. We used a Gaussian function for time segments with low values and low variance. In contrast, for time series segments with high values or high variance (e.g., inflection points or narrow vegetation peaks), we used a logistic function to weight window values, which skewed weights toward window tails, depending on the steepness and slope direction of the time series window segment. To do this, each day j (see Eq. 8) in the window of interest was first assigned a weight using where k∆ = |(d +k /2) − d |. This scaling forced weights preceding the center position to be negative and weights after the window center position to be positive. We estimated the curve steepness with a simple check of the proportional difference between the window endpoints. If the slope along the window curve was relatively flat, the window weights were calculated with a Gaussian function. Specifically, if was satisfied, then we used a Gaussian function to compute the time weight as where g σ = L(x c d+k max /2 , g 0 , g r ) and G is the Gaussian function defined as G(y, σ) = e −y 2 /2σ 2 . If, on the other hand, Eq. 11 was not satisfied, the time weight was calculated using a logistic function. To account for high rates of increasing or decreasing slopes, we emphasized window tails by applying To simplify the illustration, we only highlight two dynamic temporal smoothing calculations. The first, k , is the adjusted window size and is shown here in the dashed green line. The window size increases close to the maximum of 61 over the time series valleys, and decreases toward 15 over the time series peaks and transitions. The second calculation highlighted here is the window weight, w t j , and is shown in the solid light blue line. The filled grey areas show the width of the modified window size k . Note that we only highlight a subset of windows for illustration purposes. Over the time series peaks w t j is calculated by a Gaussian function, whereas over the inflection points the time weight is calculated by a logistic function. where and After the time weight was calculated, we weighted brightness differences (i.e., whether window values were higher or lower than the window center value) by and window color differences (i.e., radiometric distance from the window center value) by and the vector of radiometric weights RW was created by Iteratively adjusting reconstructed time series to the upper envelope of the time series curve is a common practice (Fisher et al., 2006;Yang et al., 2015) because values below the main time series curve (i.e., the upper envelope) tend to be signal noise, or outliers. Therefore, after calculating and applying the weighted average (i.e., Eq. 8), we re-adjusted the smoothed value toward the upper envelope at position d + k /2 by The weight adjustment was applied at each window iteration, thereby influencing the subsequent iteration at each step. Fig. 3 illustrates the DTS smoothing process over a short time series. In this example, the original values (grey points) were smoothed to a daily time-step (solid purple line). The dashed green line and solid cyan line illustrate changes in the adjusted window size (k ) and adjusted window weight (w t j ).

Generation of weekly time series and spatial-temporal filtering
We smoothed the daily time series with the following DTS parameters: k max = 61, k min = 15, k 0 = 0.5, k r = −10, g 0 = 0.5, g r = −10, t 0 = 0.5, t r = 15, and υ σ = 0.1. We applied three iterations for each pixel, using the smoothed values of each iteration as input to the subsequent iteration. We then sub-sampled the smoothed daily time series to weekly values, starting on the first day of each month, and reset the indexing on the first day of each month (e.g., for July, we extracted values for July at days 1, 7, 14, 21, and 28, and then restarted the process on August 1). We applied this process across two-year time steps, which resulted in 60 smoothed values per year across the time series at each pixel.
The one-dimensional DTS addresses signal noise over time, but it does not account for noise across space. Therefore, to generate more spatially-consistent data, we applied a three-dimensional spatial-temporal smoothing filter, similar to a two-dimensional bilateral filter (Tomasi et al., 1998) with an additional time weight. This step was applied to weekly data after the DTS has been applied. Specifically, the value at time p, row m, and column n in image I was adjusted as: where the brightness (bw) and radiometric (rw) weights were calculated using Eq. 16 and 17, respectively, with the brightness σ parameter set as (s b σ = 0.1). For the two-dimensional weights, the key difference between this spatial-temporal smoother and the bilateral filter is the addition of the spatial distance weight (kw), which is a Gaussian-weighted euclidean distance within a moving two-dimensional window of k x k size over temporal window t = 5 as Finally, the time weight (tw) was calculated by G(pt[−1, 1], s t σ ), with position pt normalized between a range of [−1, 1] around the window center.

Reconstruction comparison
We illustrate the efficacy of the DTS time series reconstruction in two ways: by 1) introducing artificial gaps over space and 2) assessing reconstructed values from smoothing methods against reconstructed values identified by manual interpretation of a time series.
In the first test, we assessed the reconstructed time series by introducing synthetic gaps regularly spaced across an image and selected on random dates . Then, we applied DTS to the time series with artificial gaps and calculated the bias between the reconstructed values and the original, pre-gap values.
The goal of the second test was to assess how well DTS, Cubic splines, Lowess, Savitsky-Golay, and Whittaker smoothers reconstructed EVI2 time series by comparing the smoothed values generated by the different methods against manually interpreted reconstructed values. To do this, we first generated a set of sample time series to manually interpret using a 165 x 165 row/column (nearly 5 km x 5 km) Landsat image as the reference for eleven sites across the Southern Cone. Then, we extracted a two-year time series (2017-July-1 to 2019-July-1) at every twentieth pixel along the northwest to southeast diagonal, resulting in 8 sample locations per site. To manually interpret each sample time series, we then plotted the raw EVI2 time series along with date markers at every tenth day of the time series, resulting in 18 dates for each sample. The interpreter then identified the EVI2 value they determined should be the reconstructed value at each date marker, using only the raw time series as a reference. We then reconstructed the time series for each sample with DTS, Cubic spline, Lowess, Savitsky-Golay, and Whittaker, using static smoothing parameters for the non-DTS methods. Finally, we calculated the mean squared error (MSE) of the reconstructed values and the manually interpreted values.

Time series reconstruction
The DTS adapts well to vegetation cycles of variable length and amplitude over time series. For example, individual time series for pixels composed of crops and tree cover are shown in Fig. 4 and Fig. 5, respectively. The weekly DTS results follow the upper curve of the original Landsat EVI2 values. Outliers, particularly low values (e.g., Fig. 4(b)), are largely detected and eliminated. Importantly, the DTS is also adaptive to land change occurring during the time series. Fig. 6 shows an example of Dry Chaco forest in northern Argentina that was cleared in 2005 and replaced by agriculture. The reconstructed seasonal cycles during the first four years when tree cover is still present are consistent along the upper curve, despite substantial signal noise. In 2002, for example, the cluster of low values below the upper envelope was likely caused by clouds that were not detected in Landsat Collection 2, but had only minor impact on DTS results because of the iterative upper envelope adjustments included in the algorithm.

Spatial-temporal filtering
The DTS produces 60 gap-filled values during each twelvemonth cycle, which is sufficient to capture phenological dynamics in vegetation. To illustrate, Fig. 7 shows an annual time series of EVI2 images for an area in Central Argentina dominated by croplands, and hence with pronounced phenology. In this example, the overall pattern in phenology is obvious and DTS results clearly identify field-level differences in management and crop types. A few fields with cereal grains show peak EVI2 in September, but the principle crops in this region reach maturity in February. Landsat 7 scan-line errors are not evident.

Reconstruction comparison
To provide a more comprehensive assessment of DTS results, we performed two analyses. First, we performed a visual assessment    of the quality of DTS gap-filling, which is presented in Fig. 8. The upper left panel ( Fig. 8(a)) shows the EVI2 image computed from a Landsat 7 ETM+ data acquired on February of 2019, which includes gaps from the ETM+ SLC-off. In the upper right panel (Fig. 8(b)), we systematically introduced gaps (white squares), and then used these artificial missing data to assess the DTS gap-filling. The DTS gap-filled image is shown in the lower left panel (Fig. 8(c)), and nominally corresponds to one day after the original image was acquired. As is evident, the gap-filled image is largely free of artifacts from gap locations. To further illustrate, the lower right panel (Fig. 8(d)) shows a difference image (DTS -original). Modest DTS errors along field edges are apparent. Overall, however, differences are very modest.
The results from the time series manual interpretation exercise are shown in Figs. 9 and 10. All smoothing methods performed relatively well. The DTS and Lowess smoothers had the lowest

Discussion and conclusion
We introduced a novel, dynamic time series smoother developed for remotely sensed land surface observations designed to have utility for a wide array of remote sensing applications. In general, it is most suitable for applications that require smooth and gap free reconstructions of historical Landsat imagery. The pre-reconstruction gap-filling, outlier tests, and iterative adjustments to the upper envelope included in the algorithm help the smoother avoid signal noise and results in smooth cycles without assuming or imposing a prescribed functional form (e.g., double logistic).
A key attribute of the DTS is that it is dynamic, non-linear, and adaptive. Hence, it can flexibly and freely adjust to a wide range of both inter-and intra-season variability. It therefore provides a flexible basis in support of diverse applications that range from reconstructing agricultural rotations, to estimating the phenology of natural vegetation, to creating features used for image classification. Having said this, it is important to note that the utility of the DTS is ultimately constrained by the density and depth of the image archive at any location. Stated another way, the DTS cannot overcome the lack of Landsat data prior to 1999 in places like Alaska, West Africa, and Siberia if spatial gap-filling does not result in sufficiently enhanced time series.
We demonstrated the effectiveness of the DTS with Landsat time series of EVI2 images. However, the method is not constrained to a specific sensor, wavelength, or vegetation index. The method described above includes some prescribed thresholds that are specific to the EVI2. However, it would be easy to estimate corresponding thresholds for other sensors, bands, or vegetation indices with a similar dynamic scaling. For example,  time series that combine Landsat with other sensors, such as the Harmonized Landsat Sentinel-2 (HLS) analysis ready data (Claverie et al., 2018), are a prime example of potential DTS use.

Code availability
Open satellite archives (such as Landsat and Sentinel) are already large and are rapidly growing. Hence, algorithm efficiency is important and was explicitly considered in developing the DTS method. More specifically, we intend the method to be sufficiently general for use in diverse remote sensing applications, but also computationally efficient, allowing it to scale from local to regional to global study areas. To address this, the DTS was written in Cython with parallelization over individual pixels using OpenMP to manage the shared-memory multiprocessing. The DTS is organized as a Python library available at https://github.com/jgrss/satsmooth.