Sampling Method to Estimate Fresh Fruit Bunch Production in Oil Palm Using Geostatistical Techniques

In oil palm cultivation, harvesting is a task that demands the most attention, resources, and accuracy. The logistics of harvesting activities, from bunch cutting through delivery to processing, are organized as a function of yield per hectare, and this planning depends on production samples collected from plantations to estimate the number of fresh fruit bunches that will be harvested per crop management unit. Conventionally, production is estimated through sampling grids, most commonly every five rows and every five palms (5 × 5), 7 × 7, or 10 × 10. Sampling density, however, does not typically result from a statistically rigorous method. Therefore, the objective of the present study was to define, using geostatistics, adequate sampling grids for estimating production. Specifically, seven lots in different Colombian oil palm-growing regions with planting ages between 5 and 9 years were analyzed, and the number of reproductive structures of all oil palms was recorded. From the semivariogram of each lot, different models were proposed, fitted, and cross-validated. For each sampling grid, all possible samples were simulated in the respective lots to determine the degree of accuracy in estimating oil palm production. The optimal sampling grids differed between lots, ranging from 3 × 3 to 15 × 15.

Introduction

In the last century, oil palm has become one of the most important oil crops worldwide. Among its key characteristics, productivity stands out because between 6 and 10 times more oil is produced per hectare from oil palm than from other oils [1], [2]. Colombia is the South American leader in palm oil production and ranks fourth worldwide, following Indonesia, Malaysia, and Thailand, accounting for approximately 2% of the global production [3]. Between 2016 and 2020, the Colombian area under oil palm plantations increased by 17%, from 505,965 to 590,189 ha. In 2020, the national yield averaged 3.26 tons of crude palm oil per hectare (t CPO·ha⁻¹), accounting for 6.2% of the Colombian agricultural gross domestic product (GDP) and 9.1% of the Colombian crop GDP [1]. The expansion of oil palm cultivation in Colombia can be attributed to increased global demand for vegetable oil and biofuels [4]–[6].

The intensification of oil palm cultivation in Colombia has become the key to achieving more sustainable production because sustainable farming practices can be applied to produce more palm oil and better use of both natural resources and the most efficient production factors [4], [6]. Specifically, recent efforts aimed at crop intensification have focused on narrowing the gap between the potential and actual crop yield [7].

Accordingly, a variable that should be assessed to monitor crop intensification initiatives is crop yield, expressed as tons of fresh fruit bunches (FFBs) per hectare, and this factor varies as a function of bunch count, weight, and oil content [6], [8]. Consequently, estimating FFB production is important because this information indicates the potential income that can be derived from fruits, thus enabling the allocation of resources in future activities that account for most of the production cost, such as harvesting, due to its workforce demand in cultivation, and fertilization, due to funds required for purchasing fertilizers. Together, harvesting and fertilization account for 70% of the production cost [9]. Similarly, productivity affects FFB transport logistics for processing and the use of processing plant capacity [10].

Oil palm productivity, expressed as the number of bunches or as their oil potential, is difficult to estimate given the variability within and between oil palms [11], [12]. Simulation models have been designed to characterize oil palm growth and estimate the potential crop yield; however, these models are highly demanding in terms of data requirement for their implementation and are, hence, not practical tools [2]. In Colombia, FFB production is estimated by counting female inflorescences and developing bunches—a method that was proposed in the 1950s in Malaysian plantations [10]. In this method, 5% of the oil palms in a lot are sampled twice a year. Colombian farmers use this method indistinctly, disregarding crop management, cultivar, and planting age. In other words, there is little information on the accuracy of production estimates.

Estimating oil palm production is challenging due to variability in the number of reproductive structures, which depends on the soil and climate conditions of the plantation and, consequently, in the number of bunches produced per oil palm per lot area. Accordingly, oil palm production should be estimated using spatial autocorrelation to estimate the number of palms that must be sampled for this purpose.

As mentioned above, oil palm FFB count is estimated based on systematic sampling using tools commonly known as sampling grids. These grids vary with the size of the region of interest, and the most frequently used sampling grids are every five palms and every five rows 5 × 5, 7 × 7, or 10 × 10. This sampling approach assumes that the samples are independent and that the error is random; in other words, this technique is based on the principles of classical statistics. However, these assumptions are not met if the variable of interest varies in space, leading to inadequate sampling intensity and a high degree of uncertainty [13].

In this context, other statistical methods that consider the intrinsic characteristics of FFB production must be implemented. Geostatistics is a tool for analyzing variables that vary in space and is a viable alternative to formulate accurate sampling plans [13]. This is a powerful tool to explore the underlying spatial variability of different variables. It has been used for variables such as the yield of a crop, the physical and chemical properties of the soil, and the presence of weeds, pests, and diseases, among others [13]. The study of these variables emerged in the sixties, initially focused on the mining field; Giraldo Henao [14] mentions that among the most relevant antecedents is Sichel 1947–1949, who observed the sharing of gold content in South African mines and developed basic formulas to understand its distribution. In 1951, the geologist and mining engineer Daniel Krige indicated that spatial dependence had to be taken into account [13], [15]. This theory was gradually introduced into studies of earth sciences, soil chemistry and physics, hydrology, and agronomy, among others.

In stages prior to the implementation of geostatistics, the classic statistical concepts of randomness and stratification were used to perform soil sampling [16]. For example, in the United States and Brazil, composite samples were suggested depending on the size of the area to be sampled [16], [17]. These methods have been used regardless of the management system, the type of crop (annual or perennial), or the variability of the soil [13], [18]. Classical statistics are based on the premise that samples are independent [19].

Geostatistics has had a great impact, especially in soil sampling, where the main objective is to accurately understand the state of the soil, using a minimum of samples to estimate the characteristics of interest. Understanding these characteristics allows improving crop management and can be applied to irrigation management, soil fertility and crop yield estimation [19].

Lawrence et al. [16] report research in which agricultural variables related to soil present spatial correlation. In this context, geostatistics is the tool that allows one to characterize, analyze, and explain this behavior [20]–[23]. In this order of ideas, case studies are reported on cereal crops [24]–[26], sugar cane [27], [28], mango [18], beans [29], soybeans [30], [31], and in the cultivation of oil palms [32], [33].

Rodrigues et al. [13] have proven the significant impact of geostatistics in agriculture, particularly in determining the adequate sample size and intensity for designing valid sampling plans. Knowledge of the degree of spatial variability is a parameter (known as range) allows for defining a sampling radius and, consequently, achieving an adequate sampling intensity [34]. Therefore, the present study aims to propose a method for designing a sampling plan that enables accurate estimation of oil palm crop yield, considering spatial variations, using geostatistics. Note that oil palm trees are planted at a distance of 9 meters. Strictly speaking, the variable of interest is not considered continuous in space due to the crop planting distance. Therefore, geostatistics is a practical and effective tool to capture the spatial correlation of the variable of interest, even in situations where there is no spatial continuity.

Materials and Methods

Data and Sampling Sites

The present study was based on data collected from seven lots located in four experimental fields of the Oil Palm Research Center (Centro de investigación en Palma de Aceite, Cenipalma), each of which is located in Colombian oil palm-growing regions. Lot A is located in the Palmar de La Sierra Experimental Field (Campo Experimental Palmar de La Sierra, CEPS), in Magdalena, Zona Norte [North Zone]; lots B and C are located in the Palmar de Las Corocoras Experimental Field (Campo Experimental Palmar de Las Corocoras, CEPC), in Cundinamarca, Zona Oriental [East Zone]; lots D, E, and F are located in the La Providencia Experimental Farm (Finca Experimental La Providencia, FELP), in Nariño, Zona Suroccidental [Southwestern Zone]; and lot G is located in the Palmar de La Vizcaína Experimental Field (Campo Experimental Palmar de La Vizcaína, CEPV), in Santander, Zona Central [Central Zone] (see Fig. 1). Table I outlines the most relevant characteristics of the lots. One lot was planted with Elaeis guineensis, whereas the remaining six lots were planted with oil palm O × G interspecific hybrids (E. oleífera × E. guineensis). During data collection, oil palms were growing in lots, with planting ages ranging from 4 to 9 years. Data were collected only once in 2020 by counting visible female inflorescences and bunches of all oil palms in a lot through production surveys. Additionally, the location of each palm was recorded as the East-West ($x$) and North-South $(y)$ coordinates.

Fig. 1. Location of study lots in Colombia’s oil palm-growing regions.

Table I. Details of the Studied Lots and Descriptive Statistics for Number of Bunches from Individual Palms

Lot	A	B	C	D	E	F	G
Colombian region	Northern	Eastern		Southwestern			Central
Planting age	4 years	5 years	6 years	7 years			9 years
Crops	O × G hybrid	O × G hybrid	Elaeis guineensis	O × G hybrid			O × G hybrid
Hectares (ha)	5.3	6.3	15.6	4.7	1.9	4.3	3.4
Number of palms	751	905	2235	677	285	611	475
Mean	3.6	10.1	10.5	9.4	10.2	8.8	7.3
Median	4	10	10	10	11	9	8
Minimum	0	3	3	0	0	0	0
Maximum	10	16	20	17	16	16	17
SD	2.3	2.3	4.4	2.5	2.4	3	3
CV (%)	62.9	22.8	41.8	27	23.3	33.8	41.1

Descriptive Data Analysis

Descriptive data included summary statistics (mean, median, minimum and maximum values, standard deviation, and coefficient of variation) and exploratory plots. The objective of this analysis was to characterize variations in the variable of interest, disregarding spatial variation. Variations in reproductive structures (bunches and female inflorescences) in each lot were identified, carefully searching for possible atypical data. Scatterplots of variations in the variable of interest as a function of the corresponding coordinates were created to assess the possible correlation between the variable and the location of the oil palm in the lot. If correlated, this variation would require modeling because the variable of interest should be homogeneous across the study area.

Structural Data Analysis

Structural data analysis allowed for identifying the spatial correlation structure of the characteristics of interest. Through geostatistics, the semivariance function was estimated to assess whether the variable of interest varied in space and to quantify and model this variable because semivariance indicates how similar or dissimilar the values of the variable are as they move away from a given point in any direction [35]. Semivariance was estimated using the method of moments, which estimates the variance of the characteristic of interest between all pairs of individuals located at different distances [36], for instance, at 9, 18, 27, 36, 45 m, or longer, depending on the size of the study area (note that these distances are multiples of nine because the oil palms are planted in equilateral triangles 9 m apart). The variance of the number of reproductive structures of all pairs of individuals (palms) at a distance h is known as the empirical semivariance function, and its classical estimation is shown in (1) below:

where $N\left(h\right)$ is the set of all pairs of individuals separated by distance $h$; $|N\left(h\right)|$ is the cardinal number of $N\left(h\right)$,¹ 1 The cardinal number of a set indicates the number of elements that belong to that set. and $Z_x$ and $Z_{x+h}$ are the observed values at positions $x$ and $x+h$, respectively [37]. The graphical representation of these estimated semivariance values as a function of distance is known as the semivariogram, and it is shown in Fig. 2.

(1) $$\begin{array}{r}\hat{\gamma }\left(h\right)=\frac{1}{2|N\left(h\right)|}{\sum }_{N\left(h\right)}{[Z_x-Z_{x+h}]}^2\end{array}$$

Fig. 2. Empirical semivariogram and its parameters: nugget, sill, and range.

The empirical semivariogram has three parameters, namely nugget, sill, and range [36]: (i) The nugget ($c_0$) is the variance unexplained by the model and is presented as a point discontinuity of the semivariance value at the origin (zero distance), which may result from measurement errors in the variable of interest or its measurement scale. Occasionally, a nugget is caused by the microscale effect, which occurs because the spatial structure of the variable is partly concentrated at distances smaller than the working distances [38]–[40]; (ii) The sill corresponds to the asymptote at which semivariance stabilizes and represents the total variation of the variable [41]; (iii) Finally, the range corresponds to the distance at which semivariance stabilizes and determines the point from which no spatial autocorrelation occurs. That is, the recorded values of the variable of interest are spatially independent from this distance [42].

Because the classical empirical semivariogram is sensitive to outliers, a robust empirical semivariogram was also used to determine which of the two semivariograms better fitted the estimated semivariance data [36]. The maximum distance for reliable representation of the empirical semivariogram was determined, ensuring special care when identifying the range. Since spatial autocorrelation must be repeated in different directions in space [40], the empirical semivariance was estimated in the following directions: 0°, 45°, 90°, and 135°. If a similar behavior is evidenced in shape and scale, the spatial autocorrelation is said to be isotropic. That is, it does not depend on the estimation direction. Once isotropy was confirmed, the theoretical model that best fits the empirical semivariance was determined for each lot based on the semivariograms.

Theoretical Data Modeling

As shown above, empirical semivariances were estimated at some locations in a lot, thus requiring fitting models that generalize the variation observed in the semivariogram to any possible distance, considering the area and shape of the lot [14]. Several theoretical models can be used to quantify spatial variability from an empirical semivariogram. Some of the best-known theoretical models are Exponential, Gaussian, and Circular models, among others [14], [36]. The precision and validity of this analysis depend on the appropriate choice of model. The Gaussian model is based on the normal distribution and is useful for modeling isotropic spatial dependence. This model has been implemented in ecology studies and soil studies [43], [44]. The Exponential model assumes that spatial autocorrelation decreases exponentially. In the literature, applications of the Exponential model are found in variables associated with soils, such as permeability and pH [45]. The Gneiting model is a flexible semivariance model that allows modeling spatio-temporal processes. The Wave model (or sine hole effect) allows incorporating wave patterns to model the semivariance and applications in geology and geophysics are reported. A theoretical comparison of the models is found in [41], [46], [47]. Once the initial values of the parameters of the semivariance function were estimated, the theoretical models were fitted using five estimation methods, namely ordinary least squares (OLS), ordinary weighted least squares (OWLS), and robust weighted least squares (RWLS), maximum likelihood (ML), and restricted maximum likelihood (RML).

Cross-validation was used to assess the goodness-of-fit of the theoretical models proposed for semivariance [48]. This method involves sequentially omitting the value of the study characteristic at a location to predict its value at the excluded location using other spatial locations and the selected semivariance model [14]. If the estimated semivariance model correctly fits the spatial dependence structure, the difference between the observed and predicted values should be as small as possible. This value is known as prediction error. The model with the lowest prediction error was selected, and its parameters were interpreted as the best data extracted from reality.

Based on the spatial dependence structure of all reproductive structures of each palm, the sample size and adequate sampling intensity were determined. This sampling plan was determined using the range value derived by modeling semivariance, as this parameter allows for defining a sampling distance that ensures sample independence [34]. Range is the distance from which the number of female reproductive structures of any selected palm in a lot is spatially independent of the number of female reproductive structures of other palms. In other words, range represents the distance at which bunches and female inflorescences should be sampled to estimate the total number of future bunches in the lot (i.e., FFB production). Since distance should be maintained in all directions, a sampling grid is generated with the selected distance. Note that the sampling grids are set departing from each possible palm at the lot, so that all possible samplings are considered when changing the starting point. For each simulated sample, the FFB production per palm is estimated. The goal is to determine the degree of precision when the estimated value is compared to the actual value of FFB production per palm. Data were analyzed using RStudio [49].

Results and Discussion

Descriptive Data Analysis

Table I summarizes the descriptive statistics for the number of reproductive structures per palm and the most relevant characteristics of the lots. Variation was quantified based on the coefficient of variation (CV). Lot A presented high variation (63%), followed by lot C (42%) and G (41%); the remaining lots showed moderate variation, with CV ranging from 23% to 34%. The median number of bunches per palm was 10, which is typical of O × G hybrid cultivars [50] except in lot A, where the median was 4 bunches per palm (see Fig. 3). The minima and maxima showed an amplitude ranging from 0 to 20 bunches per palm. This range may lead to the over- or underestimation of the actual value of FBB production from a lot without adequate sampling.

Fig. 3. Behavior of the number of bunches/palm along the lot.

Variations in the number of FFB per palm in each lot are presented in Fig. 3. Overall, lot A differed from the other lots and was the least productive because this lot was subjected to water deficit—a factor that directly affects the production of female inflorescences in oil palm [10], [51]. Each lot showed different soil, climate, and management characteristics, which affect the variability of FFB/palm [52]. Accordingly, each lot must be studied as a specific case by determining the appropriate sampling density to ensure sample representativeness.CV represents the degree of homogeneity of FFB/palm in each lot, with a low value indicating low heterogeneity in the area of interest. CV allows for analyzing the variation within the lot but not the spatial variation in the production of reproductive structures.

Structural Data Analysis

To characterize the FFB production behavior per palm in each lot, the data were graphically represented as a function of the coordinates in scatter diagrams (Fig. 4). A non-parametric stationarity test was carried out in each lot [53], and it was found in lot D that the assumption was rejected (p-value = 0.006) with a significance level α = 0.05, indicating a trend of the production of FFB throughout the lot. To correct this tendency, the median polishing method was used [54]. In the rest of the lots, no trend was identified, indicating that the production of bunches per palm was homogeneous throughout the study area and allowing for a direct assessment of whether this variable showed spatial correlation.

Fig. 4. Spatial distribution of bunches/palm along the lot.

Semivariograms were plotted in different directions: 0°, 45°, 90°, and 135°. It was found that lot E does not meet the isotropy assumption. Given this result and the independence mentioned above, this lot (i.e., lot E) was not taken into account for the spatial modeling of the semivariance. The rest of the lots showed a similar shape and scale, regardless of the angle at which they were estimated (see Fig. 5). The semivariance increased in the first 130 m and then stabilized, indicating that the production of female inflorescences and bunches in the oil palm varies in space. Furthermore, theoretical models that fit the empirical semivariance data were estimated to explain the spatial dependence structure of the data.

Fig. 5. Estimation of the empirical semivariance in different directions to testing the isotropy assumption.

To check whether FFB production presents spatial autocorrelation, confidence bands were estimated through permutations (9999 simulations). These results are presented in Fig. 6. If the estimated semivariance has at least one point outside the confidence bands (blue bands), the spatial process presents autocorrelation. To support this analysis, the independence hypothesis proposed by [55] was tested, and it was found that for lot E, the independence hypothesis was not rejected (p-value = 0.25) with a significance level $\alpha =0.05$.

Fig. 6. Testing of a hypothesis for spatial independence. If the estimated semivariance is contained in the blue band, there is no statistical evidence of spatial autocorrelation.

Fig. 7 shows the empirical semivariance and estimation of theoretical models for each lot. OWLS and RWLS were the estimation methods used to adequately fit the semivariance data, whereas the ML and RML methods did not produce good results. Table II outlines the fitted models for each lot, as determined using the cross-validation method. Lots C and G showed the largest distances of spatial dependence (range < 70 m), whereas the shortest autocorrelation distance was recorded in lot A (range < 25 m). In the other lots, the spatial autocorrelation distance ranged from 43 m to 52 m.

Fig. 7. Fitted models of semivariance for lots.

Table II. Details of the Fitted Models for Semivariance in Studied Lots

Lot	A	B	C	D	F	G
Model	Gneiting	Circular	Exponential	Wave	Gaussian	Wave
	4 years	5 years	6 years	7 years		9 years
Sill	4.59	4.29	21.06	4.71	7.5	8.67
Range	21.46	49.39	139.56	43.33	44.53	77.81
Nugget	3.5	2.43	14.89	2.59	4.63	6.21
Nugget-to-sill ratio (%)	76	57	70	55	62	72

In all cases, a nugget different from zero was estimated; this is mainly due to the fact that the minimum distance at which the palm is planted is 9 m, and the nugget captures the behavior of the spatial microstructure. However, the distribution of model variability and spatial autocorrelation in each lot were assessed using the association between the nugget and sill effects. If the ratio is below 25%, the nugget is smaller than the sill, and almost all variance is explained by the model. That is, there is strong spatial dependence. In contrast, the larger the nugget, the greater the variance that is no longer explained by the model and the weaker the spatial dependence. According to a ratio between 25% and 75% is considered to represent moderately strong spatial dependence, while a ratio exceeding 75% is considered to represent weak spatial dependence [56]. Except in lot A, the fitted models showed nugget-to-sill ratios between 55% and 75%, indicating moderately strong spatial dependence. Although the correlation forces found in this study are moderate to low, this structural behavior has not been recorded in the literature to date. This result offers valuable information not only for sampling designs but also for modeling FFB production in the field.

The number of samples required to estimate FFB palm can be expressed as a function of the estimated range; as such, the range is inversely proportional to the number of samples because the shorter the distance of spatial autocorrelation, the greater the variation in the area and the number of field samples required to ensure its representativeness. To ensure the spatial independence of samples, the variable of interest should be measured at a distance longer than the range. Therefore, for lot A, samples should be collected at a distance of 22 m, which is equivalent to 11% of all palms in this lot; therefore, lot A required the highest number of samples. Lots B, D, and F with O × G hybrid cultivars at planting ages of 5 and 7 years required a sampling distance of >43 m, that is, 4% of all palms in the lot. Finally, lot G, with 9-year-old palms, required a sampling distance of 78 m, and lot C, the only lot with E. guineensis, required a sampling distance of 140 m. These sampling distances correspond to 1% and 0.6% of palms in the lot, respectively.

Based on the estimated sampling distances and planting density in the lots, a sampling grid was proposed to gather all spatial data on FFB palm⁻¹. The sampling grids ranged from 3 × 3 to 15 × 15 and mainly depended on variability in the lot, area, and crop age. Note that the sampling grids are set departing from each possible palm at the lot, so that all possible samplings are considered when changing the starting point. For each simulated sample, the FFB production per palm is estimated. As shown in Table III, the mean number of estimated bunches was close to the actual number. To measure the precision of the proposed sampling grid, a hypothesis test is performed in each simulation. The purpose of this test is to compare whether the estimated average FFB per palm is equal to the actual average. Non-parametric methods are used for this test, with a significance level of 5%. It is observed that the sampling meshes yield errors lower than 10% (see Table III). Therefore, these results highlight the need for sampling grids that gather the spatial variability of all data from the lot and generate a representative sample of FFB production. The sampling design proposed here ensures that essential statistical assumptions are met, such as independence and randomness, thereby avoiding under- or overestimation of FFB production.

Table III. Simulation of Sampling Grids by Lot

Lot	A	B	C	D	F	G
Region of Colombia	Northern	Eastern		Southwestern		Central
Planting age	4 years	5 years	6 years	7 years		9 years
Crops	O × G hybrid	O × G hybrid	Elaeis guineensis	O × G hybrid		O × G hybrid
Range	21.46	49.39	139.56	43.33	44.53	77.81
Sample grid	3 × 3	5 × 5	15 × 15	5 × 5	5 × 5	9 × 9
Bunches per palm-actual	3.6	10.12	10.54	9.43	8.77	7.27
Bunches per palm-estimated	3.59	10.14	10.56	9.37	8.72	7.28
Confidence intervals (95 %)	3.10–4.09	9.38–10.9	7.16–13.95	8.48–10.25	7.51–9.92	6.25–8.31
Error	0.08	0.08	0.08	0.04	0.04	0.07

Conclusions

This document provides a statistical tool to define sampling grids. The results found are valuable because they show that it is not possible to define a universal sampling mesh and that the specific conditions of each lot must be taken into account. The benefits of doing so are the savings in terms of harvest logistics and the determination of the amounts of fertilizer that must be replaced into the soil. Note that both activities (harvest and nutrition) account for 70% of FFB production costs [9].

In oil palm-growing regions, the method used to estimate FFB production is applied regardless of the type of crop, planting year, management characteristics, and climatic variables without ensuring the representativeness and accuracy of the estimate and, hence, inadequately estimating FFB production. In contrast, the method proposed here showed satisfactory results, with estimation error of <10%, allowing for designing sampling plans based on variation in production in the lot. Therefore, specific analysis by area is warranted to determine the optimal sampling distance that ensures accurate estimation of the variable of interest.

The results of the present study demonstrate that the production of bunches and female inflorescences per oil palm shows intrinsic spatial variability that must be considered in future analyses. The study lots showed moderate-to-weak spatial dependence. However, this spatial dependence may be stronger if other variables that help characterize the variability of production, such as agronomic, climatic, or lot management variables, are included in the analysis.

Overall, by analyzing the spatial variation in FBB production, sampling grids were identified for each study lot. Most lots with O × G hybrid cultivars required 5 × 5 sampling grids. The lot with the youngest plants required a larger sample size and a denser sampling grid (3 × 3 sampling grid). Based on these results, further studies should be conducted with lots under the same conditions to better understand their variations. Finally, using the proposed method, the lot with E. guineensis required the longest sampling distance (10 × 10 grid), although the number of palms that must be sampled was nevertheless 88% lower than that required using the conventional method.