## 2 Materials and Methods
### N-Way Partial Least Squares  
In general, the PLS approach is particularly useful when one or a set of dependent variables needs to be predicted by a (very) large set of predictor variables (or time-series) that are strongly cross-correlated [25]. The strength of N-PLS is that it summarizes all latent information from a large 3-way dataset of object variables (X) and relates it to a dependent variable (y) using a relatively low number of parameters, which makes the prediction more robust [26]. This situation leads to an analysis that is able to extract the maximum information possible from samples measured at different times based on cubic structure. For these remote sensing data, the 3-PLS1 method [17] was considered here to generate a three-way array \(\textbf{X}=\mid\mid X_{i,j,k}\mid\mid\). The first of the three dimensions of cube X corresponds to the objects (\(I\)), the second to time (\(J\)) and the third to spectral bands (\(K\)). Thus, the data were organized in a three-way array of independent variables \(\textbf{X}\) (\(I\times J\times K\)) and a response vector y of size (\(I\times I\)) that is defined by the objects dimensions (Figure 2). The N-PLS regression was used to relate the three-way array \(\underline{\textbf{X}}\) (remote sensing data) to the response vector **y** (ground truth data). As [PERSON] [25] specified, PLS regression performs a simultaneous decomposition of \(\underline{\textbf{X}}\) and **y** by means of a set of latent variables that explain as much as possible of the covariance between \(\underline{\textbf{X}}\) and **y**. As [PERSON] et al. [27] detailed, if **X** (\(I\times I\)K) is a properly unfolded two-way form of a three-way array \(\underline{\textbf{X}}\) (\(I\times I\times K\)), the 3-PLS1 method can be written in mathematical form (Equation 1), where **X** and **y** are cantered along the first dimension, \(I\). \[\textbf{S}=\textbf{X}\textbf{W}\;\;y=\textbf{S}\textbf{b}+\textbf{r}. \tag{1}\]  
Estimating the weight array **W**, allowing for the scores in **S** to be expressed directly in terms of the **X**-columns, is the essential part of the method [28]. The regression coefficients **b** (Equation (1)) can be estimated afterwards by a least square approach [27]. The vector **r** in Equation (1) presents the part of **y** not explained by the model. For the estimation of **W**, the algorithm for three-way \(\underline{\textbf{X}}\) and a single response **y** proposed by [PERSON] [28] was considered. If the initial values of \(\textbf{y}_{a}\) (\(a\) is the latent variable counter) are set to the original values \(\textbf{y}\) (\(a=0\) and \(\textbf{y}_{0}=\textbf{y}\)), the algorithm can be summarized as follows [27]:  
1. Compute the reshaped covariance matrix \(\tilde{\textbf{Z}}=y_{a}^{T}\textbf{X}\) (\(I\times I\)K). Figure 1: Production of three-dimensional array from satellite imagery time-series data. Figure 2: (**a**) Representation of three-way array (\(\underline{\textbf{X}}\)) and (**b**) response vector **y** to compute a tri-linear PLS1 procedure. Define the first singular weight vectors \(\mathbf{w}_{a}^{I}\) and \(\mathbf{w}_{a}^{K}\) from \(\mathbf{Z}\): \([\mathbf{w}_{a}^{I},\mathbf{w}_{a}^{K}]=\text{svd}(\mathbf{Z},1)\). From hence, store them as additional columns in separate weight arrays \(\mathbf{W}^{I}=[\mathbf{w}_{1}^{I}\,\dots\,\mathbf{w}_{a}^{I}\,]\) and \(\mathbf{W}^{K}\,[\mathbf{w}_{1}^{K}\,\dots\,\mathbf{w}_{a}^{K}]\). 3. Compute the weight array by a Khatri-Rao product [27; 29] so that \(\mathbf{W}=\mathbf{W}^{K}\otimes\mathbf{W}^{I}\). The Khatri-Rao product defined as a column- wise Kronecker product of the two matrices [27] was used to overcome constraints related to multiway arrays decomposition. 4. Calculate \(\mathbf{S}=\mathbf{\chi}\mathbf{W}\). 5. Calculate the regression coefficients regressing \(\mathbf{y}\) on \(\mathbf{S}\) as \(\mathbf{b}=(\mathbf{S}^{T}\mathbf{S})^{-1}\,\mathbf{S}^{T}\mathbf{y}\). 6. Calculate the residuals \(\mathbf{r}=\mathbf{y}-\mathbf{S}\mathbf{b}\). 7.