Decision-Tree Application to Predict and Spatialize the Wood Productivity Probabilities of Eucalyptus Plantations

Abstract

Brazil is one of the world’s wood short-fiber producers, cultivating 7.5 million hectares of eucalypt trees. Foresters and resource managers often face difficulties in surveying reliable Eucalyptus productivity levels for the purpose of purchasing and prospecting lands. Spatial data science (DS) and machine learning (ML) provide powerful approaches to make the best use of the large datasets available today. Agriculture has made great use of these approaches, and in this paper, we explore how forestry can benefit as well. We hypothesized that both DS and ML techniques can be used to improve Eucalyptus productivity zoning based on multiple operational datasets of tree growth and environment. Based on more than 12,000 permanent forest inventory plots of commercial Eucalyptus plantations and the climate, soil, and altitude variables associated with them, a supervised ML approach was adjusted to model the forest plantation productivity. A multi-tuning of the decision-tree (DT) algorithm hyperparameters was prepared to yield 450 DT models, with a better one delivering an RMSE of 53.5 m³ ha⁻¹, split in 35 terminal nodes, here interpreted as Eucalyptus productivity zones. The DT model showed an optimum performance index of 0.83, a coefficient of determination of 0.91, a root mean squared error of 12.3 m³ ha⁻¹, and a mean absolute percentage error only of 3.1% in predicting the testing dataset throughout the study area. The DT rule set was interpreted in a user-friendly table and was prepared to classify any location within the study area in each one of the 35 productivity zones based on the required environment variables of the DT algorithm. The high quality of the model obtained made it possible to spatialize the DT rules, providing a reliable cartographic visualization of the probability levels of true Eucalyptus productivity for a huge region of forest-based industries in Brazil. These data-science techniques also provided a yield gap analysis using a very down-to-earth approach. We estimated a yield gap by an amount of 4.2 × 10⁷ m³, representing a few more than 113,000 ha, or 15% of the current forest base. This is the amount of avoided area expansion to accumulate the same wood stock in case the productivity is raised to the attainable level in each zone. This present study provided deeper analysis and reproducible tools to manage forest assets sustainably.

1. Introduction

Growing global demand for food and fiber [1], which added to the impacts of climate change [2], will require more efficient and resilient actions to consolidate the bioeconomy as a global economic system [3]. At the end of the day, the forest scientists and resource managers will provide more wood volume on less land.

Brazil is one of the world’s short-fiber wood baskets, home to some of the most productive and important multipurpose Eucalyptus plantations on the globe [4]. The country holds 7.5 million hectares of Eucalyptus plantations managed for the production of charcoal, panels, laminate flooring, paper, and mainly cellulose. The sector industries produced 21 million tons of pulp in 2020, seeing a strong increase of 90% in pulp exportation in the last 10 years worldwide (https://iba.org/eng, (accessed on 10 January 2022).

Historically, intensive research and development into forest growth in Southeastern Brazil outstandingly resulted in a fourfold increase in Eucalyptus productivity in the country from 1970 to 2015 [5]. However, most recently, the annual average Eucalyptus productivity growth rates dropped to zero when new forest expansions migrated to more climatically challenging regions of the country, which were still poorly understood ecophysiologically. New pests and issues related to variability and climate change have also exerted pressure to flatten the curve and fading productivity of Eucalyptus plantations in Brazil [6,7]. The impacts on Eucalyptus survival and productivity were well characterized during the dramatically long major droughts that struck Central–Southern Brazil between 2014 and 2015, between 2015 and 2016 in the Northern region, and between 2019 and 2021 again in Central–Southern Brazil, causing losses of hundreds of thousand hectares and a couple of billion dollars in forest assets [8,9,10].

Forest management mistakes in tropical regions generally result in larger yield gaps than in subtropical regions [11]. In-depth knowledge of production environments, their productivity levels, and the nature of their abiotic and biotic limitations are crucial points for the success of forest plantations in tropical regions. Most foresters know that the diverse physiographic and edaphoclimatic gradients found in tropical regions yield equally in diverse and complex production environments. Nonetheless, they often face difficulties in surveying their land’s reliable levels of productivity, basing their strategies on very site-specific models derived from purely empirical knowledge. Empirical models have their functionality in very specific situations where the rule set was created, as, for example, shown by Gava and Gonçalves [12]. Therefore, caution must be taken to avoid misunderstanding in any attempt of upscaling an empirical model.

The strategies used on productivity zoning should span multiple spatial scales and require a sound mechanistic understanding of the interactions between tree resource use dynamics over space and time. In Eucalyptus plantations, the use of ecophysiological and statistical models that incorporate climate variables such as rainfall or soil water deficit, adjusted with appropriate regionwide trials, are great examples to predict forest yield levels in target regions [13,14]. These studies have provided satisfactory results, but once based on experimental data generally scarcer, such approaches cannot perform distribution and probability modeling of the variable to be predicted.

Currently, most forest-based companies hold huge datasets of forest inventory and environmental monitoring with meteorological series and increasingly accurate soil and relief maps. However, in most cases, those tree growth and environment datasets are collected and used for purely operational purposes, despite foresters knowing such data remain latent for deeper applications to generate insights and other more powerful predictions. Hence, given the massive amount of data, available computational power, and balanced empirical and scientific knowledge, an opportunity arises to use spatial data-science (DS) and machine-learning (ML) concepts to improve the Eucalyptus productivity zoning.

DS is an interdisciplinary field that uses data collection, data curation, extraction, cleaning, analysis, model building, result visualization, and communication to abstract information and insights from large datasets [15]. ML refers to a set of available techniques in the toolbox of data scientists to produce predictions. ML is grouped into three categories, which include supervised, unsupervised, and reinforcement learning. Decision tree (DT) is a popular supervised classification that offers support in decision making, and it defines the logical structure, associated uncertainties, and constructive results of the decisions [15]. DT is a multifaceted machine-learning algorithm known as a glass-box model. After the DT model has determined the patterns in the training dataset, even readers with little or no machine-learning experience, after a brief explanation, might exactly and intuitively see what decisions will be made for a dataset that we want to predict. When it comes to ML, it is crucial to carefully weigh the trade off between estimation power and model interpretation. The DT algorithm excels at providing simple and straightforward results and interpretation models, showcased via elegant decision trees [16]. However, other ML methods that have a strong predictive ability often face challenges in terms of interpretation, reproducibility, and sharing capabilities [17,18,19].

In the tide of the digital world we live in today, there are many IT technicians working in large forest companies supporting digitizing processes and building big data. Nevertheless, they often take on roles beyond their professional skills, generally with a background in the basic sciences, and they rarely wet their feet in the field. Machine-learning algorithms used in forest ecology and management are prone to misinterpretation, which requires application with great caution. To avoid these interpretation problems, data scientists and foresters should go through the careful discovery phase. This includes the finding of relevant data, their formats from numerous sources to accomplish basic necessities, priorities, and objectives of the task, and also the major ecophysiological key cause–effect relationships involved in the studied case.

In a survey of machine-learning techniques for climate-based crop yield prediction, Elavarasan et al. [20] compared some existing supervised and unsupervised machine-learning models, and they concluded DT approach was faster and more precise, providing easily interpretable rules and assisting in identifying the significant fields. Pant et al. [21] made a comparison appraisal with a couple of machine-learning models to make maize, potatoes, rice, and wheat yield predictions and found that the DT regression resulted in the best accuracy for crop predictions compared to other algorithms, which linked to its reproducibility making the DT algorithm very powerful for solving problems in many applications on cash crops. Furthermore, the selection of predictor variables is another key point of a successful machine-learning implementation. Variables that are both strong predictors to target variables and simple to obtain or understand should be used. For example, for Eucalyptus productivity zoning, predictor variables such as temperature, rainfall, soil type, and altitude should be considered. Indeed, these major predictor variables fully converge with a literature review that showed these environmental features were the most applied in machine learning for crop yield prediction worldwide [22].

In the last 4 years, machine-learning usage for farm risk management has tripled [23]. Agriculture has made great use of these approaches [24], and in this paper, we explore how forestry can benefit as well. Numerous recent studies have employed various machine-learning techniques to assess the productivity of planted forests [25,26,27,28,29]. However, a notable limitation across these studies is the lack of comprehensive exploration regarding the tools utilized, deeper data visualizations, and the production of applicable and reproducible results and conclusions for end users. Furthermore, some of these studies have neglected multicollinearity in modeling analyses. This gap in the literature leaves a significant opportunity for further research and investigation to address these shortcomings and provide comprehensive and actionable insights for stakeholders in the field of planted forest productivity evaluation.

Therefore, our hypothesis suggests that the decision-tree approach can be utilized to enhance Eucalyptus productivity zoning by generating results that are transparent, reproducible, and based on multiple sets of operational tree growth data and environmental factors. The integration of decision-tree modeling allows for a more interpretable and explainable approach, shedding light on the key factors influencing Eucalyptus productivity. Additionally, this methodology facilitates the replication of results and fosters greater confidence in the zoning outcomes. Considering that, the objectives of this paper were fourfold: (1) to implement the data science and its components as discovery, data preparation, model building, model planning, model evaluation, result visualization, and communication on huge operational forest datasets; (2) to implement and assess a broad multi-tuning DT modeling approach to multiple datasets of tree growth and environment; (3) to perform zoning geographically explicit of the Eucalyptus productivity based on the adjusted DT model; (4) to determinate the Eucalyptus yield gaps using a machine-learning approach on a region-wide scale.

2. Materials and Methods

2.1. Study Area

The study area consists of approximately 11.4 million hectares located in Southeastern Brazil, covering nearly 46% of the State of São Paulo (Figure 1). According to the Spatial Database of Planted Trees [30], the State of São Paulo had in 2013/2014 approximately 823,619 hectares of Eucalyptus plantations belonging to companies in the pulp and paper, fiberboard, and wood panel industries, small growers, and participants in outgrower schemes and in state-owned forests. Of this total area, almost 95% of the Eucalyptus plantations (781,654 ha) were within the study area (Figure 1), i.e., the other 5% (41,965 ha) are well sparsely distributed in the 54% of the midwestern region of the state.

Covering different production environments and a large latitudinal gradient, the study area shows both tropical and subtropical features. Average annual temperature varies from 14 to 23 °C (on higher mountains, above 2000 m, average annual temperature can be as low as 14 °C), rainfall ranging 1200 to 2000 mm year⁻¹, and the frost days hits are frequent in the southern part are main climatological traits of study area [31,32].

The study area has complex relief, with altitude ranging from 400 to 1500 m above mean sea level and a relevant geological diversity with parent sedimentary (sandstone, claystone, and siltstone), metamorphic (phyllites, shales, gneiss) and igneous (diabase and basalt) materials, resulting in large range of deep and shallow tropical soils distributed across the study area on erratic way on part regions [33,34,35].

Such geographic distribution of Eucalyptus plantations shows that the study area is one of the most traditional forest regions in the country with almost 100 years of experience in Eucalyptus silviculture [6,36] because it is one of the most favorable in world for silvicultural success due to its high suitability for dozens of Eucalyptus species [37].

2.2. Forest Datasets

Approximately 12,000 permanent forest inventory plots of commercial Eucalyptus plantations of the Suzano SA company (www.suzano.com.br/en, (accessed on 15 June 2023) at harvest time between 2013 and 2016 (with planting year between 2008 and 2011) were available for this study. Circular permanent inventory plots had an area of 400 m² and were systematically distributed with a density of one plot per ten hectares. During field measurements, diameter at breast height (DBH, 1.3 m above the ground) of all trees, total height of the central ten trees, and the total height of the four largest trees in terms of DBH (defining dominant trees) were surveyed in the inventory plot. The mean height of these dominant trees defined the dominant height of the plot. Dominant height, basal area, and age on the inventory date were then used in local volume equations to estimate the total wood volume with bark for each plot. A total of 9953 forest inventory plots composed of 517,000 live trees (by a total of 528,000 trees) were selected for this study after tree data consistency, which included filtering them to have plots with (i) live trees between 800 and 1900 trees ha⁻¹, (ii) age between 5.5 and 6.5 years (age class center of 6 years old) at harvest time, (iii) standard size plots, (iv) only first rotation forest (not coppice), (v) free pests and diseases plots, (vi) no record of plots faced strong weed competition and fire scars, and (vii) double checking its geographic location. Finally, a visual dispersion analysis of the permanent forest inventory plots database was required to ensure reliable analysis and modeling, with the final dataset clean-up dropping any strange relationship between basal area and total wood volume (Figure S1).

The forest inventory database was randomly divided into two datasets using a geographically weighted approach. The first, referred to as the training dataset, with 80% of all plots, was used in DT modeling exercise in order to build the productivity zones for planted Eucalyptus forests. Another group, called testing dataset, with 20% remaining plots (unseen data by the model), was used to perform the validation of the productivity zones output of the DT. The forest inventory dataset consisted of 70 pure or hybrid Eucalyptus clones of the species E. grandis, E. urophylla, E. smithii, and E. dunnii.

2.3. Environmental Datasets

For each forest inventory plot, information such as geographic location, altitude obtained in the field, average terrain slope, planting date, measurement date, genetic material, planting spacing, wood volume, rotation, soil order, and soil texture was provided.

Soil information of each plot is from company-owned detailed soil maps, surveyed at scale of 1:50,000 and elaborated by professional pedologists specialized in tropical soils. Those company soil maps cover only Eucalyptus stands. Well-developed and deep soils accounted for 57.3% of the plots, as the order Latossolos (L, Oxisols) with 43.1%, Argissolos (A, Ultisols) with 13.9%, and Nitossolos (T, Oxisols) with 0.4% remaining. On the other hand, less developed and/or shallow soils totaled for 42.5% of the plots, as Neossolos (N, Entisols) with 29.6%, and Cambissolos (C, Inceptisols) with 12.9% of them. Finally, 0.2% of the plots were classified as hydromorphic soil order Gleissolos (G, Entisols). Soil order names in Portuguese are according to the Brazilian Soil Classification System [38]. The symbol letter for soil order used on DT modeling and an approximate identification based on the USDA soil taxonomy are in parentheses [39].

Climatological data such as average annual temperature, average annual rainfall, average annual potential evapotranspiration, and aridity index were associated with each forest inventory plot. Meteorological data from a network of 48 automatic weather stations of the Brazilian National Institute of Meteorology (INMET, https://mapas.inmet.gov.br, (accessed on 15 May 2021)), Integrated Center of Agrometeorological Information of Campinas Institute of Agronomy (CIIAGRO/IAC, www.ciiagro.sp.gov.br/rede.html, (accessed on 15 May 2021)), and from the forest company, covering the entire study area, were used to develop the climatological database (Figure 2). We considered a climatic period of the last 26 years (1995 to 2020) to build maps of the monthly climatological normals for both air temperature and rainfall. The weather station datasets presented some data gaps in rainfall and air temperature records. Those missing data were filled with monthly gridded data published by Xavier et al. [40] at in spatial resolution of 0.25° and from NASA/POWER—NASA Langley Research Center POWER Project, which has recently been released at 0.5° × 0.5° grid cell [41].

Average annual rainfall dataset was analyzed according to geostatistical procedures proposed by Alvares et al. [31]. The omnidirectional experimental semivariogram for the annual rainfall data was prepared, whereas only 10 lag classes of equal intervals in a geometric field of 50% range dataset (latitude and longitude) was prepared. The experimental semivariogram was tested by adjusting them to the theoretical spherical, exponential, Gaussian, and linear models since they usually covered the general dispersion of environmental spatial events [42]. The best fit was based on the smallest reduced sum of squares and on the highest determination coefficient (R²). Using these structural parameters of the theoretical semivariogram extracted from experimental semivariogram, average annual rainfall map was composed with a geographic information system. A punctual ordinary kriging estimator was used for geostatistical interpolation. With these procedures, the annual rainfall map had a spatial resolution of 100 m.

The air temperature dataset was used to establish the relationship between the dependent (average air temperature) and independent (altitude, latitude, longitude, and their combinations) variables based on a general multivariate nonlinear regression model proposed by Alvares et al. [43] as follows (Equation (1)):

$${\mathrm{T}}_{\mathrm{i}}=\mathrm{a}0+\mathrm{a}1\mathsf{\phi }+\mathrm{a}2\mathsf{\lambda }+\mathrm{a}3\mathrm{h}+\mathrm{a}4\mathsf{\phi }\mathsf{\lambda }+\mathrm{a}5\mathsf{\phi }\mathrm{h}+\mathrm{a}6\mathsf{\lambda }\mathrm{h}+\mathrm{a}7{\mathsf{\phi }}^2+\mathrm{a}8{\mathsf{\lambda }}^2+\mathrm{a}9{\mathrm{h}}^2$$

(1)

where T_i is the average temperature, monthly (i = 1, 2, …, 12) or annual (i = 13); φ is latitude, in decimal degrees (negative in the southern hemisphere); λ is longitude, in decimal degrees (negative values for west of Greenwich meridian); h is altitude in meters; and a0 to a9 are the coefficients of the multivariate regression equation.

Using map algebra techniques [42,44], the results from the multivariate regression models for average monthly air temperature were converted into maps using the software ArcGIS [45], processing the independent variables as raster layer. The altitude layer, in meters, was obtained from the digital elevation model (DEM) provided by the Shuttle Radar Topography Mission (SRTM) [46] in its fourth version [47]. Subsequently, ArcGIS was used to build the DEM for stud area at 100 m resolution (Figure 2). Latitude and longitude layers were obtained in decimal degrees using the central coordinates of each pixel corresponding to the DEM. Using geoprocessing techniques [48,49,50], monthly air temperature models were programmed and run in ArcGIS.

After building the monthly temperature maps, we proceeded with the potential evapotranspiration calculations. Average annual potential evapotranspiration was estimated via the method of Thornthwaite [51]. Such method is very practical and easy to reproduce as it only requires the average monthly temperature observed as an input. The observed climatological normals are used to calculate Thornthwaite’s monthly and annual heat indices (Equations (3) and (4)). Latitude is used to calculate the photoperiod over the course of the year. Lastly, we must consider the number of days to adjust the potential evapotranspiration of each month (Equation (2)). After calculations, the monthly potential evapotranspiration maps were summed to obtain the annual potential evapotranspiration map.

$${\mathrm{PET}}_{\mathrm{i}}=16\left(\frac{\mathrm{N}}{30}\right)\left(\frac{\mathrm{P}}{12}\right){\left(\frac{10{\mathrm{T}}_{\mathrm{i}}}{\mathrm{I}}\right)}^{\mathrm{a}}$$

(2)

$$\mathrm{I}={\displaystyle \sum _{\mathrm{i}=1}^{12}}{\left(0.2{\mathrm{T}}_{\mathrm{i}}\right)}^2$$

(3)

$$\mathrm{a}=\left(6.75\times {10}^{-7}\right){\mathrm{I}}^3-\left(7.71\times {10}^{-5}\right){\mathrm{I}}^2+\left(1.7912\times {10}^{-2}\right)\mathrm{I}+0.49239$$

(4)

where PET_i is the estimated potential evapotranspiration of month i (mm); N is the number of days in month i; P is the average daily photoperiod of month i (hours); T_i is the average daily temperature of month i (°C); and a and I are the Thornthwaite’s heat indices, which depend on the average daily temperature of month i.

Based on annual potential evapotranspiration and annual rainfall maps, the aridity index was calculated (Equation (5)). Aridity index is an indicator for the degree of site dryness [52] and is used as the surrogate of water stress on Eucalyptus plantations [8,53]. After the spatializations of annual air temperature, annual rainfall, annual potential evapotranspiration, and aridity index, extractions of these data to populate the forest inventory plot datasets were performed.

$$\mathrm{Aridity}\mathrm{index}=\left(\frac{\mathrm{PET}}{\mathrm{R}}\right)$$

(5)

where PET is the estimated annual potential evapotranspiration (mm), and R is average annual rainfall (mm).

2.4. Decision-Tree Modeling

In order to guarantee the feasibility and reproducibility of the productivity zones for planted Eucalyptus forests by the foresters, small rural producers, forest beekeepers, participants interested in forest outgrower schemes, extension technicians, and researchers, a range of easily obtainable variables of climate, soil, and terrain were used to model and train the decision tree. One object of this study is to create a method that the end users can use in a straightforward way and apply the results of the productivity zoning model of planted Eucalyptus forests, thus maintaining the practicability of the system.

After the consistency and processing of tree growth and environmental data, each dataset was joined into a data frame composed of 9953 observations and 11 variables, which are as follows: forest inventory plot ID, latitude, longitude, training or testing dataset, genetic material, tree age, total wood volume with bark as the response variable; and the predictor variables as terrain altitude, soil order, soil texture, and average annual aridity index. However, some adaptations for DT algorithm were made as the smooth discretization of aridity index on the 0.015 scale and of ground altitude on the 20 m scale in order to have adequate aesthetics in the DT and easy interpretation of the thresholds of these predictive variables.

Based on the objectives of this study, the rpart package (Recursive Partitioning and Regression Trees [54]) and its routines for R software [55] were used to predict the productivity zones for planted Eucalyptus forests. The rpart’s main outcome is a recursive partitioning model displayed as a decision tree.

Before making our DT grow, we prepared a tuning using rpart hyperparameters that control aspects of the model fit. Basically, there are hyperparameters that control the growth and pruning of the DT as well as its complexity. Minbucket parameter is the minimum number of observations in any leaf (terminal) node. Minsplit parameter is the minimum number of observations that must exist in a node for a split to be attempted. Cost-complexity factor (cp parameter) is the metric that stops splits that are not deemed important enough to tree grow up and explain data.

A broad multi-tuning of the DT hyperparameters was planned. We wrote functions to control the hyperparameters during tuning approach that ranged as follows: minbucket parameter from 79 to 1000, by 100; minsplit parameter from 158 to 1000, by 100; and cp parameter tested the complexity measures of 0.01, 0.001, 0.0005, 0.0001, and 0.00001. Once DT modeling requires caution in terms of growth, pruning and trim, and tree complexity, we wrote functions to display the accuracy of each one DT outputted using the multi-tuning approach. The same function set also produced a visualization of ranging upper hyperparameters of decision-tree tuning process, supporting the better-chosen DT setup. Before departure of hyperparameter tuning modeling, we already assume a premise that minbucket parameters start at 79, i.e., the number of observations in the leaf node must be at least 1% of training dataset, ensuring an acceptable minimum forest and environmental representations of the data entries. Hence, we can notice minsplit parameter was set to start with 158 observations.

Once the DT model was defined and its best-fitting hyperparameters were set up, the next step was to assess the performance of the model making predictions of independent data (the 20% unseen data). For that, the testing set was used to validate the classification and regression model via the following relative errors indices: mean absolute error (MAE), root mean squared error (RMSE), and mean absolute percentage error (MAPE), respectively, as shown in Equations (6)–(8). The performance of the DT model was also tested using the Performance index “Pi” (Alvares et al. 2013b [43]), which is the product of the coefficient of correlation “r” (Equation (9)), which defines models’ precision, and the refined agreement index “dr” (Equation (10)) [56], which defines the models’ accuracy. The criteria for interpreting “Pi” (Equation (11)) include the following: Pi ≥ 0.75 as optimum performance; 0.6 ≤ Pi < 0.75 as very good performance; 0.45 ≤ Pi < 0.6 as good performance; 0.3 ≤ Pi < 0.45 as tolerable performance; 0.15 ≤ Pi < 0.3 as poor performance; 0 ≤ Pi < 0.15 as bad performance; Pi < 0 as very bad performance [43].

$$\mathrm{MAE}=\frac{1}{\mathrm{N}}\sum \left|{\mathrm{y}}_{\mathrm{i}}-{\mathrm{x}}_{\mathrm{i}}\right|$$

(6)

$$\mathrm{RMSE}=\sqrt{\frac{1}{\mathrm{N}}\sum {\left({\mathrm{y}}_{\mathrm{i}}-{\mathrm{x}}_{\mathrm{i}}\right)}^2}$$

(7)

$$\mathrm{MAPE}=\frac{100}{\mathrm{N}}\sum \left|\frac{{\mathrm{y}}_{\mathrm{i}}-{\mathrm{x}}_{\mathrm{i}}}{{\mathrm{y}}_{\mathrm{i}}}\right|$$

(8)

$$\mathrm{r}=\frac{\left({\mathrm{x}}_{\mathrm{i}}-\overline{\mathrm{x}}\right)\left({\mathrm{y}}_{\mathrm{i}}-\overline{\mathrm{y}}\right)}{\sqrt{\left[{\left({\mathrm{x}}_{\mathrm{i}}-\overline{\mathrm{x}}\right)}^2\right]\left[{\left({\mathrm{y}}_{\mathrm{i}}-\overline{\mathrm{y}}\right)}^2\right]}}$$

(9)

$$\mathrm{dr}=\left\{\begin{array}{c}1-\frac{\sum \left|{\mathrm{y}}_{\mathrm{i}}-{\mathrm{x}}_{\mathrm{i}}\right|}{\mathrm{c}\sum \left|{\mathrm{y}}_{\mathrm{i}}-\overline{\mathrm{y}}\right|}, \mathrm{when} \sum \left|{\mathrm{y}}_{\mathrm{i}}-{\mathrm{x}}_{\mathrm{i}}\right|\le \mathrm{c} \sum \left|{\mathrm{y}}_{\mathrm{i}}-\overline{\mathrm{y}}\right|\\ \frac{\mathrm{c}\sum \left|{\mathrm{y}}_{\mathrm{i}}-\overline{\mathrm{y}}\right|}{\sum \left|{\mathrm{y}}_{\mathrm{i}}-{\mathrm{x}}_{\mathrm{i}}\right|}-1, \mathrm{when} \sum \left|{\mathrm{y}}_{\mathrm{i}}-{\mathrm{x}}_{\mathrm{i}}\right|>\mathrm{c} \sum \left|{\mathrm{y}}_{\mathrm{i}}-\overline{\mathrm{y}}\right|\end{array}\right.$$

(10)

$${\mathrm{P}}_{\mathrm{i}}=\mathrm{r}\times \mathrm{dr}$$

(11)

where x_i is the i-th observed value, y_i is the i-th predicted value, N is the number of samples considered, $\overline{\mathrm{x}}$ is mean of all observed values, $\overline{\mathrm{y}}$ is mean of all predicted values, r is the Pearson’s correlation coefficient, dr is refined agreement index, c is 2, and ${\mathrm{P}}_{\mathrm{i}}$ is performance index.

Finally, after obtaining the final DT model, we wrote R functions that extracted the environmental predictor variable rules for each leaf node (Eucalyptus productivity zones) and put them in organized and straightforward way in an interpretive and illustrated table to be used by the end users.

2.5. Eucalyptus Forest Productivity Zoning

In this section, we deal with the geographic visualization of the DT outputs. Each DT leaf node was interpreted as a productivity zone for planted Eucalyptus forests, i.e., a production region’s average of Eucalyptus plantations productivity, observed at training dataset, is explained via the same environmental predictor variables and its thresholds.

We already reported the spatialized environmental predictor variables such as altitude and aridity index used to populate the training and testing datasets as independent data, as obtained in Section 2.3. Finally, the semi-detailed soil map of the State of São Paulo published by Rossi [57] was used to complete the required environmental information to apply the Eucalyptus forest productivity zoning outputted via DT. In total, this soil map has 247 soil mapping units, and it presented 20 soil classes (order and suborder) in the study area according to Figure 3. In order to use the soil map on DT modeling, each soil order found at Rossi [57] in study area, as already described in Section 2.3 for the soil information of the forest inventory plots, was coded as follows: A represents Argissolos; L represents Latossolos; N represents Neossolos; C represents Cambissolos; T represents Nitossolos; G represents Gleissolos; M represents Chernossolos. Likewise, the subsurface horizon soil texture identified by Rossi [57] received the code “a” for clayey texture (clay content > 35%), “m” for medium texture (15% < clay content < 35%), and “r” for sandy texture (15% < clay content). Soil order and soil texture vector maps were converted to raster in the spatial resolution of 100 m. Lastly, a prediction was made by transforming and spatializing each leaf node on productivity zones Eucalyptus for the whole study area based on DT rules and the predictor environmental variables.

3. Results and Discussion

3.1. Climate Modeling

The average annual rainfall dataset from 48 weather stations was analyzed according to the concepts of the regionalized variables theory [58]. Among the theoretical models tested, the spherical model best described the annual rainfall experimental semivariogram (nugget effect = 13,600; partial sill = 102,980; and range of spatial dependence = 7.3 degrees). The estimated residues were low and had a high coefficient of determination (R² = 0.97). The average annual rainfall of the study area (Figure 4) has spatial variability similar to the current and official climatological normal map presented by Diniz et al. [59]. Average annual rainfall ranged from 1250 to 2050 mm, very similar to what was found in other studies [31,59].

All monthly and annual average air temperature multivariate regression models were significant (p < 0.001), and the independent variables obtained showed a significance of 5% (Table S1). Different from other studies [43,60], the combined predictor variables that used the square of altitude (h2) were significant at 5%, thus contributing to predicting the average air temperatures, both in monthly and annual time scales (Figure S2).

The average annual PET range (540–1220 mm) in the study area (Figure 4) was quite similar to that presented by Dias et al. [61]. Thereby, the ratio between annual PET and annual rainfall was applied to obtain the annual aridity index (Figure 5). The higher index identifies more arid environments. The average annual aridity index map showed the uplands and mountain regions, with altitudes above 900 m, mainly on the Serra do Mar range and the Mantiqueira region, where atmospheric water demands are very low, with an aridity index of less than 0.5. These upland regions have a water surplus of approximately 50% of the amount of annual rainfall, indicating high suitability for Eucalyptus growth [37]. Conversely, the aridity index exceeded 0.76 in most areas in the central–west region of the study area, denoting water stress throughout the year. Paraiba Valley (southeast of the study area) also has a high aridity index which certainly penalizes Eucalyptus productivity, as evidenced by many recent studies [5,8,14,62]. Our aridity index map ranges likewise to the one presented by Hubbard et al. [53]. These authors showed at tropical environments, the Eucalyptus trees grow according to soil water availability and that both high- and low-productivity Eucalyptus clones investigated transpired on average 60, 46, and 34% annual precipitation on dry (annual aridity index of 2.2), mesic (1.0), and wet (0.6) sites, respectively.

3.2. Decision-Tree Modeling

The hyperparameter tuning approach yielded 450 DT models, with RMSE ranging from 60.1 m³ ha⁻¹ for the worst model to 53.3 m³ ha⁻¹ for the best model. A scatter plot with all tuning modeling outputs allowed us to visualize a global exponential decay pattern between RMSE and the number of leaf nodes in each DT model, showing also that increasing the number of leaf nodes linearly reduces the RMSE for each one of the five cost-complexity factors applied (Figure S3). While Figure S3 provided an overview of tuning modeling, Figure S4 gave a zoom-in visualization of how all RMSE and the number of leaf nodes outputs are distributed between minsplit and minbucket hyperparameters for the cost-complexity factor = 0.0005. Among selecting the best DT by math tools like the minimum curve radius method, there is an important trade off with the empirical foresters’ knowledge. Therefore, the best and final DT was chosen based on those accuracy plot panels of the machine-learning algorithm modeling and on the practicality of the results and tools that were objectives of this study. Consequently, as expected and described previously, our premise considered that all leaf nodes must have at least 1% of the training dataset revealed via combining empirical and computational approaches because the better DT setup selected has the following hyperparameters: minbucket = 79, minsplit = 158, and cp = 0.0005, leaf node = 35, and RMSE = 53.5 m³ ha⁻¹ (Figures S3 and S4). Rossit et al. [63] carried out a big data study on forestry harvest productivity and also determined that the smallest leaf node must contain at least 1% of the dataset on DT modeling since obtaining decision rules for sets representing less than 1% is meaningless for the problem addressed. The DT selected in our study has only half of the leaf nodes of the better statistical model, but an RMSE quite similar was only 0.39% lower, guaranteeing the intelligibility and simplicity of the final machine-learning algorithm while avoiding jeopardizing both feasibility and reproducibility of the forest productivity zones so desired for the final users. Although the hyperparameter tuning process is a recognized procedure and fundamental step in machine-learning modeling, van Klompenburg et al. [22] performed a broad systematic literature review on crop yield prediction using machine learning, evidencing that only in one study had authors used hyperparameter tuning for sugarcane yield predictions.

Based on the hyperparameter setup of the selected DT, the recursive partitioning model delivered the Eucalyptus forest zone productivity classification displayed as an elegant decision-tree scheme (Figure 6). The rpart algorithm produced a DT with a maximum of nine branch levels until leaf (terminal) nodes were identified at the bottom of the tree. The rpart algorithm departed from the root node with 100% of the data (7975 observations from the training dataset) and with an average Eucalyptus total wood volume of 308 m³ ha⁻¹. The first splitting criterion was “Is the aridity index equal to or greater than 0.74?”; when the stated condition is true, we must move to the left branch. The left node with a higher aridity index (≥0.74) represents 24% of the observations, averaging 263 m³ ha⁻¹, while the right node with a lower aridity index accounts for 76% of the remnants of data, averaging 323 m³ ha⁻¹ (Figure 6). After this brief explanation, even readers with little experience in machine learning can read the DT following the Boolean choice questions until reaching the leaf node. As noted previously, the aridity index had significant importance for the recursive partitioning of our tree growth data. Soil water availability is quite well recognized for controlling the Eucalyptus wood productivity in tropical and subtropical conditions sites worldwide [64,65,66]. Soil order arose on the second splitting criterion of tree slicing the wettest zones on the following two subsets: the less productive zones with soil orders A, M, or N; the more productive zones with soil orders C, G, L, or T. Although the interpretation of the relevance of the order of soil in the DT is not so evident, we note that deeper soils are much more important for the survival and productivity of Eucalyptus in regions with higher water deficit, as demonstrated by Pinheiro et al. [67] and Silva et al. [68] on similar soils in the study area. Altitude strata were present on most of the leaf nodes, positively correlating with wood productivity. The highest zones are generally cold, which, in turn, has a lower PET rate resulting in a milder aridity index. Lastly, the soil texture was less relevant but significant to split early nodes of the drier zones and for end branches of the wetter and productive zones (Figure 6). Quite generally, sandy soils were related to low productive zones, while medium texture and clayey soils have linked to more productive zones. Many researchers found a noteworthy relationship between Eucalyptus productivity and several attributes of soil, mainly clay content, which was directly related to the amount of soil’s plant-available water content, and these findings are well documented in Gonçalves et al. [6].

The validation of the DT model showed that our algorithm had an optimum performance (Pi = 0.83) in predicting the existing testing dataset throughout the study area (top plot of Figure 7). Moreover, statistical goodness-of-fit validation tests revealed a high coefficient of determination (R² = 0.91), an RMSE of 12.3 m³ ha⁻¹, and an MAPE of only 3.1%. The validation procedure comprised making a comparison between the average Eucalyptus total wood volume of each leaf node against the Eucalyptus total wood volume predicted via the regression algorithm of the unseen data. It is necessary to highlight that the validation was performed between total wood volume averages, which basically represent a 50% probability of reaching this percentile for both DT regression and predicted results. Therefore, our initial concern for performing an adequate data division between training and testing datasets was important before starting the machine-learning modeling but, above all, having defined a minimum threshold of observations in each leaf node, pruning the tree growth’s every decision where a node had up to 1% of the training dataset, thus ensuring an appropriate total wood volume distribution for any climate, soil, and altitude conditions detected in the study area. Eucalyptus productivity zoning studies have shown that validation is a major part of modeling credibility, whether it is an empirical or a process-based model, from extremely localized applications, as presented by Almeida et al. [69]; at the regional level by Lemos et al. [70] and Attia et al. [71]; and at continental level, as shown by Caldeira et al. [72] and Elli et al. [73]. Although these studies performed a simple validation of their methods, none of them presented uncertainties in the estimates. Furthermore, even in applications validated under very localized or very experimental conditions, care must be taken with overfitting problems and site-specificity and genetic-specificity biases when planning any extrapolation [74,75,76].

Figure 7. Validation of decision-tree model via comparison of the simulated average Eucalyptus total wood volume at 6 years old of each leaf node against testing dataset independent data (top plot). See Section 2.4 for statistical analysis application. Each leaf node had equivalent data proportion between training and testing datasets, showing random division was complete (bottom plot). For example, for the leaf node 267 (Zone 29), 6.9% of the training set is 549 forest inventory plots, and 6.8% of the testing set is 132 forest inventory plots (see the gray arrow on the bottom plot). The accumulated percentages add up to 100% for both training and testing datasets.

Figure 7. Validation of decision-tree model via comparison of the simulated average Eucalyptus total wood volume at 6 years old of each leaf node against testing dataset independent data (top plot). See Section 2.4 for statistical analysis application. Each leaf node had equivalent data proportion between training and testing datasets, showing random division was complete (bottom plot). For example, for the leaf node 267 (Zone 29), 6.9% of the training set is 549 forest inventory plots, and 6.8% of the testing set is 132 forest inventory plots (see the gray arrow on the bottom plot). The accumulated percentages add up to 100% for both training and testing datasets.

3.3. Innovations in Decision-Tree Use

Our study also allows for exploring the production capabilities of each productivity zone once the data and approaches support these breakthrough deliveries. By retro-applying the DT rules to the training data and classifying it at each productivity zone, we were able to prospect the distribution and productivity probabilities in each of them. The observed productivities of each zone were adjusted to the normal distribution and the percentiles 15th, 50th, and 85th were identified in the sophisticated ridgeline plot (Figure 8). We noted that the productivity outcomes from the DT identifying the leaf node is an average value generated via the recursive partitioning and regression trees method and, therefore, slightly differs from the median values (50th percentile) found for each productivity zone in Table 1. While Zone 01 has 85% (15th percentile) and 15% (85th percentile) probabilities of productivity exceeding 308 and 462 m³ ha⁻¹, respectively, for Zone 35, the expected productivities in these same probability levels are 187 and 268 m³ ha⁻¹ (Figure 8 and Table 1). Such results and their presentation format represent a powerful tool that empowers readers to know the true productivity levels of the forest zone of interest, largely reducing uncertainty and leaning on ground-true forest data in decision making.

Figure 8. Distribution of Eucalyptus total wood volume at 6 years old from the training dataset for each leaf node (Eucalyptus forest zone). Bars within density plots identify the percentiles 15th, 50th, and 85th. The color ramp corresponds with the previous figures, but the ridgeline plot’s aesthetic required an opacity of 0.8. The percentiles values for each Eucalyptus forest zone are in Table 1.

Figure 8. Distribution of Eucalyptus total wood volume at 6 years old from the training dataset for each leaf node (Eucalyptus forest zone). Bars within density plots identify the percentiles 15th, 50th, and 85th. The color ramp corresponds with the previous figures, but the ridgeline plot’s aesthetic required an opacity of 0.8. The percentiles values for each Eucalyptus forest zone are in Table 1.

There are many other important and never-before-answered questions in the literature that we can still respond to using the productivity distributions in Figure 8. For example, a stakeholder prospecting for lands to lease a growing eucalypt forest with a baseline production planned to reach 300 m³ ha⁻¹ in 6 years (50 m3 ha⁻¹ y⁻¹) has the following question: What is the probability that the threshold will be reached at each zone? The answers with the probabilities are outlined in Figure 9, and this is one of the significant leap forwards provided by our study. Zones 01, 02, 03, and 04 have more than an 80% probability of productivity reaching the threshold of 300 m³ ha⁻¹, while for Zones 32, 33, 34, and 35, the probability dropped to less than 15% (Figure 9). These distinguished data visualizations are a major contribution to both forest science and the sector when providing unedited insights and consistent results and reducing the uncertainty of Eucalyptus productivity in the study area. Once most scientists use many tools in creating ideas and communicating new insights to peers and the public, the communication tools are fundamental and strategic because, without clear, understandable communication, brilliant insights remain weak and ineffective [77].

Figure 9. Probability (%) of total wood volume to reach 300 m³ ha⁻¹ at 6 years old for each forest zone. The probabilities were determined based on distribution of Eucalyptus total wood volume from the training dataset for each leaf node (Figure 8).

Figure 9. Probability (%) of total wood volume to reach 300 m³ ha⁻¹ at 6 years old for each forest zone. The probabilities were determined based on distribution of Eucalyptus total wood volume from the training dataset for each leaf node (Figure 8).

Table 1. Environmental rule interpretation and reading of the predictor variables forming each leaf node (Eucalyptus forest zone) outcome using the decision-tree model. Leaf nodes were ranked in decreasing order (from highest to lowest productivity), named Zone 01 to 35. The score percentiles 15th, 50th, and 85th are also provided for each leaf node as a proxy for Eucalyptus production probabilities of each forest zone.

Forest Zone	Leaf Node (m3 ha−1)		Aridity (PET/R)					Altitude (m)					Soil Order			Soil Texture		15th (m3 ha−1)	50th (m3 ha−1)	85th (m3 ha−1)
Zone 01	386	when	is	0.59	to	0.64	and	>=			900	and	is	C, G, L, or T	and	is	a	308	395	462
Zone 02	369	when	is	0.62	to	0.65	and	is	720	to	760	and	is	C, G, L, or T				307	368	419
Zone 03	358	when	is	0.59	to	0.64	and	is	800	to	900	and	is	L				300	364	412
Zone 04	354	when	is	0.59	to	0.64	and	>=			900	and	is	C, G, L, or T	and	is	m	291	351	418
Zone 05	353	when	<	0.46			and	>=			980	and	is	C, G, L, or T				275	362	420
Zone 06	349	when	is	0.61	to	0.73	and	<	720			and	is	A, M, or N	and	is	m	292	354	410
Zone 07	346	when	is	0.62	to	0.70	and	is	780	to	800	and	is	C, G, L, or T				302	344	397
Zone 08	345	when	is	0.67	to	0.74	and	<	720			and	is	C, G, L, or T	and	is	a	290	348	397
Zone 09	341	when	is	0.65	to	0.74	and	is	720	to	760	and	is	C, G, L, or T				272	342	396
Zone 10	338	when	is	0.67	to	0.74	and	<	620			and	is	C, G, L, or T	and	is	m or r	288	338	389
Zone 11	337	when	is	0.64	to	0.74	and	is	800	to	900	and	is	C, G, L, or T				270	346	393
Zone 12	334	when	is	0.46	to	0.59	and	>=			980	and	is	C, G, L, or T				261	349	405
Zone 13	333	when	is	0.59	to	0.64	and	is	800	to	900	and	is	C				266	335	402
Zone 14	331	when	<	0.47			and	<	980			and	is	C, G, L, or T				256	336	403
Zone 15	328	when	<	0.74			and	is	720	to	800	and	is	A, M, or N	and	is	a or m	246	324	404
Zone 16	327	when	is	0.61	to	0.73	and	<	720			and	is	A, M, or N	and	is	a or r	277	324	377
Zone 17	326	when	is	0.67	to	0.74	and	is	620	to	720	and	is	C, G, L, or T	and	is	m or r	281	324	370
Zone 18	324	when	is	0.62	to	0.70	and	is	760	to	780	and	is	C, G, L, or T				264	322	373
Zone 19	315	when	is	0.62	to	0.67	and	<	720			and	is	C, G, L, or T				229	332	420
Zone 20	308	when	is	0.64	to	0.74	and	>=			900	and	is	C, G, L, or T				260	308	351
Zone 21	299	when	is	0.50	to	0.59	and	is	800	to	980	and	is	C, G, L, or T				245	298	359
Zone 22	297	when	is	0.73	to	0.74	and	<	720			and	is	A, M, or N				227	285	341
Zone 23	294	when	is	0.74	to	0.77									and	is	m	254	298	336
Zone 24	290	when	<	0.74			and	>=			800	and	is	A, M, or N	and	is	a or m	235	295	358
Zone 25	277	when	is	0.50	to	0.59	and	<	800			and	is	C, G, L, or T				215	286	340
Zone 26	276	when	is	0.77	to	0.82									and	is	m	232	274	321
Zone 27	273	when	is	0.70	to	0.74	and	is	760	to	800	and	is	C, G, L, or T				238	277	312
Zone 28	269	when	<	0.74			and	>=			720	and	is	A, M, or N	and	is	r	204	272	326
Zone 29	267	when	is	0.74	to	0.82	and	<	640						and	is	a or r	230	264	307
Zone 30	265	when	<	0.61			and	<	720			and	is	A, M, or N				206	275	334
Zone 31	263	when	is	0.47	to	0.50	and	<	980			and	is	C, G, L, or T				200	264	326
Zone 32	245	when	>=			0.82	and	<	520									210	241	284
Zone 33	242	when	is	0.74	to	0.82	and	>=			640				and	is	a or r	199	244	288
Zone 34	240	when	is	0.59	to	0.62	and	<	800			and	is	C, G, L, or T				194	238	289
Zone 35	228	when	>=			0.82	and	>=			520							187	226	268

Table 1. Environmental rule interpretation and reading of the predictor variables forming each leaf node (Eucalyptus forest zone) outcome using the decision-tree model. Leaf nodes were ranked in decreasing order (from highest to lowest productivity), named Zone 01 to 35. The score percentiles 15th, 50th, and 85th are also provided for each leaf node as a proxy for Eucalyptus production probabilities of each forest zone.

Forest Zone	Leaf Node (m3 ha−1)		Aridity (PET/R)					Altitude (m)					Soil Order			Soil Texture		15th (m3 ha−1)	50th (m3 ha−1)	85th (m3 ha−1)
Zone 01	386	when	is	0.59	to	0.64	and	>=			900	and	is	C, G, L, or T	and	is	a	308	395	462
Zone 02	369	when	is	0.62	to	0.65	and	is	720	to	760	and	is	C, G, L, or T				307	368	419
Zone 03	358	when	is	0.59	to	0.64	and	is	800	to	900	and	is	L				300	364	412
Zone 04	354	when	is	0.59	to	0.64	and	>=			900	and	is	C, G, L, or T	and	is	m	291	351	418
Zone 05	353	when	<	0.46			and	>=			980	and	is	C, G, L, or T				275	362	420
Zone 06	349	when	is	0.61	to	0.73	and	<	720			and	is	A, M, or N	and	is	m	292	354	410
Zone 07	346	when	is	0.62	to	0.70	and	is	780	to	800	and	is	C, G, L, or T				302	344	397
Zone 08	345	when	is	0.67	to	0.74	and	<	720			and	is	C, G, L, or T	and	is	a	290	348	397
Zone 09	341	when	is	0.65	to	0.74	and	is	720	to	760	and	is	C, G, L, or T				272	342	396
Zone 10	338	when	is	0.67	to	0.74	and	<	620			and	is	C, G, L, or T	and	is	m or r	288	338	389
Zone 11	337	when	is	0.64	to	0.74	and	is	800	to	900	and	is	C, G, L, or T				270	346	393
Zone 12	334	when	is	0.46	to	0.59	and	>=			980	and	is	C, G, L, or T				261	349	405
Zone 13	333	when	is	0.59	to	0.64	and	is	800	to	900	and	is	C				266	335	402
Zone 14	331	when	<	0.47			and	<	980			and	is	C, G, L, or T				256	336	403
Zone 15	328	when	<	0.74			and	is	720	to	800	and	is	A, M, or N	and	is	a or m	246	324	404
Zone 16	327	when	is	0.61	to	0.73	and	<	720			and	is	A, M, or N	and	is	a or r	277	324	377
Zone 17	326	when	is	0.67	to	0.74	and	is	620	to	720	and	is	C, G, L, or T	and	is	m or r	281	324	370
Zone 18	324	when	is	0.62	to	0.70	and	is	760	to	780	and	is	C, G, L, or T				264	322	373
Zone 19	315	when	is	0.62	to	0.67	and	<	720			and	is	C, G, L, or T				229	332	420
Zone 20	308	when	is	0.64	to	0.74	and	>=			900	and	is	C, G, L, or T				260	308	351
Zone 21	299	when	is	0.50	to	0.59	and	is	800	to	980	and	is	C, G, L, or T				245	298	359
Zone 22	297	when	is	0.73	to	0.74	and	<	720			and	is	A, M, or N				227	285	341
Zone 23	294	when	is	0.74	to	0.77									and	is	m	254	298	336
Zone 24	290	when	<	0.74			and	>=			800	and	is	A, M, or N	and	is	a or m	235	295	358
Zone 25	277	when	is	0.50	to	0.59	and	<	800			and	is	C, G, L, or T				215	286	340
Zone 26	276	when	is	0.77	to	0.82									and	is	m	232	274	321
Zone 27	273	when	is	0.70	to	0.74	and	is	760	to	800	and	is	C, G, L, or T				238	277	312
Zone 28	269	when	<	0.74			and	>=			720	and	is	A, M, or N	and	is	r	204	272	326
Zone 29	267	when	is	0.74	to	0.82	and	<	640						and	is	a or r	230	264	307
Zone 30	265	when	<	0.61			and	<	720			and	is	A, M, or N				206	275	334
Zone 31	263	when	is	0.47	to	0.50	and	<	980			and	is	C, G, L, or T				200	264	326
Zone 32	245	when	>=			0.82	and	<	520									210	241	284
Zone 33	242	when	is	0.74	to	0.82	and	>=			640				and	is	a or r	199	244	288
Zone 34	240	when	is	0.59	to	0.62	and	<	800			and	is	C, G, L, or T				194	238	289
Zone 35	228	when	>=			0.82	and	>=			520							187	226	268

Although the DT design is of low complexity for reading, a user-friendly table was prepared to classify any location within the study area in each 1 of the 35 productivity zones based on required environment variables and rules defined via the machine-learning algorithm (Table 1). As described above, the total wood volume indicated in the leaf node is an average value outcome from the DT regression, close to 50% probability (50th percentile or median). Thus, it should not assume that the production of a forest zone will be the value read in the DT. The interpretation must always observe the probability of the estimates. Then, we will introduce how to read the DT and interpret the results. Suppose a reader is interested in inferring the productivity levels at location 22.7168 South latitude and 47.6113 West longitude (within Agricultural College Luiz de Queiroz—ESALQ’s facilities). Previously, the reader should have a minimum knowledge of the area of interest but enough to have, on hand, the aridity index, altitude, soil order, and soil texture of the zone to interpret Table 1. The location of interest has an aridity index of 0.75, and this PET/R ratio meets the requirements of Zones 23, 29, and 33. Considering only this input, the reader already has important information, including which the area of interest would, at most, be in position 23 of the productivity zones ranking and that the expected average production is between 242 and 294 m³ ha⁻¹ at 6 years of age. However, the reader knows that the average altitude of the interest location is 590 m, and this rule eliminates Zone 33. For the remaining two zones, and especially in this case, only the soil texture is necessary. However, it was informed that it is a medium texture Latossolo (Oxisol), thus completing the DT interpretation and converging to Zone 23. For Zone 23, the reader will find an average production of 294 m³ ha⁻¹ and 254 m³ ha⁻¹, 298 m³ ha⁻¹, and 336 m³ ha⁻¹ for the percentiles 15th (85% probability of reaching production), 50th (50% probability), and 85th (15% probability), respectively (Table 1), and also 44.5% probability of production reaching 300 m³ ha⁻¹ (Figure 9). We noted that, so far, these approaches are very innovative, incorporating data science and machine learning, dosed with a degree of forester empiricism, and delivering tools with probabilistic results to improve the Eucalyptus productivity zoning.

3.4. Decision-Tree Spatialization

The outputs from the DT were spatialized for the study area providing a reliable cartographic visualization of probability levels of the truly Eucalyptus productivity for a huge region of forest-based industries in Brazil (Figure 10). The ten most productive zones are summarized at 2,380,000 ha and are located in the wetter and higher altitude landscapes of the study area, such as a few mountain ranges in both the northern and southern coastal region and most widely represented in the southwest region of the study area. The ten lowest productivity zones totaled 4,289,007 ha and are mostly spread in the center and northern regions of the study area, where the aridity index is higher due to the combination of the lower rain amount and the higher annual temperature. Low productivity zones are well characterized in the west of the study area, mainly in the Cerrado biome, which also is observed crossing the Paraiba Valley region (Figure 10). Based on more realistic productivity numbers, we hope that our results can foster and improve fair business between the parties, providing knowledge to small rural producers, forest beekeepers, and participants interested in forest outgrower schemes. On the other hand, our study provided deeper analysis and reproducible tools for extension technicians, forestry consultants, crop and forestry companies’ engineers, and researchers to manage their forest assets in a sustainable way.

Figure 10. Decision-tree outputs spatialized for the study area. To learn how environmental variables were used on this map, see Table 1, which shows a straightforward form of how interpretation of decision tree must be performed. Leaf nodes were ranked in decreasing order (from highest to lowest productivity), named Zone 01 to 35, and their values represent the average Eucalyptus total wood volume (m³ ha⁻¹) at 6 years old. The horizontal bar chart shows the area (km²) of each Eucalyptus forest zone spatialized in study area. Blank polygons refer to urban areas, rivers, and dams.

Figure 10. Decision-tree outputs spatialized for the study area. To learn how environmental variables were used on this map, see Table 1, which shows a straightforward form of how interpretation of decision tree must be performed. Leaf nodes were ranked in decreasing order (from highest to lowest productivity), named Zone 01 to 35, and their values represent the average Eucalyptus total wood volume (m³ ha⁻¹) at 6 years old. The horizontal bar chart shows the area (km²) of each Eucalyptus forest zone spatialized in study area. Blank polygons refer to urban areas, rivers, and dams.

3.5. Yield Gap Approach

A reliable estimate of the Eucalyptus wood stock in most of the State of São Paulo was another significant leap forward supported by the approaches presented. We assume that all current 781,654 ha of Eucalyptus plantations (Figure 1) within the study area were 6 years old. The trade secret here was applying the DT for every hectare to predict by classifying it as each one of the productivity zones. We estimated the Eucalyptus wood stock in the study region was approximately 2.44 × 10⁸ m³ for a 50% probability, i.e., considering the median productivity of each zone. Based on the 15th percentiles, there was an 85% probability that the wood stock was greater than 1.99 × 10⁸ m³. Considering the 85th percentile threshold as a proxy for attainable Eucalyptus productivity at each zone, we calculated that current Eucalyptus plantations across the study area would have a potential wood stock at this productivity level by amounting to 2.86 × 10⁸ m³. Other studies have already quantified the Eucalyptus yield gap in Brazil, such as Elli et al. ([11]), but were based on an ecophysiological yield gap perspective, lacking a probabilistic and realistic view of higher productivity levels. The difference between the attainable and average productivities of wood stocks of current Eucalyptus plantations can be considered the yield gap, as that huge amount seems like 4.2 × 10⁷ m³. The sum of each zone’s yield gap ratios by its attainable productivity (85th percentile) is a few more than 113,000 ha or 15% of the current forest base. This is the amount of avoided area expansion in the current Eucalyptus forest base in the study area to accumulate the same wood stock in case the productivity of the zones was raised to the attainable level in each zone. The yield gaps were quantified for each zone, but the results of each one were not shown here. Obviously, there is a great interest in knowing and exploring opportunities for genetics and management to increase productivity in each region of the State of São Paulo. This is a topic that will continue to be addressed in the future in other studies dedicated to these purposes.

Finally, the developed tools and fully applicable results are extremely dependent on large environmental and tree growth databases. It is important to address the potential constraints of our study by other researchers and end users. Thus, it is clear that data-science insights and machine-learning predictions depend on well-collected data, sampling methodological reliability, and databases in reliable quantity and quality. Limitations can arise due to insufficient data, imbalanced data, noise in the data, or a lack of representativeness in certain subgroups. For example, if we use another climatic period, we will have a result of productivity zones different from the one presented here. ML algorithms can reflect and amplify biases present in the training data. For example, in this study, we used representative field data from high-tech forestry; thus, we have high levels of productivity that could not be achieved using forest outgrower schemes or small rural tree growers. The ML models may struggle to generalize new data that falls outside the scope of the training data. Therefore, it should be mentioned that the tools presented here should not be applied outside the study area. Above all, these integrated approaches will mainly require forest professionals with the skills to put it all together and produce insights, applied and reproducible tools, and conclusive results supporting the sustainability of the forestry business in an increasingly rapidly changing world.