Introduction 

The Sierra Madre Occidental (SMO) mountain range is of great ecological interest due to its environmental heterogeneity, which is attributed to the broad physiographical and climatic diversity in the area. Moreover, the SMO shelters uneven-aged and species-rich pine oak forests that are economically important ecosystems for Mexico and other parts of the world. The SMO crosses several states in western Mexico, including the state of Durango (the SMO occupies 71.5% of the surface area of the state). The state covers a total area of 12.3 million ha, of which 9.1 million ha (74.35% of the land in the state) is forestland managed by 11 Regional Forest Management Units (UMAFORES). A large part of the forestland (4.9 million ha) is occupied by temperate forest and is subjected to precipitation levels of between 800 and 1200 mm per year, with frost occurring in winter as a result of the combination of low temperatures and humid winds from the Pacific Ocean; a smaller area of this land (0.5 million), affected by warmer climate, is occupied by forest classified as rainforest ([42]). The state of Durango generates between 25% and 30% of the national timber production, producing a total of 1.5 million m³ of roundwood per year, and promotes forest reserves as important sources of environmental services ([30]).

Estimation of diverse dasometric variables in forest stands is a fundamental aspect of forest inventories. Commercial volume is particularly important in forest management and is valuable for the conservation of diverse forest ecosystems and for improving the economic productivity of forests ([19]). Other stand variables are also important in forest management, such as diameter at breast height, total height, crown cover, basal area and plant biomass. Stand variables are usually measured on the ground by forest inventory teams in field surveys covering large areas ([43]); however, such surveys are both expensive and time-consuming ([21]).

Use of remote sensors is a cost-effective method of collecting data in large inaccessible areas and provides an alternative source of data for extrapolating and estimating forest variables ([16], [37]). The data obtained can be rapidly updated and compared with existing data, and the method can be easily integrated with Geographical Information Systems. These methods have been successfully combined in several studies, and statistically significant correlations between the remotely-sensed data (captured by Landsat satellite sensors) and field-measured data have been reported for diameter at breast height, age and total height ([11]), volume ([23], [37]) and biomass ([31], [27], [50]). The aim of the present study was to analyze the relationship between a set of field-measured dasometric variables and spectral data obtained by the Landsat-8 OLI (Operational Land Imager) sensor and data derived from a digital elevation model, aimed at estimating the total number of trees per plot, crown cover, mean diameter at base (0.30 m), mean diameter at breast height (1.30 m), mean height, mean volume and total volume per plot in a mixed uneven-aged forests in the Sierra Madre Occidental.

  Material and methods 

Study area

The study was carried out in the Ejido La Victoria (surface area: 10.810 ha), located in the Sierra Madre Occidental in the municipality of Pueblo Nuevo, state of Durango, Mexico (Fig. 1). Site elevation ranges between 2245 and 2870 m above sea level, according to the data obtained from the digital elevation model (DEM) developed by the National Institute of Geography and Statistics of Mexico ([26]). The forests in the study area show a rich biodiversity and include at least 12 coniferous tree species (of which 10 are Pinus species) and 5 species of Quercus. The dominant types of vegetation in the area are mixed conifer forest and pine-oak forest ([46], [22]). The main species present are: Pinus cooperi, P. durangensis, P. engelmannii, P. leiophylla, P. strobiformis, P. teocote, P. herrerae, Juniperus deppeana, Quercus sideroxil and Q. crassifolia ([33]).

Fig. 1 - Maps showing the location of the study area and sample plots used in the study: (right, top to bottom) the State of Durango in Mexico; the municipality of Pueblo Nuevo in the State of Durango; the Ejido La Victoria in the municipality of Pueblo Nuevo. The area covered by the Ejido La Victoria is outlined in the large inset.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

Field sampling

A systematic sampling method (1 × 1 km between plots) was used to locate the sampling plots in the study area. In total, 83 circular sampling plots each of surface area 1000 m² were established. The plots were georeferenced relative to a reference system (WGS 84 Datum / UTM zone 13N). The plot data were obtained with the aid of a submeter GNSS receiver, tree calipers (Haglof, Sweden), hypsometer (Vertex IV^®, Haglof, Sweden) and metric tape. The following variables were measured in each plot, taking into account all the trees with diameter at breast height (DBH) larger than 5 cm: tree number (N), diameter at base (DB in cm, measured at 0.30 m above the ground), DBH (in cm, measured at 1.30 m above the ground), total height (HT in m) and crown diameter (CD in m). The volume (V in m³) per plot was estimated using specific allometric equations developed by Corral-Rivas et al. ([8]) for the same study area; the descriptive statistics corresponding to the response variables are shown in Tab. 1.

Datasets

Source of spectral data

The image used in the study was captured by the Landsat-8 OLI sensor (LC80310442014018LGN00) and supplied by the USGS Global Visualization Viewer Server (⇒ http:/­/­glovis.­usgs.­gov/­). The image was captured at about the same time as the field sampling was carried out (January, 2014). The images were obtained as a L1T level product (i.e., they were geometrically corrected with ground control points and the DEM). The satellite image was digitally preprocessed by radiometric correction and atmospheric and topographic techniques, following the procedures established in the ATCOR3^® module (Atcor for Erdas Imagine - [13]) developed for mountainous areas and implemented with ERDAS^® IMAGINE^® 2013 software ([13]). The digital levels (DL) of each of the bands (Tab. 2) were converted to surface reflectance values by reading the model parameters and the orbit of the geometrical satellites included in the image metadata. Bands 1 to 7 (B1, B2, B3, B4, B5, B6 and B7) of Landsat-8 OLI were used in the present study.

A number of vegetation indexes and other derived parameters (normalized difference vegetation index, NDVI; soil adjusted vegetation index, SAVI; leaf area index, LAI; fraction of photosynthetically active radiation, FPAR; albedo, ALB; absorbed shortwave solar radiation, ASR) were computed from the atmospherically and topographically corrected image bands and then included in the stand variables estimation models for evaluation as possible regressor variables (Tab. S1 in Supplementary material).

Terrain variables

These variables have influence on forest species composition, tree height growth, and other forest stand variables ([32], [40]). First order variables included: elevation (DEM), slope (SLOPE), transformed aspect (TRASP), profile curvature (PROFCURV), plan curvature (PLANCURV) and curvature (CURV); and second order variables were: terrain shape index (TERRSHPIN2) and wetness index (WETIND). Terrain parameters were generated in the software ArcGIS 10^® ([14]), with a 5×5 low pass filtered DEM ([26]) and included as candidate variables in the models (Tab. S2 in Supplementary material). These parameters are potentially related to key features affecting forest stand development, such as overall climate characteristics, insolation, evapotranspiration, run-off, infiltration, wind exposure and site productivity. Some of these terrain features are widely used in hydrological, geomorphological and ecological studies ([47]), whereas others are used more specifically for vegetation and forest assessment ([32], [40]).

Dataset integration

The sample plots were geopositioned with the aim of extracting the pixel value average with an associated buffer of 18-meter radius for geolocalization of the plots. The values were extracted using the raster package available in the R statistical software ([39]). Finally, a database was constructed with the mean stand variables values for each plot: the corrected bands of the Landsat-8 OLI sensor (7 bands), the vegetation indexes (6 indexes) and the terrain variables derived from the DEM (8 variables).

Statistical analysis

As the Shapiro-Wilks test revealed that the data were not normally distributed, the multivariate adaptive regression splines (MARS) technique was used for model construction. The MARS technique is an innovative, flexible method that constructs models for continuous, categorical or binary-type dependent and/or independent variables. It is a non-parametric regression technique based on constructing flexible models by fitting data to linear regressions in segments, combination of which generates local models where the relationship between the predictive and response variables may be linear or non-linear and which also incorporate interactions between variables. Each segment of the regression line represents a stretch of linear response to the variable and is denominated the basic function; the final model consists of the entire set of basic functions. The model can be represented so that it separately identifies additive contributions and those associated with different multivariate interactions ([18]). The general form of the MARS non parametric regression model, formulated on the dependent y variable and the x predictors, is as follows (eqn. 1):

where ε is the error, and f (x) is the unknown regression function, derived as follows (eqn. 2):

where β₀ is the intercept of the model, B_m (x) is the m^th basis function, β_m is the coefficient of the m^th basis function, and M is the number of basis functions in the model. Each basis function B_m(x) takes one of the following two forms: (i) a hinge function of the form max(0, x-k) or max(0, k-x), where k is a constant value; (ii) a product of two or more hinge functions, which can therefore model the interaction between two or more predictors (x).

In this study, MARS analysis was performed using the “earth” package in R ([39]). To improve predictions and prevent over-fitting, MARS uses a modified form of the “cross-validation method” (GCV - [45]). GCV is an approximation of the error that would be determined by leave-one-out validation and considers both the residual error and the model complexity. Therefore, the optimal model is that yielding the lowest GCV (eqn. 3):

where y_i represents the observed values of the dependent variable, n is the number of observations, f_m(x_i) is the MARS model with M basis functions, x_i represents the observed values of the predictors included in the MARS model and p_M is the number of parameters of the MARS model.

The importance of the measurement variable was also considered, to take into account whether model information was used or not. Use of the previously mentioned GCV parameter, together with the parameters N of the subsets (N_sub) and the residual sum of squares (RSS), yields reliable results. The N_sub criterion considers the number of model subsets that include the variable, i.e., one subset for each model size from one up to the selected model size (in the final pruned model). The wrapper model utilizes the performance of the algorithm MARS to determine which features should be selected. The subset that produces the highest overall accuracy is deemed the best. For the RSS criterion, the decrease in the RSS is first calculated for each subset relative to the previous subset. For each variable, the decreases are then summed over all subsets that include that variable ([3]).

Three criteria were used to evaluate the model performance: (i) the determination coefficient (R²) as a measure of how well the model fits the training data; (ii) the generalized determination coefficient (GR²), which is based on the GCV statistic and the root mean square error (RMSE) for each model to provide a generalization of the predictive power of the model; and (iii) the common method of k-fold cross-validation (for k = 2 to 10 folds), in which error statistics are calculated across all k trials. Additional information on MARS can be found in Friedman ([18]) and Hastie et al. ([24])

Cross validation was used with the aim of defining a dataset to test the model in the training phase, to minimize problems such as overfitting and to predict how the model will perform with an independent dataset. This provides a good indication of how well the classifier will perform with unseen data and also indicates the degree of variability in the k-fold cross-validation results.

For the purpose of comparison, multiple linear regression (MLR) was also carried out to compare how it performs relative to the MARS technique for modeling stand variables in the study area. MLR is the technique most commonly used in this kind of studies, as it is easy to understand and widely used in most scientific disciplines ([15]). Stepwise selection (forward and backward) was performed to select the most informative variables that were included in the model using the stepAIC function and the exact Akaike’s information criterion (AIC) in the MASS package ([44]). The aim of this procedure was to identify the model with the smallest AIC by removing or adding predictors.

The equations obtained with the MARS technique (R² and RMSE) were finally used to generate a map of the stand variables by means of the map algebra and conditional tools of the GIS package ArcGIS 10^® ([14]).

  Results 

Tab. 3 shows the results produced by the MLR and MARS techniques for each of the dasometric variables considered in the study. Although MARS generally performed best, MLR yielded better prediction of the H mean variable.

Tab. 4 shows the results obtained with the MARS models that performed best for estimating the dasometric variables considered in each plot of mixed and uneven aged forest plot in the study area.

According to the selection criteria, the model used to estimate the mean diameter at breast height (1.30 m), mean volume and the mean crown diameter were those that performed best with the training data set (i.e., highest R² = 0.73, 0.57 and 0.55; GR² = 0.26, 0.34 and 0.30; and lowest RMSE = 3.01 cm, 0.20 m³ and 0.61 m; 16.79%, 47.39% and 15.83% of the mean variables, respectively). The model used to estimate the total volume per plot also performed reasonably well (R² = 0.41, GR² = 0.31, and RMSE = 9.00 m³; 36.02% of the mean total volume per plot). The model used to estimate the number of trees per plot yielded the poorest fit (R² = 0.22, GR² = 0.04, with RMSE = 51.01; 64.14% of the mean number of trees per plot) and the model used to estimate the mean height also provided a poor fit to the training data set (R² = 0.25, GR² = 0.21 and RMSE = 2.73 m; 24.96% of the mean height). Thus a large difference between GR² and R² indicates a severely over-parametrized model, i.e., a model that may show a good fit to the training dataset, but would not be generally applicable.

The predictor variables (spectral bands and indices and topographic variables) used to construct the basic functions can be evaluated after fitting a MARS model to the different dasometric variables. The parameter estimates of the MARS models fitted to the complete database are shown in Tab. 5.

The number of the predictor variables and the percentage contribution to the predictive power of the models varied widely in the models fitted (Fig. 2). For instance, the following variables were the most important (100% of GCV) in some models: slope (SLOPE), for estimating number of trees per plot; soil adjusted vegetation index (SAVI), for estimating mean diameter at base per plot, mean total tree height, mean crown diameter and mean volume; and the green band (B3) and short wave infrared 2 band (B7) for estimating the mean diameter at breast height and the total volume per plot. Likewise, the following variables contributed least to the predictive ability of the models (<25% of GCV): fraction of photosynthetically active radiation (FPAR), for estimating number of trees per plot; albedo (ALB) for estimating the mean diameter at base and the total volume per plot; the blue band (B2), for estimating the mean diameter at breast height, and DEM, wetness index (WETIND) and the red band (B4), for estimating the mean volume.

Fig. 2 - Importance and selection of predictor variables in the multivariate adaptive regression splines (MARS) models for each radiometric correction algorithm considered.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

Fig. 3 shows the variation in k-fold cross-validation (k = 2 to 10), which is used to visualize problems such as overfitting. And finally, the spatial distribution of the dasometric variables in the study area obtained by application of the MARS models (Tab. 5) is shown in Fig. 4.

Fig. 3 - Variations in the goodness of fit statistics obtained by applying the cross-validation technique (nfolds = 2-10) to the MARS models.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

Fig. 4 - Spatial distribution of the stand variables in the study area, obtained by application of the MARS technique.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

  Discussion 

Predictive modeling is not an exact science ([17]), and model performance can be influenced by the method of statistical analysis, topographical factors (slope, exposure and elevation) and forest structure (proportion of young and mature trees, species and type of silvicultural management). The results obtained in this and other studies ([34], [51], [1]) were based on nonparametric techniques. In this case, the MARS technique yielded significantly more accurate predictions than the MLR excepting the variable H mean as a consequence of the smaller number of predictors (1 of 21) in the MARS model.

However, the use of very high resolution sensors can improve prediction of forest structural parameters by linear regression analysis, as reported by Ozdemir & Karnieli ([36]). These authors used multispectral data (WordView-2) and a stepwise regression technique to predict the number of trees per hectare (R² = 0.38, RMSE = 109.56 tree ha^-1; 28.45%) and total volume per plot (R² = 0.42, RMSE = 27.18 m³ ha^-1; 44.20%) in a plantation forest in Israel.

Trotter et al. ([43]) obtained high goodness-of-fit values (R² ≥ 0.85 and RMSE = 39.00 m³ ha^-1) for band 4, and for bands 3 and 4 of Landsat TM satellite images for estimating timber volume by simple and multiple linear regression methods, respectively. In a study monitoring growth of a young forest plantation in UK, they used generalized linear models (GLM) with logarithmic links and normal errors, and obtained higher values (R² = 0.85 and RMSE = 0.88 m; 28.38%) for the mean height, and lower values for the diameter at breast height (R² = 0.36 and RMSE = 2.15 cm; 27.24%) and for density (R² = 0.03 and RMSE = 3.260.59 tree ha^-1; 98.62%) than those obtained in the present study. The differences in the results can be explained by the more heterogeneous conditions of the native vegetation in the present study area and, in the case of height, the higher R² values obtained in UK may also be related to the larger latitude in this country, which implies longer shadows at the time of the Landsat overpass, making it easier to estimate tree heights ([11]). In a study modeling forest stand structure attributes by multivariate regression and discriminant function analysis, Hall et al. ([23]) obtained higher error values for height and crown cover (R² = 0.57 and 0.65, and RMSE = 12.0 and 2.8 m, respectively). Coulston et al. ([9]) estimated crown cover in forests in Utah and Michigan (USA) by using a beta regression model and the random forest algorithm, and obtained also high error values (R² = 0.65 and 0.69 and RMSE = 0.87 and 0.89 m, respectively).

Cruz-Leyva et al. ([10]) modeled even-aged Pinus patula and P. teocote forest stand variables in Hidalgo, Mexico by GIS variables and reflectance data derived from a multispectral SPOT 5 image and multiple linear regression. The volume was significantly related to tree crown cover, leaf area index (LAI), slope, elevation, mean annual temperature, maximum temperature of warmest month, mean annual precipitation, and precipitation of the wettest quarter (R² = 0.72 and 0.88, p < 0.001).

Li et al. ([29]) used stepwise linear regression and regression tree analysis, together with cross validation methods, to estimate height in young forests and obtained the following values for the first and second modeling techniques, respectively: R² = 0.15 and 0.19 and RMSD (root mean square of the difference) = 10.32 and 6.07 (m) for the first group of variables by combining bands 1-5 and 7 to obtain the integrated forest index (IFI), a measure of the probability that a pixel corresponds to forest, the normalized difference vegetation index (NDVI) and the normalized burn ratio index (NBRI). The same authors obtained R² = 0.88 and 0.89 and RMSD = 2.42 and 2.13 (m) for the second group of variables (age since disturbance, derived from Vegetation Change Tracker [VCT] and Landsat time series stacks [LTSS] disturbance year map), and R² = 0.69 and 0.68 with RMSD =3.85 and 3.73 (m) successively for the third group of variables (combining accumulated indexes: IFI, NDVI and NBRI).

The null contribution of the NDVI (Fig. 2) can be attributed to the significant association between this index and both the leaf area index and the fraction of photosynthetically active radiation intercepted ([28]). The values of these variables are low in winter, which is when the satellite image was captured and the field sampling was carried out for the present study.

Some estimates are difficult to compare, because of differences in the sources of data, estimation methods used, period of data compilation and large differences between the study areas ([50]).

The top performing model (R²) was that used to estimate the mean diameter at breast height, probably because this model included the largest number of independent variables (Tab. 4), indicating problems of overfitting (GR²) in the training data set (Fig. 3). The remaining unexplained variance in the spatial patterns may be related to the high heterogeneity of the dependent variables ([20]), due to the uneven-aged and tree species-rich structure of the Pueblo Nuevo forests in Durango. We therefore recommend increasing the number of observations in order to decrease the possibility of occurrence of this effect.

Considering the advantages and limitations of different remote sensing images, the medium-resolution (pixel size, 30 m) Landsat series is one of the most widely used for estimating dasometric variables ([2], [38], [12], [52]). The advantages of using the Landsat series are that numerous historical spatiotemporal archives are available and that the sensor is cheaper and more accessible than high resolution sensors, particularly for analysis of large areas ([48]).

  Conclusions 

We conclude that the spectral bands of the Lansat-8 OLI sensor can be used to estimate four of the seven stand variables (mean diameter at breast height, mean crown diameter, mean volume and total volume per plot) in the study region. Furthermore, the inclusion of aggregated values of vegetation indexes and terrain variables derived from the DEM contributed significantly to increasing the the goodness-of-fit of the models.

Multivariate adaptive regression splines (MARS) models are more flexible than linear regression models. Application of the MARS technique facilitated handling and interpretation of large amounts of data representing complex structures. The predictors thus obtained for the study area were accurate and the method proved suitable for modeling a large number of predictive variables.

  References

  Supplementary Material

Tab. S1 - Vegetation indexes and other derived parameters generated from the spectral bands of the sensor.
Tab. S2 - Terrain attributes derived from the digital elevation model (DEM).