## Materials and Methods
### Methodology
The temporal gradient of a given spectral band, TGSB, at the \(i\)th (\(t_{i}\)) and \(j\)th (\(t_{j}\)) times are then generated to obtain the temporal changes of crops as follows:  
\[\text{TGSB }(t_{i},t_{j})=\text{SB }\big{(}t_{j}\big{)}-\text{SB }(t_{i}) \tag{1}\]  
The temporal gradients of VIs, TGVI, are also computed and considered as new features, which means for a given VI at the \(i\)th (\(t_{i}\)) and \(j\)th (\(t_{j}\)) times, a set of VI gradients in time sequence is generated as follows:  
\[\text{TGVI }(t_{i},t_{j})=\text{VI }\big{(}t_{j}\big{)}-\text{VI }(t_{i}) \tag{2}\]  
Figure 5: The missing data retrieval method by the K-nearest method of cloudless pixels in the feature space. \begin{table}
\begin{tabular}{c c c c} \hline \hline
**VIs. No.** & **Index** & **Formula** & **Reference** \\ \hline
**1** & Normalized Difference Vegetation Index & \(\text{NDVI}=\frac{(\text{NIR-R})}{(\text{NIR-R})}\) & [34] \\ \hline
**2** & Simple Ratio & \(\text{SR}=\frac{\text{NIR}}{\text{NIR}}\) & [35] \\ \hline
**3** & Stress-related Vegetation Index 1 & \(\text{STVI}=\text{MIR}\times\Big{(}\frac{\text{R}}{\text{NIR}}\Big{)}\) & [36] \\ \hline
**4** & Stress-related Vegetation Index 3 & \(\text{STVI}3=\frac{\text{NIR}}{(\text{R+MIR})}\) & [37] \\ \hline
**5** & Stress-related Vegetation Index 4 & \(\text{STVI}4=\text{NIR }-\Big{(}\frac{\text{R+MIR}}{\text{NIR+MIR}}\Big{)}\) & [37,38] \\ \hline
**6** & Enhanced Vegetation Index & \(\text{EVI}=25\times\Big{(}\frac{\text{NIR-R}}{\text{NIR+R+259+1}}\Big{)}\) & [9,39] \\ \hline
**7** & Modified Soil Adjusted Vegetation Index & \(\text{MSAVI}=\frac{(\text{NIR+1-v}^{\prime}(\text{NIR+1-v}^{\prime}-8(\text{NIR -R}))}{2}\) & [40] \\ \hline
**8** & Modified Normalized Difference Vegetation Index & \(\text{MNDVI}=\text{NDVI}\times\Big{(}\frac{\text{NIR}_{\text{NIR}_{\text{NIR}_{\text{NIR}}}}}{(\text{ SWIK}_{\text{NIR}_{\text{NIR}}}})\) & [41] \\ \hline
**9** & Soil Adjusted Vegetation Index & \(\text{SAVI}=\frac{(\text{NIR-R})}{(\text{NIR+R+1})}\times(1+L)\) & [42] \\ \hline \hline \end{tabular}
\end{table} 
Table 3: Remote sensing vegetation indices used in this study. **Step 5**: The RF algorithm is implemented to explore the achieved features, including spectral bands at different times, VIs, and their temporal gradients, to obtain an ensemble of decision trees. As mentioned, eight images from path 162 and two images from path 163, having less than 20% cloud cover, are selected. Therefore, for each pixel, a feature vector is obtained that includes spectral bands and VIs, along with their temporal gradients. To investigate the effect of multi-temporal image classification issues by the RF classification, different scenarios are considered according to Table 4.  
The RF algorithm can generally handle high dimensional data and use a large number of trees in the ensemble. In RF, each tree in the forest is built and tested independently from other trees. During training, each tree receives a new bootstrapped training set generated from the original training set by subsampling with a replacement [43]. Those samples which are not included during the training of a tree are called out-of-bag (OOB) samples of that tree. These samples can be used to compute the out-of-bag-error (OOBE) of the tree, as well as the ensemble, which is an unbiased estimate of the generalization error [44; 45; 46]. In the experiments, replacement and OOB are set to 15 times and 33%, respectively. First, the optimum number of training sample size along with the number of trees are investigated. Then, the effect of each other items in Table 4 issue is assessed to obtain an appropriate model for the multi-temporal classification. Decision-makers aim to measure the crop area from remotely-sensed data in earlier seasons with appropriate accuracy. Therefore, according to Table 5, based on the desired classes, two cases of classification are determined by the decision-makers.