## 2 Data and Methods
### Simulation Experiment
#### 2.2.3 Building Regression Model  
We used Python as a programming language and performed the calculations on Amazon Web Services (AWS, Seattle, WA, USA). We processed 20,438 Landsat scenes for this study. Apart from satellite imagery processing, calculations associated with the feature selection process (from a few thousand features) were among the most computationally intensive. Figure 2 demonstrates the data flow for the dataset creation process. The formulas for calculating weighted parameter values and a regression model are as follows:  
\[F^{j}_{weighted}=\sum_{i=1}^{N}w_{i}*F^{j}_{i} \tag{7}\]  
\[NEE=H(F^{1}_{weighted},F^{2}_{weighted},\cdots,F^{M}_{weighted}) \tag{8}\]  
where \(w_{i}\)--weight value for the \(i\)-th grid cell (see Figure 3); \(F^{j}_{i}\)-- value of \(j\)-th predictive feature for \(i\)-th grid cell; \(F^{j}_{weighted}\)--weighted sum for \(j\)-th feature; \(NEE\)--eddy covariance data; \(N\)--number of grid cells; \(M\)--number of features; and \(H\)--non-linear operator. Figure 3 shows an example of original features calculated from Landsat data and results of LPDM calculations of the footprint for the up-scaling process. The weighted sum for each parameter may be interpreted as the dot product of the feature and footprint map. To train the model, we match the resulting input dataset (Table 1) with eddy covariance data (Table 2). We selected an algorithm based on the kernel method kernel ridge regression (KRR) [46; 47]. This method is similar to the more widely known class of support vector machine (SVM) methods [48]. However, it uses a different loss function. The KRR method uses linear least squares with l2-norm regularization, meaning it has scale factor for squared feature coefficients as an additional term in the loss function. This allows for creating non-linear dependencies between the features according to the selected kernel and then searching for the optimal coefficients of these dependencies with control of the regularization process. The result of solving such a problem is a hyperplane in the multidimensional feature space that describes the dependencies of the target values on the selected features. This method works well with small datasets, which is the case in this study and is also faster than the original support vector regression (SVR). KRR has three main hyperparameters: \(\alpha\)--the regularization parameter of the model (the higher the value, the stronger the regularization); kernel--kernel interaction between the features included in the model (see Equations (9) and (10)); \(\gamma\)--internal kernel parameter (see Equations (9) and (10)). \[k(x,y) =\mathrm{e}^{-\gamma\|x-y\|_{1}} \tag{9}\] \[k(x,y) =\mathrm{e}^{-\gamma\|x-y\|_{2}} \tag{10}\]  
where \(e\)--exponent; \(x\), \(y\)--input vectors; \(\|x-y\|_{1}\)--Manhattan distance between the input vectors; and \(\|x-y\|_{2}\)--Euclidean distance between the input vectors. Another critical aspect of preparation is the feature selection process. As mentioned in Section 2.2.2, in addition to the primary remote sensing and weather parameters, 87 additional features in total were prepared. In addition, we used products between all meteorological and remote sensing features that led to a few thousand features. For each IGBP class, a small set of features has to be selected to allow optimal prediction of NEE at the stations. The exact solution for finding the optimal combination of features can be obtained by a complete enumeration of all possible variants. The complexity of this approach is \(2^{n}\), where \(n\) is a number of features, so the approach becomes inapplicable for a large number of features. The utilities in the package machine learning extensions (Mlxtend) [49] were used to select the features. The selection process was divided into two stages. The first stage estimates mutual information for each feature in the list. A standard method from the scikit-learn library was used [50], which allows us to reduce the number of features to select.