## 4 Discussion  
Overall, we found that PX and TM were the best algorithms for phenology estimation, with the lowest errors (Table S1 and Figure 5). The algorithms had different performances when predicting SD or HD; we found that the TM algorithm was more suitable to predict sowing while PX was to predict the harvest. Packages based on the threshold method (such as PX and TM) showed higher consistency, except for GB (Figures S3 and S4). Large discrepancies were observed between the algorithms based on the derivative method (PF, GB, and DT using faster growth and senescence rate). Due to diverse methodologies, these inconsistencies are a source of uncertainty in detecting LSP [4; 61]. Thus, even the best-performing algorithms lack the precision to drive in-season decisions (best performances with RMSE between 10 and 20 days--Table S1). Still, these methods can provide trends and assess regional practice changes much quicker than surveys. For example, detecting sowing practices changes with the introduction of new crop cultivars, i.e., a cultivar that could be sown deeper and, hence, earlier in dry lands to capture early rains [62] or crop models that could be informed to evaluate large regions [63]. Although small differences were found, the remote sensing data sources that showed the best results were Sentinel-2 and MCD13 for most packages (Table S1 and Figure 5). Studies that have analyzed different data sources have observed that estimates of phenological metrics could be affected when the temporal [64] and spatial resolution [65] of the sensor changed. However, here, the TM package was the most robust and showed similar behavior and results regardless of the data sources (Table S1 and Figure 7a,b). On the other side, GB and DT seem to be more sensitive to the data source used (Figure 7c,d). \begin{table}
\begin{tabular}{c c c c c c} \hline \hline  & **CP** & **DT** & **GB** & **PF** & **PX** & **TM** \\ \hline Time (s) & 0.343 & 6.879 & 0.130 & 0.146 & 10.278 & 76.544 \\ \hline \hline \end{tabular}
\end{table} 
Table 2: Runtime for 100 NDVI time series for each package, with CP (CropPhenology), DT (DEA tools), GB (greenbrown), PF (phenofit), PX (phenex), and TM (TIMESAT). Figure 8: Dispersion between observed field data and sowing (**a**) and harvest (**b**) estimates for different numbers of seasons. Only the TM package and MCD13 data source were used as they showed the best fit (see Figure 5). The black line represents the 1:1 line. The use of different VIs had a stronger impact on the phenology extraction performance than the changing spatial and temporal resolution (Figure 4). This is probably related to the fact that EVI and EVI2 were developed to overcome problems with saturation, minimize atmospheric contamination, and reduce the bare-ground signal present in NDVI [44; 66; 67]. Estimating phenological metrics in agricultural regions with more than one cycle makes this process even more challenging [68; 69]. Double cropping, a common practice in the studied regions (soy and corn in the south region [40] and soy, corn, cotton, or non-commercial crops in the north [70]), generates confusion in the algorithms (Figure 8) because of the spectral pattern similarity of these crops [70]. Although TM has features that allow to identify the number of the crop cycle and set the desired cycle [71], it often ends up selecting the wrong cycle regardless--particularly harvest dates (Figure 8b). This problem could be solved by inserting the time series corresponding exactly to the period of soybeans grown in the region. However, this information is not available, and it would require a classification step that would make estimates even more challenging for large regions because soybeans could be the first or the second crop in the season. In addition, the validation of phenological metrics estimates for agricultural areas is still difficult because of the gap between remotely sensed data and those observed in the field [17; 72].