A computational framework for processing time-series of earth observation data based on discrete convolution: global-scale historical Landsat cloud-free aggregates at 30 m spatial resolution

A computational framework for time-series processing

The implementation of several tools used to preprocess Earth Observation (EO) time-series involves the usage of a discrete convolution. Generally speaking, convolution is the mathematical operation that takes two functions, slides them one over the other, multiply them point by point, and integrates their product. In fact, the methods that involve a sliding window in time (or also in space) are often based on convolution, and many of them are used for EO data processing. Combining the convolution operation with the Hadamard (element-wise) product and division, TSIRF can also be used to perform time-series reconstruction and temporal data aggregation. In addition, the implementation of the framework is provided using three different computational back-ends. The selection of the most convenient back-end depends on the application and used computational infrastructure. Implementation details are provided after the definitions and notation.

Theoretical background and notation

Consider two real-valued functions $v(t)$ and $w(t)$, where $v(t)$ represents the variable of interest in function of time (e.g., blue band value for a certain pixel in EO time-series), and $w(t)$ is the convolution kernel, input that is defined depending on the target processing task. The uni-dimensional convolution between the two variables is defined as

$$(v\ast w)(t)={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}v(\tau )w(t-\tau )d\tau ,$$where the integrated variable $\tau$ is used to slide one of the function on the other one (Oppenheim & Schafer, 1975). The time-discrete version of convolution can be defined for uniformly sampled versions of the two variables as

$$(v\ast w)[n]=\sum _{m=-\mathrm{\infty }}^{\mathrm{\infty }}v[m]w[n-m],$$where $n$ is the index of a time sample, with sampling period $T_s$ and frequency $f=\frac{1}{T_s}$. As detailed in the following, in order to optimize computational performance in case of long time-series, it is convenient to perform the discrete convolution using Fast Fourier Transform (FFT). In order to do that, we extend here the notation to the also define the circular convolution. Noting that the observed variable $v[n]$ is only available for a limited number of samples $N_s$, relative to the observation time-frame $(T_0,T_0+N_s{T}_s)$, and supposing to have a convolution kernel with limited support in the sample range $(-N_p,N_f)$, where $N_p$ is the number of non zero samples in the “relative past” and $N_f$ is the number of non zero samples in the “relative future”, it is possible to define a periodic extensions of the two time-discrete variables, $v_e[n]$ and $w_e[n]$, with periodicity $N_e{T}_s$ where

$$N_e=N_s+max(N_p,N_f).$$

In order to facilitate the reading of the equations, the symbology is summarized in Table 1. Following this notation, we define their circular convolution (Bamieh, 2018) as

(1) $$(v_e⊛w_e)[n]=\sum _{m=0}^{N_e-1}{v}_e[m]w_e[n-m].$$

Table 1:

Summary of defined symbols.

Symbol	Description
$T_s$	Sampling period
$N_s$	Number of samples of the time series
$N_p$	Number of elements of the convolution kernel in the relative past
$N_f$	Number of elements of the convolution kernel in the relative future
$N_w$	Total number of elements of the convolution kernel
$N_e$	Number of the periodic extension of the time-series
w	Vector containing the convolution kernel samples
${\mathit{w}}_p$	Sub-vector of w with elements from the relative past
${\mathit{w}}_f$	Sub-vector of w with elements from the relative future
$N_r$	Number of time-series to which a the method is applied
$N_a$	Aggregation factor for temporal aggregation
$\epsilon$	Floating point machine precision

DOI: 10.7717/peerj.18585/table-1

Note that from previously defined range of the support of $w[n]$, $(-N_p,N_f)$, the total number of element in $w[n]$ support is $N_w=N_p+N_f+1$. Assuming that $N_p<N_s>N_f$, periodic extension of $v_e[n]$ and $w_e[n]$ can be performed with zero-padding for newly introduced elements such that $(v\ast w)[n]=(v_e⊛w_e)[n]$. A visual interpretation of the circulant convolution operation, inspired by the one in Guo et al. (2023), is shown in Fig. 1.

In order to describe the different computational backends for the convolution, we introduce the vectorial notation of the previously defined quantities. The values of the time-discrete variables $v[n]$ and $w[n]$ are respectively stored in the row vectors.

$$\mathit{v}={\left[\begin{array}{c}{v}_0\\ v_1\\ ⋮\\ v_{N_s-2}\\ v_{N_s-1}\end{array}\right]}^{\mathrm{\top }}={\left[\begin{array}{c}v[0]\\ v[1]\\ ⋮\\ v[N_s-2]\\ v[N_s-1]\end{array}\right]}^{\mathrm{\top }}and\mathit{w}={\left[\begin{array}{c}{w}_0\\ w_1\\ ⋮\\ w_{N_f-1}\\ w_{N_w-N_p}\\ ⋮\\ w_{N_w-2}\\ w_{N_w-1}\end{array}\right]}^{\mathrm{\top }}={\left[\begin{array}{c}w[0]\\ w[1]\\ ⋮\\ w[N_f-1]\\ w[-(N_p-1)]\\ ⋮\\ w[-2]\\ w[-1]\end{array}\right]}^{\mathrm{\top }}.$$

Note that the indexing notation concatenates the values of the relative past in the convolution kernel after those of the relative future, from the least recent, $w[-(N_p-1)]$, to the most recent, $w[-1]$ (Fig. 1). Similarly, the vectors storing their periodic extensions are defined as

$${\mathit{v}}_e={\left[\begin{array}{c}{v}_{e;0}\\ v_{e;1}\\ ⋮\\ v_{e;N_s-1}\\ v_{e;N_e-max(N_p,N_f)}\\ ⋮\\ v_{e;N_e-2}\\ v_{e;N_e-1}\end{array}\right]}^{\mathrm{\top }}={\left[\begin{array}{c}v[0]\\ v[1]\\ ⋮\\ v[N_s-1]\\ 0\\ ⋮\\ 0\\ 0\end{array}\right]}^{\mathrm{\top }},{\mathit{w}}_e={\left[\begin{array}{c}{w}_{e;0}\\ w_{e;1}\\ ⋮\\ w_{e;N_f-1}\\ w_{e;N_f}\\ ⋮\\ w_{e;N_f-N_p-1}\\ w_{e;N_e-N_p}\\ ⋮\\ w_{e;N_e-2}\\ w_{e;N_e-1}\end{array}\right]}^{\mathrm{\top }}={\left[\begin{array}{c}w[0]\\ w[1]\\ ⋮\\ w[N_f-1]\\ 0\\ ⋮\\ 0\\ w[-(N_p-1)]\\ ⋮\\ w[-2]\\ w[-1]\end{array}\right]}^{\mathrm{\top }}.$$

As highlighted in the notation, the (zero) padding elements are inserted at the end of the input time-series and in the center of the convolution kernel, in order to correctly align the two vectors.

Following this indexing convention, the convolution between $v[n]$ and $w[n]$ can be computed as

(2) $$(v_e⊛w_e)[n]=\sum _{m=0}^{N_e-1}{v}_m{w}_{e;(N_e+m-n)\mathrm{\%}{N}_e}.$$

We now define the vector with elements containing the result of the convolution and of the circulant convolution respectively as

$$\mathit{v}\ast \mathit{w}={\left[\begin{array}{c}(v\ast w{)}_0\\ (v\ast w{)}_1\\ ⋮\\ (v\ast w{)}_{N_s-2}\\ (v\ast w{)}_{N_s-1}\end{array}\right]}^{\mathrm{\top }}={\left[\begin{array}{c}(v\ast w)[0]\\ (v\ast w)[1]\\ ⋮\\ (v\ast w)[N_s-2]\\ (v\ast w)[N_s-1]\end{array}\right]}^{\mathrm{\top }},{\mathit{v}}_e⊛{\mathit{w}}_e={\left[\begin{array}{c}(v_e⊛w_e{)}_0\\ (v_e⊛w_e{)}_1\\ ⋮\\ (v_e⊛w_e{)}_{N_e-2}\\ (v_e⊛w_e{)}_{N_e-1}\end{array}\right]}^{\mathrm{\top }}={\left[\begin{array}{c}(v_e⊛w_e)[0]\\ (v_e⊛w_e)[1]\\ ⋮\\ (v_e⊛w_e)[N_e-2]\\ (v_e⊛w_e)[N_e-1]\end{array}\right]}^{\mathrm{\top }}.$$

From these definitions, we have that the elements of vector $\mathit{v}\ast \mathit{w}$ coincide with the first $N_s$ elements of ${\mathit{v}}_e⊛{\mathit{w}}_e$. Finally, we introduce the circulant matrix with columns derived by rolling the vector ${\mathit{w}}_e$ as

$${\mathit{W}}_e=\left[\begin{array}{ccccc}{w}_{0,0}& w_{0,1}& \dots & w_{0,N_e-2}& w_{0,N_e-1}\\ w_{1,0}& w_{1,1}& \dots & w_{1,N_e-2}& w_{1,N_e-1}\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ w_{N_e-2,0}& w_{N_e-2,1}& \dots & w_{N_e-2,N_e-2}& w_{N_e-2,N_e-1}\\ w_{N_e-1,0}& w_{N_e-1,1}& \dots & w_{N_e-1,N_e-2}& w_{N_e-1,N_e-1}\end{array}\right]=\left[\begin{array}{ccccc}{w}_{e;0}& w_{e;N_e-1}& \dots & w_{e;2}& w_{e;1}\\ w_{e;1}& w_{e;0}& \dots & w_{e;3}& w_{e;2}\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ w_{e;N_e-2}& w_{e;N_e-3}& \dots & w_{e;0}& w_{e;N_e-1}\\ w_{e;N_e-1}& w_{e;N_e-2}& \dots & w_{e;1}& w_{e;0}\end{array}\right].$$

The same indexing notation is used for all the other matrices. Clearly, the circulant matrix ${\mathit{W}}_e$ can be used to compute the elements of ${\mathit{v}}_e⊛{\mathit{w}}_e$ performing the vector-matrix product

$${\mathit{v}}_e⊛{\mathit{w}}_e={\mathit{v}}_e{\mathit{W}}_e.$$

From the top left quarter of the matrix ${\mathit{W}}_e$, we also define the sub-matrix

$$\mathit{W}=\left[\begin{array}{ccccc}{w}_{e;0}& \underset{\_}{w_{e;N_e-1}}& \dots & \underset{\_}{w_{e;N_s+1}}& \underset{\_}{w_{e;N_s}}\\ \overline{w_{e;1}}& w_{e;0}& \dots & \underset{\_}{w_{e;N_s+2}}& \underset{\_}{w_{e;N_s+1}}\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ \overline{w_{e;N_s-2}}& \overline{w_{e;N_s-3}}& \dots & w_{e;0}& \underset{\_}{w_{e;N_e-1}}\\ \overline{w_{e;N_s-1}}& \overline{w_{e;N_s-2}}& \dots & \overline{w_{e;1}}& w_{e;0}\end{array}\right]$$where elements of the convolution kernel associated with the relative past, below the main diagonal, are under-lined, while the ones associated with the relative future, above the main diagonal, are upper-lined to emphasize the interpretability of the defined operation. In the case of zero-padding, this results in

(3) $$\mathit{v}⊛\mathit{w}=\mathit{v}\mathit{W}.$$

Note that matrices ${\mathit{W}}_e$ and $\mathit{W}$, respectively, are of circulant and Toeplitz type. It is well known that matrix-vector products involving circulant or Toeplitz matrices can be efficiently computed using FFT (Strang, 1986). For several applications, this represents a convenient solution in terms of computational cost and time, as shown in the results section.

Targeting large-scale applications, we consider the case in which the same processing method is applied to multiple time-series. We so define the matrix $\mathit{V}$, which contains on each row one of the $N_r$ input time-series. Clearly, the matrix resulting from the matrix-matrix product $\mathit{V}\mathbf{W}$ will contain on each row the convolution between the corresponding time-series and the convolution kernel $w[n]$. However, the same result can be obtained by explicit iterative summation, as in Eq. (1), or as detailed in the following sections, using FFT. Finally, in order to deal with missing values in the original time-series, we define the matrix

$$\mathit{M}=\left[\begin{array}{ccccc}{m}_{0,0}& m_{0,1}& \dots & m_{0,N_s-2}& m_{0,N_s-1}\\ m_{1,0}& m_{1,1}& \dots & m_{1,N_s-2}& m_{1,N_s-1}\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ m_{N_r-2,0}& m_{N_r-2,1}& \dots & m_{N_r-2,N_s-2}& m_{N_r-2,N_s-1}\\ m_{N_r-1,0}& m_{N_r-1,1}& \dots & m_{N_r-1,N_s-2}& m_{N_r-1,N_s-1},\end{array}\right]$$containing the validity-mask of the matrix V, with elements value defined as

$$m_{i,j}=\{\begin{array}{cc}1& \mathrm{i}\mathrm{f}{v}_{i,j}\text{is a valid sample};\\ 0& \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}.\hfill \end{array}$$

Guidelines for convolution kernel definition

As anticipated, different time-series processing tasks can be accomplished by properly defining the convolution kernel in TSIRF. This section provides some details on how to set the convolution kernel based on the target application. All convolution kernel elements are supposed to be in the range $\epsilon <w[n]\le 1$, where $\epsilon$ is the machine accuracy previously mentioned (e.g., for double float $\epsilon =2.22\cdot {10}^{-16}$). Using an higher upper limit instead of 1 would have no impact on the final result, but it is a preferable choice to avoid potential overflows. For all the following application the value of $w[0]$ is set to 1 and, for sake of brevity, this will not repeated. As a general note, in case the application requires to respect the causality principle, it is necessary to do not use nonzero weights for the values associated with the relative future of the convolution kernel. A property of TSIRF is to be an interpolant, so the reconstructed values will be numerically bounded between the minimum and the maximum of the available ones.

As first application, we describe how to perform temporal aggregation by averaging a sliding window of images. If we want an aggregation factor of $N_a$ simply by setting $N_a$ consecutive elements of the convolution kernel to 1. For instance, these elements could be centered in $w[0]$ (included in the 1 element count). However, after the application of TSIRF, the output $\tilde{\mathit{V}}$ would still have $N_s$ columns as input data. To obtain the aggregated time-series it is necessary to keep one column every $N_a$. More conveniently, when using the Matrix backend, it is possible to remove the unused columns directly from W before performing the matrix-matrix, in order to save unnecessary computation. A graphical example and summary is reported in the first row of Table 2.

Table 2:

Summary and graphical representation of the convolution kernel setup guidelines for different applications.

Method	Convolution kernel shape	Description
Temporal aggregation by averaging		Set to 1 a number of consecutive elements of the convolution kernel equivalent to the aggregation factor $N_a$. Only keep one column out of $N_a$ in $\tilde{\mathit{V}}$ (or more conveniently of $\mathit{W}$).
Approximated linear interpolation		Using the convolution kernel described by Eq. (7) results in a linear interpolant of the available values, with properties similar to the standard piecewise linear interpolation.
Approximated most recent image propagation		Using the convolution kernel described by Eq. (8), when the machine precision is not too coarse nor $N_s$ too high, results in a numerical approximation of the most recently available value propagation.
Seasonally weighted average (linear)		The function is described by Eq. (9) and can be used to reconstruct time series with known seasonality of period $T_{seas}$ (in this case one year).
Seasonally weighted average (exponential)		Similar to the previous row, but also including an envelope attenuation to reduce land-cover changes propagation. The function is described by Eq. (10).

DOI: 10.7717/peerj.18585/table-2

We now describe the implementation different methods for time series reconstruction. Some of the propose methods do not exactly match their standard counterpart but are approximations of them. For instance, piecewise linear interpolation, a popular choice for simple time-series reconstruction (Potapov et al., 2021; Yin et al., 2017), can be approximated applying TSIRF with a convolution kernel described by the function

(7) $$w(t)=1-\left|\frac{t}{N_s{T}_s}\right|+\epsilon$$and represented in the second row of Table 2. The operation is indeed interpolant and linear, however not piecewise as usually implemented in literature. Compared to the standard piecewise version, this implementation does not require to iterate over the time-series to fine the two nearest available values. The vectorization of the operation generally reduce the computational cost and as shown in the results section also lead to improved reconstruction performance compared to the standard version.

Another widely used strategy is the most recently available value propagation to reconstruct missing values. Similarly to the linear interpolation, this method can be approximated with TSIRF. In particular, using the convolution kernel

(8) $$w(t)={10}^{\left|\frac{t}{N_s{T}_s}\right|lo{g}_{10}\epsilon },$$the weight associated with each available image decay exponentially with relative time distance, as shown in the third column of Table 2. If the numerical precision is not too coarse nor the number of samples too high, each weight can differ of order of magnitudes, propagating numerically the most recent available image.

An additional strategy for EO time-series reconstruction implemented within the TSIRF framework relies on the prior that most of the global land surface reflectance shows an yearly seasonal pattern (typically with one year periodicity). In particular, a higher weight is given to images with relative time distance nearer to integer multiples of one season period $T_{seas}$. From this comes its name Seasonally Weighted Average (SWA). An example of such a convolution kernel is represented in the forth row of Table 2. This can be implemented using a triangular wave function of period $T_{seas}$, that reads like

(9) $$w(t)=1-2|\frac{t}{T_{seas}}-\lfloor \frac{t}{T_{seas}}+\frac{1}{2}\rfloor |+\epsilon ,$$where $\lfloor \rfloor$ is the floor function.

Furthermore, to avoid potential propagation land-cover change along the time-series, it is preferable to also give higher priority to images that are more recent. This is done by subtracting a scaled modulo function from Eq. (9). To enhance the numerical impact of the weighting, like in the case of most recent image propagation, the function is used as an exponent of 10. Two input parameters are introduced to define the seasonal attenuation $A_{seas}$, and the envelope attenuation $A_{env}$, with the latter used to limit land-cover changes propagation. Both attenuations are interpreted in decibels (dB), leading to the final version of the weighting function

(10) $$w(t)={10}^{-2\frac{A_{seas}}{10}\left|\frac{t}{T_{seas}}-\lfloor \frac{t}{T_{seas}}+\frac{1}{2}\rfloor \right|-\frac{A_{env}}{10}\left|\frac{t}{N_s{T}_s}\right|}.$$

Note that in consideration of the numerical problems previously described, it is recommendable to respect the condition ${\log }_{10}\epsilon <-\left(\frac{A_{seas}}{10}+\frac{A_{env}}{10}\right)$. A graphical example of the convolution kernel is represented in the bottom row of Table 2, while a detailed interpretation of the method is shown in Fig. 2. On top-left the input time-series with missing samples that need to be reconstructed. In the bottom left of the scheme, the graph shows the weights of the convolution kernel $w[t]$ associated with ${\mathit{w}}_e$. In the central column, for each of the three missing values, the convolution kernel is shifted and overlapped with the original time series to highlight the method functioning. The right most column we find the reconstructed time-series.

Analysis ready and cloud-free Landsat bimonthly aggregates

The Landsat ARD provided by the GLAD team at the University of Maryland is one of the few globally consistent archives for historical time-series of normalized surface reflectance harmonizing the entire Landsat satellite collections (Potapov et al., 2020). However, the raw GLAD Landsat ARD images include cloud contamination and image artifacts that could propagate in derived product if used as direct input for modeling land-use (LU) or land-cover (LC). In addition, the storage size of the whole collection amounts to about 1.4 PB, a size that can make prohibitive its usage due to download time and required hardware. For these reasons, we produced a derived collection removing contaminated part of images, performing weighted temporal aggregation, and applying SWA for reconstructing the times-series associated with each pixel. The so-produced images can be used as direct input for modeling or to produce global time-series of biophysical indices. Stored as cloud optimized GeoTIFF (COG) on S3 systems, they result in an ARCO product that can also be straightforwardly used for geo-spatial analysis. The production pipeline, schematized in Fig. 3 and described in the following, represents a first large-scale application of TSIRF.

The input dataset GLAD Landsat ARD is a tiled multi-band 16-days interval composite of images with 23 images per year from 1997 to 2022, for a total of 598 GeoTIFF images per tile saved as unsigned 16-bit integer, and a total storage size of about 1.4 PB (Potapov et al., 2020). The GLAD Landsat ARD harmonizes all available images from the different Landsat satellites (Landsat 5 TM, Landsat 7 ETM+ and Landsat 8 OLI/TIRS) from 1997 to present, providing one of the best sources of historical 30 m resolution surface reflectance time-series. In addition to the seven reflectance bands, blue, green, red, near infrared (NIR), short-wave infrared 1 (SWIR1) and short-wave infrared 2 (SWIR2) and the thermal band, the GLAD Landsat ARD includes a detailed quality flag (QA) that classifies each pixel as: land, water, cloud, cloud shadow, topographic shadow, hill shade, snow, haze, cloud proximity, shadow proximity, other shadows or buffered proximity of the previously mentioned ones. From this flag a space-time not clear-sky mask was derived, including the points that are classified as cloud, cloud shadow, haze, cloud proximity, shadow proximity or other shadows. Those points are considered as gaps and are imputed with the following procedure.

From the 16-days time interval time-series, a bimonthly product was derived by performing weighted temporal aggregation. The aggregation was obtained using four images from the original product to produce one image in the aggregated product. Since each year has 23 images, the last image of the year, associated to the November-December time-frame, was produced aggregating the last three images of the relative year and the first images of the following year of the GLAD Landsat ARD (excluding November–December for the year 2022 for which only three images were used). To minimize cloud-induced artifacts in the time-series, a heuristic approach was employed, prioritizing images with lower cloud cover. For each image in the bi-month time-frame, a clear-sky fraction of the relative image was computed (e.g., in an image where one-fourth of the pixels are considered as gap the clear-sky fraction is 0.75). The aggregation was then performed as a weighted average between the available values, where each pixel was weighted by the clear-sky fraction of the associated image. If no valid observation was found, the pixel was considered as “no-data” value in the output. The computation was performed using the previously described temporal aggregation of TSIRF with $N_a=4$ and element-wise scaling of the values based on the clear-sky fraction. Finally, each time-frame/tile/band combination was saved in a separate GeoTIFF as unsigned 8-bit (byte) integer, with values ranging from 0 to 250 and associating the no-data values to 255. After compressing per chunks of 1,024 by 1,024 pixels with Deflate (Deutsch, 1996), the total storage size resulted in about 60 TB.

However, the aggregated product still contains several missing values in the images time-series. In particular, this is due to: (i) no available images during winter months in northern areas, (ii) cloud cover during the rain season in tropical areas, and (iii) a high revisit period during years in which fewer Landsat satellites are active, in particular when only Landsat 7 producing images affected by striping due to sensor malfunctioning was active (Hermosilla et al., 2015). To produce complete time-series, SWA was applied on top of the bimonthly aggregated product, for each band and each pixel. In particular, the $A_{seas}$ and $A_{env}$ parameters were optimized by grid-search optimization on the benchmark dataset (Consoli et al., 2023), consisting of 16-days time-series of NDVI images derived from the harmonized Landsat–Sentinel (HLS) collection from 2015 to 2022. The $A_{env}$ parameter was scaled to be consistent with the different time range of the two products. As a result of the optimization, we used $A_{seas}=45dB$ and $A_{env}=46dB$. To avoid reconstructions of missing values with future images, ${\mathit{w}}_f$ was left empty. In addition to respect causality, this strategy allows one to consistently update the product when images when more recent ones are available applying the time-series reconstruction only on them, without modifying the previous images. However, this implies that missing values could still be present in case no images are available at the beginning of the time-series. Since this generally impacts only few regions and for a few time frames, the choice on how to perform imputation is left to the final user. For instance, it is possible to apply SWA again using images from the future only for the imputation of the remaining missing values. In addition to the seven Landsat reconstructed bands, we also produced a clear sky flag specifying if the value was reconstructed or already available in the aggregated product. In case the value was already available in the aggregated product the flag is set to 250, while 255 represent again the no-data values. All other values of the flag are associated with reconstructed pixels. Also in this case, each time-frame/tile/band combination was saved and compressed, as for the aggregated product, leading to a total storage size resulting in about 100 TB of data.

As shown in the bottom of Fig. 3, the Landsat bimonthly reconstructed images can be used as input for LU/LC modeling, or alternatively, to produce remote-sensed biophysical indices. This was the case for the soil health monitoring data produced in Tian et al. (2024). The production pipeline of each index used the reconstructed images to calculate the index values. Then each bimonthly time frame is landmasked to remove sea water pixel, and finally aggregated in global (or regional) mosaics saved as cloud-optimized GeoTIFF (COG). Storing data in S3 cloud systems allows fast access and inspection of the data for users, providing an ARCO product.

Comparison of discrete convolution computational backends

The three strategies to perform the discrete convolution presented in “Efficient Computation of Circular Convolutions” were compared in terms of computational time as a function of the three parameters that could influence it. In particular, the involved parameters are the time-series length $N_s$, the number of processed time-series $N_r$, corresponding to the number of rows in the matrix V, and the size of the convolution kernel. For the latter, considering the assumption $N_p<N_s>N_f$, implying $N_w<2N_s+1$, and that the ratio between $N_w$ and $N_s$ impacts the sparsity of the matrix ${\mathit{W}}_e$, we introduce the sparsity parameter S. This parameter is used to create a benchmark experiment in which random data are used to create $\mathit{V}$, ${\mathit{w}}_p$, ${\mathit{w}}_f$ and $w_0$, where the size of ${\mathit{w}}_p$ and ${\mathit{w}}_f$ is $\lceil S{N}_s\rceil$, with S in the range $(0,1)$ and $\lceil \rceil$ the ceil function. With $S=1$, ${\mathit{W}}_e$ is a dense matrix. This experiment was run with each backend, Matrix, Summation and FFT, for each combination of $N_s=\left[1-2,200\right]$, $N_r=\left[1-40\right]$, and $S=\left[0.05,0.1,0.5,1.0\right]$, repeating for three iterations each and averaging the total computational time for each combination of parameters and backend (for a total of about three million runs). The results were used to find the fastest backend at each point of the parameter space. The computation was performed on a workstation HP^® Z840 with Intel^® Xeon^® CPU E5-2650 v4 @ 2.20 GHz $\times$ 48 and 128 GB of Hynix^® DDR4 Synchronous Registered (Buffered) 2,666 MHz RAM.

Since the results heavily depend on the used computational infrastructure and the used library, we combined the asymptotic behavior of each method and manual tuning of constant to fit functions that can describe the observed results. Those functions were used to produce Fig. 5 that helps users predict which backend works better depending on the application. The figure shows curves delimiting regions where each backend performs better, but without specifying the axis value, since this would be nonsensical for different computational infrastructures. However, since the asymptotic behavior should instead hold, shape and behavior of the regions should generally hold, and serve as baseline to decide which backend to use. The parameters $N_s$ and $N_r$ are on the $x$-axis and $y$-axis, respectively, and the three sets of curves are plot for S equal to $0.05$, $0.1$ and $1.0$. Region I is the region in which the summation backend is the fastest and shrinks toward the origin of the axis (low $N_s$ and low $N_r$) with increasing S. Region II, associated with the FFT backend, is located on the bottom right (high $N_s$ and low $N_r$) and expands with increasing S. Finally, for intermediate values of $N_s$ and high $N_r$ we generally fall in Region III, associated with the matrix backend. The selection of the best backend depending on the application is discussed in “Feasibility of Petabyte Scale Remote Sensing Time-Series Processing”.

Time-series reconstruction methods comparison

The benchmark dataset previously described was used to assess the effectiveness of time-series reconstruction methods under variable gap conditions by introducing artificial gaps into the extracted Landsat bimonthly time-series. These artificial gaps, also represented by “no data” values, were randomly inserted into the time-series to simulate real-world data loss. By varying the number of artificial gaps, the performance of different reconstruction methods under increasing gap density was evaluated, to test the robustness and project the error to the original gap fraction. The evaluation was done by comparing the imputed values and the original values corresponding to the artificial gaps. This comparison enabled the determination of the ability of each method to recover the original land cover dynamics despite the presence of both cloud-induced and artificial gaps. Note that, to have a fair comparison, the optimization of the method’s parameters was performed on a different dataset (Consoli et al., 2023).

We compared the performance of five different methods: (i) piecewise linear interpolation between available samples, in the following labeled Piecewise LI, (ii) TSIRF version of linear interpolation between available samples, as defined in Table 2, labeled TSIRF LI, (iii) TSIRF version of most recent images propagation smoothed with Savitzky–Golay filters, labeled MR-SG, (iv) SWA with the parameters used for production, labeled SWA, and (v) SG smoothed version of SWA, labeled as SWA-SG. The Piecewise LI was implemented imputing the value of the linear interpolation between the most recently available samples in the past and in the future of each missing value. For the other TSIRF-based methods the convolution kernel was set using the function described in the guidelines. For all the reconstruction methods, no weights for the future were used. However, for the SG a polynomial of order two was used to smooth a window of five samples centered on the third sample with zero-padding, making it non-causal (but only locally). An alternative could be to use non-symmetric SG by shifting the SG coefficients toward the ${\mathit{w}}_p$ vector. The SG filters were applied on top of SWA and MR results to evaluate the impact of smoothing.

We computed the Root Mean Square Error (RMSE) for each method in 10 experiments adding a varying fraction of artificial gaps. The results are reported for the Landsat NIR and red bands and for NDVI. For the NDVI, we also compared the impact of computing the NDVI first and then applying the reconstruction methods with the case in which the raw bands are first reconstructed and then used to compute the NDVI. The results are shown in Fig. 6. In the supplementary computational notebooks the same results are reported for all the bands. Projecting the curves toward the original gap fraction, for all the four cases, SWA and its smoothed version show a RMSE between 10% and 14% lower then the other methods. However, since the comparison only involves TSIRF-based or simple methods like piecewise linear interpolation, this comparison is not intended to comprehensive, but mainly used to compare alternatives for the Landsat aggregates reconstruction.

Furthermore, we compared the coefficient of determination ( $R^2$) for different methods and different bands, similarly to what is done in Hermosilla et al. (2015). Figure 7 shows the density scatter plot of the original NDVI against its reconstructed version using SWA on a the benchmark dataset with 10% of artificial gaps. The same plot is available for the other methods in the supplementary computational notebooks. In particular, the $R^2$ of each method in this experiment was 0.87 for Piecewise LI, 0.88 for TSIRF LI, 0.89 for MR-SG, 0.91 for SWA and 0.91 for SWA-SG. Notably, the ranking of the $R^2$ values is similar to the one of RMSE.

Land-cover and spectral dependency of time-series reconstruction error

To analyze the performance of the the time-series reconstruction for different LC, the RMSE was computed for each strata from the aggregated ESA CCI classes and the LC change samples for each band. In this case, only the results for a $10\mathrm{\%}$ of additional gaps are considered. Figure 8 shows the result for SWA with the previously described input parameters. The rows, representing the different LC classes, are ordered based on the total fraction of gap present in the input dataset, to see the correlation between RMSE and the number of gaps (except the LC changes class that is placed on the bottom of the graph). The error in each band is represented by a different symbol/color. The scatter plots and $R^2$ values are reported in Fig. 9. Finally, different error metrics in are computed for each band and summarized in Table 3, where also normalized versions of the RMSE and Concordance Correlation Coefficient (CCC) are included. The results were produced using the data available at https://zenodo.org/records/11150343 processed with the open-source code available at https://github.com/openlandmap/scikit-map/blob/feat_tsirf/. Note that the notebooks uses a pure Python implementation of TSIRF in order to facilitate reproducibility of the results. The production code based on hybrid Python/C++ code is also available on other branches of the same folders, but requires to be compiled before usage.

Landsat bimonthly cloud-free global ARCO

This section reports some results of the GLAD Landsat ARD processing and reconstruction. Figure 10 shows a zoomed-in RGB composite of the passage from four cloud screened images to a bimonthly aggregated and reconstructed images of a tile in Papua New Guinea, summarizing the passages for the production of a single time-frame of one tile.

Figure 11, on top, shows a region of the global mosaic of the NDVI cloud-free ARCO product relative to the bimonthly timeframe January–February 2022. From the global mosaics, three zoomed-in regions, corresponding to one 16th of the original tiles, are analyzed to show details and artifacts. In particular, site A is located in Turkey (30E, 37N), site B in Papua New Guinea (144E, 7S) and site C in Norway (13E, 61N). For each site, in the bottom of the figure, we find from left to right the Landsat bimonthly aggregated NDVI (with no-data in gray), the Landsat bimonthly reconstructed NDVI, the QA layer associated to the reconstruction and the high-resolution RGB composites from Google Hybrid, used to compare with some reference data. In these last images, we also locate with a red cross the central location used as query point to analyze the associated time-series reported in Fig. 12. The values already available in the Landsat bimonthly aggregated NDVI are plotted with orange crosses in function of time from the year 1997 to the year 2022, for a total of 156 samples per time-series. The Landsat bimonthly reconstructed NDVI time-series is instead plotted with a blue line. Due to the usage of only images from the past, it is possible to observe missing values in the reconstructed time-series in case of complete lack of images at the beginning of it. A green box highlights the time frame associated with the images shown in Fig. 11.

In Fig. 13 we compared the reconstructed Landsat images against Sentinel-2 L2A. We selected three locations in Europe with different latitudes and landforms. To each location we associated a different spectral band, in particular, the blue band with a location in Norway, the NIR band with a location in France and the SWIR2 band with a location in Italy (top-left and bottom-right corner coordinates are reported in the figure). The first location is $1^{\circ }\times 1^{\circ }$ in size, while the other two ${0.5}^{\circ }\times {0.5}^{\circ }$, in order to show additional details and artifact. Aggregated Sentinel-2 images to the time-frames (Q1) 2019/12/02-2020/03/20, (Q2) 2020/02/21-2020/06/24 (Q3) 2020/06/21-2020/09/24 (Q4) 2020/09/12-2020/12/01. In order to have most overall with the Sentinel-2 time-frames, Landsat images are referred to the bimonthly time-frames Jan–Feb, May–June, Jul–Aug, Sep–Oct of 2020. To test reconstruction performance, for the corresponding time-frame of the visualization was considered as gap for all pixels. Since NIR and SWIR2 values are quite skewed non linear color-maps are used for visualization. It is evident that the quality of the blue band of Landsat is not comparable with the one of Sentinel-2 images, while NIR and SWIR2 bands mostly match. In the Landsat NIR image for the Jan–Feb time-frame it is possible to noticed some stripes artifacts propagated during time-series reconstruction. The Sentinel-2 images used for comparison are available at https://stac.ecodatacube.eu (Witjes et al., 2023).

Finally, two complete global mosaics, corresponding to the January–February and July–August timeframes of 2022 for NDVI are available at https://github.com/openlandmap/scikit-map/blob/feat_tsirf/CATALOG.md. The values range from 0 to 250, where 0 corresponds to −1 NDVI values and 250 to +1 NDVI value. In “Discussion” we discuss the cause of the artifacts present in the data and possible solution to attenuate them that will be developed in future works. That the land-mask is derived from the global land cover maps produced form the GLAD group in Potapov et al. (2020), considering as not-land any pixel associated to a water-related class at least one time in the time-series from 2000 to 2020. Note that since the LC map is derived from the first version of the GLAD Landsat ARD the arctic regions and other smaller ones are not included in the masked and mosaiced data, but still available in the tiled products.

Feasibility of petabyte scale remote sensing time-series processing

The storage size of several collections of satellite images today reaches the petabyte scale. In addition to the Landsat collection mentioned above, as of August 2023, the Sentinel-2 archive amounted to a storage size of 20 PiB (Bauer-Marschallinger & Falkner, 2023; Radeloff et al., 2024). In addition, both missions keep generating about 0.5 PiB and 2 PiB of data per year, respectively (Frantz, 2019), without considering the distribution of different user-levels of the same product. Similar storage sizes are also required for satellite images used to monitor atmospheric conditions and solar activity, such as the ones produced by the Geostationary Operational Environmental Satellites (GOES) program. Indeed, as of July 2022, the GOES-16 and GOES-17 archives require more than 4.7 PB of storage size (Willett et al., 2023). While this does not represent a problem for most academic institution or big IT companies with high performance computing (HPC) infrastructure and programming know-how, the usage of such data is prohibitive for organizations with limited resources. Even more limited is the hosting of the data since it requires distributed storage systems and high-performance interconnect solutions. Platforms like Amazon Web Services (AWS), Microsoft Azure (Azure), and Google Cloud Platform (GCP) can help overcome hardware limitations. However, cutting-edge-derived products also require the usage of highly optimized and parallel software. Unfortunately, few open-source libraries provide a computational framework that can be used for such large-scale geospatial data applications. For these reasons, we developed $\mathsf{scikit}\text{-}\mathsf{map}$ (https://github.com/openlandmap/scikit-map), a Python library that uses ad-hoc developed C++ functions for a computationally intensive task by directly exposing the data structures between the two languages, zero-copy operations and performing most of the operations with high parallelism. The library includes the described implementations of TSIRF, both with Python and C++ backends. In summary, scikit-map tires to combine computational performance with flexibility for large scale EO data processing, focusing in particular to time-series reconstruction.

In the result section, we compared the computational performance of the three algorithmic backends to perform the circulant convolution operation. In the following sections, we discuss the reasons and implications of this result, following the Bachmann–Landau notation for the asymptotic analysis. The computational complexity of the matrix backend is dominated by the matrix-matrix product(s). Most of the available libraries for matrix operations are based on BLAS (Basic Linear Algebra Subprograms) that chose the best algorithm based on the computational architecture and sizes of the involved matrices for each operation. However, for most cases V is a rectangular matrix for which fast methods for the matrix-matrix product are convenient only for very large-scale and specific cases (Knight, 1995), the computational complexity of the selected algorithms is likely asymptotically bounded by $\mathrm{\Theta }(N_r{N}_s^2)$. For the FFT backend instead, we know that for $N_s\to \mathrm{\infty }$ the complexity of the convolution of each time-series scales asymptotically as $\mathcal{O}(N_s{\log }_2{N}_s)$ (and the same goes for $\lceil S\rceil$) (Flannery et al., 1992), while for $N_r\to \mathrm{\infty }$ as $\mathcal{O}(N_r)$. Finally, the summation backend is probably asymptotically bounded by $\mathrm{\Theta }(N_r{\lceil S{N}_s\rceil }^2)$. This analysis helps to understand the results presented in Fig. 5. For example, for $N_s\to \mathrm{\infty }$, the FFT backend scales as $\mathcal{O}(N_s{\log }_2{N}_s)$, while the matrix and summation ones as $\mathcal{O}(N_s^2)$, making the first the best choice for processing of time-series with many samples in time (Region II in the figure), as in the case of historical data from GOES. On the other hand, for $N_r\to \mathrm{\infty }$, all methods behave asymptotically as $\mathcal{O}(N_r)$. However, from empirical results, the matrix backend is generally the best option (Region III in the figure). In this case, the dominant factor comes not from computational complexity, but from the fact that the BLAS code is highly optimized, and probably caching optimization allows it to reuse the content of W for all time-series without reloading it from memory. As already mentioned, a scenario that falls into Region III is when the same processing is applied to several time-series, as for the case of a moderate and high spatial resolution EO dataset (MODIS, Copernicus Sentinel, and Landsat collections). Finally, for $S\to 0$, we clearly expect the summation backend to perform better, since its complexity reduces quadratically with S, while the computational complexity of the FFT backend only linearly and that of the matrix backend is not impacted at all. This defines the boundary of region Region I, which could be of interest when using Savitzky–Golay or other digital filters with small window length compared to time-series length.

The Python library SciPy also offer different backends to perform the convolution (https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.convolve.html). In particular, its backend full matches the backend summation in TSIRF, while the FFT one intuitively matches the FFT in TSIRF. However, there is no equivalent for the matrix backend, which was instead the most convenient option for our processing of the GLAD Landsat ARD archive. In addition, with the used software setup and computing infrastructure, even if using the same number of CPUs, we noticed that the computational efficiency of the SciPy library was not comparable to the one achieved by combining the C++ libraries Eigen (https://eigen.tuxfamily.org/) for algebraic operations and OpenMP (https://www.openmp.org/) for multi-threading. This is probably due to limitations imposed by the Python’s Global Interpreter Lock (GIL) and different implementation strategies.

Time-series reconstruction for EO data imputation

Even if in literature several works used purely spatial approaches for EO data imputation (Desai & Ganatra, 2012), in this work, we only focused on time-series reconstruction methods (including hybrid approaches). This is because purely spatial methods risk creating patches in the imputed images, reducing the actual spatial resolution of the final product. Yu et al. (2021) categorizes time-series reconstruction methods into three categories: local window, global window, and others. Differently, Siabi, Sanaeinejad & Ghahraman (2020) and others (e.g., Desai & Ganatra (2012), Gerber et al. (2018), Zeng, Shen & Zhang (2013)) divided the methods between spatial, temporal, and spatio-temporal. Following this categorizations, our framework would fall under the categories temporal and global window, since it only uses the temporal information but consider the whole time-series to reconstruct each missing sample.

The mentioned articles, also provides a comprehensive review of imputation methods for EO data. However, only few of them represents a good option for reconstruction of petabyte size historical EO time-series, and most of them are not openly available. Similarly to what is done in Julien & Sobrino (2018), to assent the performance of SWA, we decided to create a benchmark dataset based on the Landsat product that we presented in this work, since different time-series reconstruction methods could perform quite differently depending on the data set to which they are applied. Indeed, an intrinsic limitation of SWA, is that applying it to time-series where there is no seasonality (or known periodicity in general) would be actually equivalent to the numerical approximation of most recent available value propagation. As already mentioned, the benchmark dataset is openly available and we invite other groups to use it to compare the result with other methods.

From Fig. 6, projecting the error towards the original gaps fraction, we can observe that for all bands/indices, SWA, or its smoothed version, are the methods with best performance. In addition, we also wanted to study whether computing the NDVI before performing the reconstruction or doing the vice versa would have an impact on the final error. However, comparing the two graphs in the bottom, we can observe only a slight increase in error for the smoothed version of SWA. We are aware that this comparison is not comprehensive of space-time or data-fusion approaches that could potentially outperform SWA, however, since we did not find openly available code to compare such methods and we currently did not develop other methods that can be applied to large-scale dataset, we limited the analysis to what is reported. More comprehensive comparisons are available in Kandasamy et al. (2013) and Julien & Sobrino (2010).

Figure 8 instead analyses the reconstruction-error on different LC classes for each of the 8 Landsat spectral bands. In this case, only the result relative to SWA with a 10% of additional artificial gaps is considered. The aggregated ESA CCI classes of the benchmark dataset are listed on the left. To help taking into account that different land covers could present a larger fraction of missing values, we used this parameter to sort the classes, in addition to plot it on the right of the figure. Observing the figure, it can be noticed that for most of the LCs, there is an inverse correlation between the RMSE of each band end their wavelength. This is probably due to the fact that bands with shorter wavelengths are more sensitive to atmospheric contamination, introducing more noise in the time-series. Indeed, considering the spectral absorption profile of snow/ice shows, it is evident that more blue light is reflected compared to other LC classes, reducing the relative noise in the satellite sensor reading. It is also interesting to note that the error in samples with changes in land cover is not particularly high. This is probably due to the fact that because of the deforestation of the last years, many LC changes involve the transformation of forest pixels into crops or grassland pixels, for which the error is generally in the same order of magnitude. This also validates the efficacy of using a envelop attenuation for the weights of SWA to track LC changes.

Analysis-readiness of EO data and Landsat processing

The analysis-readiness of the EO data is a trending topic of debate (Baumann, 2024; Truckenbrodt et al., 2019). The definition of the term is per se controversial. In fact, anyone could argue: “Ready for what?”. Following the definition of the USGS (https://www.usgs.gov/landsat-missions/landsat-us-analysis-ready-data), their Landsat ARD products drastically reduce the magnitude of data processing required for direct use in monitoring and evaluating landscape change (Dwyer et al., 2018). The product includes a Quality Assessment (QA) layer that allows, for instance, to filter out cloudy pixels from the analysis. Frantz (2019) presents the FORCE (Framework for Operational Radiometric Correction for Environmental monitoring), designed to produce ARD Landsat and Sentinel data. The work highlights that the USGS Landsat ARD still amounts to prohibitive size data volume, in addition to suffering from poor spatial and temporal consistency, and attempts to address these issues. However, FORCE is only provided as a tool, and not as a global and historical consistent dataset. For this reason, in Potapov et al. (2020), the GLAD group presents the already mentioned GLAD Landsat ARD, which represents a unique product in terms of consistency and usability of the data, to the best of our knowledge.

Evidently, the cloud presence could impact for the production of biophysical indices maps such as normalized differences, albedo, gross and net primary productivity (GPP and NPP), or water vapor. In addition, for applications like LU/LC classification, most machine learning modeling techniques (e.g., random forests, artificial neural networks, support vector machines) could produce artifacts in the output map if the contaminated pixels are provided in input. An option is to label such pixels as no-data and eventually propagate the no-data value to the output. However, considering that, at the global level from 1997 to 2022, in the bimonthly aggregated Landsat product presented in this work, about 40.1% of the pixels are no data, propagating this to the output would mean to do not classify large areas. Yearly statistics of valid pixels fraction for each year is reported in Table 4, where is interesting to observe the increasing trend of images availability. A per-pixel count of the presence of data gaps in the Landsat bimonthly aggregated time-series is shown in Fig. 14, for which source COG data at 30 m resolution can be downloaded or visualized at https://github.com/openlandmap/scikit-map/blob/feat_tsirf/CATALOG.md. In addition to higher gap presence in cloudy regions, it is possible to observe patterns determined by scenes overlap. However, all these statistics also depends on the aggregation factor and on the used flags criteria to define a pixel valid or not. For instance, to consider snow a valid pixel or to aggregate the product quarterly would reduce the amount of gaps. In any case, the percentage of gaps are not negligible, especially in some regions. In fact, to produce some products derived from Landsat data, the GLAD group applies time-series reconstruction techniques such as piecewise linear regression (Potapov et al., 2021). This approach still requires to load and preprocess the time-series to use these images, again limiting the usage of the data for some user, in addition to potentially achieving poorer results in the time-series reconstruction when compared to other methods. Figure 15 compare reconstructed Landsat summer and winter images with averaged annual Sentinel-2B images and shows that our product still suffer from artifact presence, in particular in northern area or rain forests. In addition, it is possible to notice the Landsat 7 strips and the shapes of different scenes from which each image is produced (Fig. 11). Similar artifacts can be observed in Fig. 13. In future versions more efforts will be invested in reducing these limitations. Note that all the produced mosaics were landmasked using pixels classified as land-related classes, for at least one time frame, in the land-cover product described in Potapov et al. (2022). In the boundaries between water bodies and land, depending also on the timeframe, some not masked pixels could actually be derived from water reflectance and show heavy contrasts compared to land pixels.

An important first achievement in processing is the reduction in the storage space of the derived products. In particular, compared to the 1.4 PiB of the input collection, aggregated and reconstructed bimonthly only require 60 and 100 TB of storage, respectively. Derived products of global time-series with the same temporal and spatial resolution saved as COG require approximately 20 TB each. However, to achieve such compression, the temporal interval increased from 16-days to bimonthly, for some applications such detection of multiple corp cycle during one year, this could be a limiting factor. Nevertheless, considering that globally for several regions in some part of each year there are anyway no available images, this does not necessarily impact the analysis at global scale. Another factor impacting the compression is the storage of data as 8-bit unsigned integers instead of 16-bit ones. This could be relevant, for example, for the mapping of different vegetation species distributions or soil moisture. With these considerations we would expect a total compression factor of about 7.7 times, while we observe a compression factor of about 14 times. Part of the compression comes indeed from the storage of the images as 1,024 by 1,024 compressed chunks instead of 1 by 4,004 chunks. This emphasizes the importance of storing and compression strategies. In addition, we released an yearly aggregated version using the percentile 50 for each pixel, also available at https://stac.openlandmap.org. For applications that do not require a temporal resolution superior to one year, this represents the most compact way to access the data, in addition to being less affected by artifacts.

Another problem connected to the time-series reconstruction strategy can be observed in the time-series of the point queried at site C shown at the bottom of Fig. 12. In this case, it can be noticed that the complete absence of available images during winter months in northern regions brings to the propagation of values from summer into them. This result in the shown case in an NDVI time-series for those regions that resembles the one of a tropical area as the one shown for cite B. To overcome this issue, it could be possible to do not consider as gap pixels flagged as snow and ice. Finally, as anticipated, the application of SWA without using values from the past implies that if no values are available at the beginning of the time-series, missing values will still be present. This could be solved, for example, by using future images only for such samples and assigning them to ${\mathit{w}}_f$ lower values compared to ${\mathit{w}}_p$. All these consideration will be taken into account for the production of the next version.

Additional applications and future works

As for TSIRF, the Savitzky–Golay (SG) filters are a flexible tool that can be used for different applications by simply adjusting the input parameters. Through local polynomial fitting, SG filters smooth out noise while keeping the analytical signal features, thus, is widely applied for data smoothing and time-series reconstruction (Chen et al., 2004; Savitzky & Golay, 1964). By accurately calculating the elements of the convolution kernel of TSIRF, it is possible to combine several time-series reconstruction techniques with SG filters without the need to apply them in a second step, saving additional computational resources. As future development we plan to automatize the convolution kernel definition to integrate SG filter on top of other convolution kernels for arbitrary parameters.

Another approach to perform time-series reconstruction with TSIRF could be to set all the values of the convolution kernel using pure data-driven approaches. For instance, it could be possible to tune the convolution kernel using the gradient descent algorithm to minimize the RMSE on artificially generated data-gaps. Other loss functions or minimization algorithms could also apply. However, such an approach would lead to less interpretability about the application of the time-series reconstruction.

The flexibility of the framework, can also allow to extend it from purely temporal to a spatio-temporal approach (Liu et al., 2019). For instance, for each pixel time-series, it could be possible to create another time-series containing the value of the spatially nearest available pixel for each timeframe. These time-series can be concatenated to the extended ones in Eq. (4) (while maintaining the proper padding distance), and the convolution kernel defined to also include those values, possibly with lower weights for them compared to the original pixel time-series. A similar approach could also be used to perform data fusion. Indeed, concatenating to the original one a time-series from an harmonized product and then performing the convolution could allow the use of information from other sensors while keeping the same spatial resolution.

An option to improve quality of the Landsat ARD bimonthly product could indeed be the data fusion with the MODIS collection, which has lower spatial resolution but shorter revisit period (8-day). However, this approach requires the usage of time-series harmonized in space, time and variable spectrum. Then, a modified convolution kernel based on SWA can be used, in combination with the concatenated and padded time-series. In particular, the part of the convolution kernel associated with MODIS data can be defined with same seasonal end envelope attenuations as the Landsat one, but scaled by a factor $f<1$, such that higher priority will be given to Landsat values with the same seasonality and temporal distance from the value under reconstruction. This could help for instance to reduce the reconstruction artifacts in winter months in Northern regions. Since the result would still be a weighted average along the time-series that still involves 30 m resolution values from the Landsat data, by properly tuning the value of f, we expect that the presence of artifacts due to double resolution will be marginal. A combination of Landsat and MODIS satellite collections for data fusion was proposed in Luo, Guan & Peng (2018). However, their product is unfortunately not openly available, so it is also difficult to validate and reproduce.

Another aspect that will be explored in future version is the usage of different projection strategies to further reduce the storage size. As demonstrated in Bauer-Marschallinger & Falkner (2023), the projection used to store Sentinel-2 images increases about 33% the amount of information that they actually provide. Bauer-Marschallinger, Sabel & Wagner (2014) proposed a solution to the oversampling issues consisting of a girding system of seven different continental projections in order to minimize area deformation, from which the name “Equi7 Grid”, reducing the global mean oversampling to about 2%. In addition, since the area deformation is minimal, data store in Equi7 can be directly used as input for application where the a match between graph-distance and geodesic-distance is required (e.g., convolutional neural networks or hydrological modeling). However, this strategy has some drawbacks when working on a global scale. Since different projection systems are used, to visualize global mosaics requires different re-projection on each continent, with consequent computational over-head. In addition, applying spatial convolutions could be tedious in transitional regions, and padding strategies could be required. An alternative is to work with “Discrete Global Grid System” (DGGS) (Kmoch et al., 2022), where the spatial distribution of the data is obtained by associating each location with an equal-area cell from Earth surface tessellation. Nevertheless, this approach requires the usage of proper data-structure to store and process the data. In summary, both Equi7 and DGGS are elegant solutions for storing global EO data, but both require infrastructure and additional efforts to process the data.

Recent works targeting production of analysis ready geospatial images focus on usage multi-sensor data. This is done using different approaches, such as deep-features extraction (Zhou et al., 2023), multi-spectral data-fusion (Moreno-Martínez et al., 2020) and the creation of geospatial foundation models (Jakubik et al., 2023; Han et al., 2024). Deep-learning (DL) with convolutional neural networks demonstrated effectiveness in extracting features from raw surface reflectance data suitable for LU/LC classification (Dou et al., 2024). Applying unsupervised or transfer learning techniques on such models, it is possible to extract abstract features from a set of EO images that can vary in acquisition time, spectral band or even satellite sensor and use them for different modeling tasks. The advantages of such approaches consist in the possibility to combine several raw EO data sources and compress their useful information in few general purpose abstract features. On the other side, the usage abstract features lead to lack of interpretability of the features importance ranking of a model. In addition, the abstract features can not be used to directly compute well established biophysical indices (Montero et al., 2023). This implies the usage of modeling techniques to produce such indices, while using the reconstructed bands produced in this work they can be computed from simple arithmetic formulas. Overall, deep-features extraction could be a promising enabling approach for several applications in EO data classification. Similarly, foundation models such the one proposed in Jakubik et al. (2023) are innovative and flexible options for several geospatial applications, including ARD generation. In future works could be interesting to compare complete satellite images produced by foundation models with the ones obtained with classical imputation methods. Another front that could be explored is the training of a foundation model based on a set of reconstructed images selected by expert based visual interpretation, taking advantage of both approaches and potentially avoiding some model biases. As a final note, we want to highlight that the boundary in the definition of ML (or DL) and other data driven approaches can be blurred. Indeed, if as previously proposed the weights of TSIRF are all obtained by minimization of a loss function, TSIRF could be also defined as ML model.