Use of supervised and unsupervised approaches to make zonal application maps for variable-rate application of crop growth regulators in commercial cotton fields

Abstract

Background

Zonal application maps are designed to represent field variability using key variables that can be translated into tailored management practices. For cotton, zonal maps for crop growth regulator (CGR) applications under variable-rate (VR) strategies are commonly based exclusively on vegetation indices (VIs) variability. However, VIs often saturate in dense crop vegetation areas, limiting their effectiveness in distinguishing variability in crop growth. This study aimed to compare unsupervised framework (UF) and supervised framework (SUF) approaches for generating zonal application maps for CGR under VR conditions. During 2022–2023 agricultural seasons, an UF was employed to generate zonal maps based on locally collected field data on plant height of cotton, satellite imagery, soil texture, and phenology data. Subsequently, a SUF (based on historical data between 2020–2021 to 2022–2023 agricultural seasons) was developed to predict plant height using remote sensing and phenology data, aiming to replicate same zonal maps but without relying on direct field measurements of plant height. Both approaches were tested in three fields and on two different dates per field.

Results

The predictive model for plant height of SUF performed well, as indicated by the model metrics. However, when comparing zonal application maps for specific field-date combinations, the predicted plant height exhibited lower variability compared with field measurements. This led to variable compatibility between SUF maps, which utilized the model predictions, and the UF maps, which were based on the real field data. Fields characterized by much pronounced soil texture variability yielded the highest compatibility between the zonal application maps produced by both SUF and UF approaches. This was predominantly due to the greater consistency in estimating plant development patterns within these heterogeneous field environments. While VR application approach can facilitate product savings during the application operation, other key factors must be considered. These include the availability of specialized machinery required for this type of applications, as well as the inherent operational costs associated with applying a single CGR product which differs from the typical uniform rate applications that often integrate multiple inputs.

Conclusion

Predictive modeling shows promise for assisting in the creation of zonal application maps for VR of CGR applications. However, the degree of agreement with the actual variability in crop growth found in the field should be evaluated on a field-by-field basis. The SUF approach, which is based on plant heigh prediction, demonstrated potential for supporting the development of zonal application maps for VR of CGR applications. However, the degree to which this approach aligns itself with the actual variability in crop growth observed in the field may vary, necessitating field-by-field evaluation.

Introduction

Brazil is one of the world's main cotton producers. In terms of volume, it ranks behind China and India and is very close to the United States production (Organisation for Economic Co-operation and Development and Food and Agriculture Organization of the United Nations (OECD-FAO) 2023). In the 2022–2023 agricultural season, Brazil had a cultivated area of 1.66 million hectares of cotton. During the same agricultural season, the seedcotton yield was 4 257 kg·ha⁻¹ (with an average lint yield of 41%). At the national level, the Midwest region stands out: during 2022–2023, this macro-region had a cultivated area of 1.291 million hectares and an average yield of 4 588 kg·ha⁻¹. The largest national producing state, Mato Grosso, had a cultivated area of 1.190 million hectares (92% of the Midwest region and 72% of the national territory) and an average yield of 4 580 kg·ha⁻¹ (Companhia Nacional de Abastecimento (CONAB) 2023). In Mato Grosso and most of the Midwest region, where much of the country's agricultural activity is located, the predominant biome is the Cerrado (Brazilian savanna).

Agricultural areas within the Cerrado are predominantly characterized by weathered soils, high acidity, and low natural fertility. These regions are typically vast and flat, allowing for extensive mechanization and the establishment of large production units. Cotton production in the Cerrado is considered a high-value but high-risk activity due to the significant input costs and risks associated with environmental conditions. Production is largely rainfed, with the rainy season generally spanning from October to next April, delivering average annual precipitation totals between 1 200 and 2 000 mm (Souza et al. 2013), with the highest values occurring closer to the southern edge of the Amazon. Cotton is a unique crop in terms of its cultivation. Although it is a perennial plant, it is grown as an annual crop. Achieving consistently high yields and enabling efficient mechanization require careful management of the plant's vegetative and reproductive growth, which must be well balanced (Silvertooth et al. 1996). Environmental factors, such as temperature, cloud cover, and rainfall, alongside management practices, including nitrogen application, directly influence the growth dynamics of cotton. Frequent monitoring is necessary, and the application of crop growth regulator (CGR) has become a common practice in managing the plant’s development. CGRs help reduce plant height, control the height-to-node ratio, limit excessive vegetative growth, and improve the balance of photosynthate partitioning. Field assessments, particularly of plant height and structures like the number of nodes, are widely used for crop monitoring to guide CGR application (Echer et al. 2022). In the Brazilian Cerrado, the first CGR application is typically recommended at approximately 35 days after emergence (DAE), with the final application occurring around 115 DAE, which halts vegetative growth completely (Chiavegato et al. 2012). The specifics of CGR application, such as the number of treatments, dosage, and choice of commercial products, are determined based on local conditions and farm management practices (Berger et al. 2019).

As with many inputs in the production cycle, it is recommended that CGRs be applied according to site-specific needs within the field, rather than relying on an average value, within the context of precision agriculture (PA) (Plant 2001). This approach is particularly important for large commercial production areas—fields exceeding 100 hectares—which invariably exhibit spatial variability in environmental conditions like soil texture, ultimately influencing plant development. Site-specific management of crop inputs, enabled by variable-rate (VR) application techniques, is crucial for achieving more sustainable production, optimizing product use, and maximizing crop profitability—the key objectives of PA (Plant 2001; Taylor et al. 2010). VR applications in PA typically follow two primary methods: map-based or sensor-based systems. Map-based systems tend to be slower, as data must first be collected, analyzed through algorithms, and then used to generate zonal application maps for users (Taylor et al. 2003). Sensor-based systems, on the other hand, offer faster, real-time measurement and processing. However, their widespread adoption is currently limited by high costs, especially for small producers, and even large producers face significant expenses when considering the machinery needed for large-scale applications. This cost factor makes map-based systems a more feasible option for incorporating PA techniques in large commercial fields.

Many sensor- and/or map-based systems rely on the plant's reflectance response, commonly assessed through vegetation indices (VIs) like the normalized difference vegetation index (NDVI). VIs are algorithms that provide insights into the spectral signatures of vegetation, reflecting its structure, physiology, and biochemistry (Janse et al. 2018). These methods are widely used in PA and crop modelling studies to link vegetation characteristics with radiation captured via remote sensing. However, it's essential to note that the relationship between VIs and plant traits is indirect, as VIs act as indicators for these traits. Despite the broad application of NDVI throughout plant development, it tends to saturate under medium to high biomass conditions. This saturation phenomenon becomes apparent in row crops during stages of denser vegetation, limiting its effectiveness in monitoring further plant growth at those stages (Gitelson 2004; Payero et al. 2004; Gutierrez et al. 2012). In cotton cultivation, plant height is a critical biophysical characteristic closely linked to the need for CGR applications (Taylor et al. 2010). Although relationships between remote sensing bands, VIs, and plant height are well-documented (Payero et al. 2004), the linear association between VIs and plant height only holds at the beginning of the crop cycle. As the plant reaches a certain phenological stage, the saturation of VIs limits their ability to provide meaningful data for continued growth monitoring, underscoring the limitations of relying solely on VIs for CGR management (Taylor et al. 2010).

In PA, management zones (MZ) are commonly used to guide the site-specific application of crop inputs. By definition, an MZ is an area where crop development is relatively homogeneous and is influenced by manageable limiting factors. These zones are typically delineated based on a range of factors, which may be quantitative and stable (e.g., topography, soil texture, organic matter, pH), quantitative and dynamic (e.g., productivity, soil moisture, salinity, soil nitrogen status), qualitative and stable (e.g., soil survey maps, immobile nutrients like phosphorus (P)), or intuitive and historical (e.g., farmers' knowledge, past practices, and crop rotation). Inputs guided by MZ often include nitrogen (N), P, potassium (K), lime, and seeds (Doerge 2000). MZ are generally used to represent the productive potential of areas. These zones may change significantly from year to year (Plant 2001), or remain relatively stable when based on a longer-term understanding of the field (Santi et al. 2016; Damian et al. 2018). Intra- and inter-annual variability within fields, which influence site-specific input management, can be attributed to climate fluctuations (Alesso et al. 2023), as well as dynamic field factors like pests and diseases (Méndez-Vazquez et al. 2019). A challenge with traditional PA methods is delineating MZ solely through temporally stable factors, which may not account for dynamic occurrences. Incorporating dynamic information layers, such as those related to the biometric properties of plants, could improve zone accuracy and adaptability throughout the season (Pierce et al. 1999). Moving beyond VIs and addressing their tendency to saturate in high-biomass conditions offer further potential for enhancing MZ-based approaches.

A critical aspect of developing zonal application maps is the choice and adaptation of algorithms. While ready-to-use algorithms from research institutions are available, local adaptation is crucial to capture the specific conditions of the region (Liu et al. 2021; Vaz et al. 2023). Most studies rely on unsupervised multivariate techniques, such as clustering and principal component analysis (PCA), to delineate management zones (Fridgen et al. 2004; Córdoba et al. 2016; Damian et al. 2018; Ouazaa et al. 2022; Vaz et al. 2023). PCA's primary objective is to reduce dimensionality while retaining the maximum variance from the variables, creating linear combinations of conceptually relevant variables. The application of PCA in agriculture has been widely recognized in both statistical fields (Jolliffe 2002) and PA (Córdoba et al. 2016; Ortuani et al. 2019; Ouazaa et al. 2022; Vaz et al. 2023) across various crops and management aspects. This study builds on the findings of Andrea et al. (2023), where a supervised approach was tested to predict plant height based on remote sensing data. That study evaluated the predictive capabilities of different algorithms. The present study extends this by incorporating predicted plant height data and comparing two approaches for constructing zonal application maps that will be used for VR growth regulator management. Key points of the present study include:

• Large commercial fields exhibit significant variability in soil texture and crop development.

• Many CGR prescription maps are based solely on VIs, but VI saturation concerning crop height occurs during the reproductive phase.

• An unsupervised framework (UF) was developed to generate zonal application maps during the agricultural planting season, using crop height measurements, soil texture, and VIs.

• A supervised framework (SUF) was used as an alternative to field measurements of crop height, followed by unsupervised analysis to create zonal application maps.

• The zonal application maps generated by both UF and SUF were compared in terms of the compatibility and their relationship with the input variables.

The objective of this study was to compare zonal application maps designed for VR applications of CGR under two approaches: (i) using field-measured crop height in an unsupervised framework and (ii) using a combination of supervised and unsupervised frameworks without relying on field measurements of crop height.

Materials and methods

Study area

The data used in this study were collected from three commercial farms located in Mato Grosso, a state in the Central-West region of Brazil. These farms, owned by the same company, span more than 100 000 hectares, primarily dedicated to cotton and soybean production. Management practices across the farms are standardized, particularly in terms of nutrition and plant health, with adjustments made to account for local variations in weather, soil texture, and pest and disease pressures.

Input data

Field data: plant height, soil texture, days after emergence

In this study, the predictive algorithm was developed using data collected from more than 200 fields (Table 1) across the three previously mentioned farms. This extensive dataset was derived from routine field monitoring activities carried out on the farms over three growing seasons (2020–2021, 2021–2022, and 2022–2023). The dataset is representative of numerous fields, each monitored with varying frequency during the harvest, as part of the farms' commercial operations. The number and distribution of samples differed among fields and even within the same fields over time, primarily due to logistical factors. Monitoring activities are typically influenced by labor availability, coordination with other farm tasks, and environmental conditions, such as maintaining safety intervals after product applications or adjusting for rainfall and other weather events. The data presented in Table 1, combined with plant height measurements, resulted in approximately 30 000 georeferenced observations over the three seasons.

Table 1 General information about farms and fields that contributed to the predictive model training dataset

For the development of zonal application maps for crop growth regulator, three fields were selected (Table 2) in collaboration with local agronomists. These fields were chosen to enable more frequent and evenly distributed sampling throughout the agricultural seasons. The fields varied in size, ranging from 100 to 250 hectares, and were located in areas with diverse terrain characteristics. This included flat regions with up to 3% slope and medium-clayey soil texture, as well as more marginal areas exhibiting greater variability in both soil texture and slope. Detailed characteristics of these fields are provided in Table 2.

Table 2 Data on cultivation in the fields that were closely monitored and for which application maps were developed

In the fields that were closely monitored (Table 2), there was an average of 30 measurement points per monitoring event, with a frequency of approximately every ten days. However, this frequency could fluctuate across farms due to factors such as heavy rainfall, pesticide applications, and operational planning dynamics. Plant height data were gathered manually using a calibrated ruler, with each recorded measurement representing the average height of at least three adjacent plants. The coordinates of each plant height measurement were logged. Soil texture data were sourced from the farm's existing database and were gathered according to the farm's regular soil sampling protocols. The soil dataset was represented by a grid with a spatial resolution of 10 m. The variable DAE was used to represent plant phenology, with data coming from the farm's records on seedling emergence dates for each field.

Remote sensing data and dataset

Satellite multispectral imagery from Planet’s Planetscope satellite constellation, which includes over 180 CubeSats, was utilized in this study. The imagery provided a 4-band multispectral product covering the visible and near-infrared spectra, with the following spectral ranges: blue (455–515 nm), green (500–590 nm), red (590–670 nm), and near-infrared (780–860 nm). The Planetscope products offer surface reflectance with 16-bit depth and a 3-m pixel size, and are corrected for geometric and atmospheric effects using the 6SV2.1 radiative transfer code (Planet 2022). While the data is typically available with daily temporal resolution, this frequency may be reduced between October and next March due to rain and heavy cloud cover, which can obstruct image capture. The images used in this study were specifically acquired for cotton fields and were provided to the end-client with precise field boundaries, minimizing interference from other vegetation types. Details on the remote sensing data, including bands and VIs, are listed in Table 3.

Table 3 Summary of the remote sensing bands and VIs used in the present study

The selection of specific bands and VIs is based on previous research (Andrea et al. 2023). The variability in these products allows for a comprehensive assessment of crop traits, phenological phases (Gutierrez et al. 2012), and environmental conditions (e.g., soil and atmospheric effects) (Payero et al. 2004; Xue et al. 2017). Most VIs utilize the near-infrared (NIR) and red (R) bands, except for the Green NDVI (GNDVI), due to their well-documented sensitivity to photosynthetically active biomass (Tucker 1979; Gutman et al. 2021). On the days when plant height measurements were taken in the field, corresponding satellite images for that field—acquired within a maximum of 5 days away from the field sampling—were used to compile the dataset for the respective field.

General framework

The main steps involved in both approaches of this study are illustrated in Fig. 1. Although the input data for both UF and SUF approaches are nearly identical, their application differs based on the availability of field data on plant height and the use of predictive modelling. Further details on the procedures for UF and SUF approaches are provided in the following subsections.

Fig. 1

Flowchart summarizing the methodological approaches used in this study. SUF: supervised framework; UF: unsupervised framework; PCA: Principal component analysis

Unsupervised framework: PCA-based approach

As presented in Fig. 1, the UF approach began with the assessment of fields for the 2022–2023 agricultural season. Data from these fields comprising soil texture, VIs, and plant height samples underwent manipulation, which included interpolating the plant height measurements and constructing the dataset using a spatial join algorithm. Plant height sample data, which had a relatively low spatial resolution (typically 30 points per field), were interpolated using a radial basis function (RBF) algorithm with a thin plate spline kernel. This algorithm is renowned for its effectiveness in producing meshes, particularly in environmental sciences with scattered data (Sastry et al. 2015; Keller et al. 2019). To ensure consistent spatial resolution, all data were adjusted to a 10-m grid, aligning with the resolution of the soil texture data. For integrating all datasets, a core set KDTree (or cKDTree) algorithm was used, which creates an index for the data and performs nearest neighbor queries. More details on this algorithm can be found in Narasimhulu et al. (2021). Interpolation and spatial join processes were executed using the SciPy library (https://zenodo.org/records/6940349), and all analyses were conducted in Python. For each dataset, an unsupervised approach based on PCA was applied, culminating in the creation of zonal application maps. This PCA approach was central to the delineation of the zonal application maps. The goal of PCA was to summarize and reduce the information in the data by describing it with fewer concepts than originally present and deriving factor scores that replace the original values. The PCA was chosen due to the presumed existence of an underlying structure in the variables (Jolliffe 2002). The methodology began with centering the data in a coordinate system and identifying the first principal component (PC), which represents the direction of maximum variance among the observations. Each PC consists of two essential elements: the loading, which is a K-dimensional direction vector defining the best-fit line, and the score, which represents the projection value for each observation. The process then continues with the determination of the second PC, which accounts for the maximum remaining variance in the data. The combination of these two PCs creates a plane representing the latent variable model, providing optimal data representation (Jolliffe et al. 2016).

More information about the procedure for achieving these results is outlined in Table 4. The variable representing the majority of the variance in the data, weighted by the selected PC scores, was classified (three classes) using percentiles. This process resulted in a map with three distinct classes for the zonal application maps.

Table 4 Summary of the main steps taken in the PCA-based approach

Supervised framework: predictive modelling approach

The SUF approach, as illustrated in Fig. 1, began with the development of a dataset from a large number of fields over three agricultural seasons, as shown in Table 1. The dataset was then manipulated to include field-measured plant height, associated with respective DAE, and matched with remote sensing information based on proximity to crop height measurements. This manipulation created the modelling dataset, which was used for predictive modelling. In this dataset, plant height was the dependent variable, while remote sensing information and DAE were the independent variables. The data were divided into training and validation sets, with 85% allocated for this purpose (70% for training and 15% for validation) and the remaining data reserved for testing. For model optimization, a hyperparameter tuning process was performed using Bayesian optimization. This technique was chosen for its efficiency and precision advantages over other methods, such as grid search and random search. Bayesian optimization, implemented through a Python library (Nogueira 2014), uses a probabilistic model (often Gaussian processes) to approximate the true performance of the objective function, such as mean squared error (MSE), which was utilized in this study. The optimization process employs an acquisition function (or exploration strategy), like the upper confidence bound or expected improvement, to balance exploration and exploitation, emphasizing more promising regions. Further details on Bayesian optimization algorithms are provided by Wu et al. (2019) and Snoek et al. (2012). Random forest (RF), a robust machine learning technique widely used in agriculture (Jeong et al. 2016; Geetha et al. 2020; Carneiro et al. 2023), was the algorithm used in SUF. The supervised algorithm was applied to datasets from the fields evaluated in this study (Table 2) to predict plant height at dates beyond the modelling dataset's coverage. The predictions were assessed using key error metrics: R-squared (Chicco et al. 2021), mean absolute error (MAE) (Botchkarev 2019), and maximum absolute error (MAEr) (Botchkarev 2019). Soil texture data was incorporated into this resulting dataset, and the same unsupervised process previously described in the UF approach (PCA-based) was used to generate the zonal application maps.

Zonal application map comparison analysis

This section aims to present the techniques used for the comparative analysis of the maps generated by each framework. The maps were evaluated using classification metrics (Eqs. 1, 2, and 3), with UF serving as the reference. The performance metrics indicate the degree of compatibility between the UF and SUF approaches in defining classes.

$$Precision=\frac{TP}{TP+FP}$$

(1)

$$Recall=\frac{TP}{TP+FN}$$

(2)

$$F1-Score=\frac{2\ast\left(Precision\ast Recall\right)}{Precision+Recall}$$

(3)

Where, TP represents the true positives; FP represents the false positives; FN represents the false negatives. Precision measures how accurate a model’s positive predictions are, while recall measures how well a model identifies all positive instances. A high precision indicates fewer false positives, while a high recall indicates fewer false negatives. The F1-score balances precision and recall, providing a comprehensive evaluation of a model’s performance, especially in cases of class imbalance (Grandini et al. 2020). The analysis also presentes the mismatch—the inconsistencies in class assignments (classes 0, 1, or 2) between the two frameworks (UF and SUF). This provides a spatial distribution view of how well the approaches align in creating the zonal application maps.

An economic analysis was conducted to compare the UF and SUF approaches using a uniform-rate (UR) scenario as a reference, specifically considering the consumption of CGR.The input consumption was calculated using Eqs. 4 and 5.

$$Uniform\;rate\;consumption=\left(\sum areas\;of\;classes\right)\ast UR$$

(4)

$$Variable\;rate\;consumption={\textstyle\sum_{class=0}^2}\;\left(Area\ast Rate\;of\;class\right)$$

(5)

where the area of each class was calculated based on the map of application and the respective spatial resolution.

Results

Unsupervised framework: analysis and zonal application maps

The PCA-based approach was carried out to delineate the zonal application maps. Results of the main steps of this process can be seen in Table 5.

Table 5 Main parameters of the PCA-based approach of the unsupervised framework (UF) for fields-dates to create the zonal application maps

The PCA-based approach, with its main parameters detailed in Table 5, was primarily utilized within the UF. In this context, zonal application maps were created during the agricultural planting season based on farm demands and available field measurements. The correlations between plant height and other variables (NDVI and soil texture), presented in Table 5, showed significant variation across different field-date combinations. This variation indicates that the relationships among these variables are not static but change throughout the plant development cycle. Such dynamic relationships highlight the complexity of modelling plant growth parameters and underscore the challenges of evaluating these relationships during the growing season. The Bartlett test, conducted during the development of all zonal application maps, showed significance (P-value < 0.05) for all fields and dates in the UF, indicating the suitability of PCA. To determine the number of components to represent the data, the percentage of variance criterion was adopted after the exploratory analysis. The first two PCs were chosen, as they explained between 71% and 85% of the variance for all fields and dates, which was considered satisfactory. Communality scores, representing the proportion of each variable’s variance explained by the underlying latent variables, ranged from 0.553 to 0.993 in this study. Following the PCA-based approach in UF, the resulting zonal application maps are presented in Figs. 2, 3, and 4.

Fig. 2

Plant height and UF zonal application maps from field A on the first (upper row) and second dates (bottom row). a, c Plant height; b, d UF zonal application maps

Fig. 3

Plant height and UF zonal application maps from field B on the first (upper row) and second dates (bottom row). a, c Plant height; b, d UF zonal application maps

Fig. 4

Plant height and UF zonal application maps from field C on the first (upper row) and second dates (bottom row). a, c Plant height; b, d UF zonal application maps

In Figs. 2, 3, and 4, it is possible to see the spatial distribution of plant height data collected in the field alongside the zonal application maps for CGR management. The size of the points on the zonal application maps represents the location and relative classification of the sampled plant sizes. Application zones 0, 1, and 2 correspond to relatively homogeneous regions classified by the average plant height: low, medium, and high, respectively. Overall, the zonal application maps from the UF method visually align with the field-collected plant height data, as indicated by the points and their sizes. Regions in the field with taller plants corresponded to class number 2 in the zonal application maps. Further information regarding the spatial distribution and the variability in soil texture and NDVI can be found in the Supplementary Fig. S1.

Predictive modelling for supervised framework

As noted earlier, detecting the saturation of the NDVI is crucial, depending on the local data and research questions. During the predictive modelling phase, the data was first assessed for NDVI saturation in relation to plant height and phenology (Fig. 5). The results indicated that, under local conditions, NDVI saturation related to plant height occurs around 70 DAE. Although the NDVI shows a linear relationship with plant height early in the growth cycle, this linearity ends during the flowering phase, where continued vegetative growth is no longer reflected by the NDVI values.

Fig. 5

Relationships between DAE and NDVI (left figure) and between the plant height and NDVI across agricultural seasons (right figure) of the fields assessed in the this study

The modelling dataset was divided into training (70%), validation (15%), and test (15%) sets. The selection of variables was guided by previous studies and aimed to maintain simplicity for users who wish to replicate or extend the research. During the modelling phase, exploratory analyses were conducted to understand the relationships being evaluated. The global correlations between variables (correlation in the modelling dataset), shown in Fig. 6, include observations from all fields and seasons. The strongest correlation with plant height was observed with DAE, indicating that plant height tends to increase significantly over time, although not linearly. Following this, strong positive correlations were also found with VIs and the NIR band, suggesting that these variables can be used as indicators of plant height. On the other hand, the red, green, and blue bands exhibited medium-intensity negative correlations with plant height, also suggesting importance in the predictive process.

Fig. 6

Global correlation between variables obtained during 2020–2023 agricultural planting seasons

Model tuning was performed using a Bayesian optimization approach, which identified the optimal hyperparameter configuration by iteratively updating the probabilistic model to estimate the objective function. For the RF regressor, the optimization process focused on the following hyperparameters: the number of estimators, minimum samples split, minimum samples leaf, and maximum depth, with MSE as the target metric. The values corresponding to the best model performance are presented in Table 6. Further details about the hyperparameters can be found in the scikit-learn (Python library) documentation (https://scikit-learn.org/stable/modules/generated/sklearn.ensemble.RandomForestRegressor.html). The estimator's performance was evaluated using fivefold cross-validation, and the Bayesian optimization process was set to run for 30 iterations.

Table 6 Results of the best parameter settings and targets from Bayesian optimization in hyperparameter tuning

Regarding the algorithm performance, a comparison of the predicted versus true values and the performance metrics of the regressor on the test set have been presented in Fig. 7 and Table 7.

Fig. 7

Relationships between true and predicted values for plant height in the test set determined using the hyperparameter-tuned random forest regressor

Table 7 Metrics of the test set using the hyperparameter-tuned random forest regressor

The relationship between the predicted and observed values of plant height in the test set was illustrated in Fig. 7. This test set was exclusively used after the training, hyperparameter optimization, and validation phases to ensure the model's performance closely reflects real-world conditions. Although there is a general alignment between the predicted and observed values, some variability is present, as indicated by the points deviating from the 1:1 line in the left figure and the discrepancies between the paired bars representing plant height in the right figure. The algorithm’s performance, as quantified in Table 7, demonstrates satisfactory results across most metrics. However, the maximum error metric, which reflects the inherent variability in the data, shows some limitations, reflecting the inherent variability in the data.

Supervised framework: analysis and zonal application maps

The same procedure (PCA-based) was carried out after using the algorithm (Table 6) and predicting plant height for each specific field and date. The main parameters of this analysis in SUF were observed (Table 8). The correlations of the estimated plant height with the NDVI and soil texture for each assessed field-date were slightly different from the patterns found in the UF. In SUF, higher and positive correlations of plant height with NDVI and intermediate to null correlations with soil texture were found for all field-dates. The Bartlett's test resulted in large differences compared with UF, since for all field-date combinations in SUF, the P-values were not significant.

Table 8 Main parameters of the PCA-based approach of the supervised framework (SUF) for fields-dates to create the zonal application maps

The accumulated variance explained by the first two components was relatively high, around 90%, slightly higher than those found in the UF model. The communality values in the SUF model were also slightly higher than those previously found in the UF model, potentially due to the stronger correlations observed in this approach. Following the PCA-based approach, the zonal application maps for CGR management generated from the SUF are presented in Figs. 8, 9, and 10.

Fig. 8

Plant height and SUF zonal application maps from field A on the first (upper row) and second dates (bottom row). a, c Plant height; b, d SUF zonal application maps

Fig. 9

Plant height and SUF zonal application maps from field B on the first (upper row) and second dates (bottom row). a, c Plant height; b, d SUF zonal application maps

Fig. 10

Plant height and SUF zonal application maps from field C on the first (upper row) and second dates (bottom row). a, c Plant height; b, d SUF zonal application maps

In Figs. 8, 9, and 10, it is possible to see the spatial distribution of plant height predicted by the model, along with the zonal application maps in the fields resulting from the SUF model. Similar to the presentation in UF (Figs. 2, 3, and 4), the plots on the right show the sample points of plant height data collected in the field on the corresponding date. This visual comparison aimed to assess how well the application zones agrees with the observed plant height data. In this approach, the variability in plant height (left plots) was relatively low, and the general agreement between the application zones and the plant height sample points was also generally lower than that observed in UF, varying by field. For instance, in Field C, there was a stronger alignment, with larger points primarily located in Zone 2 and smaller points in Zone 0, compared with other fields.

Figure 11 provides additional insights into the variability and distribution of plant height, whether interpolated from field samples or predicted, as used in the development of zonal application maps through the different approaches: UF and SUF. The variability in plant height values across fields and dates, as shown in Figs. 2, 3, and 4 and 8, 9, and 10, is further highlighted in Fig. 11 for both the UF and SUF approaches. The UF approach exhibited greater variability in plant height values (especially in field C) when compared with the SUF approach across all dates and fields. Additionally, the first date of each field in the SUF approach showed slightly lower variability than the second date. These results suggested that while the predictive model performed well, the predicted data may underestimate the true variability observed in field conditions.

Fig. 11

Plant height variability across Fields A, B, and C, with respective dates (upper and lower rows) displayed for UF (left-side figures) and SUF (right-side figures) for each field. The orange line represents the median plant height values

Zonal application maps comparison

The comparative analysis of the maps based on the spatial distribution of all points revealed class mismatches between the UF and SUF approaches, in date 1 (left column) and date 2 (right column), regardless of the designated class (0, 1, or 2), as presented in Fig. 12. In Field A, which exhibited the lowest variability in soil texture, NDVI, and plant height, the lowest compatibility between the UF and SUF maps was observed (indicated by a higher number of red points) compared with other fields. Conversely, Field C, with the highest relative variability in soil texture, NDVI, and plant height, showed greater compatibility between the zonal application maps from both approaches. Field B displayed intermediate results.

Fig. 12

Mismatches in the classification of application zones between the UF and SUF approaches across Fields A (a), B (b), and C (c). For each field, the first and second dates correspond to the initial and subsequent map dates. Red points indicate mismatches between any of the classes 0, 1, and 2

The spatial distribution of plant height, whether based on field-collected or model-predicted data, is a critical aspect of this analysis. Field-sampled and interpolated data exhibited greater variability, aligning with observations made by field agronomists during the agricultural seasons. In contrast, the predictive model results were less variable, potentially underestimating real-world conditions. In order to quantify the consistency of zone delineation under both frameworks, a confusion matrix was developed using the classes from UF as a reference (Supplementary Fig. S2). They do not represent an absolute truth but rather the results of the approach considered as reference in this study (UF) in comparison with an alternative (SUF) for developing the zonal application maps for CGR management. The precision for the different classes showed wide variability: values between 0.20 and 0.30 (field A), between 0.05 and 0.85 (field B), and between 0.45 and 0.85 (field C). The greatest variability of precision was found for field B on the first date. Consistent higher values were found in field C and were generally lower in field A. In the present study, the precision and recall results were relatively similar across all fields and dates, thus the F1-scores also did not differ from the precision and recall. In conclusion, an overall superior compatibility between the UF and SUF maps of field C was found (Fig. 12).

Cost analysis

Cost analysis comparing the VR approaches for growth regulator application under both UF and SUF approaches compared with a UR strategy was presented (which uses a single dose of growth regulator) across the entire field (Table 9). This analysis focuses solely on the cost of the growth regulator product.

Table 9 Summary of the cost analysis of applying the growth regulator in the field and dates evaluated under the uniform and variable rates techniques

The results indicate that both VR strategies (UF and SUF) resulted in lower input costs compared with the uniform rate approach, with reductions between 32% and 43%. The cost reduction was observed in all VR scenarios. The values were based on rates applied in the field and were kept constant within each class (0, 1, or 2) to simplify the comparisons. This analysis demonstrates that the VR strategies from both UF and SUF provide significant cost savings over the UR approach.

Discussion

The relationship between plant height and VIs in our modelling dataset revealed complex dynamics. While global correlations were relatively high (r> 0.70) (Fig. 6), field-date-specific assessments demonstrated significant variability throughout crop development. This variability reflects the dynamic interactions between plant height and VIs, influenced by factors such as species, planting management, and phenological stages. The NIR band showed strong positive correlations with plant height, whereas blue, green, and red bands exhibited medium to strong negative correlations. This pattern likely stems from increased canopy cover and biomass, which reduce reflectance in visible bands while enhancing NIR reflectance. Our findings align with previous research. Hoffmeister et al. (2016) using terrestrial laser scanning found poor correlations between NDVI and plant height (r < 0.4), while Payero et al. (2004) and Gutierrez et al. (2012) highlighted that VI-biomass relationships vary significantly across phenological stages. Recent studies by Herr et al. (2023) emphasized the potential of unmanned aerial vehicles (UAV) for crop phenotyping, though scalability remains challenging in large commercial production areas.

Selecting the number of components in PCA involves inherent methodological challenges. Various approaches—including total cumulative variance, principal component variance magnitude, scree plots, and log-eigenvalue plots—provide guidance, yet retain a significant degree of subjectivity (Jolliffe 2002). These methods fundamentally depend on the researcher's analytical interpretation. Communalities—a critical metric in PCA—conventionally should exceed 0.4, with ideal values above 0.8 (Mittal 2020). However, Costello et al. (2019) noted that real-world data typically exhibited communalities between 0.4 and 0.7, consistent with our findings, particularly in the SUF approach where values were slightly elevated. Statistical testing revealed important methodological considerations. While Bartlett’s test confirmed data suitability in the UF approach, the SUF approach showed contrasting results, raising critical questions about method applicability (Single 1986; Jolliffe 2002).

Plant height estimation using multispectral satellite imagery remains relatively uncommon, with most remote sensing studies focusing on biomass and yield estimation (Kaplan et al. 2023; Herr et al. 2023). Weiss et al. (2020) noted that while plant height is integral to radiative transfer processes, its modelling complexity is lower compared with other agronomic variables, with more advanced estimation techniques emerging in LiDAR and photogrammetry. Our model's performance, validated through hyperparameter tuning, demonstrates promising potential. Recent studies illustrate diverse approaches to plant height prediction: Kaplan et al. (2023) achieved R² values of 0.70–0.95 using Sentinel data, Teodoro et al. (2021) obtained r > 0.77 with machine learning on UAV imagery, and Jamali et al. (2023) successfully predicted wheat plant height with an R² of 0.82 using Landsat-7 images and deep learning. Despite the growing popularity of UAV-based approaches, orbital satellite imagery remains crucial for large-scale agricultural monitoring. However, significant challenges persist in predicting plant height during later growth stages, where VI saturation and reduced model linearity complicate accurate estimations. Key insights include the potential of relatively simple modelling techniques—as demonstrated by Guo et al. (2022), who achieved an R² of 0.83 using basic linear regression—and the importance of understanding model limitations across crop development stages.

In Figs. 2, 3, and 4 (UF) and 8, 9, and 10 (SUF), it is possible to see the plant height distributions, providing critical insights for variable-rate growth regulator application. By mapping average plant heights across different zones, the analysis enables targeted management strategies that can enhance crop management efficiency. The PCA-based zoning approach demonstrated effectiveness in creating relatively homogeneous zones, classifying plant height into three levels (0, 1, 2). This technique proves particularly valuable when VIs are limited by saturation, as multiple variables can compensate for information variability. However, significant challenges persist. Satellite data availability is often constrained by weather conditions, especially during rainy seasons (Shafi et al. 2019; Prudente et al. 2020). Moreover, the SUF method's reliance on labor-intensive field-collected plant height data poses scalability limitations for large-scale commercial agricultural operations (Tanriverdi 2006; Ampatzidis et al. 2020). Despite these constraints, the approach demonstrates potential for precision agriculture, offering a nuanced method for understanding and managing spatial crop variability.

Comparative analysis of SUF and UF approaches revealed distinct spatial variability characteristics across three fields. SUF generally showed lower agreement with field samples, with notable misalignments in plant height zones, particularly in fields A and B. In contrast, field C demonstrated more consistent agreement between field data and defined classes across both approaches. SUF zones demonstrated greater stability, aligning more closely with soil texture, especially in areas of significant textural variability. This stability stems from low plant height variability, which allowed the PCA approach to capture more variance from soil texture. Field C, characterized by high soil textural variability and significant sand content, exhibited the most consistent zone delineations. Soil texture's critical role in crop growth was evident, influencing water retention, nutrient availability, and organic matter distribution (Basso et al. 2001; Jiang et al. 2020; Feng et al. 2022). Feng et al. (2022) notably highlighted seasonal variations in soil properties' impact, particularly under water stress, with sand content negatively correlating with soil moisture. SUF-derived plant height estimates consistently showed lower variability compared with UF across all fields, with field C presenting the most stable growth patterns despite typical management practices.

The PCA-based zonal mapping approach introduces a dynamic methodology that contrasts with traditional static zoning techniques relying on historical multi-year data (Plant 2001; Nawar et al. 2017; Damian et al. 2018). This adaptive approach aligns with precision agriculture principles, enabling seasonal re-delineation and improved spatial response understanding as computational capabilities evolve. Méndez-Vázquez et al. (2019) and Bökle et al. (2023) underscored the limited comparative research on zonal mapping methods, particularly for managing CGR applications. While VR technologies are widely used for inputs like limestone and nitrogen, CGR application remains challenging due to complex interactions between environmental factors, yield response functions, and marginal input benefits (Pedersen et al. 2017). Despite these complexities, zone delineation remains crucial for input management (Plant 2001; Taylor et al. 2003; Nawar et al. 2017; Ouazaa et al. 2022), with dynamic zoning offering a promising approach to navigate agricultural variability.

Comparing zonal maps reveals significant methodological challenges, primarily due to the absence of an absolute reference layer for plant height. The study found that zone compatibility between UF and SUF approaches varied inversely with textural homogeneity. Field A, with the least textural variability, showed the lowest compatibility, while field C, characterized by high sand content and significant textural variations, demonstrated the most consistent class compatibility. Class differentiation improves as inter-zone variability becomes more pronounced, underscoring the importance of soil texture in crop development. However, these effects were not static, being subject to inter-annual climate variability and crop management interactions (Guo et al. 2022). Bökle et al. (2023) emphasized the need to define robust zonal application map criteria, recognizing that optimal approaches may differ across field profiles. This highlights the complexity of precision agriculture and the necessity for tailored, adaptive strategies.

Cost analysis of VR for CGR applications reveals multiple critical considerations. Practical application expenses extend beyond input costs, incorporating machinery, labor, and operational factors. UR applications often combine multiple treatments, while VR operations are exclusive, potentially altering overall cost structures. The economic viability of VR technology remains context-dependent. While all assessed fields in this study showed similar input cost reductions, research results suggest VR provides the greatest benefits in heterogeneous fields with significant soil type variations (Pannell et al. 2019; Späti et al. 2021). Lawes et al. (2011) noted that VR could add value, but financial gains were typically modest. Challenges persist in VR implementation, including limited machinery availability, variable dosage requirements, and potential economic uncertainty. Studies like Boyer et al. (2011) have found no significant economic advantages, particularly in small-scale farming contexts. Consequently, VR profitability must be rigorously evaluated within specific local agricultural systems (Plant 2001; Fabiani et al. 2020; Tenreiro et al. 2023).

Conclusions

This study compared frameworks for generating zonal application maps for crop growth regulator management using a comprehensive dataset of cotton plant height from hundreds of commercial fields. Both approaches employed a PCA-based methodology, differing primarily in plant height data acquisition—field measurements versus model predictions. Field-measured data provided the most accurate representation of local variability, generating unique maps for each reading. The predictive model approach showed promise but exhibited variable effectiveness in capturing crop growth patterns. Soil texture emerged as a critical factor: fields with high textural variability demonstrated more consistent crop development and stable zonal maps, while fields with low textural variability showed less predictable growth patterns. Despite the potential of variable-rate strategies to enhance input efficiency, practical implementation remains challenging. Operational constraints, including machinery availability and cost considerations, particularly impact large agricultural properties, underscoring the complexity of precision agriculture adoption. Future studies could delve deeper into crop management practices, climate conditions, and their interactions to better understand the impact on crop height evolution throughout the agricultural season. Incorporating yield maps in future analyses could also offer valuable insights. Additionally, exploring alternative sampling methods and arrangements may enhance the spatial distribution layer of crop height data used as input.