Mapping Spatial Variability of Sugarcane Foliar Nitrogen, Phosphorus, Potassium and Chlorophyll Concentrations Using Remote Sensing

Abstract

Various nutrients are needed during the sugarcane growing season for plant development and productivity. However, traditional methods for assessing nutritional status are often costly and time consuming. This study aimed to determine the level of nitrogen (N), phosphorus (P), potassium (K) and chlorophyll of sugarcane plants using remote sensing. Remotely sensed images were obtained using a MicaSense RedEdge-P camera attached to a drone. Leaf chlorophyll content was measured in the field using an N-Tester chlorophyll meter, and leaf samples were collected and analyzed in the laboratory for N, P and K. The highest correlation between field samples and predictor variables (spectral bands, selected vegetation indices, and plant height from Light Detection and Ranging (LiDAR)), were noted.The spatial distribution of chlorophyll, N, P, and K maps achieved 60%, 75%, 96% and 50% accuracies, respectively. The spectral profiles helped to identify areas with visual differences. Spatial variability of nutrient maps confirmed that moisture presence leads to nitrogen and potassium deficiencies, excess phosphorus, and a reduction in vegetation density (93.82%) and height (2.09 m), compared to green, healthy vegetation (97.64% density and 3.11 m in height). This robust method of assessing foliar nutrients is repeatable for the same sugarcane variety at certain conditions and leads to sustainable agricultural practices in Costa Rica.

1. Introduction

Sugarcane (Saccharum officinarum) is a large, tropical grass belonging to the Graminaceae family [1]. It is one of the most economically important crops due to its potential in the agricultural and industrial sectors in Costa Rica [2]. Costa Rica is a major sugarcane producer, with Guanacaste province having the largest cultivation area at around 34,991.9 ha. This is due to the region’s favorable climate conditions, including an ideal average temperature of 27 °C, which promotes optimal crop growth [3]. Flat lands with loam and/or clay loam texture also facilitate the plant’s root expansion, enhancing its ability to extract nutrients [3].

The sugarcane phenological cycle consists of four stages: (i) germination from 0 to 90 days, (ii) vegetative growth from 90 to 180 days, (iii) vegetative development from 180 to 270 days, and (iv) maturation from 270 to 360 days; and the vegetative development stage shows the highest concentrations of nutrients within the plantation [4]. During all stages of growth, the concentration of nutrients must be maintained at an optimal level, with specific ranges for the key elements: (i) nitrogen (N): 2.0–2.6%, (ii) phosphorus (P): 0.20–0.30%, and (iv) potassium (K): 1.0–1.6% [5]. N, P and K are the elements with the highest demand and extraction by the plant, as they contribute to stimulating crop growth, stem development, sugar transport and photosynthesis [6]. Chlorophyll is a photoreceptor pigment which is crucial for light harvesting and converting solar energy into chemical energy [7,8]. It is also proportionally related to the nitrogen content in the plant, making it an indirect indicator of nitrogen levels [7,8]. Hence, a level of chlorophyll greater than 550 μmoles is considered optimal [9]. Additionally, K influences stem growth and development, while P is essential for chlorophyll synthesis and is related to the formation of sucrose and other sugars [4]. A deficiency in N, P, K and chlorophyll causes a decrease in stem diameter and a reduction in tillering, while an excess leads to dark green coloration of the leaves, abundant foliage, delayed maturation, and limited root system growth [10].

In Costa Rica, the primary method for analyzing and measuring chlorophyll and the nourishing elements (N, P, K) in crops is through foliar analysis, which involves collecting leaf samples from the study area and their subsequent analysis in the laboratory; this method provides detailed knowledge of the crop conditions throughout its cycle but requires extensive manual labor [11]. This process is costly, time consuming, and does not provide real-time information, highlighting the need for a more reliable and efficient method. However, Costa Rica needs a repeatable, time- and cost-effective alternative method that is simple and uses limited resources to assess crop nutrients. In this context, remote sensing has emerged as a potential tool for nutrient evaluation, which can significantly improve agricultural practices [12]. Remote sensors allow detecting and monitoring of specific characteristics of an area by measuring the radiation reflected or emitted by those areas [13]. Remote sensing methods for plant nutrient estimation vary from simple linear regression models using spectral bands and vegetation indices to complex machine learning algorithms [14,15,16]. The selection is mainly based on the correlation between nutrients and remote sensing predictors. For instance, a linear regression model will provide higher accuracies if the linear correlation is high (greater than 0.5) [17]. For example, Martins [14] and Lin [15] used leaf samples, spectral curves, and vegetation indices to develop predictive models for the nitrogen and chlorophyll content in sugarcane through simple linear regression. Both studies measured leaf spectral responses in the field using a portable hyperspectral sensor, achieving high accuracy with a determination coefficient (R²) of 0.89, root mean square error (RMSE) of 0.42, and relative prediction errors (RE) of 6.9%. Similarly, Narmilan et al. [18] used multispectral images captured by a drone to calculate vegetation indices, analyze the variable correlation using Pearson’s correlation coefficient, and apply machine learning algorithms to develop a chlorophyll predictive model. This model demonstrated robust performance, with R² values of 90% and an RMSE of 0.96, demonstrating the effectiveness of remote sensing methods.

Vegetation indices are numerical measures generated from combining multispectral bands and used to evaluate crop cover, vigor, and growth dynamics [12]. The commonly used vegetation indices are Response Index (RI), Green Chlorophyll Index (GCI), Visible Atmospherically Resistant Index (VARI), Excess Green Index (ExG), Normalized Blue-Green Difference Index (NGBDI), Triangular Greenness Index (TGI), Structure Insensitive Pigment Index (SIPI) and Leaf Chlorophyll Index (LCI), all of which have demonstrated strong correlations (Pearson’s coefficient > 0.5) with various nutrient levels, including N, P, K, and chlorophyll [16,18,19,20,21,22,23,24]. Light Detection and Ranging (LiDAR) technology also provides high spatial resolution three-dimensional data, allowing estimates of plant height, density, and topography [25]. These aspects are fundamental for establishing more accurate relationships between nutritional conditions and spectral information, enabling the identification of vegetation losses or height reduction due to nutrient deficiency or excess [25]. A Digital Terrain Model (DTM), Digital Surface Model (DSM), and Canopy Height Model (CHM) provide accurate information for plant height and, thus, relationships between nutritional conditions, vegetation indices, and the topography of the area [26].

When using remote sensing, isolating the vegetated areas before applying any algorithm is essential. Duro et al. [27] and Sumesh et al. [22] highlighted that the accuracy of object-based classification is higher, allowing for a more specific representation of the crop. Some studies showed lower accuracy when using Pixel-Based Image Analysis (PBIA) to separate sugarcane from elements unrelated to sugarcane (noise, soil, and grass); however, the Maximum Likelihood method is known to provide reliable results for classifying high-density sugarcane crops, making it a crucial tool for its analysis [28].

Multiple linear regression (MLR) is widely used in agricultural studies to model the relationship between dependent and independent variables [29]. This simple approach can be implemented at the field level with limited resources [29]. However, when considering the MLR for these types of applications, the data should satisfy the main assumption of MLR, which are the residuals or error terms must be normally distributed [30]. Hence, an exploratory data analysis is necessary.

Nutrient management at the farm level plays a crucial role in ensuring agricultural sustainability by optimizing the use of soil nutrients to enhance crop productivity while minimizing environmental impacts [1]. It involves strategic planning, efficient nutrient application, and adopting sustainable farming practices to maintain soil health, improve farm profitability, and contribute to broader environmental goals [6]. Effective nutrient management is fundamental for the sustainability of farming systems, as it boosts crop yields and protects ecological resources such as water, soil, and biodiversity [2]. By adopting practices such as precision agriculture, organic amendments, crop rotation, and conservation tillage, farmers can enhance the sustainability of their operations [2]. Hence, this study is a novel approach to sustainable sugarcane production in Costa Rica by implementing precision agriculture techniques at the farm level. The study developed a robust method that uses remote sensing data and an MLR algorithm to assess foliar nutrients, filling the existing research gap. The complete data analysis workflow, including developed codes for exploratory data analysis, can be repeated at the farm level.

This research aimed to develop a robust method that can be implemented at the farm level to assess sugarcane foliar nutrients using remote sensing. The attention was given to N, P, K and leaf chlorophyll concentrations, which are crucial during the maturation stage of sugarcane. The specific objectives were to (a) collect leaf samples and analyze them in the laboratory for nutrient content (N, P and K as calibration and validation data), (b) measure leaf chlorophyll in the field using a chlorophyll meter, (c) acquire multispectral images and LiDAR over the study area, (d) develop regression models between sugarcane nutrients and remote sensing data, and (e) map the spatial distribution of nutrients over the study area.

2. Materials and Methods

2.1. Study Area and Data

This study was conducted on a nine-hectare sugarcane plantation located in Cañas, Guanacaste, Costa Rica (center coordinates: 10°22′21.2″ N, 85°07′57.4″ W) (Figure 1). The farm belongs to an independent producer and is planted with the “TW08110” variety, with a furrow spacing of 1.7 m [31]. The vegetative cycle of sugarcane in the Chorotega region, where the study area is located, has an average duration of 12 months [2]. When sampling, the sugarcane was 165 days old (5.5 months) and was in the full vegetative growth stage [32].

According to field data, the soil in the study area is characterized by a clay loam texture, with a particle distribution of 45% sand, 27% silt, and 28% clay [33]. Similarly, the volumetric moisture content was 51%, at the time of sampling, corresponding to a value consistent with field capacity, given the present texture [34].

Twenty-five sampling plots were selected based on the plant’s homogeneity, spatial distribution, and access to the plots. Each sampling plot was georeferenced using a Global Navigation Satellite System (GNSS) with Real Time Kinematic (RTK) fix solution. In each plot, 30 sugarcane leaves were randomly collected and stored for N, P and K laboratory analysis. The nutritional analysis was conducted in the laboratories of DISAGRO, Abonos del Pacífico Costa Rica, S.A., following the methodologies established by AGQ Labs & Technological Services, S.L. [35]. Additionally, the chlorophyll content of the same 30 leaves was measured in situ using the N-Tester chlorophyll meter, which determines the chlorophyll present by measuring light absorbance at two wavelengths [36]. All samples obtained within each sampling plot were averaged at the end, resulting in twenty-five sample values.

A DJI Matrice 300 drone [37] equipped with the MicaSense RedEdge P multispectral camera [38] and the DJI Zenmuse L1 LiDAR sensor [39] were used to acquire multispectral images and LiDAR data. The flying height was 50 m, and the resultant spatial resolution of images was 3 cm. The positional accuracy of the system was increased to the centimeter level using the DJI D-RTK2 differential correction base station [40]. Multispectral images were radiometrically calibrated using the calibration panel provided for the MicaSense RedEdge P camera (AgEagle, Kansas, USA) [38], correcting reflectance values to reduce the impact of variable lighting conditions and ensure image consistency. The DN values of raw data ranged from 0 to 65,535 (16-bit unsigned). Once calibrated, the spectral reflectance values ranged from 0 to 0.25 for the blue band, from 0.01 to 0.51 for the green band, 0 to 0.33 for the red band, 0.04 to 0.98 for near the infrared band, and 0.02 to 0.69 for the RedEdge band. An orthomosaic was generated using the Pix4DMapper software (Version 4.9, Pix4D S.A., Prilly, Switzerland), and the positional accuracy (image alignment accuracy) of the project was 0.01 m.

2.2. Method

The LiDAR data were processed in DJI Terra software (Version 4.3.0, DJI Technology Co., Ltd., Guangdong, China) [41], removing all points classified as noise. Ground and above-canopy points were separated through iterative adjustments available in the DJI Terra software, optimizing the distinction between ground and non-ground points based on elevation and spatial distribution. The multispectral camera has five bands ranging from blue to the near-infrared region of the electromagnetic spectrum and one panchromatic band. The multispectral images were processed in Pix4D software (Version 4.9, Pix4D S.A., Prilly, Switzerland) to generate the orthomosaic with spectral reflectance values [42]. Figure 2 shows the complete flowchart of this study.

The Digital Terrain Model (DTM) and the Digital Surface Model (DSM) were generated from the LiDAR data. The DTM represents the topography of the study area, while the DSM includes both ground and non-ground points [43]. The Canopy Height Model (CHM) was generated in ArcGIS Pro software (Version 3.3.0, Esri, California, USA) [44], with a resolution of 0.05 m using Equation (1).

CHM = DSM − DTM,

(1)

A supervised classification was carried out to separate sugarcane and other features on the image using the Maximum Likelihood classification algorithm in ENVI software (Version 6.1, NV5 Geospatial, Long Beach, CA, USA) [45]. Over two hundred training samples were collected manually for soil, other vegetation, and sugarcane. The accuracy of the classification was tested based on validation samples.

According to

Table 1. Selected vegetation indices.

Index	Equation	Purpose	Reference
RI	$\mathrm{R}\mathrm{I}={\displaystyle \frac{\mathrm{R}-\mathrm{G}}{\mathrm{R}+\mathrm{G}}}$	To estimate the response of the plants to the nutrient concentrations present in them.	[20]
GCI	$\mathrm{G}\mathrm{C}\mathrm{I}={\displaystyle \frac{\mathrm{N}\mathrm{I}\mathrm{R}}{\mathrm{G}}}-1$	To estimate the leaf area index and green leaf biomass.	[12]
VARI	$\mathrm{V}\mathrm{A}\mathrm{R}\mathrm{I}={\displaystyle \frac{(\mathrm{G}-\mathrm{R})}{(\mathrm{G}+\mathrm{R}-\mathrm{B})}}$	To indicate the green level of a specific area in an image and it can also be used to detect stress in areas.	[21]
ExG	$\mathrm{E}\mathrm{x}\mathrm{G}={\displaystyle \frac{2\times \mathrm{G}-\mathrm{R}-\mathrm{B}}{\mathrm{R}+\mathrm{G}+\mathrm{B}}}$ 1	To estimate biomass and used to extract sugarcane pixels in the sample plot.	[22]
NGBDI	$\mathrm{N}\mathrm{G}\mathrm{B}\mathrm{N}\mathrm{I}={\displaystyle \frac{\mathrm{G}-\mathrm{B}}{\mathrm{G}+\mathrm{B}}}$	To predict canopy nitrogen concentration.	[16]
TGI	$\mathrm{T}\mathrm{G}\mathrm{I}=-0.5\left[\right({\mathsf{\lambda }}_{\mathrm{r}}-{\mathsf{\lambda }}_{\mathrm{b}}\left)\right(\mathrm{R}-\mathrm{G})-({\mathsf{\lambda }}_{\mathrm{r}}-{\mathsf{\lambda }}_{\mathrm{g}}\left)\right(\mathrm{R}-\mathrm{B}\left)\right]$ 2	To estimate nitrogen requirements and chlorophyll content in the leaves.	[23]
LCI	$\mathrm{L}\mathrm{C}\mathrm{I}={\displaystyle \frac{\mathrm{N}\mathrm{I}\mathrm{R}-\mathrm{R}\mathrm{E}}{\mathrm{N}\mathrm{I}\mathrm{R}+\mathrm{R}}}$	To relate the concentration of nutrients in the plant and determine the synthesis of chlorophyll.	[12]
SIPI	$\mathrm{S}\mathrm{I}\mathrm{P}\mathrm{I}={\displaystyle \frac{(\mathrm{N}\mathrm{I}\mathrm{R}-\mathrm{B})}{(\mathrm{N}\mathrm{I}\mathrm{R}-\mathrm{R})}}$	To estimate the concentration of pigments and nutrients in leaves related to chlorophyll.	[24]

¹ In ExG, factor 2 is to amplify the green band. ² λ_r, λ_b and λ_g represent the center wavelengths for the red, blue and green band, respectively. In addition, the factor 0.5 adjusts the scale of the formula.

, different indices were generated for further analysis within sugarcane coverage only. These indices were selected based on their capacity to isolate the plant vigor and stress and, thus, various nutrients. For example, TGI, which uses red, green and blue bands, estimates plants’ nitrogen and chlorophyll concentration. The attention will be given to the deviations of absorption near the red and blue wavelength regions and the reflectance in the green region from the healthy vegetation spectral profile. Table 1 summarizes the purposes of each index. Finally, each sampling plot’s average spectral values, vegetation height, and vegetation indices were calculated (only the sugarcane area within each plot was considered).

Exploratory data analysis assessed the normality of N, P, K, and chlorophyll samples using index plots, histograms, box plots, and normal Q–Q plots. Outliers were identified and removed. After that, four Pearson correlation matrices were calculated to analyze the relationship between each dependent variable (N, P, K, and chlorophyll) and independent variables (spectral bands, vegetation indices, and plant height). Independent variables with a correlation coefficient of 0.6 or higher were selected for further analysis. Based on this selection, a multiple linear regression (MLR) was tested using the “lm” package in R Studio software (Version 2024.09.0, Posit, Bostont, USA) [46], obtaining the intercept and the coefficients for each independent variable in the model [46]. Equation (2) shows an example of the MLR model.

$$Y={\mathsf{\beta }}_0+{\mathsf{\beta }}_1{\mathrm{X}}_1+{\mathsf{\beta }}_2{\mathrm{X}}_2+\dots +{\mathsf{\beta }}_{\mathrm{n}}{\mathrm{X}}_{\mathrm{n}}+\mathsf{ϵ}$$

(2)

where Y is the dependent variable to be predicted; X₁, X₂ and X_n are the independent variables; β₀ is the intercept; β₁, β₂, and β_n are the coefficients of each independent variable; and ϵ is the error [47].

The four generated models were evaluated using the coefficient of determination (R²), the adjusted R² and the p-value. The root mean square error (RMSE) and the mean absolute error (MAE) were calculated by comparing the predicted values from the models with the field collected data. The evaluated MLR models were used to generate the spatial distribution of each predictor, that is sugarcane nutrients. The accuracy of each predicted raster was assessed.

The visual investigation of the orthomosaic showed a distinct pattern in some areas of the sugarcane field. For instance, when displaying the image in true color combination, they were more reddish than other areas, and the density of the vegetation was lower. Several areas were identified with these differences, and 1000 random points were generated within those areas. Also, another 1000 points were generated outside these areas. The spectral information was extracted for both sets of points; spectral profiles were created and compared. After that, the elevation, plant height, and density of the representative areas were also extracted and tested against the optimal nutrient levels established for sugarcane in Costa Rica [5].

3. Results

3.1. Separating Sugarcane Coverage

The classification performed using the Maximum Likelihood method allowed the exclusion of areas not corresponding to sugarcane, such as soil, shadow and other vegetation. An accuracy of 99.76% and a kappa coefficient of 0.99 were obtained. Figure 3a shows one of the sampling plots with the unclassified orthomosaic, and Figure 3b shows the corresponding classified image. This classification integrates the various shades of sugarcane plants identified within the orthomosaic.

3.2. Exploratory Data Analysis

Four graphs were generated to identify the presence of outliers and evaluate the distribution of the dependent variables (N, P, K and chlorophyll). Figure 4 shows the data distribution for each variable, highlighting the number and sequence of collected samples on the Y-axis. According to the index plot of residuals (Figure 4a), no patterns are associated with residuals. The boxplots confirm the absence of outliers, and normal Q–Q plots indicate a normal distribution of each variable.

Additionally, nutrient levels were within optimal ranges, suggesting that at the time of sampling, the vegetation in the analyzed plots had an adequate nutritional content corresponding to its phenological stage.

3.3. Mapping Sugarcane Nutrients

Table 2. Correlation between crop nutrients and spectral bands, vegetation indices and plant height (RI = Response Index; ExG = Excess Green Index; NGBDI = Normalized Blue-Green Difference Index; TGI = Triangular Greenness Index; G = Green Band; RE = RedEdge Band; NIR = Near-Infrared Band; GCI = Green Chlorophyll Index; LCI = Leaf Chlorophyll Index; B = Blue Band).

Chlorophyll		Nitrogen		Phosphorus		Potassium
Independent Variable	Correlation Coefficient	Independent Variable	Correlation Coefficient	Independent Variable	Correlation Coefficient	Independent Variable	Correlation Coefficient
RI	0.616	G	−0.665	Plant Height	−0.898	G	−0.545
ExG	−0.671	RE	−0.600	GCI	−0.655	B	−0.557
NGBDI	−0.629	NIR	−0.611	NIR	−0.664	RE	−0.517
TGI	−0.651	TGI	−0.729	LCI	−0.779	TGI	−0.520

shows the most significant correlation between nutritional conditions (chlorophyll, nitrogen, phosphorus, and potassium) and independent variables (plant height, spectral bands, and vegetation indices). Independent variables with a correlation higher than 0.6 for chlorophyll, nitrogen, and phosphorus and higher than 0.5 for potassium were selected for further analysis.

For chlorophyll, the best correlation was with the ExG index (−0.671). Nitrogen showed the highest correlation with the TGI (−0.729). For phosphorus, the highest correlation was with plant height, with a value of −0.898. As for potassium, although no correlation coefficient exceeded 0.6, the blue band showed the best correlation coefficient, being −0.557.

Equations (3)–(6) show the prediction models for each nutrient. The accuracies of these models were evaluated based on the determination coefficient (R²), adjusted determination coefficient (R² adjusted), p-value, root mean square error (RMSE), and mean absolute error (MAE) (

Table 3. The regression model validation.

Value	Chlorophyll	Nitrogen	Phosphorus	Potassium
Coefficient of Determination (R2)	0.61	0.71	0.96	0.44
Adjusted coefficient of determination (R2 adjust)	0.50	0.60	0.89	0.28
p-value	0.0106	0.005	0.07	0.41
RMSE	13.95	0.05	0.002	0.07
MAE	11.64	0.04	0.002	0.06

). These metrics allowed for verifying the models’ reliability when using 70% of the data for prediction and 30% for validation.

Finally, these models were used to map the spatial distribution of different nutrients in the sugarcane plantation in the study area. The comparison between the measured and predicted average values for sample plots is shown in Figure 5 and

Table 4. Comparison between average measured and predicted nutrient values for sampling plots.

Variable	Average Values for Sample Plots
	Measured	Predicted
Chlorophyll (μmoles)	605.88	619.31
Nitrogen (%)	2.02	2.13
Phosphorus (%)	0.22	0.29
Potassium (%)	1.49	1.66

.

Chl = 5933.083 + 60557.184 × RI + 70389.000 × ExG − 60622.683 × NGBDI − 28.013 × TGI

(3)

N = 2.186 − 33.503 × G + 48.657 × RE + 0.681 × NIR − 1.432 × TGI

(4)

P = 0.3825 + 0.2241 × NIR − 0.0683 × GCI − 0.1546 × CHM + 0.7639 × LCI

(5)

K = 3.275 + 147.433 × G − 173.258 × B − 0.583 × RE − 1.763 × TGI

(6)

The accuracy of spatial variability of each nutrient map was determined using the validation data. The results showed 60%, 75%, 96% and 50% accuracy for chlorophyll, nitrogen, phosphorus and potassium, respectively. The highest difference was obtained for potassium; this may reflect the lowest accuracy of the potassium predictive model (Table 3). In general, these values demonstrate the representativeness of these models for the study area and highlight the applicability of these models for real-world field applications.

3.4. Analyzing Spectral Variations Within the Study Area

Figure 6 shows spectral profiles from two sets of points generated within selected areas. They are slightly different, which confirms the visual observations. Sugarcane areas where healthy sugarcane appears (green) replicate the spectral profile of healthy green areas. In contrast, the other areas with light green or brownish leaves show slightly lower values in the NIR region and higher values in the other areas.

Figure 7 compares the nutritional composition of two distinct sugarcane areas. According to

Table 5. Summary of the characteristics of the analyzed areas.

Average Nutrients and Other Characteristics	Healthy Green Areas	Distinct Areas with Less Green Leaves
Chlorophyll	603.4 μmoles	578.6 μmoles
N	1.82%	1.27%
P	0.32%	0.37%
K	1.21%	0.88%
Elevation	29.0–30.8 m	28.8–30.6 m
Plant height	3.11 m	2.09 m
Density	97.64%	93.82%

, there are differences in terms of average nutrient levels, elevation, plant heights and sugarcane density. The spatial distribution of these nutrients in Figure 5A–E highlights the differences. According to Figure 7, the left and right columns show different value ranges. For example, in Figure 7B, the right column has more leaves with lower nitrogen levels than the left. However, the area with less green/brownish leaves shows high phosphorus content (Figure 7C and Table 5). Compared to the elevation, phosphorus is high in low-elevation areas.

4. Discussion

The study aimed to map the spatial distribution of nutrients in a sugarcane field using remote sensing. A multiple linear regression algorithm was tested since the correlation between different nutrients, spectral bands and various vegetation indices was high. Later, the areas with distinct spectral information displayed were further investigated for nutrient variability.

The Maximum Likelihood classification algorithm performed well when isolating sugarcane coverage from other features on the image. There was a high agreement between the validation samples and the final classification (kappa = 0.99). As noted by Chikezie and Eme [48], when a kappa coefficient is greater than 0.74, it indicates an excellent agreement between the validation samples and results. Therefore, it can be concluded that the Maximum Likelihood classification is satisfactory for separating sugarcane coverage from other features.

An exploratory data analysis was conducted for field-collected data (N, P, K and chlorophyll) to check their data distributions. The visual inspection confirms their normality (Figure 4). For example, index plots show a random distribution of residual points without a pattern, and normal Q–Q plots resemble straight lines. Although the histograms show some skewness or minor irregularities, it can be resolved by changing the bin values. The distribution still follows the expected trend, with values concentrated in a specific range while maintaining symmetry. The existing patterns are not severe enough to suggest a significant departure from the normality. Further, the numerical analysis (p-value) confirmed those patterns. Hence, the perfectly normally distributed data ensures the testing of MLR models for further analysis [30]. The study did not consider implementing advanced machine learning algorithms because the new method should be simple enough to implement in the field environment.

When selecting predictor variables for MLR models, all variables that exhibited correlations above 0.5 were selected. This aligns with the robustness standards suggested by Hernández et al. [19]. The marginal correlation value (0.5) specified in this study is also consistent with those reported in similar research. Zhang et al. [49] reported correlation values of 0.63, 0.54, and 0.69 in predictive nitrogen models based on vegetation indices. Further, referring to several other studies [14,15,18,50], it can be confirmed that the marginal correlation value is suitable for selecting predictor variables for this study.

Significant correlations between chlorophyll, nitrogen, phosphorus, potassium, and vegetation indices such as RI, ExG, NGBDI, TGI, GCI and LCI highlighted the usefulness of these indices for estimating nutrient concentrations (Table 2). These findings align with previous studies that emphasize the importance of spectral bands, both directly and indirectly, in analyzing absorption and photosynthetic activity in the plantation [16,23,51]. For instance, RI uses red and green bands that are sensitive to chlorophyll absorption. Healthy green sugarcane typically shows reduced reflectance in these bands due to high chlorophyll absorption for photosynthesis [20]. The ExG index, which involves red, green, and blue bands, highlights green vegetation by capturing higher reflectance in the green band [22]. The NGBDI index is sensitive to chlorophyll content and leaf structure, showing increased green band reflectance [16]. GCI and LCI indices, using NIR and RE band, respectively, also reflect well in healthy vegetation and help detect soil moisture, as higher moisture increases absorption in these bands [18].

Apart from these indices, the study also considered the biophysical characteristics of sugarcane. For example, phosphorus is highly correlated with vegetation height. When there is a phosphorus deficiency, the vegetation is lower in height. This is mainly because the nutrient levels are related to plant growth and, thus, plant height [52]. However, further analysis throughout the crop’s phenological cycle is required to confirm this pattern.

As explained by Table 3 and Figure 5, the MLR models and the accuracies of the spatial distribution of different nutrients are acceptable. For instance, when developing the MLR model, the study considered the null hypothesis that the predictors do not correlate with the dependent variable. Since the p-values are statistically significant (less than 0.05) for chlorophyll, N and P, it was confirmed that they were correlated with the predictors. These results are consistent with Niedbala [53], whose study highlighted the usefulness of MLR in the agricultural field. The MLR model developed for potassium was not statistically significant, indicating no strong relation between potassium field samples and predictor variables used. Hence, the potassium model exhibited moderate results with different predictor variables. This led to the conclusion that a dense sampling strategy should be used to represent field variations accurately [54,55,56].

The robustness of the models is supported by determination coefficients (R²) greater than 0.6, as observed in this study. An R² of 0.6 indicated that 60% of the variability in the dependent variable is explained by the predictors, demonstrating strong predictive capacity. The RMSE and MAE values, which quantify the average magnitude of prediction errors, were below 0.3 for nitrogen, phosphorus, and potassium, meaning the prediction error was less than 0.3 units relative to the observed values for these nutrients [14]. For chlorophyll, the RMSE value of 14 represents a small error in absolute terms compared to field data, further supporting the model’s accuracy [15]. These metrics align with previous studies that validate RMSE and MAE as reliable model performance indicators [14,15]. Ultimately, the accuracy of the spatial distribution of nutrient maps was assessed based on field samples, confirming the models’ stability with acceptable accuracies (60% for chlorophyll, 75% for nitrogen, 96% for phosphorus, and 50% for potassium).

The visual analysis of the orthomosaic of sugarcane exhibited a pattern—patches of a mix of brownish or reddish, light green leaves in specific areas (upper left of Figure 5A). The average spectral profile of these areas was compared with the healthy sugarcane spectral profile (Figure 6). The average nutrients and biophysical variables were also noted (Table 5). Healthy green vegetation exhibited lower reflectance in the blue band due to its high absorption of light than distinct areas with light or brownish green leaves [57]. In the green band, healthy green vegetation reflected more light than areas with light or brownish green leaves, which reflected less light due to the plant’s structure and chlorophyll. In the red band, healthy green vegetation showed lower reflectance, as expected due to its greater light absorption for photosynthesis [20]. The RedEdge band better reflected the cellular structure in healthy green vegetation, allowing for the detection of water stress in areas with light or brownish leaves [20]. Finally, in the near-infrared band, healthy green vegetation reflected more due to cellular integrity, while areas with light or brownish green leaves affected by cellular damage reflected less. Then, it was found that these areas have higher moisture content. It can be concluded that areas with high moisture generated red “patchiness” in the orthomosaic (Figure 5A), and this is supported by the findings of Senyurek et al. [58], who emphasized that red and near-infrared bands are key to predicting soil and crop moisture. It is important to note that the humidity in the low-elevation areas of the terrain is due to excess water in the drainage of the plantation, originating from nearby flooded rice fields. This directly impacts sugarcane. Additionally, the natural saturation of the terrain and the excess water causes accumulation in the drainage areas, preventing adequate drainage and leading to physiological weakening of the crop.

Based on the above and analyzing the information from Table 5, the nutritional conditions in green, healthy areas (Figure 7), with a ground elevation between 29.04 m and 30.78 m, were optimal. According to Yara [5], the vegetation density (97.36%) and the average height (3.11 m) also confirmed the optimal growth of these green areas. These areas were supported by adequate moisture (good drainage based on the change of elevation), which enhances nutrient absorption by the plant, promotes vigorous growth, and increases crop density [59]. The spectral profile also shows optimal growth (Figure 6). In the distinct areas with light or brownish leaves (Figure 7), with heights between 28.75 m and 30.63 m, nitrogen and potassium deficiencies and an excess of phosphorus were detected. Yara [5] further confirmed this. Leaching and denitrification in saturated soils hinder nutrient absorption by the plant [60]. Also, mineralization in moist soils leads to the rapid release of organic phosphorus that is quickly absorbed, coupled with the reduction of phosphorus-fixing compounds that increase its availability to the plant [61]. These deficiencies resulted in lower vegetation coverage (93.82%), a decrease in plant height in these areas (2.09 m), stunted growth, and yellowing of the lower leaves, with symptoms of water stress [62]. The above information confirms that, although the correlation between nutritional conditions and vegetation height is not strong, a reduction of approximately one meter in plant height is observed in areas with nutritional deficiencies. Due to the sampling strategy of this study (limited to healthy, green and accessible areas only), there were no samples from areas where patchiness was exhibited; however, the final spatial variability maps differentiated those areas based on spectral values. A representative sampling strategy is recommended for the future. Terrain topography is analyzed using remote sensing techniques, as moisture influences sucrose production and maturation, while topography plays a fundamental role in controlling water accumulation in the soil [3,4].

For future research, conducting a denser, representative sampling in the study area is recommended to establish relationships between nutrition and vegetation height. Additionally, it is essential to analyze these variations throughout the crop’s phenological cycle to confirm the relation between soil moisture, growth, and crop nutrition. This will allow the development of regression models and relationships at different stages of the sugarcane phenological cycle. Advanced machine learning algorithms can also be tested on this application, requiring more resources (computers and software) and technological knowledge to implement at the farm level.

Finally, this innovative approach enables more precise crop development and production optimization monitoring at the farm level. The developed method is time- and cost-effective and technologically simple. Further, nutrient management at the farm level plays a crucial role in ensuring agricultural sustainability while minimizing environmental impacts. It involves strategic planning, efficient nutrient application, and adopting sustainable farming practices to maintain soil health, improve farm profitability, and contribute to broader environmental goals. Therefore, this study is a novel approach to sustainable sugarcane production in Costa Rica that implements precision agriculture techniques at the farm level.

5. Conclusions

Sugarcane is a crop of great global importance, particularly in Costa Rica, due to its impact on the economic, social, energy, and technological sectors. The edaphic and climatic conditions of the Guanacaste province favor its development, positioning Costa Rica as a leader in sugar production. Adequate nutrition of this crop is crucial for maximizing its growth and productivity. For example, nitrogen, phosphorus, potassium, and chlorophyll play a key role in stem development, root growth, and sugar transport.

Remotely sensed images were acquired using a MicaSense RedEdge-P camera, leaf samples were collected to assess the foliar nutrients (N, P and K) in the laboratory, and leaf chlorophyll concentration was measured in the field. Field samples were used for calibration and validation purposes. Remotely sensed images were radiometrically and atmospherically corrected, and an orthomosaic was generated. The sugarcane area was isolated using a Maximum Likelihood classification algorithm. The classification received 99% of the accuracy compared to the validation samples with a kappa coefficient of 0.99. Light Detection And Ranging (LiDAR) data was acquired to generate a Digital Terrain Model (DTM), a Digital Surface Model (DSM), and thus a Canopy Height Model (CHM). The CHM represents the heights of sugarcane plants. A selected set of vegetation indices (Response Index (RI), Excess Green Index (ExG), Normalized Blue-Green Difference Index (NGBDI), Triangular Greenness Index (TGI), Green Chlorophyll Index (GCI), and Leaf Chlorophyll Index (LCI)) were calculated for the isolated sugarcane area. After that, the correlation between field samples and the vegetation indices was checked. Chlorophyll, nitrogen and phosphorus obtained a moderate correlation (greater than 0.6), and the correlation for the potassium was close to 0.5. Phosphorus exhibited a negative correlation (−0.898) with vegetation height. The predictive models achieved an acceptable accuracy–coefficient of determination (R²) of 0.61 for chlorophyll, 0.71 for nitrogen, 0.96 for phosphorus, and 0.44 for potassium, with low RMSE and MAE errors. The accuracies of the maps of spatial distribution of chlorophyll, nitrogen, phosphorus and potassium were 60%, 75%, 96% and 50%, respectively.

The results indicated a relationship between nutritional conditions (nitrogen, phosphorus, and potassium), soil moisture and plant health. For example, healthy green areas showed adequate nutrient levels, while regions with light or brownish green leaves showed nitrogen and potassium deficiencies but higher phosphorus concentrations. This was also reflected in vegetation height, with healthy green areas reaching above 3 m (average height was 3.11 m) and areas with lighter green/brownish leaves reaching about 2 m (average height was 2.09 m). Further, spatial variability of these nutrient maps indicated a relation in differences in nutritional conditions and crop density, with 97.64% coverage in healthy green areas compared to 93.82% in areas with light or brownish green leaves, indicating the direct impact of moisture on plant nutrition due to processes such as leaching, denitrification, and mineralization. Excessive moisture was linked to water drainage issues from nearby flooded rice fields and natural terrain saturation, affecting the drainage in the study area.

In conclusion, the study developed a time- and cost-effective, technically acceptable remote sensing method to estimate the spatial variability of sugarcane nutrients. The models were tested for the sugarcane variety TW08110 with a clay loam soil texture during the irrigation season. Further, the study will help to evaluate the site-specific nutrient levels and their contribution to plant growth. For example, it was found that soil moisture, detected through the plant’s spectral information, directly impacts the reduction of nutrient content in the plant, which in turn decreases vegetation height. It can be suspected that there is a direct relationship between vegetation height and nutrients; however, this needs more investigation to confirm. It is recommended that future research include a denser leaf sampling to cover various areas of the field, perhaps a stratified random sampling scheme. The strata should be generated based on the orthometric heights of the field to identify the relationship between plant growth, especially height, nutrients and soil moisture. Also, the method should be repeated throughout the growing season to ensure consistency and investigate the impact on yield during the harvest.