High-Resolution Yield Map** for Eucalyptus grandis—A Case Study

Abstract

Yield data represent a valuable layer for supporting decision-making as they reflect crop management results. Forestry decision-makers often rely on coarse spatial resolution data (e.g., forest inventory plots) despite the availability of modern harvesters that can provide high-resolution forestry yield data. The objectives of this study were to present a method for generating high-resolution Eucalyptus grandis yield data (individual tree-level) and explore their applications, such as correlation analysis with soil attributes to aid nutrient recommendations. Two evaluations were conducted at two sites in Brazil: (a) assessing the positioning accuracy of the global navigation satellite system (GNSS) receiver positioning, and (b) analyzing the yield data and their correlation with the soil attributes. The results indicated that positioning the GNSS receiver at the harvesting head provided higher accuracy than placement at the top of the harvester cabin for individual tree-level data. Reliable yield data were generated despite the GNSS receiver’s increased susceptibility to damage when mounted on a harvest head. The linear correlation analysis between the Eucalyptus grandis yield data and soil attributes showed both negative (Clay, B, S, coarse sand, and potential acidity − H + Al) and positive correlations (K, Mg, pH-SMP, Ca, sum of bases, pH, base saturation, fine sand, total sand, and silt content). This study demonstrates the feasibility of obtaining high-resolution yield data at the individual tree-level and their correlation with soil attributes, providing valuable insights for improving forestry decision-making.

1. Introduction

Yield data are considered the result of interactions among environmental and anthropogenic factors, providing valuable information for many decision-making processes, such as nutrient removal maps, profitability maps, management zone establishment, and others [1].

Yield monitoring systems are commonly found for annual crops like maize [2], rice [3], soybean [4], and wheat [5], as well as for semiperennial crops like sugarcane [6]. However, perennial crops such as Eucalyptus spp., lack these systems due to challenges associated with automated yield accountability. Unlike others crops, the yield in forestry involves volumetric measurements of the tree trunk rather than mass measurements.

The most common method for obtaining the yield measurements for Eucalyptus spp. or any forest is through forest inventory plots. The biometrics methods employed for forest plantation measurements involve assessing the tree features (diameter, total height, and shape) within rectangular or circular plots of varying sizes, usually ranging from 200 to 600 m². These plots are distributed randomly or systematically within the population, with the sampling density per one plot from 5 to 20 ha within a stand [7]. The method relies on the manual use of a tree caliper and hypsometer that are costly, time-consuming, and provide coarse spatial resolution.

In this sense, to overcome these issues, advances in forest inventory have been made towards the use of orbital, aerial and terrestrial Light Detection and Ranging (LiDAR) sensor systems, resulting in high-spatial-resolution data [8,9,10]. However, there are limitations for data generation. Orbital and aerial systems face problems related to tree density and leaf cover. Terrestrial systems have major issues associated with occlusions, terrain, and aboveground objects.

Given the significance of yield data for foresters and resource managers, improved models to understand or predict the yield and its components are in development. However, for reliable yield prediction equations, regardless of the approach (individual tree [11] or stand-level [11]), it is necessary to obtain ground-truth data in the field, which are commonly obtained by manual measurements. Even for LiDAR sensor systems, it is necessary to calibrate with manual measurements as a methodological procedure [12].

Yield data are the most common variable evaluated in field experiments [13,14,15,16,17]. However, the majority of these experiments assess the yield based on low-resolution data, which underrepresent the crop’s spatial behavior in the field. This approach contradicts the principles of precision agriculture, which aim to enhance the resource use efficiency, productivity, quality, profitability, and sustainability of agricultural production by considering temporal, spatial, and individual plant data [18].

The development of modern harvesters and forwarders offers a potential solution for acquiring high-resolution yield data, specifically at the individual tree-level. These machines can record data under a standard known as the StanForD standard, using a specific file format (“.stm”) [19]. However, in the forestry industry’s routine operations, instead of utilizing the individual data of harvested trees, summarized data from each stand are often used. This presents an opportunity for the researchers engaged in big data analysis to gain new insights based on overlooked data, such as the timber log diameter and length information from the mechanical harvesting process.

Considering that the methods to obtain high-spatial-resolution data are still in development and acknowledging the pressing need to enhance the crop production efficiency to align with the sustainable development goals [20], it becomes crucial to minimize the errors related to crop management practices. Particularly for Eucalyptus spp., a long-term crop [21], the common approach to obtaining yield data is still based on a coarse spatial resolution (e.g., inventory level with about one sample for every 400 m²). There is a gap in improving the efficiency of the managerial practices supported by the use of yield data layers.

In this context, exploring the data produced by harvesting systems is promising. This study explores the hypothesis that it is feasible to obtain high-resolution yield data using modern harvesters. The objectives of this study were to demonstrate a methodology for generating high-resolution yield data for Eucalyptus grandis and its possible applications, such as correlation analysis with soil attributes to aid nutrient recommendation.

2. Materials and Methods

2.1. Description of the Study Site

The study was conducted at two experimental sites cultivated with Eucalyptus grandis, located at the State of Sao Paulo, Brazil, that belong to Suzano SA Company. Site 1 coordinates are 23°42′31′′ S and 48°22′21′′ W and site 2 coordinates are 22º54′57′′ S and 48º14′33′′ W (Figure 1).

First forest site has an average altitude of 650 m above mean sea level and is covered by Red Oxisol with a clayey texture [22]. This site has an average annual rainfall of 1350 mm and an average annual temperature of 19 °C, characteristic of a humid subtropical climate with hot summer and without dry season [23]. Second forest site has an average altitude of 580 m and is covered by quartzipsamments with a sandy texture. This site has an average annual rainfall of 1200 mm and an average annual temperature of 21.5 °C, characteristic of a tropical climate with dry winter [23].

Both forest sites studied are within one of the most traditional forestry regions in Brazil, renowned for being one of the most favorable areas in the world for successful forestry due to its high suitability for dozens of Eucalyptus species [24].

2.2. Experimental Design (GNSS Placement) and Harvesting Procedures

Site 1 was planted in 2010 across an area of 11.07 ha, with a row and plant spacing of 3.30 m and 2.27 m, respectively, and mechanically harvested in 2018. Site 2, spanning 18.56 ha, was planted in 2011 with a row and plant spacing of 3.80 m and 2.10 m, respectively, and was mechanically harvested in 2019.

The harvesting system adopted by eucalypt-based pulp companies is typically cut-to-length logging. In this system, mechanical harvesters are used for felling trees, delimbing (removing branches), bucking (cutting the tree into logs of specific lengths), and arranging these timber logs into bunks for transportation. In this study, yield data were obtained using a harvesting head (model H77euca—Ponsse, Vierema, Finland) coupled with a hydraulic excavator (model PC200, Kotmatsu, Sao Paulo, Brazil). The harvester head, equipped with a sensor system integrated into a multifunctional data logger and a global navigation satellite system (GNSS) receiver (model Opti4G 4.715—Ponsse, Vierema, Finland), initially measures the tree’s diameter at breast height (DBH). Subsequently, during delimbing and bucking, it measures the diameter and length of each cut wood timber log to calculate the tree volume (m³). The minimum diameter processed (measured) by the harvester head is 0.03 m, meaning that harvest obtains the commercial tree volume.

For georeferenced data, testing the placement of the GNSS receiver on the harvester was necessary to ensure accurate data collection by tree. At site 1, the receiver was installed at the top of the harvester cabin as per the default protocol (Figure 2A), which could potentially compromise the accuracy of tree distribution data. At site 2, to improve the precision and accuracy in capturing individual tree data, the GNSS receiver was placed at the harvesting head (Figure 2B), accounting for radial movement and accurately representing the real position of individual harvested trees.

A tree caliper (Ponsse, Vierema, Finland) was used to calibrate the harvesting sensor system. The tree caliper provided manual measurements and records of the diameters and lengths of the timber logs from sampled and harvested trees. After the measurement, it was connected to the harvester yield monitor, feeding the sensor system in a self-calibration mode.

2.3. Sampling Procedures

In the second site, aiming to explore the linear correlation between tree yield and soil physical and chemical attributes, 54 sampling points (about three sampling points per ha) were collected at a depth of 0.00–0.20 m. Each sampling point comprised four subsamples: two within the row spacing and two within the plant spacing. These samples were subjected to physical (clay, sand, and silt contents) and chemical (Al: aluminum, B: boron, Ca: calcium, CEC: cation exchange capacity, Cu: copper, Fe: iron, H + Al: potential acidity, K: potassium, Mg: magnesium, Mn: manganese, OM: organic matter, pH: potential of hydrogen—calcium chloride (CaCl₂), pH-SMP: potential of hydrogen—Shoemaker–McLean–Pratt, S: sulfur, and Zn: zinc) analyses and coincided with the inventory points. Given the higher density of yield points compared to the soil sampling points, a buffer with a 20 m radius was applied using the soil sampling as a reference. This buffer aimed to obtain an average yield value for each soil sampling and to approximate the overlap between the inventory and buffer (Figure 3).

2.4. Data Analysis

The evaluation of the yield estimation values obtained by the harvesting system involved two approaches. The first method, inherent to the harvester system as described earlier, involves summing the volumes of the timber logs to calculate the individual volume of the harvested tree. The second method is based on the Schumacher and Hall (1933) [22] equation (Equation (1)), which estimates the individual volume of the felled tree using measurements of DBH and commercial height obtained at the harvester.

ln(V) = β₀ + β₁ × ln(DBH) + β₂ × ln(HT) + ε_i,

(1)

where ln = natural logarithm, V = tree volume, β_i = parameters to be estimated, DBH = diameter at the breast height, HT = total height, and ε_i = random error.

The equation was adjusted based on precise measurements from the same trees felled and evaluated for self-calibration of the machine. Both yield estimation methods were compared with manual measurements conducted on forest inventory plots. This step aimed to evaluate the best method for generating yield maps from harvester data. The comparison was conducted using the forest inventory as the reference, which involved measuring DBH, commercial tree height, and volume of all the trees within an 11.28 m radius but, in this case, at a high density of two plots per hectare. Analysis of variance followed by the Tukey multiple comparison test at a significance level of 5% were applied to compare the tree volume.

The best model was considered the one that presented similar yield data to the forest inventory and the lowest average error. After choosing the model, a data interpolation and filtering process was applied. Due to the absence of consensus among researchers on the yield data filtering process, values considered unreal to the scenario were removed. In this case, values below 70 m³ ha⁻¹ and above 421 m³ ha⁻¹ were eliminated, with a standard deviation (σ) equal to 0 (indicating that the observed and predicted datasets are equal, which would be unlikely because it is a continuous variable).

The Pearson correlation (r) analysis was conducted at a significance level of 5% for the soil sampling results and yield data. Variables that presented significant correlations were used to generate surface maps with a spatial resolution of 3 m × 3 m. Data extraction from the StanForD and analysis were carried out in RStudio [25], geostatistical analysis was conducted in Vesper [26], and surface maps were generated in a geographic information system (QGIS 3.16) [27]. In the presence of spatial dependence, ordinary kriging was applied; otherwise, inverse distance weighting was used [28].

3. Results and Discussion

Figure 4 shows the yield data obtained by different GNSS receiver placements. It demonstrates that the GNSS receiver positioned at the top of the harvester cabin mainly captured data related to the harvester’s path, which jeopardizes the generation of accurate yield maps as the points scarcely represent the tree’s location (Figure 4A). On the other hand, positioning the GNSS receiver at the harvesting head provided the most accurate tree-specific yield data (Figure 4B). However, due to the GNSS receiver’s error (C/A code, port L1) of approximately 5 m, it also affected the data’s position concerning the planting line. Enhancing the accuracy would necessitate a more accurate and precise receiver. Despite these challenges, generating yield maps with the GNSS receiver at the harvesting head is feasible. Nevertheless, caution is warranted as the device is located in a region susceptible to greater physical impact during the harvesting operation.

Considering the data from site 2, the volumetric yield measurements obtained from the forest inventory were compared to the volume estimated from two harvesting system approaches: sum of timber log volumes and Schumacher and Hall’s equation. The analysis of variance and Tukey’s test did not show statistically different differences (p > 0.263). However, Schumacher and Hall’s equation presented a lower average error (−0.001 m³) compared to the data from the yield monitoring system (0.162 m³). Therefore, the yield data were generated by applying Schumacher and Hall’s equation with data from the harvester as input variables (DBH and commercial height).

Figure 5 shows the yield map obtained from the harvester after the ordinary kriging and filtering process. Figure 5A shows the spatial distribution of the yield data, potentially at the individual tree-level. After the interpolation and filtering process (Figure 5B), the surface maps facilitated the visualization of the yield data, enabling decision-makers to identify patterns.

Based on the proposed setup, crop decision management should utilize high-spatial-resolution yield data (individual tree-level—plant and row spacing, e.g., 3.80 m × 2.10 m ~ 8 m²) and accept lower errors instead of the higher errors derived from low-spatial-resolution data (e.g., inventory level –3.14 × 11.28 m × 11.28 m~400 m²) [21]. Furthermore, using high-spatial-resolution yield data (e.g., individual tree-level) can enhance and support the estimates of wood stock produced by individual tree models, which have been under development since the 1970s [11,29,30].

Moreover, precise and accurate yield data can mitigate the uncertainties arising from the discrepancies between the traditional field measurements conducted during inventory assessments and the actual volume of wood entering the mill. The traditional methods rely on coarse spatial resolution estimations, which can negatively impact the representation of reality. For instance, even with high-resolution yield map** (Figure 6), gaps between data points can be observed. After interpolating the data, it is expected that areas with closer data points will have lower standard deviation errors, while areas with more dispersed data points will exhibit higher standard deviation errors. Consequently, lower-spatial-resolution data have a higher probability of generating larger measurement errors compared to higher-spatial-resolution data.

Based on Figure 6, the soil attribute variables that were positively significant at 5% with yields include K, Mg, pH-SMP, Ca, sum of bases, pH, base saturation, fine sand, total sand, and silt content. The negatively significant variables include clay, B, S, coarse sand, and potential acidity (H + Al). To avoid collinearity, the variables with strong correlations (r > 0.82) were simplified. For example, sum of bases, Ca, pH (CaCl2), and base saturation were strongly correlated. Hence, only Ca was used for the analysis because sum of bases and base saturation are directly related to Ca, K, and Mg. Although total sand content did not show a significant correlation with fine and coarse sand content, it was disregarded in the study for creating surface maps due to its strong negative correlation with clay content (r = −0.96). Appendix B contains the correlogram with the r values (Figure A1).

Considering the chemical soil attributes, most of the significant variables are commonly explored in forestry studies, especially with Eucalyptus spp. [13,14,15,16]. For example, calcium presented a linear and positive relationship with yield (r = 0.26). This correlation is expected because highly productive Eucalyptus plantations require large amounts of Ca [14]. This accumulation of Ca in leaves increases over time until leaf abscission, which differs from other cations. After canopy closure, the concentration of other cations tends to reduce due to internal translocation to growth zones [31].

Despite knowing that Eucalyptus growth is variable and influenced by soil properties along the field or sites, or both [17], the use of high-resolution yield data is still minimal, thus jeopardizing crop management efficiency. This might occur because producers consider larger errors acceptable under the premise that Eucalyptus spp. is a long-term crop [21] or they lack access to high-resolution yield map**, or both.

In the literature, a consistent and positive correlation between regional Eucalyptus spp. yield and clay content is observed, with values ranging from 100 to 650 g 1000 g⁻¹ [32,33]. However, our study reveals a contrasting negative correlation. It is crucial to note that the clay content in our investigation ranged from 64.10 to 92.63 g 1000 g⁻¹ (Figure 7A), falling outside the range reported by Alvares (2011) [32] and Gonçalves et al. (2012) [33]. This result underscores the importance of more comprehensive studies, encompassing a broader clay content range, ideally spanning soils with textures ranging from sandy to very clayey. Such nuanced investigations would significantly contribute to the understanding of the intricate relationship between clay content and Eucalyptus spp. yield.

The range of values of the soil chemical attributes found in this study can be classified as “adequate” only for Ca⁺² and Mg⁺² and “low” for K⁺, according to the class of interpretation from Neto et al. 2020 [13], which established a critical level of nutrients for eucalypts in soil and plants. Thus, aware of the variability in the soil physical (Figure 7) and chemical attributes (Figure 8) within the field and the capacity to obtain high-resolution yield data, crop decision-making must be improved, such as with enhanced nutrient recommendation considering both the yield and soil attributes content.

The attributes that presented spatial dependence (clay—Figure 7A, coarse sand—Figure 7B, silt content—Figure 7C, fine sand—Figure 7D, B—Figure 8C, K—Figure 8E, and pH-SMP—Figure 8F) were submitted to ordinary kriging. As a result, the surface maps are smoother in comparison to those without spatial dependence (Ca, Mg, and S), which were submitted to the inverse distance weighting. Appendix A (

Table A1. Semivariogram parameters obtained for ordinary kriging.

Variable	n	C0	C1	A1	Variogram Model	RMSE	AIC
Boron	54	7.40 × 10−04	5.91 × 10−04	145.60	Gaussian	6.49 × 10−05	−163.80
Clay	54	84.23	252.70	462.50	Gaussian	28.07	308.10
Coarse Sand	54	0.00 × 1000	1.00 × 10+04	185.00	Gaussian	1326.60	427.40
Fine Sand	54	0.00 × 1004	1.00 × 10+04	184.70	Gaussian	1154.40	420.70
Potassium	54	4.30 × 10−04	2.59 × 10−04	294.90	Spherical	5.69 × 10−05	−198.80
pH-SMP	54	0.05	0.06	380.60	Gaussian	0.01	−104.10
Silt	54	5.33	7.55	198.20	Gaussian	1.44	75.67

n = number of samples, C0 = nugget effect, C1 = structural variance, A1 = range, RMSE = root mean square error, and AIC = Akaike Information Criterion.

) shows the semivariogram parameters used for ordinary kriging.

Clay (Figure 7A) and coarse sand content (Figure 7B) presented an inverse pattern in comparison to the yield (Figure 5B), indicating that areas with higher values of clay or coarse sand content tend to negatively affect the yield. Moreover, the opposite effect was found for silt content (Figure 7C) and fine sand content (Figure 7D).

Looking at the chemical surface maps (Figure 8), they show the spatial variability in the attributes (Ca, Mg, B, S, K, and pH-SMP). Aware of the variability and its significance with yield, it can guide or support decision-makers, or both, regarding the nutrition recommendation to improve the nutritional management practices for the area, avoiding over- and/or under-nutrient application to reach high-yielding conditions. For example, using the layer of Ca and pH, it can improve the efficiency of the lime application to meet the crop requirements.

In sites where the yield is not responsive to fertilizer applications but presents field variability, the use of high-resolution yield data can guide fertilizer recommendations based on the nutrient removal maps. This approach helps to prevent soil nutrient depletion [34] and plant deficiency [34]. For example, the Brazil Eucalyptus Potential Productivity (BEPP) project found that the yield was not responsive to high rates of fertilization (e.g., 500 to 1000 kg of N ha⁻¹) compared to the current rates (15 to 70 kg of N ha⁻¹). This finding underscores the importance of evaluating the field variability to optimize the fertilization recommendations [35].

Although soil nutrient depletion is generally not considered a risk in Eucalyptus plantations [36], there are still nutrient-related issues to consider, such as the overapplication of nutrients and the creation of nutrient pools. Evaluating nutrient pools requires multiyear measurements from the soil and forest inventory [37], which are laborious and expensive. However, by leveraging high-resolution yield data, it is possible to undertake more informed and efficient fertilizer applications, thereby reducing the risk of overapplication and mitigating the potential environmental impacts.

This study presented a method to obtain high-resolution yield data from a harvester and demonstrated their applicability through pair-wise evaluations with soil attributes. The high-resolution data enhanced the understanding of the crop’s spatial behavior and its possible relationship with soil attributes, which can contribute to achieving the sustainable development goals [20]. However, the correlation analysis was limited to the results from soil cores taken at a depth of 0.00–0.20 m. Future studies could improve this by including deeper sampling depths and higher sampling densities for further evaluation, as found in [17].

Since this study did not focus on the effects of nutrient treatments, the conclusions regarding the relationship between soil attributes and yield are limited to the studied area. Nevertheless, the dataset obtained can provide valuable insights and guidance for experimental designs at the farm level, such as on-farm experimentation.

Our case study highlights the common challenges faced by the planted forest sector. Primarily, there is a prevailing mindset among technicians and managers that high-resolution information is valuable, yet there is minimal effort undertaken to implement it. Deploying this study to other contexts must consider challenges such as the difficult operating conditions in forest harvesting environments, the risk of equipment damage from residues, and the high cost of precision agriculture equipment. Additionally, there is a shortage of qualified forestry technicians familiar with precision agriculture concepts, and the lack of standardization among the equipment and operational procedures further limits widespread adoption.

Future research efforts should prioritize the evaluation of the multiyear measurements acquired from the yield harvester system across diverse locations. These efforts should focus on proposing advanced decision-making frameworks that leverage high-resolution data layers, such as the development of management zones or the implementation of site-specific application strategies. Moreover, incorporating yield data layers can significantly enhance both individual and stand-level modeling approaches.

4. Conclusions

Our case study demonstrates the feasibility of generating high-resolution yield data for Eucalyptus grandis using modern harvesting technology. By positioning the GNSS receiver at the harvesting head, we obtained accurate yield data at the individual tree-level, overcoming the challenges associated with the traditional forest inventory methods. Additionally, our comparison of the yield estimation approaches revealed the effectiveness of Schumacher and Hall’s equation in estimating the tree volume using the harvester data as the input.

The availability of high-resolution yield data reveals new possibilities for enhancing the forestry decision-making processes significantly. By correlating yield data with soil attributes, we identified significant relationships that can inform site-specific management practices and improve resource allocation. Furthermore, the use of high-resolution yield maps can aid in optimizing fertilizer application, mitigating nutrient-related risks, and promoting sustainable forestry practices.