Developing a farm-level productivity indicator for coffee and cocoa via NDVI time series and National Agricultural Survey data in Peru

Developing a farm-level productivity indicator for coffee and cocoa via NDVI time series and National Agricultural Survey data in Peru | Research Square window.SnipcartSettings = { analytics: { enabled: false } }; (function() { var accessVector = localStorage.getItem('access_vector') || ''; window.dataLayer = window.dataLayer || []; if (accessVector) { window.dataLayer.push({ user: { profile: { profileInfo: { snid: accessVector } } } }); } })(); (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0],j=d.createElement(s),dl=l!='dataLayer'?'&l='+l:'';j.async=true;j.src='https://www.googletagmanager.com/gtm.js?id='+i+dl;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-K279D39R'); Browse Preprints In Review Journals COVID-19 Preprints AJE Video Bytes Research Tools Research Promotion AJE Professional Editing AJE Rubriq About Preprint Platform In Review Editorial Policies Our Team Advisory Board Help Center Sign In Submit a Preprint Cite Share Download PDF Research Article Developing a farm-level productivity indicator for coffee and cocoa via NDVI time series and National Agricultural Survey data in Peru Robinson López-Monzon, José Antonio Otoya-Barrenechea, Juan Leonardi Romero-Vasquez, and 1 more This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-9464917/v1 This work is licensed under a CC BY 4.0 License Status: Posted Version 1 posted You are reading this latest preprint version Abstract Context: Farm-level productivity indicators are still limited in Peruvian agricultural statistics, particularly for smallholder perennial crops such as coffee and cocoa, where strong within-region heterogeneity and complex management conditions make productivity difficult to measure consistently. Moreover, increasing access to satellite imagery creates new opportunities to complement survey-based agricultural information with spatially explicit indicators. Aims: This study aimed to develop and validate a farm-level indicator of productive capacity for smallholder coffee and cocoa systems in Peru by integrating National Agricultural Survey (ENA) microdata with Sentinel-2 normalized difference vegetation index (NDVI) time series. Methods: ENA 2023–2024 microdata were linked to georeferenced farm polygons and cloud-filtered 10 m Sentinel-2 normalized difference vegetation index (NDVI) imagery aggregated into 12 bimonthly composites. NDVI snapshots and temporal features were derived and used first in ridge and random forest regression models for continuous yield prediction and then in gradient boosting ensemble models for ordinal classification into low-, medium-, and high-productivity classes. A simplified decision-tree indicator was subsequently derived for operational use. Key Results: Continuous-yield regression performed poorly (R² ≈ 0; RMSE ≈ 340 kg ha⁻¹), indicating that the NDVI alone was insufficient for reliable parcel-level yield estimation. In contrast, ordinal classification performed well, reaching 0.73 validation accuracy, 0.85 evaluation accuracy, and Cohen’s κ = 0.75. The simplified indicator achieved lower but still useful performance (56.7% accuracy; κ = 0.31). Conclusion: Survey microdata and NDVI time series can be combined to generate an operational, farm-level indicator of productive capacity for smallholder perennial systems, even when continuous yield prediction remains weak. Implications and Impacts: The proposed indicator provides a practical way to translate satellite time series and survey microdata into an operational tool for official statistics. In coffee and cocoa systems, it can support stratified sampling, improve the review of anomalous survey records, and guide more spatially explicit interventions, with potential relevance for other contexts after local calibration. remote sensing yield mapping multispectral time series ensemble models official statistics tropical perennial crops Figures Figure 1 Figure 2 Figure 3 Figure 4 Figure 5 Figure 6 Figure 7 Figure 8 KEYPOINTS/HIGHLIGHTS This study develops and validates a farm-level indicator of productive capacity for smallholder coffee and cocoa systems in Peru by integrating National Agricultural Survey (ENA) microdata with Sentinel-2 normalized difference vegetation index (NDVI) time series data. The results show that NDVI-based models are insufficient for reliable, continuous yield prediction at the parcel level, as regression performance remained weak (R² ≈ 0; RMSE ≈ 340 kg ha⁻¹). Reframing the problem as an ordinal classification into low-, medium-, and high-productivity categories substantially improved model performance, reaching 0.73 validation accuracy, 0.85 evaluation accuracy, and Cohen’s κ = 0.75. Although less accurate than ensemble models, the simplified decision tree indicator offers a more transparent and operationally feasible tool for implementation in agricultural statistical systems. The proposed indicator is a practical tool for supporting agricultural decision-making, particularly in improving sampling strategies, strengthening data quality control, and helping target extension and technical assistance more effectively. IMPACT: This study demonstrates that Sentinel-2 NDVI time series and official agricultural survey microdata can be combined to produce an operational indicator of productive capacity for smallholder coffee and cocoa systems in Peru. While developed with the microdata of the National Agricultural Survey (ENA) of Peru, the methodology also offers a potentially useful framework for other contexts seeking to make spatial variability actionable for sampling design, record review, and targeted support, subject to local calibration and validation. INTRODUCTION Agriculture remains a strategic sector for Peru, both for its contribution to the gross domestic product (GDP) and for its role in employment, territorial cohesion and rural livelihoods (World Bank, 2023 , 2024 ). However, the lack of objective and comparable indicators of productivity at the farm level limits the capacity to target public policies, technical assistance and financial services to the areas and producers where they are most needed (World Bank, 2017 ; Gil et al., 2019 ). In practice, most decisions are still based on averages calculated at the regional or domain level, which masks the strong microterritorial heterogeneity observed within the same crop and year. Peru has made significant progress in agricultural statistics through the National Agricultural Survey (ENA), which provides key information on area, production, yields and inputs for major crops (MIDAGRI, 2022). Nevertheless, estimating productive capacity at the farm level remains challenging, especially in perennial and smallholder-dominated systems such as coffee and cocoa, where plant density, age structure, management and environmental conditions vary widely even within small areas. In parallel, the availability of freely accessible satellite imagery and cloud-computing platforms has opened new possibilities for the use of remote sensing to monitor crop vigor and yield with high spatial and temporal resolutions (Flores-Anderson et al., 2023 ; Wuepper et al., 2025 ). Vegetation indices derived from multispectral images, particularly the normalized difference vegetation index (NDVI), have been consistently associated with crop development and yield in a wide range of systems (Groten, 1993 ; Labus et al., 2002 ; Wang et al., 2005 ; Fuentes et al., 2025). In coffee and cocoa, recent studies have shown that the NDVI and related indices can help characterize yield variability, stress patterns and climatic gradients (Bernardes et al., 2012 ; Nogueira et al., 2018 ; Martello et al., 2022 ; Atalaya-Marín et al., 2024; Íñiguez Freites et al., 2025). However, the application of NDVI-based models in official operational statistics poses specific challenges. First, yield information in surveys such as ENA is collected through farmer-reported production and area, which may introduce measurement error and bias (FAO, 2017 ). Second, perennial crops exhibit complex phenological patterns and alternate bearing and management cycles that are not easily summarized by single-date NDVI snapshots (Anyimah et al., 2021 ; Manoel et al., 2024 ). Third, statistical offices need simple and interpretable indicators that can be integrated into existing workflows rather than black-box models that are difficult to communicate and maintain (Wuepper et al., 2025 ). The objective of this study is to develop and validate an annual farm-level indicator of productive capacity for smallholder coffee and cocoa systems in Peru through the integration of ENA microdata and Sentinel-2 NDVI time series. To achieve this goal, a harmonized database linking individual farms with satellite-derived NDVI metrics is first constructed; the capacity of these time series to predict continuous yield is then evaluated; then, ordinal classification models are developed to group farms into low-, medium-, and high-productivity categories; and finally, a simple and interpretable decision-rule indicator is derived for potential use within the Peruvian agricultural statistical system. MATERIALS AND METHODS Study area and data sources This study focused on smallholder coffee and cocoa systems located in the eastern Andean slopes of Peru, particularly in the regions of San Martín, Amazonas and Ucayali, which constitute a large share of national production. Three main data sources were integrated. First, microdata from the National Agricultural Survey (ENA 2023–2024) was used, as this survey provides farm-level information on cultivated area, production, management practices, and structural characteristics (INEI, 2025 ). For coffee and cocoa, plot-level data on harvested area (ha) and production (kg) were extracted, and yields (kg ha⁻¹) were calculated as the ratio of production to area. Records with missing or inconsistent values for area or production were excluded in accordance with ENA quality-control criteria. Second, georeferenced polygons corresponding to each sampled farm were obtained from ENA and used to aggregate satellite-derived NDVI data at the farm level. Third, NDVI time series were generated from Sentinel-2 Level-2A surface reflectance products at a 10 m spatial resolution, using bands 4 (red) and 8 (near-infrared). These data were processed for the 2023–2024 agricultural season over a 12-month period and aligned with ENA yield reporting. Cloud- and cloud shadow-contaminated pixels were removed via quality flags and conservative thresholds, and the remaining observations were summarized into 12 bimonthly normalized difference vegetation index (NDVI) composites by calculating the median number of cloud-free pixels (Flores-Anderson et al., 2023 ). Data integration and analytical sample For each sampled farm, the farm polygon was overlaid on the series of NDVI composites, and zonal statistics were calculated via the median NDVI value at each bimonthly interval. Farms with very small coffee or cocoa cultivation areas relative to the total polygon, as well as those with insufficient cloud-free observations, were excluded from the analysis. After applying all the consistency checks, the final analytical sample included several hundred coffee and cocoa farms. In the case of coffee, the yield distribution was strongly skewed, with most observations concentrated below 500 kg ha⁻¹ and relatively few above 1,000 kg ha⁻¹, which is consistent with previous evidence on smallholder productivity in the Peruvian Amazon (World Bank, 2017 ; Morales et al., 2022 ). For cocoa, yields were more concentrated at approximately 600 kg ha⁻¹, with very few observations exceeding 1,500 kg ha⁻¹. These distributions are illustrated by histograms of coffee (Fig. 1 ) and cocoa yields (Fig. 2 ), which show a marked right-skewness and a small number of high-yield outliers. The planting density distributions derived from the ENA confirmed substantial deviations from the theoretical or recommended densities; thus, the plant density was excluded as a direct predictor. Feature Engineering Two families of NDVI-based features were constructed. The first consisted of snapshot metrics , defined as the NDVI values observed during specific bimonthly periods corresponding to key phenological stages of coffee and cocoa, such as vegetative growth, flowering, and grain filling. For each crop, these periods were identified on the basis of agronomic calendars and local expert knowledge, and the median NDVI value of each farm polygon was retained for the selected dates. The second family consisted of derived temporal features calculated from the full 12-step NDVI series. These included the annual maximum and minimum NDVI, the maximum minus minimum (NDVI) amplitude, the coefficient of variation, the interquartile range, the selected percentiles, and the cumulative NDVI, which were used as approximations of the area under the NDVI curve over the agricultural season. The full set of features was standardized and screened to avoid redundant predictors and reduce collinearity. A correlation matrix of NDVI-derived features was used to guide feature selection and aggregation, and a feature importance ranking from preliminary models helped identify the most informative predictors. Table 1 provides a detailed description of all the dependent, explanatory and control variables used in the analysis, including their roles, sources (ENA vs. Sentinel-2) and units. Table 1 Definitions, roles and sources of variables used in the models. Variable Role Operational definition/Source ENA Fields/Sentinel Unit Scale Procedure REND Dependent Total production/area P117_SUP_ha; totales/destinos Cap. 200AB (P220_*, en kg) kg/ha Ratio Regression and binning, classification RENDPLANT Dependent REND/PLANTHA P126 (plantas), P117_SUP_ha kg/plant Ratio Regression; sensitivity to biases PLANTHA Control/Explanatory P126/P117_SUP_ha ENA plants/ha Ratio Control; quality filters REGION Control Geographical Region REGION/NOMBREDD — Nominal Fixed effects/blocks STARTMES Control Month when harvest starts P124_MES month Ordinal Phenological alignment ndvi_cv_all Explanatory (global) std(NDVI)/mean(NDVI) 2023–2024 Sentinel-2 (zonal stats) adim. Ratio Key feature; permutation importance ndvi_mean_all Explanatory(global) Mean NDVI 2023–2024 Sentinel-2 adim. Interval Feature; control for baseline status ndvi_amp_all Explanatory(global) max − min NDVI 2023–2024 Sentinel-2 adim. Interval Seasonality signal AUC_main_23/24 Explanatory(annual) Area under the curve (main window) Sentinel-2 NDVI· bimonthly Interval Within-year phenology auc_ratio_main_minor_23/24 Explanatory(annual) AUC_main/AUC_minor Sentinel-2 adim. Ratio Balance of peaks rise_pre_23/24, fall_harv_23/24, rebound_23/24 Explanatory(annual) Slopes by window Sentinel-2 ΔNDVI/Δbimonthly Interval Vigor and recovery skew_23/24, kurt_23/24 Explanatory(annual) Shape of NDVI distribution Sentinel-2 adim. Interval Cycle stability delta_ Explanatory(interannual) 2024 − 2023 for each feature Sentinel-2 according to metric — Interannual change NDVI* (selectas) Explanatory (snapshot) NDVI in critical bimonthly periods Sentinel-2 adim. Interval Only if phenologically justified To document collinearity patterns and guide feature selection, a correlation heatmap of all NDVI-derived variables was computed, which is presented in Fig. 3 . Continuous yield modeling NDVI time series were evaluated to validate whether continuous yields (kg ha⁻¹) can be predicted at the holding level. Two supervised regression models were trained: (i) ridge regression , a linear model with L2 regularization that handles correlated predictors while preserving interpretability; and (ii) random forest regression , an ensemble of decision trees that are able to capture nonlinear relationships and interactions. The dataset was randomly split into training (60%), validation (20%) and evaluation (20%) subsets, which were stratified by yield quantiles. The hyperparameters (e.g., the regularization strength for ridges and the number and depth of trees for the random forest) were tuned via cross-validation on the training and validation sets. Model performance was evaluated on the held-out evaluation set via the root mean squared error (RMSE), mean absolute error (MAE) and coefficient of determination (R²). Table 2 summarizes the performance of the continuous regression models, reporting the RMSE, MAE and R² for both ridge and random forest. Table 2 Performance metrics (RMSE, MAE, R²) of the ridge and random forest regression models for predicting continuous yields Método RMSE MAE R2 Ridge 339.4397 243.262 -0.0524 RandomForest 346.477 235.234 -0.0965 Ordinal classification of productive capacity Given the weak performance of the continuous models, the problem was reformulated as an ordinal classification task with three mutually exclusive categories of productive capacity: low, medium, and high, as illustrated in Fig. 4 . Cutoff points were established via ENA yield data and agronomic criteria so that each category would represent meaningful productivity thresholds for policy and extension purposes, including very low yields well below regional averages, intermediate yields, and top-performing farms. Gradient boosting ensemble models were then trained via CatBoost and LightGBM, both of which are well suited to structured tabular data. Multiclass and ordinal classification approaches were tested. The same training, validation, and evaluation strategies were applied, and model performance was assessed via overall accuracy, class-specific precision and recall, and Cohen’s κ to account for agreement occurring by chance. Confusion matrices were also examined to identify misclassification patterns, particularly between the medium- and high-productivity categories. Table 3 presents the main performance indicators of the best ensemble models and the simplified indicator, including overall accuracy, balanced accuracy and κ. Table 3 Performance of the ensemble ordinal classifiers and simplified indicator (overall accuracy, balanced accuracy and Cohen’s κ). # model score_val eval_metric 0 CatBoost_r177_BAG_L2 0.727 accuracy 1 WeightedEnsemble_L3 0.723 accuracy 2 LightGBM_r131_BAG_L2 0.723 accuracy 3 CatBoost_r9_BAG_L2 0.722 accuracy 4 LightGBM_BAG_L2 0.722 accuracy 5 CatBoost_BAG_L2 0.771 accuracy 6 LightGBMXT_BAG_L2 0.712 accuracy 7 LightGBMLarge_BAG_L2 0.707 accuracy The confusion matrix (Fig. 4 ) for the selected three-class ensemble classifier illustrates how often low, medium and high categories are correctly or incorrectly predicted. Feature dimensionality reduction. While ensemble models achieve high predictive performance, their complexity poses challenges for operational use within national statistical offices. Therefore, feature correlations were analyzed to select variables with low mutual dependencies, as illustrated in Fig. 5 . This procedure reduced the original 54-variable feature space to the 20 most informative predictors. Gradient boosting ensemble models were trained via CatBoost and LightGBM, both of which are well suited to structured tabular data. Both multiclass and ordinal classification settings were tested. The same training, validation, and evaluation strategies were applied throughout, and model performance was assessed via overall accuracy, class-specific precision and recall, and Cohen’s κ to account for agreement by chance. In addition, confusion matrices were examined to identify misclassification patterns, particularly between the medium- and high-productivity categories. RESULTS Yield distributions of coffee and cocoa Coffee yields derived from the ENAl sample showed a markedly right-skewed distribution, with most smallholder farms concentrated below 500 kg ha⁻¹ and relatively few observations above 1,000 kg ha⁻¹. Cocoa yields displayed a similar shape but were more concentrated at approximately 600 kg ha⁻¹, with very few farms exceeding 1,500 kg ha⁻¹ (Figs. 1 – 2 ). These patterns confirm large productivity gaps within the same regions and crops and motivate the construction of an indicator capable of stratifying farms by productive capacity. Relationship between the NDVI and continuous yield NDVI dynamics were compared across the agricultural season for two contrasting groups of farms: those with the lowest yields and those with the highest yields in the sample. As shown in Fig. 6 , both groups displayed the typical seasonal pattern of coffee and cocoa systems, characterized by an early-season decline followed by recovery; however, differences were observed in both average NDVI levels and temporal stability. Farms with higher yields maintained consistently higher NDVI values, experienced less pronounced troughs, and showed faster recovery, whereas lower-yield farms exhibited deeper declines and slower rebounds. Exploratory analysis revealed no strong linear relationship between the NDVI metrics and continuous yields at the farm level. Scatter plots of the annual cumulative NDVI and maximum NDVI against yield (Fig. 7 ) revealed large dispersion and negligible linear trends, suggesting that multiple unobserved factors (management, age, shade, pests) modulate the NDVI–yield relationship. The ridge regression model achieved RMSE values of 339–340 kg ha⁻¹ and a MAE of approximately 240 kg ha⁻¹, with a slightly negative R², indicating that the model did not improve upon a naive mean-based prediction. The random forest model yielded similar RMSE and MAE values and likewise failed to capture substantial variance in yields. These results suggest that, in this context, NDVI time series alone are not sufficient to provide reliable continuous yield estimates at the scale of individual smallholder farms. Performance of ordinal classification models When yields were grouped into three categories of productive capacity (low, medium and high), model performance improved substantially. Gradient boosting ensemble classifiers trained on NDVI-derived features achieved validation accuracies of approximately 0.72–0.73 across several CatBoost and LightGBM configurations, with the best CatBoost model achieving approximately 0.77. On the held-out evaluation set, the selected three-class ensemble classifier attained overall accuracies of approximately 0.73 (validation) and 0.85 (evaluation), with Cohen’s κ ≈ 0.75, indicating substantial agreement beyond chance. Table 3 presents the main performance indicators of the best ensemble models and the simplified indicator, including overall accuracy, balanced accuracy and κ. Misclassifications occurred mainly between adjacent categories (low vs. medium, medium vs. high), whereas confusion between the lowest and highest classes was rare. Importantly, the models were particularly effective at identifying the lowest-productivity farms, which is a critical segment for targeting public support and extension. The confusion matrix (Fig. 8 ) for the selected three-class ensemble classifier illustrates how often low, medium and high categories are correctly or incorrectly predicted. Performance of the simplified decision-cut indicator The simplified decision-cut tree, built from a subset of NDVI features and constrained to a shallow depth, achieved an overall accuracy of 56.7% and κ = 0.31 on the evaluation set (Table 3 ). While clearly inferior to the ensemble models, this indicator still performed substantially better than random assignment and retained reasonable sensitivity for identifying low-yield farms. The final tree contains a small number of rules that combine thresholds on the NDVI in one or two critical windows with the maximum NDVI or NDVI amplitude. For example, farms with persistently low NDVI throughout the season and low maximum NDVI are classified as “low capacity”, whereas those that reach high NDVI peaks and maintain moderate variability are classified as “high capacity”. Intermediate combinations are assigned to the “medium” category. From an operational standpoint, this indicator can be implemented using preprocessed NDVI data stored in a central database, without the need for complex machine learning pipelines. It can be recalibrated periodically as new ENA campaigns and satellite data become available, adjusting cutoff points to reflect changes in management, varieties or climate. DISCUSSION This study explored the potential of the Sentinel-2 NDVI time series to support the construction of a parcel-level indicator of productive capacity for smallholder coffee and cocoa systems in Peru in combination with official agricultural survey data. The main findings can be summarized in three points. First, the ability of NDVI-based models to predict continuous yields at the individual holding scale was very limited. Second, when yields were grouped into broad categories of productive capacity, ensemble classifiers achieved substantially better performance. Third, a simplified rule-based indicator derived from a small set of NDVI features provided a transparent and operationally feasible approximation, albeit with lower accuracy than the full ensemble models. Our results confirm that the NDVI time series is a useful but incomplete proxy for productive capacity at the farm level in perennial smallholder systems. The lack of strong predictive power for continuous yields is consistent with recent meta-analyses showing that the NDVI–yield relationship weakens at fine spatial scales and in complex systems with strong management and environmental heterogeneity (Fuentes et al., 2025; Wuepper et al., 2025 ). In coffee and cocoa, factors such as tree age structure, pruning cycles, shade composition and pest and disease incidence can significantly affect yields without necessarily being fully captured by the NDVI. Nonetheless, the capacity of NDVI-derived features to support ordinal classification into meaningful categories of productivity suggests that they encode enough information to distinguish broad performance strata. This approach fits with the way in which satellite data have been used to monitor stress, canopy vigor and broad-scale yield patterns in coffee and cocoa landscapes (Anyimah et al., 2021 ; Ashiagbor et al., 2020; Atalaya-Marín et al., 2024; Martello et al., 2022 ). For national statistical systems, the value of the proposed indicator lies not in replacing survey-based yield estimation but in complementing it. By combining ENA microdata with NDVI time series, statistical agencies can construct an external reference that helps design stratified sampling frames that explicitly incorporate productivity gradients, identify extreme or inconsistent yield reports that merit further review, and characterize the distribution of productive capacity within and across regions at a finer spatial resolution than is possible with survey data alone. This aligns with international recommendations to integrate remote sensing into agricultural statistics as a way to improve efficiency, reduce respondent burden and increase spatial resolution (FAO, 2017 ; Wuepper et al., 2025 ). In addition, the indicator can inform the design of targeted policies and programmes in priority value chains such as coffee and cocoa, where productivity gaps, climate risks, and deforestation concerns intersect. Limitations and future research This analysis has certain limitations. First, it is based on one or a few agricultural seasons; extending the series to multiple campaigns would allow a better characterization of interannual variability, alternate bearing, and climate shocks (García et al., 2022 ; Manoel et al., 2024 ). Second, the models rely solely on the NDVI and do not incorporate other potentially informative layers, such as radar backscatter (Sentinel-1), climatic variables, soil properties or topography (Lu & Weng, 2007 ; Flores-Anderson et al., 2023 ). Third, yields are derived from farmer-reported data in the ENA, which may contain measurement errors; combining objective plot-level measurements with survey data in subsamples could help calibrate and validate the indicator more robustly. Finally, the study focuses on coffee and cocoa; assessing the transferability of the approach to other crops and agroecological contexts within Peru would be valuable for scaling. Building on these constraints, future work should incorporate multiyear time series and enrich the feature set with complementary sensors and environmental data. CONCLUSION This study shows that integrating official agricultural survey microdata with Sentinel-2 normalized difference vegetation index (NDVI) time series enables the construction of a farmer-level indicator of productive capacity for smallholder coffee and cocoa systems in Peru. While NDVI-based models exhibit very limited performance for continuous yield prediction at the farm level, they support robust ordinal classification into low-, medium-, and high-productivity categories when combined with appropriate modeling strategies. Gradient boosting ensembles achieve high classification performance but are difficult to implement operationally within statistical offices. A simplified decision-cut indicator derived from a small set of NDVI features offers a practical alternative that balances interpretability and predictive capacity. It can be embedded in the national agricultural statistical system to support sampling design, quality control and the targeting of extension and support services. The proposed indicator represents a first step toward a more integrated, spatially explicit, and data-rich approach to measuring productive capacity in Peruvian agriculture. Over time, expanding and refining this approach could help close productivity gaps, strengthen the resilience of smallholder systems and support more sustainable coffee and cocoa value chains. Declarations Competing Interest The authors declare that they have no competing interests. Funding The research leading to these results was funded by the Ministerio de Desarrollo Agrario y Riego del Perú (MIDAGRI) through the project Mejoramiento del Sistema de Información Estadística Agraria y del Servicio de Información Agraria para el Desarrollo Rural del Perú (PIADER), implemented by the Unidad Ejecutora de Gestión de Proyectos Sectoriales. Funding was awarded through the public call Concurso de Investigación sobre la Encuesta Nacional Agraria (Orden de Servicio No. 2672-2025). Data Availability ENA data are managed by the National Institute of Statistics and Informatics (INEI) of Peru. ENA public-use microdata can be downloaded from its official portal. Geospatial polygons of sampled farms are not publicly released and can only be accessed under the institution’s standard procedures. Sentinel-2 imagery is available through the Copernicus Open Access Hub. ACKNOWLEDGMENTS The authors thank the teams of UEGPS (Ministerio de Desarrollo Agrario y Riego del Peru) and INEI for facilitating access to ENA microdata and geospatial information, as well as the PIADER project team for their financial and technical support in the construction of the analytical database. Any remaining errors are the sole responsibility of the authors. Author Contribution R.L.-M. contributed to conceptualization, data curation, formal analysis, funding acquisition, investigation, methodology, resources, software, supervision, validation, visualization, and writing the original draft and reviewing and editing the manuscript. J.A.O.-B. contributed to funding acquisition, supervision, and review and editing. J.L.R.-V. contributed to funding acquisition, supervision, and review and editing. D.Q.-G. contributed to conceptualization, data curation, funding acquisition, investigation, project administration, resources, and writing the original draft and reviewing and editing the manuscript. All authors reviewed the manuscript. References Anyimah FO, Osei Jnr EM, Nyamekye C, 2021. Detection of stress areas on cocoa farms using GIS and remote sensing: A case study of Offinso Municipal and Offinso North District, Ghana. Environ Challenges 4: 100087. Ashiagbor G, Forkuo EK, Asante WA, Acheampong E, Quaye-Ballard JA, Boamah P, Mohammed Y, Foli E, 2020. Pixel-based and object-oriented approaches in segregating cocoa from forest in the Juabeso-Bia landscape of Ghana. Remote Sens Appl Soc Environ 19: 100349. Atalaya-Marín N, Goñas M, Tineo D, Chuquibala-Checan B, Arce-Inga M, Tarrillo C, Alvarez-Robledo YA, Tafur-Culqui J, Cabrera-Hoyos H, Gómez-Fernández D, et al., 2024. Integrating remote sensing and in situ data to determine climate diversity and variability in cocoa systems in the provinces of Jaén and San Ignacio, Cajamarca (NW Peru). Trees People 19: 100749. Bernardes T, Moreira MA, Adami M, Giarolla A, Rudorff BFT, 2012. Monitoring biennials bearing effect on coffee yield via MODIS vegetation indices. Remote Sens 4(9): 2492-2509. FAO, 2017. Handbook on remote sensing for agricultural statistics. Global Strategy to improve Agricultural and Rural Statistics (GSARS). FAO, Rome, Italy. Flores-Anderson AI, Groom G, Brumby SP, et al., 2023. Spatial and temporal availability of cloud-free optical observations over the tropics. Sci Data 10: 574. Fuentes I, Al-Shammari D, Al-Nasrawi AKM, Wang Y, Wang J, Lebrini Y, Chen Y, Jones BG, Bishop TFA, 2025. The normalized difference vegetation index as an analytic tool for wheat crop yield prediction: A review and meta-analysis. Precis Agric 26: 55. García GM, Perez-Esteban J, Tarquis AM, 2022. Reproductive trade-offs in a perennial crop: Exploring the causes of alternate bearing in coffee. Agric For Meteorol 323: 109049. Gil JDB, Reidsma P, Giller KE, Todman LC, Whitmore AP, van Ittersum MK, 2019. Sustainable development goal 2: Improved targets and indicators for agriculture and food security. Ambio 48(7): 685-698. Groten SME, 1993. NDVI-crop monitoring and early yield assessment of Burkina Faso. Int J Remote Sens 14(8): 1495-1515. Íñiguez Freites EA, Echevarría Zamora E, Ramírez Sánchez A, 2025. Geomorphological, biophysical and photogrammetric analysis of agroforestry systems associated with cacao (Theobroma cacao L.) cultivation in the Dominican Republic. Open Access Libr J 12(7): 1-16. INEI, 2025. Encuesta Nacional Agropecuaria 2023–2024: Metodología y principales resultados. Instituto Nacional de Estadística e Informática, Lima, Peru. Labus MP, Nielsen GA, Lawrence RL, Engel R, Long DS, 2002. Wheat yield estimates using multitemporal NDVI satellite imagery. Precis Agric 3(1): 15-26. Lu D, Weng Q, 2007. A survey of image classification methods and techniques for improving classification performance. Int J Remote Sens 28(5): 823-870. Manoel MC, Rosa MR, Queiroz AP, 2024. Analysis of the biennial productivity of Arabica coffee with Google Earth Engine in the northeast region of São Paulo, Brazil. Remote Sens 16(20): 3833. Martello M, Molin JP, Wei MCF, Canal Filho R, Nicoletti JV, 2022. Coffee-yield estimation using high-resolution time-series satellite images and machine learning. AgriEngineering 4(4): 888-902. Ministerio de Desarrollo Agrario y Riego del Perú (MIDAGRI), 2022. Plan Estadístico Agrario Nacional 2022–2026 del Sistema Integrado de Estadística Agraria. Ministerio de Desarrollo Agrario y Riego, Lima, Peru. Morales LV, Robiglio V, Baca M, Bunn C, Reyes M, 2022. Planning for adaptation: A system approach to understand the value chain’s role in supporting smallholder coffee farmers’ adaptive capacity in Peru. Front Clim 4: 788369. Nogueira SMC, Moreira MA, Volpato MML, 2018. Relationships between coffee crops productivity and vegetation indices derived from OLI/Landsat-8 sensor data with and without topographic correction. Eng Agric 38(3): 387-394. Wang J, Rich PM, Kettle WD, Price KP, 2005. Relationships between the NDVI and grassland production and crop yield in the Central Great Plains. Geocarto Int 20(3): 5-11. World Bank, 2017. Gaining momentum in Peruvian agriculture: Opportunities to increase productivity and enhance competitiveness. World Bank Group, Washington, DC. World Bank, 2023. World Development Indicators: Agriculture, forestry and fishing, value added (% of GDP) – Peru. https://data.worldbank.org/indicator/SL.AGR.EMPL.ZS?locations=PE [10 October 2025]. World Bank, 2024. World Development Indicators: Agriculture, forestry and fishing, value added (% of GDP) – Peru. https://data.worldbank.org/indicator/NV.AGR.TOTL.ZS [10 October 2025]. Wuepper D, Oluoch WA, Hadi H, 2025. Satellite data in agricultural and environmental economics: Theory and practice. Agric Econ 56(3): 493-511. Additional Declarations No competing interests reported. Cite Share Download PDF Status: Posted Version 1 posted You are reading this latest preprint version Research Square lets you share your work early, gain feedback from the community, and start making changes to your manuscript prior to peer review in a journal. As a division of Research Square Company, we’re committed to making research communication faster, fairer, and more useful. We do this by developing innovative software and high quality services for the global research community. Our growing team is made up of researchers and industry professionals working together to solve the most critical problems facing scientific publishing. Also discoverable on Platform About Our Team In Review Editorial Policies Advisory Board Help Center Resources Author Services Accessibility API Access RSS feed Manage Cookie Preferences © Research Square 2026 | ISSN 2693-5015 (online) Privacy Policy Terms of Service Do Not Sell My Personal Information {"props":{"pageProps":{"initialData":{"identity":"rs-9464917","acceptedTermsAndConditions":true,"allowDirectSubmit":true,"archivedVersions":[],"articleType":"Research Article","associatedPublications":[],"authors":[{"id":629218999,"identity":"25572039-2b02-47a7-b42d-881725f693c4","order_by":0,"name":"Robinson López-Monzon","email":"","orcid":"","institution":"Ministry of Agriculture of Peru","correspondingAuthor":false,"prefix":"","firstName":"Robinson","middleName":"","lastName":"López-Monzon","suffix":""},{"id":629219000,"identity":"a9c74906-d391-42e1-8886-22110d5d7d53","order_by":1,"name":"José Antonio Otoya-Barrenechea","email":"","orcid":"","institution":"Ministry of Agriculture of Peru","correspondingAuthor":false,"prefix":"","firstName":"José","middleName":"Antonio","lastName":"Otoya-Barrenechea","suffix":""},{"id":629219001,"identity":"a348b9de-86a5-4e02-a756-2b9b8a4cdca3","order_by":2,"name":"Juan Leonardi Romero-Vasquez","email":"","orcid":"","institution":"Ministry of Agriculture of Peru","correspondingAuthor":false,"prefix":"","firstName":"Juan","middleName":"Leonardi","lastName":"Romero-Vasquez","suffix":""},{"id":629219002,"identity":"6cbc2df5-fe83-4926-86aa-05287ab9e58d","order_by":3,"name":"Diana Carolina Quintero-Giraldo","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAABE0lEQVRIiWNgGAWjYDCCw1CaDULZyPGDqIQCorQwg6g0Y8kGkBYDPFoOwFlgLYcTN4BF8GjhO8787MEPBuvEPun+g48rKpgTN59fnfjhgQGDPL/YAaxaJA+zmRv2MKQntskcZjY8c4bNeNuNt5slgA4znDk7AasWg8MMZhI8QPe0SSSzSTa28chuu3F2A0hLgsFtXFrYv0n+gWhh/9nYJsG4ecbZzT/wa+Exk4bZwtjYZqC4gb93G15bJA/zlEnLGKQbA7UAw/dMgrHEDd5tFgkGEjj9wnf++DbJNxXWsvNnJD782FDxX46//+zmmz8qbOT5pbFrgTqPGYkjAVYpgUc5GCBr4T9ASPUoGAWjYBSMMAAAr6tdHQr821kAAAAASUVORK5CYII=","orcid":"","institution":"Ministry of Agriculture of Peru","correspondingAuthor":true,"prefix":"","firstName":"Diana","middleName":"Carolina","lastName":"Quintero-Giraldo","suffix":""}],"badges":[],"createdAt":"2026-04-19 21:53:15","currentVersionCode":1,"declarations":"","doi":"10.21203/rs.3.rs-9464917/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-9464917/v1","draftVersion":[],"editorialEvents":[],"editorialNote":"","failedWorkflow":false,"files":[{"id":107964773,"identity":"7376ac37-84fa-451f-a088-292d2db63801","added_by":"auto","created_at":"2026-04-28 05:27:54","extension":"jpg","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":52903,"visible":true,"origin":"","legend":"\u003cp\u003eHistogram of coffee yields (ENA sample)\u003c/p\u003e","description":"","filename":"1.jpg","url":"https://assets-eu.researchsquare.com/files/rs-9464917/v1/9762644759589db62a3a476d.jpg"},{"id":108007329,"identity":"b315ec6e-5c5b-4dcc-bc38-39d9cd1b9d91","added_by":"auto","created_at":"2026-04-28 12:59:29","extension":"jpg","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":55477,"visible":true,"origin":"","legend":"\u003cp\u003eHistogram of cocoa yields (ENA sample)\u003c/p\u003e","description":"","filename":"2.jpg","url":"https://assets-eu.researchsquare.com/files/rs-9464917/v1/5ff563ec6c389ae0100e9e9e.jpg"},{"id":107964766,"identity":"d06ced62-ff34-4c9b-a8fe-ee549f8dc6db","added_by":"auto","created_at":"2026-04-28 05:27:53","extension":"jpg","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":138206,"visible":true,"origin":"","legend":"\u003cp\u003eCorrelation matrix of the NDVI-derived features\u003c/p\u003e","description":"","filename":"3.jpg","url":"https://assets-eu.researchsquare.com/files/rs-9464917/v1/be6779b5d883df5bbbf699d5.jpg"},{"id":108007206,"identity":"3cb697b4-fc57-4375-9250-7f8084f5a433","added_by":"auto","created_at":"2026-04-28 12:58:54","extension":"jpg","order_by":4,"title":"Figure 4","display":"","copyAsset":false,"role":"figure","size":69356,"visible":true,"origin":"","legend":"\u003cp\u003eDecision tree cuts with dotted lines\u003c/p\u003e","description":"","filename":"4.jpg","url":"https://assets-eu.researchsquare.com/files/rs-9464917/v1/7c6329b02fffddc7e9fab59b.jpg"},{"id":107964768,"identity":"194aa320-1a6c-4323-b5f5-81c0e9792cab","added_by":"auto","created_at":"2026-04-28 05:27:53","extension":"jpg","order_by":5,"title":"Figure 5","display":"","copyAsset":false,"role":"figure","size":108010,"visible":true,"origin":"","legend":"\u003cp\u003eRanking of the most important NDVI-derived features for the ensemble classifier\u003c/p\u003e","description":"","filename":"5.jpg","url":"https://assets-eu.researchsquare.com/files/rs-9464917/v1/afe8ae5ec8328f4a8a8d404d.jpg"},{"id":107964769,"identity":"18016538-6e51-458b-a96b-50746c7aa061","added_by":"auto","created_at":"2026-04-28 05:27:53","extension":"jpg","order_by":6,"title":"Figure 6","display":"","copyAsset":false,"role":"figure","size":76887,"visible":true,"origin":"","legend":"\u003cp\u003eNDVI time series for the Top 5 vs. the Bottom 5 farms over 12 bimonthly periods\u003c/p\u003e","description":"","filename":"6.jpg","url":"https://assets-eu.researchsquare.com/files/rs-9464917/v1/fd8332bb5359bc60ae79af41.jpg"},{"id":107964770,"identity":"4d680d3d-0045-4a57-8b31-3726e8ecfa48","added_by":"auto","created_at":"2026-04-28 05:27:53","extension":"jpg","order_by":7,"title":"Figure 7","display":"","copyAsset":false,"role":"figure","size":39940,"visible":true,"origin":"","legend":"\u003cp\u003eScatter plot(s) of a key NDVI metric (e.g., cumulative NDVI or annual maximum) versus yield\u003c/p\u003e","description":"","filename":"7.jpg","url":"https://assets-eu.researchsquare.com/files/rs-9464917/v1/ffc65194703b8bf05f6b1e6f.jpg"},{"id":107964771,"identity":"ffe20f32-b2f8-4a3e-ac25-e489ac37b756","added_by":"auto","created_at":"2026-04-28 05:27:53","extension":"jpg","order_by":8,"title":"Figure 8","display":"","copyAsset":false,"role":"figure","size":39143,"visible":true,"origin":"","legend":"\u003cp\u003eConfusion matrix for the selected three-category ensemble classifier (Low, Mid, High)\u003c/p\u003e","description":"","filename":"8.jpg","url":"https://assets-eu.researchsquare.com/files/rs-9464917/v1/7be0d5eadd50af2915270344.jpg"},{"id":108804333,"identity":"23917954-ef19-4f2f-8421-71d2e7eb7c39","added_by":"auto","created_at":"2026-05-08 15:19:25","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":881890,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-9464917/v1/17af24a3-4a8a-4aff-b94f-63d382026256.pdf"}],"financialInterests":"No competing interests reported.","formattedTitle":"Developing a farm-level productivity indicator for coffee and cocoa via NDVI time series and National Agricultural Survey data in Peru","fulltext":[{"header":"KEYPOINTS/HIGHLIGHTS","content":"\u003cul class=\"decimal_type\"\u003e\n \u003cli\u003eThis study develops and validates a farm-level indicator of productive capacity for smallholder coffee and cocoa systems in Peru by integrating National Agricultural Survey (ENA) microdata with Sentinel-2 normalized difference vegetation index (NDVI) time series data.\u003c/li\u003e\n \u003cli\u003eThe results show that NDVI-based models are insufficient for reliable, continuous yield prediction at the parcel level, as regression performance remained weak (R\u0026sup2; \u0026asymp; 0; RMSE \u0026asymp; 340 kg ha⁻\u0026sup1;).\u003c/li\u003e\n \u003cli\u003eReframing the problem as an ordinal classification into low-, medium-, and high-productivity categories substantially improved model performance, reaching 0.73 validation accuracy, 0.85 evaluation accuracy, and Cohen\u0026rsquo;s \u0026kappa; = 0.75.\u003c/li\u003e\n \u003cli\u003eAlthough less accurate than ensemble models, the simplified decision tree indicator offers a more transparent and operationally feasible tool for implementation in agricultural statistical systems.\u003c/li\u003e\n \u003cli\u003eThe proposed indicator is a practical tool for supporting agricultural decision-making, particularly in improving sampling strategies, strengthening data quality control, and helping target extension and technical assistance more effectively.\u003c/li\u003e\n\u003c/ul\u003e\n\u003cp\u003e\u003cstrong\u003eIMPACT:\u003c/strong\u003e This study demonstrates that Sentinel-2 NDVI time series and official agricultural survey microdata can be combined to produce an operational indicator of productive capacity for smallholder coffee and cocoa systems in Peru. While developed with the microdata of the National Agricultural Survey (ENA) of Peru, the methodology also offers a potentially useful framework for other contexts seeking to make spatial variability actionable for sampling design, record review, and targeted support, subject to local calibration and validation.\u003c/p\u003e"},{"header":"INTRODUCTION","content":"\u003cp\u003eAgriculture remains a strategic sector for Peru, both for its contribution to the gross domestic product (GDP) and for its role in employment, territorial cohesion and rural livelihoods (World Bank, \u003cspan citationid=\"CR43\" class=\"CitationRef\"\u003e2023\u003c/span\u003e, \u003cspan citationid=\"CR45\" class=\"CitationRef\"\u003e2024\u003c/span\u003e). However, the lack of objective and comparable indicators of productivity at the farm level limits the capacity to target public policies, technical assistance and financial services to the areas and producers where they are most needed (World Bank, \u003cspan citationid=\"CR41\" class=\"CitationRef\"\u003e2017\u003c/span\u003e; Gil et al., \u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e2019\u003c/span\u003e). In practice, most decisions are still based on averages calculated at the regional or domain level, which masks the strong microterritorial heterogeneity observed within the same crop and year.\u003c/p\u003e \u003cp\u003ePeru has made significant progress in agricultural statistics through the National Agricultural Survey (ENA), which provides key information on area, production, yields and inputs for major crops (MIDAGRI, 2022). Nevertheless, estimating productive capacity at the farm level remains challenging, especially in perennial and smallholder-dominated systems such as coffee and cocoa, where plant density, age structure, management and environmental conditions vary widely even within small areas.\u003c/p\u003e \u003cp\u003eIn parallel, the availability of freely accessible satellite imagery and cloud-computing platforms has opened new possibilities for the use of remote sensing to monitor crop vigor and yield with high spatial and temporal resolutions (Flores-Anderson et al., \u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e2023\u003c/span\u003e; Wuepper et al., \u003cspan citationid=\"CR47\" class=\"CitationRef\"\u003e2025\u003c/span\u003e). Vegetation indices derived from multispectral images, particularly the normalized difference vegetation index (NDVI), have been consistently associated with crop development and yield in a wide range of systems (Groten, \u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e1993\u003c/span\u003e; Labus et al., \u003cspan citationid=\"CR25\" class=\"CitationRef\"\u003e2002\u003c/span\u003e; Wang et al., \u003cspan citationid=\"CR39\" class=\"CitationRef\"\u003e2005\u003c/span\u003e; Fuentes et al., 2025). In coffee and cocoa, recent studies have shown that the NDVI and related indices can help characterize yield variability, stress patterns and climatic gradients (Bernardes et al., \u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e2012\u003c/span\u003e; Nogueira et al., \u003cspan citationid=\"CR37\" class=\"CitationRef\"\u003e2018\u003c/span\u003e; Martello et al., \u003cspan citationid=\"CR31\" class=\"CitationRef\"\u003e2022\u003c/span\u003e; Atalaya-Mar\u0026iacute;n et al., 2024; \u0026Iacute;\u0026ntilde;iguez Freites et al., 2025).\u003c/p\u003e \u003cp\u003eHowever, the application of NDVI-based models in official operational statistics poses specific challenges. First, yield information in surveys such as ENA is collected through farmer-reported production and area, which may introduce measurement error and bias (FAO, \u003cspan citationid=\"CR9\" class=\"CitationRef\"\u003e2017\u003c/span\u003e). Second, perennial crops exhibit complex phenological patterns and alternate bearing and management cycles that are not easily summarized by single-date NDVI snapshots (Anyimah et al., \u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e2021\u003c/span\u003e; Manoel et al., \u003cspan citationid=\"CR29\" class=\"CitationRef\"\u003e2024\u003c/span\u003e). Third, statistical offices need simple and interpretable indicators that can be integrated into existing workflows rather than black-box models that are difficult to communicate and maintain (Wuepper et al., \u003cspan citationid=\"CR47\" class=\"CitationRef\"\u003e2025\u003c/span\u003e).\u003c/p\u003e \u003cp\u003eThe objective of this study is to develop and validate an annual farm-level indicator of productive capacity for smallholder coffee and cocoa systems in Peru through the integration of ENA microdata and Sentinel-2 NDVI time series. To achieve this goal, a harmonized database linking individual farms with satellite-derived NDVI metrics is first constructed; the capacity of these time series to predict continuous yield is then evaluated; then, ordinal classification models are developed to group farms into low-, medium-, and high-productivity categories; and finally, a simple and interpretable decision-rule indicator is derived for potential use within the Peruvian agricultural statistical system.\u003c/p\u003e"},{"header":"MATERIALS AND METHODS","content":"\u003cdiv id=\"Sec3\" class=\"Section2\"\u003e \u003ch2\u003eStudy area and data sources\u003c/h2\u003e \u003cp\u003eThis study focused on smallholder coffee and cocoa systems located in the eastern Andean slopes of Peru, particularly in the regions of San Mart\u0026iacute;n, Amazonas and Ucayali, which constitute a large share of national production.\u003c/p\u003e \u003cp\u003eThree main data sources were integrated. First, \u003cem\u003emicrodata from the National Agricultural Survey\u003c/em\u003e (ENA 2023\u0026ndash;2024) was used, as this survey provides farm-level information on cultivated area, production, management practices, and structural characteristics (INEI, \u003cspan citationid=\"CR23\" class=\"CitationRef\"\u003e2025\u003c/span\u003e). For coffee and cocoa, plot-level data on harvested area (ha) and production (kg) were extracted, and yields (kg ha⁻\u0026sup1;) were calculated as the ratio of production to area. Records with missing or inconsistent values for area or production were excluded in accordance with ENA quality-control criteria. Second, \u003cem\u003egeoreferenced polygons\u003c/em\u003e corresponding to each sampled farm were obtained from ENA and used to aggregate satellite-derived NDVI data at the farm level. Third, \u003cem\u003eNDVI time series\u003c/em\u003e were generated from Sentinel-2 Level-2A surface reflectance products at a 10 m spatial resolution, using bands 4 (red) and 8 (near-infrared). These data were processed for the 2023\u0026ndash;2024 agricultural season over a 12-month period and aligned with ENA yield reporting. Cloud- and cloud shadow-contaminated pixels were removed via quality flags and conservative thresholds, and the remaining observations were summarized into 12 bimonthly normalized difference vegetation index (NDVI) composites by calculating the median number of cloud-free pixels (Flores-Anderson et al., \u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e2023\u003c/span\u003e).\u003c/p\u003e \u003c/div\u003e\n\u003ch3\u003eData integration and analytical sample\u003c/h3\u003e\n\u003cp\u003eFor each sampled farm, the farm polygon was overlaid on the series of NDVI composites, and zonal statistics were calculated via the median NDVI value at each bimonthly interval. Farms with very small coffee or cocoa cultivation areas relative to the total polygon, as well as those with insufficient cloud-free observations, were excluded from the analysis. After applying all the consistency checks, the final analytical sample included several hundred coffee and cocoa farms. In the case of coffee, the yield distribution was strongly skewed, with most observations concentrated below 500 kg ha⁻\u0026sup1; and relatively few above 1,000 kg ha⁻\u0026sup1;, which is consistent with previous evidence on smallholder productivity in the Peruvian Amazon (World Bank, \u003cspan citationid=\"CR41\" class=\"CitationRef\"\u003e2017\u003c/span\u003e; Morales et al., \u003cspan citationid=\"CR35\" class=\"CitationRef\"\u003e2022\u003c/span\u003e). For cocoa, yields were more concentrated at approximately 600 kg ha⁻\u0026sup1;, with very few observations exceeding 1,500 kg ha⁻\u0026sup1;.\u003c/p\u003e \u003cp\u003eThese distributions are illustrated by histograms of coffee (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003e) and cocoa yields (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003e), which show a marked right-skewness and a small number of high-yield outliers.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eThe planting density distributions derived from the ENA confirmed substantial deviations from the theoretical or recommended densities; thus, the plant density was excluded as a direct predictor.\u003c/p\u003e\n\u003ch3\u003eFeature Engineering\u003c/h3\u003e\n\u003cp\u003eTwo families of NDVI-based features were constructed. The first consisted of \u003cem\u003esnapshot metrics\u003c/em\u003e, defined as the NDVI values observed during specific bimonthly periods corresponding to key phenological stages of coffee and cocoa, such as vegetative growth, flowering, and grain filling. For each crop, these periods were identified on the basis of agronomic calendars and local expert knowledge, and the median NDVI value of each farm polygon was retained for the selected dates.\u003c/p\u003e \u003cp\u003eThe second family consisted of \u003cem\u003ederived temporal features\u003c/em\u003e calculated from the full 12-step NDVI series. These included the annual maximum and minimum NDVI, the maximum minus minimum (NDVI) amplitude, the coefficient of variation, the interquartile range, the selected percentiles, and the cumulative NDVI, which were used as approximations of the area under the NDVI curve over the agricultural season.\u003c/p\u003e \u003cp\u003eThe full set of features was standardized and screened to avoid redundant predictors and reduce collinearity. A correlation matrix of NDVI-derived features was used to guide feature selection and aggregation, and a feature importance ranking from preliminary models helped identify the most informative predictors.\u003c/p\u003e \u003cp\u003eTable\u0026nbsp;\u003cspan refid=\"Tab1\" class=\"InternalRef\"\u003e1\u003c/span\u003e provides a detailed description of all the dependent, explanatory and control variables used in the analysis, including their roles, sources (ENA vs. Sentinel-2) and units.\u003c/p\u003e \u003cp\u003e \u003cdiv class=\"gridtable\"\u003e\u003ctable float=\"Yes\" id=\"Tab1\" border=\"1\"\u003e \u003ccaption language=\"En\"\u003e \u003cdiv class=\"CaptionNumber\"\u003eTable 1\u003c/div\u003e \u003cdiv class=\"CaptionContent\"\u003e \u003cp\u003eDefinitions, roles and sources of variables used in the models.\u003c/p\u003e \u003c/div\u003e \u003c/caption\u003e \u003ccolgroup cols=\"7\"\u003e \u003cdiv align=\"left\" class=\"colspec\" colname=\"c1\" colnum=\"1\"\u003e\u003c/div\u003e \u003cdiv align=\"left\" class=\"colspec\" colname=\"c2\" colnum=\"2\"\u003e\u003c/div\u003e \u003cdiv align=\"left\" class=\"colspec\" colname=\"c3\" colnum=\"3\"\u003e\u003c/div\u003e \u003cdiv align=\"left\" class=\"colspec\" colname=\"c4\" colnum=\"4\"\u003e\u003c/div\u003e \u003cdiv align=\"left\" class=\"colspec\" colname=\"c5\" colnum=\"5\"\u003e\u003c/div\u003e \u003cdiv align=\"left\" class=\"colspec\" colname=\"c6\" colnum=\"6\"\u003e\u003c/div\u003e \u003cdiv align=\"left\" class=\"colspec\" colname=\"c7\" colnum=\"7\"\u003e\u003c/div\u003e \u003cthead\u003e \u003ctr\u003e \u003cth align=\"left\" colname=\"c1\"\u003e \u003cp\u003eVariable\u003c/p\u003e \u003c/th\u003e \u003cth align=\"left\" colname=\"c2\"\u003e \u003cp\u003eRole\u003c/p\u003e \u003c/th\u003e \u003cth align=\"left\" colname=\"c3\"\u003e \u003cp\u003eOperational definition/Source\u003c/p\u003e \u003c/th\u003e \u003cth align=\"left\" colname=\"c4\"\u003e \u003cp\u003eENA Fields/Sentinel\u003c/p\u003e \u003c/th\u003e \u003cth align=\"left\" colname=\"c5\"\u003e \u003cp\u003eUnit\u003c/p\u003e \u003c/th\u003e \u003cth align=\"left\" colname=\"c6\"\u003e \u003cp\u003eScale\u003c/p\u003e \u003c/th\u003e \u003cth align=\"left\" colname=\"c7\"\u003e \u003cp\u003eProcedure\u003c/p\u003e \u003c/th\u003e \u003c/tr\u003e \u003c/thead\u003e \u003ctbody\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003eREND\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003eDependent\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003eTotal production/area\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003eP117_SUP_ha; totales/destinos Cap. 200AB (P220_*, en kg)\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c5\"\u003e \u003cp\u003ekg/ha\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c6\"\u003e \u003cp\u003eRatio\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c7\"\u003e \u003cp\u003eRegression and binning, classification\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003eRENDPLANT\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003eDependent\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003eREND/PLANTHA\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003eP126 (plantas), P117_SUP_ha\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c5\"\u003e \u003cp\u003ekg/plant\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c6\"\u003e \u003cp\u003eRatio\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c7\"\u003e \u003cp\u003eRegression; sensitivity to biases\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003ePLANTHA\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003eControl/Explanatory\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003eP126/P117_SUP_ha\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003eENA\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c5\"\u003e \u003cp\u003eplants/ha\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c6\"\u003e \u003cp\u003eRatio\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c7\"\u003e \u003cp\u003eControl; quality filters\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003eREGION\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003eControl\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003eGeographical Region\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003eREGION/NOMBREDD\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c5\"\u003e \u003cp\u003e\u0026mdash;\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c6\"\u003e \u003cp\u003eNominal\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c7\"\u003e \u003cp\u003eFixed effects/blocks\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003eSTARTMES\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003eControl\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003eMonth when harvest starts\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003eP124_MES\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c5\"\u003e \u003cp\u003emonth\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c6\"\u003e \u003cp\u003eOrdinal\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c7\"\u003e \u003cp\u003ePhenological alignment\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003endvi_cv_all\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003eExplanatory (global)\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003estd(NDVI)/mean(NDVI) 2023\u0026ndash;2024\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003eSentinel-2 (zonal stats)\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c5\"\u003e \u003cp\u003eadim.\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c6\"\u003e \u003cp\u003eRatio\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c7\"\u003e \u003cp\u003eKey feature; permutation importance\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003endvi_mean_all\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003eExplanatory(global)\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003eMean NDVI 2023\u0026ndash;2024\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003eSentinel-2\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c5\"\u003e \u003cp\u003eadim.\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c6\"\u003e \u003cp\u003eInterval\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c7\"\u003e \u003cp\u003eFeature; control for baseline status\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003endvi_amp_all\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003eExplanatory(global)\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003emax\u0026thinsp;\u0026minus;\u0026thinsp;min NDVI 2023\u0026ndash;2024\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003eSentinel-2\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c5\"\u003e \u003cp\u003eadim.\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c6\"\u003e \u003cp\u003eInterval\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c7\"\u003e \u003cp\u003eSeasonality signal\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003eAUC_main_23/24\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003eExplanatory(annual)\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003eArea under the curve (main window)\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003eSentinel-2\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c5\"\u003e \u003cp\u003eNDVI\u0026middot; bimonthly\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c6\"\u003e \u003cp\u003eInterval\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c7\"\u003e \u003cp\u003eWithin-year phenology\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003eauc_ratio_main_minor_23/24\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003eExplanatory(annual)\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003eAUC_main/AUC_minor\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003eSentinel-2\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c5\"\u003e \u003cp\u003eadim.\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c6\"\u003e \u003cp\u003eRatio\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c7\"\u003e \u003cp\u003eBalance of peaks\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003erise_pre_23/24, fall_harv_23/24, rebound_23/24\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003eExplanatory(annual)\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003eSlopes by window\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003eSentinel-2\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c5\"\u003e \u003cp\u003eΔNDVI/Δbimonthly\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c6\"\u003e \u003cp\u003eInterval\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c7\"\u003e \u003cp\u003eVigor and recovery\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003eskew_23/24, kurt_23/24\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003eExplanatory(annual)\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003eShape of NDVI distribution\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003eSentinel-2\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c5\"\u003e \u003cp\u003eadim.\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c6\"\u003e \u003cp\u003eInterval\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c7\"\u003e \u003cp\u003eCycle stability\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003edelta_\u0026lt;m\u0026eacute;trica\u0026gt;\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003eExplanatory(interannual)\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003e2024\u0026thinsp;\u0026minus;\u0026thinsp;2023 for each feature\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003eSentinel-2\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c5\"\u003e \u003cp\u003eaccording to metric\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c6\"\u003e \u003cp\u003e\u0026mdash;\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c7\"\u003e \u003cp\u003eInterannual change\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003eNDVI* (selectas)\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003eExplanatory (snapshot)\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003eNDVI in critical bimonthly periods\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003eSentinel-2\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c5\"\u003e \u003cp\u003eadim.\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c6\"\u003e \u003cp\u003eInterval\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c7\"\u003e \u003cp\u003eOnly if phenologically justified\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003c/tbody\u003e \u003c/colgroup\u003e \u003c/table\u003e\u003c/div\u003e \u003c/p\u003e \u003cp\u003eTo document collinearity patterns and guide feature selection, a correlation heatmap of all NDVI-derived variables was computed, which is presented in Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003e.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e\n\u003ch3\u003eContinuous yield modeling\u003c/h3\u003e\n\u003cp\u003eNDVI time series were evaluated to validate whether continuous yields (kg ha⁻\u0026sup1;) can be predicted at the holding level. Two supervised regression models were trained: (i) \u003cem\u003eridge regression\u003c/em\u003e, a linear model with L2 regularization that handles correlated predictors while preserving interpretability; and (ii) \u003cem\u003erandom forest regression\u003c/em\u003e, an ensemble of decision trees that are able to capture nonlinear relationships and interactions.\u003c/p\u003e \u003cp\u003eThe dataset was randomly split into training (60%), validation (20%) and evaluation (20%) subsets, which were stratified by yield quantiles. The hyperparameters (e.g., the regularization strength for ridges and the number and depth of trees for the random forest) were tuned via cross-validation on the training and validation sets. Model performance was evaluated on the held-out evaluation set via the root mean squared error (RMSE), mean absolute error (MAE) and coefficient of determination (R\u0026sup2;).\u003c/p\u003e \u003cp\u003eTable\u0026nbsp;\u003cspan refid=\"Tab2\" class=\"InternalRef\"\u003e2\u003c/span\u003e summarizes the performance of the continuous regression models, reporting the RMSE, MAE and R\u0026sup2; for both ridge and random forest.\u003c/p\u003e \u003cp\u003e \u003cdiv class=\"gridtable\"\u003e\u003ctable float=\"Yes\" id=\"Tab2\" border=\"1\"\u003e \u003ccaption language=\"En\"\u003e \u003cdiv class=\"CaptionNumber\"\u003eTable 2\u003c/div\u003e \u003cdiv class=\"CaptionContent\"\u003e \u003cp\u003ePerformance metrics (RMSE, MAE, R\u0026sup2;) of the ridge and random forest regression models for predicting continuous yields\u003c/p\u003e \u003c/div\u003e \u003c/caption\u003e \u003ccolgroup cols=\"4\"\u003e \u003cdiv align=\"left\" class=\"colspec\" colname=\"c1\" colnum=\"1\"\u003e\u003c/div\u003e \u003cdiv align=\"char\" char=\".\" class=\"colspec\" colname=\"c2\" colnum=\"2\"\u003e\u003c/div\u003e \u003cdiv align=\"char\" char=\".\" class=\"colspec\" colname=\"c3\" colnum=\"3\"\u003e\u003c/div\u003e \u003cdiv align=\"char\" char=\".\" class=\"colspec\" colname=\"c4\" colnum=\"4\"\u003e\u003c/div\u003e \u003cthead\u003e \u003ctr\u003e \u003cth align=\"left\" colname=\"c1\"\u003e \u003cp\u003eM\u0026eacute;todo\u003c/p\u003e \u003c/th\u003e \u003cth align=\"left\" colname=\"c2\"\u003e \u003cp\u003eRMSE\u003c/p\u003e \u003c/th\u003e \u003cth align=\"left\" colname=\"c3\"\u003e \u003cp\u003eMAE\u003c/p\u003e \u003c/th\u003e \u003cth align=\"left\" colname=\"c4\"\u003e \u003cp\u003eR2\u003c/p\u003e \u003c/th\u003e \u003c/tr\u003e \u003c/thead\u003e \u003ctbody\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e\u003cb\u003eRidge\u003c/b\u003e\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c2\"\u003e \u003cp\u003e339.4397\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c3\"\u003e \u003cp\u003e243.262\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c4\"\u003e \u003cp\u003e-0.0524\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e\u003cb\u003eRandomForest\u003c/b\u003e\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c2\"\u003e \u003cp\u003e346.477\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c3\"\u003e \u003cp\u003e235.234\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c4\"\u003e \u003cp\u003e-0.0965\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003c/tbody\u003e \u003c/colgroup\u003e \u003c/table\u003e\u003c/div\u003e \u003c/p\u003e\n\u003ch3\u003eOrdinal classification of productive capacity\u003c/h3\u003e\n\u003cp\u003eGiven the weak performance of the continuous models, the problem was reformulated as an \u003cem\u003eordinal classification task\u003c/em\u003e with three mutually exclusive categories of productive capacity: low, medium, and high, as illustrated in Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003e. Cutoff points were established via ENA yield data and agronomic criteria so that each category would represent meaningful productivity thresholds for policy and extension purposes, including very low yields well below regional averages, intermediate yields, and top-performing farms.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eGradient boosting ensemble models were then trained via CatBoost and LightGBM, both of which are well suited to structured tabular data. Multiclass and ordinal classification approaches were tested. The same training, validation, and evaluation strategies were applied, and model performance was assessed via overall accuracy, class-specific precision and recall, and Cohen\u0026rsquo;s κ to account for agreement occurring by chance. Confusion matrices were also examined to identify misclassification patterns, particularly between the medium- and high-productivity categories.\u003c/p\u003e \u003cp\u003eTable\u0026nbsp;\u003cspan refid=\"Tab3\" class=\"InternalRef\"\u003e3\u003c/span\u003e presents the main performance indicators of the best ensemble models and the simplified indicator, including overall accuracy, balanced accuracy and κ.\u003c/p\u003e \u003cp\u003e \u003cdiv class=\"gridtable\"\u003e\u003ctable float=\"Yes\" id=\"Tab3\" border=\"1\"\u003e \u003ccaption language=\"En\"\u003e \u003cdiv class=\"CaptionNumber\"\u003eTable 3\u003c/div\u003e \u003cdiv class=\"CaptionContent\"\u003e \u003cp\u003ePerformance of the ensemble ordinal classifiers and simplified indicator (overall accuracy, balanced accuracy and Cohen\u0026rsquo;s κ).\u003c/p\u003e \u003c/div\u003e \u003c/caption\u003e \u003ccolgroup cols=\"4\"\u003e \u003cdiv align=\"left\" class=\"colspec\" colname=\"c1\" colnum=\"1\"\u003e\u003c/div\u003e \u003cdiv align=\"left\" class=\"colspec\" colname=\"c2\" colnum=\"2\"\u003e\u003c/div\u003e \u003cdiv align=\"char\" char=\".\" class=\"colspec\" colname=\"c3\" colnum=\"3\"\u003e\u003c/div\u003e \u003cdiv align=\"left\" class=\"colspec\" colname=\"c4\" colnum=\"4\"\u003e\u003c/div\u003e \u003cthead\u003e \u003ctr\u003e \u003cth align=\"left\" colname=\"c1\"\u003e \u003cp\u003e#\u003c/p\u003e \u003c/th\u003e \u003cth align=\"left\" colname=\"c2\"\u003e \u003cp\u003emodel\u003c/p\u003e \u003c/th\u003e \u003cth align=\"left\" colname=\"c3\"\u003e \u003cp\u003escore_val\u003c/p\u003e \u003c/th\u003e \u003cth align=\"left\" colname=\"c4\"\u003e \u003cp\u003eeval_metric\u003c/p\u003e \u003c/th\u003e \u003c/tr\u003e \u003c/thead\u003e \u003ctbody\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e0\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003eCatBoost_r177_BAG_L2\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c3\"\u003e \u003cp\u003e0.727\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003eaccuracy\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e1\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003eWeightedEnsemble_L3\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c3\"\u003e \u003cp\u003e0.723\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003eaccuracy\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e2\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003eLightGBM_r131_BAG_L2\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c3\"\u003e \u003cp\u003e0.723\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003eaccuracy\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e3\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003eCatBoost_r9_BAG_L2\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c3\"\u003e \u003cp\u003e0.722\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003eaccuracy\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e4\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003eLightGBM_BAG_L2\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c3\"\u003e \u003cp\u003e0.722\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003eaccuracy\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e5\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003eCatBoost_BAG_L2\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c3\"\u003e \u003cp\u003e0.771\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003eaccuracy\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e6\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003eLightGBMXT_BAG_L2\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c3\"\u003e \u003cp\u003e0.712\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003eaccuracy\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003e7\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003eLightGBMLarge_BAG_L2\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"char\" char=\".\" colname=\"c3\"\u003e \u003cp\u003e0.707\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003eaccuracy\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003c/tbody\u003e \u003c/colgroup\u003e \u003c/table\u003e\u003c/div\u003e \u003c/p\u003e \u003cp\u003eThe confusion matrix (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003e) for the selected three-class ensemble classifier illustrates how often low, medium and high categories are correctly or incorrectly predicted.\u003c/p\u003e \u003cp\u003e \u003cb\u003eFeature dimensionality reduction.\u003c/b\u003e \u003c/p\u003e \u003cp\u003eWhile ensemble models achieve high predictive performance, their complexity poses challenges for operational use within national statistical offices. Therefore, feature correlations were analyzed to select variables with low mutual dependencies, as illustrated in Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003e. This procedure reduced the original 54-variable feature space to the 20 most informative predictors.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eGradient boosting ensemble models were trained via CatBoost and LightGBM, both of which are well suited to structured tabular data. Both multiclass and ordinal classification settings were tested. The same training, validation, and evaluation strategies were applied throughout, and model performance was assessed via overall accuracy, class-specific precision and recall, and Cohen\u0026rsquo;s κ to account for agreement by chance. In addition, confusion matrices were examined to identify misclassification patterns, particularly between the medium- and high-productivity categories.\u003c/p\u003e"},{"header":"RESULTS","content":"\u003cdiv id=\"Sec9\" class=\"Section2\"\u003e \u003ch2\u003eYield distributions of coffee and cocoa\u003c/h2\u003e \u003cp\u003eCoffee yields derived from the ENAl sample showed a markedly right-skewed distribution, with most smallholder farms concentrated below 500 kg ha⁻\u0026sup1; and relatively few observations above 1,000 kg ha⁻\u0026sup1;. Cocoa yields displayed a similar shape but were more concentrated at approximately 600 kg ha⁻\u0026sup1;, with very few farms exceeding 1,500 kg ha⁻\u0026sup1; (Figs.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003e\u0026ndash;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003e). These patterns confirm large productivity gaps within the same regions and crops and motivate the construction of an indicator capable of stratifying farms by productive capacity.\u003c/p\u003e \u003c/div\u003e\n\u003ch3\u003eRelationship between the NDVI and continuous yield\u003c/h3\u003e\n\u003cp\u003eNDVI dynamics were compared across the agricultural season for two contrasting groups of farms: those with the lowest yields and those with the highest yields in the sample. As shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003e, both groups displayed the typical seasonal pattern of coffee and cocoa systems, characterized by an early-season decline followed by recovery; however, differences were observed in both average NDVI levels and temporal stability. Farms with higher yields maintained consistently higher NDVI values, experienced less pronounced troughs, and showed faster recovery, whereas lower-yield farms exhibited deeper declines and slower rebounds.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eExploratory analysis revealed no strong linear relationship between the NDVI metrics and continuous yields at the farm level. Scatter plots of the annual cumulative NDVI and maximum NDVI against yield (Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003e) revealed large dispersion and negligible linear trends, suggesting that multiple unobserved factors (management, age, shade, pests) modulate the NDVI\u0026ndash;yield relationship.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eThe ridge regression model achieved RMSE values of 339\u0026ndash;340 kg ha⁻\u0026sup1; and a MAE of approximately 240 kg ha⁻\u0026sup1;, with a slightly negative R\u0026sup2;, indicating that the model did not improve upon a naive mean-based prediction. The random forest model yielded similar RMSE and MAE values and likewise failed to capture substantial variance in yields. These results suggest that, in this context, NDVI time series alone are not sufficient to provide reliable continuous yield estimates at the scale of individual smallholder farms.\u003c/p\u003e \u003cdiv id=\"Sec11\" class=\"Section2\"\u003e \u003ch2\u003ePerformance of ordinal classification models\u003c/h2\u003e \u003cp\u003eWhen yields were grouped into three categories of productive capacity (low, medium and high), model performance improved substantially. Gradient boosting ensemble classifiers trained on NDVI-derived features achieved validation accuracies of approximately 0.72\u0026ndash;0.73 across several CatBoost and LightGBM configurations, with the best CatBoost model achieving approximately 0.77. On the held-out evaluation set, the selected three-class ensemble classifier attained overall accuracies of approximately 0.73 (validation) and 0.85 (evaluation), with Cohen\u0026rsquo;s κ\u0026thinsp;\u0026asymp;\u0026thinsp;0.75, indicating substantial agreement beyond chance.\u003c/p\u003e \u003cp\u003eTable\u0026nbsp;\u003cspan refid=\"Tab3\" class=\"InternalRef\"\u003e3\u003c/span\u003e presents the main performance indicators of the best ensemble models and the simplified indicator, including overall accuracy, balanced accuracy and κ. Misclassifications occurred mainly between adjacent categories (low vs. medium, medium vs. high), whereas confusion between the lowest and highest classes was rare. Importantly, the models were particularly effective at identifying the lowest-productivity farms, which is a critical segment for targeting public support and extension.\u003c/p\u003e \u003cp\u003eThe confusion matrix (Fig.\u0026nbsp;\u003cspan refid=\"Fig8\" class=\"InternalRef\"\u003e8\u003c/span\u003e) for the selected three-class ensemble classifier illustrates how often low, medium and high categories are correctly or incorrectly predicted.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec12\" class=\"Section2\"\u003e \u003ch2\u003ePerformance of the simplified decision-cut indicator\u003c/h2\u003e \u003cp\u003eThe simplified decision-cut tree, built from a subset of NDVI features and constrained to a shallow depth, achieved an overall accuracy of 56.7% and κ\u0026thinsp;=\u0026thinsp;0.31 on the evaluation set (Table\u0026nbsp;\u003cspan refid=\"Tab3\" class=\"InternalRef\"\u003e3\u003c/span\u003e). While clearly inferior to the ensemble models, this indicator still performed substantially better than random assignment and retained reasonable sensitivity for identifying low-yield farms.\u003c/p\u003e \u003cp\u003eThe final tree contains a small number of rules that combine thresholds on the NDVI in one or two critical windows with the maximum NDVI or NDVI amplitude. For example, farms with persistently low NDVI throughout the season and low maximum NDVI are classified as \u0026ldquo;low capacity\u0026rdquo;, whereas those that reach high NDVI peaks and maintain moderate variability are classified as \u0026ldquo;high capacity\u0026rdquo;. Intermediate combinations are assigned to the \u0026ldquo;medium\u0026rdquo; category.\u003c/p\u003e \u003cp\u003eFrom an operational standpoint, this indicator can be implemented using preprocessed NDVI data stored in a central database, without the need for complex machine learning pipelines. It can be recalibrated periodically as new ENA campaigns and satellite data become available, adjusting cutoff points to reflect changes in management, varieties or climate.\u003c/p\u003e \u003c/div\u003e"},{"header":"DISCUSSION","content":"\u003cp\u003eThis study explored the potential of the Sentinel-2 NDVI time series to support the construction of a parcel-level indicator of productive capacity for smallholder coffee and cocoa systems in Peru in combination with official agricultural survey data. The main findings can be summarized in three points. First, the ability of NDVI-based models to predict continuous yields at the individual holding scale was very limited. Second, when yields were grouped into broad categories of productive capacity, ensemble classifiers achieved substantially better performance. Third, a simplified rule-based indicator derived from a small set of NDVI features provided a transparent and operationally feasible approximation, albeit with lower accuracy than the full ensemble models.\u003c/p\u003e \u003cp\u003eOur results confirm that the NDVI time series is a useful but incomplete proxy for productive capacity at the farm level in perennial smallholder systems. The lack of strong predictive power for continuous yields is consistent with recent meta-analyses showing that the NDVI\u0026ndash;yield relationship weakens at fine spatial scales and in complex systems with strong management and environmental heterogeneity (Fuentes et al., 2025; Wuepper et al., \u003cspan citationid=\"CR47\" class=\"CitationRef\"\u003e2025\u003c/span\u003e). In coffee and cocoa, factors such as tree age structure, pruning cycles, shade composition and pest and disease incidence can significantly affect yields without necessarily being fully captured by the NDVI.\u003c/p\u003e \u003cp\u003eNonetheless, the capacity of NDVI-derived features to support ordinal classification into meaningful categories of productivity suggests that they encode enough information to distinguish broad performance strata. This approach fits with the way in which satellite data have been used to monitor stress, canopy vigor and broad-scale yield patterns in coffee and cocoa landscapes (Anyimah et al., \u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e2021\u003c/span\u003e; Ashiagbor et al., 2020; Atalaya-Mar\u0026iacute;n et al., 2024; Martello et al., \u003cspan citationid=\"CR31\" class=\"CitationRef\"\u003e2022\u003c/span\u003e).\u003c/p\u003e \u003cp\u003eFor national statistical systems, the value of the proposed indicator lies not in replacing survey-based yield estimation but in complementing it. By combining ENA microdata with NDVI time series, statistical agencies can construct an external reference that helps design stratified sampling frames that explicitly incorporate productivity gradients, identify extreme or inconsistent yield reports that merit further review, and characterize the distribution of productive capacity within and across regions at a finer spatial resolution than is possible with survey data alone.\u003c/p\u003e \u003cp\u003eThis aligns with international recommendations to integrate remote sensing into agricultural statistics as a way to improve efficiency, reduce respondent burden and increase spatial resolution (FAO, \u003cspan citationid=\"CR9\" class=\"CitationRef\"\u003e2017\u003c/span\u003e; Wuepper et al., \u003cspan citationid=\"CR47\" class=\"CitationRef\"\u003e2025\u003c/span\u003e). In addition, the indicator can inform the design of targeted policies and programmes in priority value chains such as coffee and cocoa, where productivity gaps, climate risks, and deforestation concerns intersect.\u003c/p\u003e \u003cdiv id=\"Sec14\" class=\"Section2\"\u003e \u003ch2\u003eLimitations and future research\u003c/h2\u003e \u003cp\u003eThis analysis has certain limitations. First, it is based on one or a few agricultural seasons; extending the series to multiple campaigns would allow a better characterization of interannual variability, alternate bearing, and climate shocks (Garc\u0026iacute;a et al., \u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e2022\u003c/span\u003e; Manoel et al., \u003cspan citationid=\"CR29\" class=\"CitationRef\"\u003e2024\u003c/span\u003e). Second, the models rely solely on the NDVI and do not incorporate other potentially informative layers, such as radar backscatter (Sentinel-1), climatic variables, soil properties or topography (Lu \u0026amp; Weng, \u003cspan citationid=\"CR27\" class=\"CitationRef\"\u003e2007\u003c/span\u003e; Flores-Anderson et al., \u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e2023\u003c/span\u003e). Third, yields are derived from farmer-reported data in the ENA, which may contain measurement errors; combining objective plot-level measurements with survey data in subsamples could help calibrate and validate the indicator more robustly. Finally, the study focuses on coffee and cocoa; assessing the transferability of the approach to other crops and agroecological contexts within Peru would be valuable for scaling. Building on these constraints, future work should incorporate multiyear time series and enrich the feature set with complementary sensors and environmental data.\u003c/p\u003e \u003c/div\u003e"},{"header":"CONCLUSION","content":"\u003cp\u003eThis study shows that integrating official agricultural survey microdata with Sentinel-2 normalized difference vegetation index (NDVI) time series enables the construction of a farmer-level indicator of productive capacity for smallholder coffee and cocoa systems in Peru. While NDVI-based models exhibit very limited performance for continuous yield prediction at the farm level, they support robust ordinal classification into low-, medium-, and high-productivity categories when combined with appropriate modeling strategies.\u003c/p\u003e \u003cp\u003eGradient boosting ensembles achieve high classification performance but are difficult to implement operationally within statistical offices. A simplified decision-cut indicator derived from a small set of NDVI features offers a practical alternative that balances interpretability and predictive capacity. It can be embedded in the national agricultural statistical system to support sampling design, quality control and the targeting of extension and support services.\u003c/p\u003e \u003cp\u003eThe proposed indicator represents a first step toward a more integrated, spatially explicit, and data-rich approach to measuring productive capacity in Peruvian agriculture. Over time, expanding and refining this approach could help close productivity gaps, strengthen the resilience of smallholder systems and support more sustainable coffee and cocoa value chains.\u003c/p\u003e"},{"header":"Declarations","content":"\u003cp\u003e\u003cstrong\u003eCompeting Interest\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe authors declare that they have no competing interests.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eFunding\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe research leading to these results was funded by the Ministerio de Desarrollo Agrario y Riego del Perú (MIDAGRI) through the project Mejoramiento del Sistema de Información Estadística Agraria y del Servicio de Información Agraria para el Desarrollo Rural del Perú (PIADER), implemented by the Unidad Ejecutora de Gestión de Proyectos Sectoriales. Funding was awarded through the public call Concurso de Investigación sobre la Encuesta Nacional Agraria (Orden de Servicio No. 2672-2025).\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eData Availability\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eENA data are managed by the National Institute of Statistics and Informatics (INEI) of Peru. ENA public-use microdata can be downloaded from its official portal. Geospatial polygons of sampled farms are not publicly released and can only be accessed under the institution’s standard procedures. Sentinel-2 imagery is available through the Copernicus Open Access Hub.\u0026nbsp;\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eACKNOWLEDGMENTS\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe authors thank the teams of UEGPS (Ministerio de Desarrollo Agrario y Riego del Peru) and INEI for facilitating access to ENA microdata and geospatial information, as well as the PIADER project team for their financial and technical support in the construction of the analytical database. Any remaining errors are the sole responsibility of the authors.\u003c/p\u003e\u003ch2\u003eAuthor Contribution\u003c/h2\u003e\u003cp\u003eR.L.-M. contributed to conceptualization, data curation, formal analysis, funding acquisition, investigation, methodology, resources, software, supervision, validation, visualization, and writing the original draft and reviewing and editing the manuscript. J.A.O.-B. contributed to funding acquisition, supervision, and review and editing. J.L.R.-V. contributed to funding acquisition, supervision, and review and editing. D.Q.-G. contributed to conceptualization, data curation, funding acquisition, investigation, project administration, resources, and writing the original draft and reviewing and editing the manuscript. All authors reviewed the manuscript.\u003c/p\u003e"},{"header":"References","content":"\u003col\u003e\n\u003cli\u003eAnyimah FO, Osei Jnr EM, Nyamekye C, 2021. Detection of stress areas on cocoa farms using GIS and remote sensing: A case study of Offinso Municipal and Offinso North District, Ghana. Environ Challenges 4: 100087.\u003c/li\u003e\n\u003cli\u003eAshiagbor G, Forkuo EK, Asante WA, Acheampong E, Quaye-Ballard JA, Boamah P, Mohammed Y, Foli E, 2020. Pixel-based and object-oriented approaches in segregating cocoa from forest in the Juabeso-Bia landscape of Ghana. Remote Sens Appl Soc Environ 19: 100349.\u003c/li\u003e\n\u003cli\u003eAtalaya-Mar\u0026iacute;n N, Go\u0026ntilde;as M, Tineo D, Chuquibala-Checan B, Arce-Inga M, Tarrillo C, Alvarez-Robledo YA, Tafur-Culqui J, Cabrera-Hoyos H, G\u0026oacute;mez-Fern\u0026aacute;ndez D, et al., 2024. Integrating remote sensing and in situ data to determine climate diversity and variability in cocoa systems in the provinces of Ja\u0026eacute;n and San Ignacio, Cajamarca (NW Peru). Trees People 19: 100749.\u003c/li\u003e\n\u003cli\u003eBernardes T, Moreira MA, Adami M, Giarolla A, Rudorff BFT, 2012. Monitoring biennials bearing effect on coffee yield via MODIS vegetation indices. Remote Sens 4(9): 2492-2509.\u003c/li\u003e\n\u003cli\u003eFAO, 2017. Handbook on remote sensing for agricultural statistics. Global Strategy to improve Agricultural and Rural Statistics (GSARS). FAO, Rome, Italy.\u003c/li\u003e\n\u003cli\u003eFlores-Anderson AI, Groom G, Brumby SP, et al., 2023. Spatial and temporal availability of cloud-free optical observations over the tropics. Sci Data 10: 574.\u003c/li\u003e\n\u003cli\u003eFuentes I, Al-Shammari D, Al-Nasrawi AKM, Wang Y, Wang J, Lebrini Y, Chen Y, Jones BG, Bishop TFA, 2025. The normalized difference vegetation index as an analytic tool for wheat crop yield prediction: A review and meta-analysis. Precis Agric 26: 55.\u003c/li\u003e\n\u003cli\u003eGarc\u0026iacute;a GM, Perez-Esteban J, Tarquis AM, 2022. Reproductive trade-offs in a perennial crop: Exploring the causes of alternate bearing in coffee. Agric For Meteorol 323: 109049.\u003c/li\u003e\n\u003cli\u003eGil JDB, Reidsma P, Giller KE, Todman LC, Whitmore AP, van Ittersum MK, 2019.\u003c/li\u003e\n\u003cli\u003eSustainable development goal 2: Improved targets and indicators for agriculture and food security. Ambio 48(7): 685-698.\u003c/li\u003e\n\u003cli\u003eGroten SME, 1993. NDVI-crop monitoring and early yield assessment of Burkina Faso. Int J Remote Sens 14(8): 1495-1515.\u003c/li\u003e\n\u003cli\u003e\u0026Iacute;\u0026ntilde;iguez Freites EA, Echevarr\u0026iacute;a Zamora E, Ram\u0026iacute;rez S\u0026aacute;nchez A, 2025. Geomorphological, biophysical and photogrammetric analysis of agroforestry systems associated with cacao (Theobroma cacao L.) cultivation in the Dominican Republic. Open Access Libr J 12(7): 1-16.\u003c/li\u003e\n\u003cli\u003eINEI, 2025. Encuesta Nacional Agropecuaria 2023\u0026ndash;2024: Metodolog\u0026iacute;a y principales resultados. Instituto Nacional de Estad\u0026iacute;stica e Inform\u0026aacute;tica, Lima, Peru.\u003c/li\u003e\n\u003cli\u003eLabus MP, Nielsen GA, Lawrence RL, Engel R, Long DS, 2002. Wheat yield estimates using multitemporal NDVI satellite imagery. Precis Agric 3(1): 15-26.\u003c/li\u003e\n\u003cli\u003eLu D, Weng Q, 2007. A survey of image classification methods and techniques for improving classification performance. Int J Remote Sens 28(5): 823-870.\u003c/li\u003e\n\u003cli\u003eManoel MC, Rosa MR, Queiroz AP, 2024. Analysis of the biennial productivity of Arabica coffee with Google Earth Engine in the northeast region of S\u0026atilde;o Paulo, Brazil. Remote Sens 16(20): 3833.\u003c/li\u003e\n\u003cli\u003eMartello M, Molin JP, Wei MCF, Canal Filho R, Nicoletti JV, 2022. Coffee-yield estimation using high-resolution time-series satellite images and machine learning. AgriEngineering 4(4): 888-902.\u003c/li\u003e\n\u003cli\u003eMinisterio de Desarrollo Agrario y Riego del Per\u0026uacute; (MIDAGRI), 2022. Plan Estad\u0026iacute;stico Agrario Nacional 2022\u0026ndash;2026 del Sistema Integrado de Estad\u0026iacute;stica Agraria. Ministerio de Desarrollo Agrario y Riego, Lima, Peru.\u003c/li\u003e\n\u003cli\u003eMorales LV, Robiglio V, Baca M, Bunn C, Reyes M, 2022. Planning for adaptation: A system approach to understand the value chain\u0026rsquo;s role in supporting smallholder coffee farmers\u0026rsquo; adaptive capacity in Peru. Front Clim 4: 788369.\u003c/li\u003e\n\u003cli\u003eNogueira SMC, Moreira MA, Volpato MML, 2018. Relationships between coffee crops productivity and vegetation indices derived from OLI/Landsat-8 sensor data with and without topographic correction. Eng Agric 38(3): 387-394.\u003c/li\u003e\n\u003cli\u003eWang J, Rich PM, Kettle WD, Price KP, 2005. Relationships between the NDVI and grassland production and crop yield in the Central Great Plains. Geocarto Int 20(3): 5-11.\u003c/li\u003e\n\u003cli\u003eWorld Bank, 2017. Gaining momentum in Peruvian agriculture: Opportunities to increase productivity and enhance competitiveness. World Bank Group, Washington, DC.\u003c/li\u003e\n\u003cli\u003eWorld Bank, 2023. World Development Indicators: Agriculture, forestry and fishing, value added (% of GDP) \u0026ndash; Peru. https://data.worldbank.org/indicator/SL.AGR.EMPL.ZS?locations=PE [10 October 2025].\u003c/li\u003e\n\u003cli\u003eWorld Bank, 2024. World Development Indicators: Agriculture, forestry and fishing, value added (% of GDP) \u0026ndash; Peru. https://data.worldbank.org/indicator/NV.AGR.TOTL.ZS [10 October 2025].\u003c/li\u003e\n\u003cli\u003eWuepper D, Oluoch WA, Hadi H, 2025. Satellite data in agricultural and environmental economics: Theory and practice. Agric Econ 56(3): 493-511.\u003c/li\u003e\n\u003c/ol\u003e"}],"fulltextSource":"","fullText":"","funders":[],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":true,"hideJournal":true,"highlight":"","institution":"","isAcceptedByJournal":false,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":false,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true},"keywords":"remote sensing, yield mapping, multispectral time series, ensemble models, official statistics, tropical perennial crops","lastPublishedDoi":"10.21203/rs.3.rs-9464917/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-9464917/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"\u003cp\u003e\u003cstrong\u003eContext: \u003c/strong\u003eFarm-level productivity indicators are still limited in Peruvian agricultural statistics, particularly for smallholder perennial crops such as coffee and cocoa, where strong within-region heterogeneity and complex management conditions make productivity difficult to measure consistently. Moreover, increasing access to satellite imagery creates new opportunities to complement survey-based agricultural information with spatially explicit indicators.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eAims:\u003c/strong\u003e This study aimed to develop and validate a farm-level indicator of productive capacity for smallholder coffee and cocoa systems in Peru by integrating National Agricultural Survey (ENA) microdata with Sentinel-2 normalized difference vegetation index (NDVI) time series.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eMethods:\u003c/strong\u003e ENA 2023–2024 microdata were linked to georeferenced farm polygons and cloud-filtered 10 m Sentinel-2 normalized difference vegetation index (NDVI) imagery aggregated into 12 bimonthly composites. NDVI snapshots and temporal features were derived and used first in ridge and random forest regression models for continuous yield prediction and then in gradient boosting ensemble models for ordinal classification into low-, medium-, and high-productivity classes. A simplified decision-tree indicator was subsequently derived for operational use.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eKey Results:\u003c/strong\u003e Continuous-yield regression performed poorly (R² ≈ 0; RMSE ≈ 340 kg ha⁻¹), indicating that the NDVI alone was insufficient for reliable parcel-level yield estimation. In contrast, ordinal classification performed well, reaching 0.73 validation accuracy, 0.85 evaluation accuracy, and Cohen’s κ = 0.75. The simplified indicator achieved lower but still useful performance (56.7% accuracy; κ = 0.31).\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eConclusion:\u003c/strong\u003e Survey microdata and NDVI time series can be combined to generate an operational, farm-level indicator of productive capacity for smallholder perennial systems, even when continuous yield prediction remains weak.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eImplications and Impacts:\u003c/strong\u003e The proposed indicator provides a practical way to translate satellite time series and survey microdata into an operational tool for official statistics. In coffee and cocoa systems, it can support stratified sampling, improve the review of anomalous survey records, and guide more spatially explicit interventions, with potential relevance for other contexts after local calibration.\u003c/p\u003e","manuscriptTitle":"Developing a farm-level productivity indicator for coffee and cocoa via NDVI time series and National Agricultural Survey data in Peru","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2026-04-28 05:27:48","doi":"10.21203/rs.3.rs-9464917/v1","editorialEvents":[{"type":"communityComments","content":0}],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true}}],"origin":"","ownerIdentity":"c2fb6d32-dab5-4a1a-9cd7-b1621e5708bd","owner":[],"postedDate":"April 28th, 2026","published":true,"recentEditorialEvents":[{"type":"decision","content":"Rejected","date":"2026-05-05T12:39:49+00:00","index":"","fulltext":""},{"type":"editorAssigned","content":"","date":"2026-05-05T12:01:01+00:00","index":"","fulltext":""},{"type":"checksComplete","content":"","date":"2026-05-05T12:00:23+00:00","index":"","fulltext":""}],"rejectedJournal":[],"revision":"","amendment":"","status":"posted","subjectAreas":[],"tags":[],"updatedAt":"2026-05-05T12:56:10+00:00","versionOfRecord":[],"versionCreatedAt":"2026-04-28 05:27:48","video":"","vorDoi":"","vorDoiUrl":"","workflowStages":[]},"version":"v1","identity":"rs-9464917","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-9464917","identity":"rs-9464917","version":["v1"]},"buildId":"XKTyCvWXoU3ODBz1xrDgd","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}