Evaluation of Radiometric and Atmospheric Correction Algorithms for Aboveground Forest Biomass Estimation Using Landsat 5 TM Data

remote sensing Article Evaluation of Radiometric and Atmospheric Correction Algorithms for Aboveground Forest Biomass Estimation Using Landsat 5 TM Data Pablito M. López-Serrano 1, José J. Corral-Rivas 2,†, Ramón A. Díaz-Varela 3,†, Juan G. Álvarez-González 4,† and Carlos A. López-Sánchez 2,*,† 1Doctorado en Ciencias Agropecuarias y Forestales, DICAF, Universidad Juárez del Estado de Durango, Boulevard del Guadiana 501, Ciudad Universitaria, Torre de Investigación, 34120 Durango, Dgo, Mexico; [email protected] 2Instituto de Silvicultura e Industria de la Madera, Universidad Juárez del Estado de Durango, Boulevard del Guadiana 501, Ciudad Universitaria, Torre de Investigación, 34120 Durango, Dgo, Mexico; [email protected] 3Departamento de Botánica–GI-1908-BIOAPLIC, Universidad de Santiago de Compostela, Escuela Politécnica Superior, 27002 Lugo, Spain; [email protected] 4Departamento de Ingeniería Agroforestal, Universidad de Santiago de Compostela, Escuela Politécnica Superior, 27002 Lugo, Spain; [email protected] *Correspondence: [email protected]; Tel.: +52-618-8251886 † These authors contributed equally to this work. Academic Editors: James Campbell and Prasad S. Thenkabail Received: 21 January 2016; Accepted: 20 April 2016; Published: 29 April 2016 Abstract: Solar radiation is affected by absorption and emission phenomena during its downward trajectory from the Sun to the Earth’s surface and during the upward trajectory detected by satellite sensors. This leads to distortion of the ground radiometric properties (reflectance) recorded by satellite images, used in this study to estimate aboveground forest biomass (AGB). Atmospherically-corrected remote sensing data can be used to estimate AGB on a global scale and with moderate effort. The objective of this study was to evaluate four atmospheric correction algorithms (for surface reflectance), ATCOR2 (Atmospheric Correction for Flat Terrain), COST (Cosine of the Sun Zenith Angle), FLAASH (Fast Line-of-sight Atmospheric Analysis of Spectral Hypercubes) and 6S (Second Simulation of Satellite Signal in the Solar), and one radiometric correction algorithm (for reflectance at the sensor) ToA (Apparent Reflectance at the Top of Atmosphere) to estimate AGB in temperate forest in the northeast of the state of Durango, Mexico. The AGB was estimated from Landsat 5 TM imagery and ancillary information from a digital elevation model (DEM) using the non-parametric multivariate adaptive regression splines (MARS) technique. Field reference data for the model training were collected by systematic sampling of 99 permanent forest growth and soil research sites (SPIFyS) established during the winter of 2011. The following predictor variables were identified in the MARS model: Band 7,Band 5,slope ( β ), Wetness Index (WI), NDVI and MSAVI2. After cross-validation, 6S was found to be the optimal model for estimating AGB (R 2 = 0.71 and RMSE = 33.5 Mg¨ha´1 ; 37.61% of the average stand biomass). We conclude that atmospheric and radiometric correction of satellite images can be used along with non-parametric techniques to estimate AGB with acceptable accuracy. Keywords: multivariate adaptive regression splines; remote sensing; radiometric correction algorithms; terrain features 1. Introduction Apart from sensor gains/offsets, solar irradiance and Sun-Earth geometry, the characteristics of electromagnetic energy detected by remote sensing optical sensors is affected by particles and other Remote Sens. 2016,8, 369; doi:10.3390/rs8050369 www.mdpi.com/journal/remotesensing Remote Sens. 2016,8, 369 2 of 19 components present in the atmosphere. Hence, as solar radiation passes through the atmosphere (Sun-surface-sensor), it is affected by absorption or scattering by particles in suspension (aerosols) and other atmospheric elements, thus creating a hazy effect that distorts the radiometric properties of satellite images [ 1 – 4 ]. Various algorithms have been developed to correct such effects in the image. These include calibrations based on sensor parameters, solar-Earth geometry, dark object subtraction and radiative transfer [ 5 – 8 ]. Correction of atmospheric effects is important in relation to improving data quality [ 9 , 10 ]. However, the use of correction algorithms to minimize atmospheric effects, especially the scattering and absorption caused by aerosols, remains challenging [ 11 ]. This particularly applies to the parameterization of algorithms for calculating the surface reflectance for its eventual use in estimating aerial or aboveground forest biomass (AGB) [12]. AGB is an important parameter that provides information about the current status of forests and therefore facilitates forest management decisions [ 13 , 14 ]. The use of optical remote sensors of medium spatial resolution, such as the Landsat Thematic Mapper (TM), has become increasingly common in the last few decades. The use of these sensors for large-scale monitoring of AGB dynamics, particularly in the case of multi-temporal/multi-scene analyses or unevenly-distributed atmospheric effects in single-image analyses, frequently involves the application of different atmospheric correction algorithms. Such algorithms correct distortions between or within images other than those related to real land cover differences [15–19]. The aim of the present study was to evaluate four different atmospheric correction algorithms (surface reflectance)—COST (Cosine of the Solar Zenith Angle), a modification of the dark object subtraction method [ 20 ] with the addition of a multiplicative correction for atmospheric transmittance [ 21 ], ATCOR2 (Atmospheric Correction for Flat Terrain), [ 22 ]), FLAASH (Fast Line-of-sight Atmospheric Analysis of Spectral Hypercubes, [ 23 , 24 ]) and 6S (Second Simulation of Satellite Signal in the Solar Spectrum, [ 25 ]), and one radiometric correction algorithm (apparent reflectance at sensor), ToA (Top of Atmosphere), for estimating AGB from spectral data captured by the Landsat 5 TM sensor and topographic variables derived from a digital elevation model (DEM) generated using data collected by sampling of the permanent forest growth and soil research sites (SPIFyS) established in a temperate forest in the northwest of the state of Durango, Mexico. Different methods could be used for remote AGB estimation in forests, including parametric and nonparametric approaches. However, since forest structure and biomass often entail nonlinear variability, variable interaction across scales and autocorrelation, nonparametric approaches often markedly outperform parametric methods [26]. In this study, we investigate remote AGB estimation of mixed and uneven-aged forests using the multivariate adaptive regression splines (MARS). To the best of our knowledge, MARS have rarely been used for remote estimation of AGB [ 27 – 29 ], but in these cases, the method has performed best for prediction than other parametric and non-parametric approaches, such us linear models, classification, regression trees and artificial neural networks [27] or hybrid tree-based algorithms [28,29]. 2. Materials and Methods 2.1. Study Area The study area is located in the northwest of the state of Durango (Figure 1). The forest vegetation comprises pine (Pinus spp.), oak (Quercus spp.), fir (Abies spp.) and plurispecific stands with different proportions of species of the genera Pinus and Quercus, in accordance with the Land and Vegetation Use Map, scale 1:250,000, Series V [ 30 ]. The weather is cold in the canyons (with temperatures ranging from ´ 20 ˝ C in winter to 20 ˝ C in summer) and mild or warm in the valleys (with temperatures ranging from 10 ˝C in winter to 40 ˝C in summer). Remote Sens. 2016,8, 369 3 of 19 Remote Sens. 2016, 8, 369 3 of 18 Figure 1. Location of the study area and of the permanent forest growth and soil research sites (SPIFyS). 2.2. Field Data The field data were AGB values taken from a database of stand variables measured in 99 permanent forest growth and soil research sitesSPIFyS established during the winter of 2011 in accordance with the method developed by Corral-Rivas et al. [31]. The SPIFyS are square plots of a size of 50 × 50 m (surface area, 2500 m2). The location of the sites was carried out by systematic sampling at the regional level using a squared grid of five kilometers, although a small squared grid (3 km) was used on very steep areas to include their variability. In each sample plot, the tree species were recorded, and the diameter at breast height (cm) and total height (m) of all standing trees were measured. Species-specific models developed for the study area were used to estimate individual tree aboveground biomass [32]. The goodness of fit for such statistical models ranged from 0.87–0.99 (R2), and the root mean square error (RMSE) varied from 22.8–95.2 kg. Once the tree aboveground biomass was estimated, all values from each sampling plot were summed and expressed on a per hectare basis. Table 1 shows the descriptive statistics of the main stand variables. Table 1. Descriptive statistics of the main stand variables estimated from the 99 permanent sample plots. The dominant stand height was calculated as an average value from the 100 thickest trees per hectare. Variable Mean Standard Deviation Minimum Value Maximum Value Number of stems per ha 655.47 322.25 224 2264 Diameter at breast height (cm) 18.44 3.46 11.69 31.12 Dominant height (m) 14.62 3.72 6.87 24.81 Stand biomass (Mg·ha−1) 89.03 43.45 2.70 234.03 2.3. Spectral Data from the Landsat 5 TM Sensor The image used in the study, which was captured by the Landsat 5 TM sensor in April 2011, is available at the U.S. Geological Service (USGS) website [33]. The sensor operates in seven bands of the electromagnetic spectrum: Band 1 (blue, 450–520 nm), Band 2 (green, 520–600 nm), Band 3 (red, Figure 1. Location of the study area and of the permanent forest growth and soil research sites (SPIFyS). 2.2. Field Data The field data were AGB values taken from a database of stand variables measured in 99 permanent forest growth and soil research sitesSPIFyS established during the winter of 2011 in accordance with the method developed by Corral-Rivas et al. [ 31 ]. The SPIFyS are square plots of a size of 50 ˆ 50 m (surface area, 2500 m 2 ). The location of the sites was carried out by systematic sampling at the regional level using a squared grid of five kilometers, although a small squared grid (3 km) was used on very steep areas to include their variability. In each sample plot, the tree species were recorded, and the diameter at breast height (cm) and total height (m) of all standing trees were measured. Species-specific models developed for the study area were used to estimate individual tree aboveground biomass [ 32 ]. The goodness of fit for such statistical models ranged from 0.87–0.99 (R 2 ), and the root mean square error (RMSE) varied from 22.8–95.2 kg. Once the tree aboveground biomass was estimated, all values from each sampling plot were summed and expressed on a per hectare basis. Table 1shows the descriptive statistics of the main stand variables. Table 1. Descriptive statistics of the main stand variables estimated from the 99 permanent sample plots. The dominant stand height was calculated as an average value from the 100 thickest trees per hectare. Variable Mean Standard Deviation Minimum Value Maximum Value Number of stems per ha 655.47 322.25 224 2264 Diameter at breast height (cm) 18.44 3.46 11.69 31.12 Dominant height (m) 14.62 3.72 6.87 24.81 Stand biomass (Mg¨ha´1)89.03 43.45 2.70 234.03 2.3. Spectral Data from the Landsat 5 TM Sensor The image used in the study, which was captured by the Landsat 5 TM sensor in April 2011, is available at the U.S. Geological Service (USGS) website [ 33 ]. The sensor operates in seven bands of the electromagnetic spectrum: Band 1 (blue, 450–520 nm), Band 2 (green, 520–600 nm), Band 3 (red, 630–690 nm ), Band 4 (near-infrared, 760–900 nm), Bands 5 and 7(mid-infrared region, 1550–2350 nm) Remote Sens. 2016,8, 369 4 of 19 and Band 6 (which provides information in the thermal infrared region and is not used in this type of study). Likewise, the spectral indices Normalized Difference Vegetation Index (NDVI, [ 34 ]) and Modified Soil Adjusted Index (MSAVI2, [ 35 ]) were used. Such indexes are potentially less sensitive than single band values to artefacts due to differences in light conditions, exposed soil or terrain slope and orientation, which might eventually affect AGB estimation. 2.4. Radiometric Correction Algorithms A useful prior step to the interpretation of satellite images is converting the digital numbers (DNs) stored in the original image into biophysical variables of standard significance (reflectance). These variables are comparable in the same sensor over time and over scenes, as well as between different sensors and between remote sensing and other methods of detecting electromagnetic energy. The correction is also advisable in the case of unevenly-distributed atmospheric effects in the image and also when vegetation indexes based on band ratios are included in the analyses [ 36 – 41 ]. In order to calculate surface reflectance, the atmospheric effects must be removed. This involves estimating the atmospheric transmissivity (both up and down welling), diffuse irradiance and atmospheric radiance due to scattering [ 42 ]. In the present study, we evaluated five different radiometric correction methods: the first four methods (namely ATCOR2, COST, FLAASH and 6S) including atmospheric corrections (surface absolute reflectance), while the fifth (ToA) without considering atmospheric effects (apparent reflectance at the sensor) for the final estimation of AGB. 2.4.1. Atmospheric Correction for Flat Terrain (ATCOR2) The aim of ATCOR2 correction [ 22 ] is to remove the atmospheric effects in order to recover the physical parameters of the Earth’s surface, including the surface reflectance, soil visibility and the temperature. The correction is carried out with the ATCOR module [ 43 ] included in ERDAS IMAGINE software, Version 2013 [44]. 2.4.2. Cosine of the Sun Zenith Angle (COST) This is a radiometric calibration method that considers the atmospheric effect and is based entirely on the characteristics of the satellite image, in contrast with other methods of atmospheric correction, like ATCOR2, FLAASH or 6S, requiring some extra parameters, such as atmospheric profiles, the aerosol models or visibility [ 21 , 45 ]. COST applies Dark Object Subtraction (DOS, [ 20 ]) to compensate for the additive components of the atmosphere, which mainly affect the shortest wavelengths. DOS does not take into account the multiplicative effect on the longer wavelengths. For initial estimation of the multiplicative effect, the value of the atmospheric transmittance along the path from the ground to the sensor (TAUz) was computed from the cosine of the solar zenith angle [ 46 ]. This correction was carried out by implementing the algorithm in the Model Maker®module of ERDAS IMAGINE software, Version 2013 [44]. 2.4.3. Fast Line-of-Sight Atmospheric Analysis of Spectral Hypercubes (FLAASH) This is an advanced atmospheric correction module based on the MODTRAN4 algorithm of radiative transfer developed by Spectral Sciences Inc. (Burlington, MA, USA) under the sponsorship of the U.S. Air Force Research Laboratory [ 23 , 24 ]. The technique is initially based on the standard equation of spectral radiance for each pixel of the sensor, which applies to the range of wavelengths of solar radiation (thermal emission is omitted) and a Lambertian and flat surface or their equivalents. It considers the radiance reflected from the Earth’s surface and scattered by the atmosphere towards the sensor. The difference between these two radiances is due to the adjacency effect (spatial mixture of radiance among nearby pixels) caused by atmospheric scattering. The correction was carried out with the FLAASH module in ENVI®software, Version 5.1 [47]. Remote Sens. 2016,8, 369 5 of 19 2.4.4. Second Simulation of Satellite Signal in the Solar Spectrum (6S) This procedure eliminates atmospheric effects on the reflectance values in images captured by sensors on board satellite or aircraft platforms. It is based on an advanced code of radiative transfer, designed for simulating the interaction between solar radiation and an atmospheric-surface system, together with a wide range of atmospheric, spectral and geometric conditions [ 25 ]. The code acts on the basis of the successive order of scattering (SOS) method and explains the polarization of radiation in the atmosphere by calculating the different components. The scenes used to belong to the National Landsat Archive Processing System (NLAPS) and correspond to the product obtained by Landsat 4–5 Thematic Mapper Level 1 of surface reflectance processed by the Standard Landsat Product Generation System (LPGS) and using the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) algorithm (available at the U.S. Geological Service website [33]). 2.4.5. Apparent Reflectance at the Top of Atmosphere (ToA) This technique enables calculation of the apparent reflectance in a satellite image and consists of converting the DNs to radiance values and then to reflectance values. The word “apparent” refers to the fact that the reflectance has not been corrected for atmospheric effects and represents an initial normalization of the image [ 48 ]. The correction was carried out with the “apparent reflectance” method implemented in the “landsat” package [49] in the R software [50]. 2.5. Parameters Derived from the Digital Elevation Model (DEM) The primary and secondary terrain parameters were calculated from the DEM (Table 2), after the application of a low bandpass filter with a moving window of 5 ˆ 5, with the aim of reducing the banding effect present in the original archive [ 51 ]. These parameters are considered as predictor variables for estimating AGB, because they are directly related to the distribution and development of forest species and can subsequently be evaluated [52,53]. Table 2. Additional field variables used to model the aboveground forest biomass (AGB). Variable Formula Reference Elevation Digital Elevation Model Slope (β)β= arctan ”pG2`H2q1{2ı Transformed aspect (Trasp)Trasp =1´cosppπ{180qpα´30qq 2[53] Terrain Shape Index (TSI)TSI =Z/R [52] Wetness Index (WI)WI = ln (As/tanβ)[54] Profile curvature (Ø)Ø“ ´2DG2`EH2`FGH G2`H2[55] Plan curvature (ω)ω“2DH2`EG2`FGH G2`H2 Curvature (x) x=ω´Ø where: Z , average elevation; R, point radio altitude units; As, drainage area specified; tan β , local slope angle; D, F,Gand Hwere derived according to the equation of [56]. 2.6. Generation of a Database The database used to estimate AGB was generated by extracting the mean values of the satellite image pixels after each of the radiometric correction techniques and of the different DEM parameters by considering a buffer of 25 m for the geolocalization of the SPIFyS plots. The values were extracted using the “raster” package [57], implemented in R statistical software [50]. 2.7. Statistical Analysis 2.7.1. Analysis of Variance In order to establish whether there were any significant differences between the radiometric correction algorithms considered, analysis of variance (ANOVA) was used to compare the values of Remote Sens. 2016,8, 369 6 of 19 each spectral band of the sensor obtained by processing the images of the SPIFyS and applying the different correction algorithms. The correction algorithm was used as a fixed factor: yi“µ`CAi`εi(1) where y i are the values of each spectral band of the sensor corrected by the algorithm i;CA i is the correction algorithm used (fixed factor) and εiare the associated errors. Prior to the application of the ANOVA, the dependent variables normality and homogeneity of variance were checked using respectively Shapiro–Wilk’s test [ 58 ] and Levene’s test [ 59 ]. When the results of these tests indicated that the ANOVA assumptions were not satisfied, the dependent variables were transformed by a natural log transformation. When the ANOVA indicated significant differences between correction algorithms ( α = 0.05), Tukey’s multiple comparison test was used to identify which algorithms were different. Likewise, box plots were elaborated for each variable (spectral band) with the aim of facilitating graphical interpretation of the results. All analyses were carried out with R statistical software [50]. 2.7.2. Multivariate Adaptive Regression Spline (MARS) The statistical analysis was carried out by the non-parametric multivariate adaptive regression splines (MARS) proposed by Friedman [ 60 ], using the “earth” package [ 61 ] implemented in the R software [ 50 ]. This method involves constructing a non-linear regression model based on a product of functions denominated “smoothed basis functions” (splines). The predictors are incorporated in their structure as part of a function generating a model for the dependent variable, which may be continuous or binary, and at the same time automatically selects the predictors in the final model. The general form of the MARS non-parametric regression model, formulated on the dependent “y” variable and the “x” predictors, is as follows: y“fpxq ` ε(2) where εis the error and f(x) is the unknown regression function, derived as follows: fpxq “ β0` M ÿ m“1 βmBmpxq(3) where β0 is the intercept of the model, B m (x) is the m-th basis function, β mis the coefficient of the m-th basis function and Mis 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 kis a constant value; (ii) a product of two or more hinge functions that, therefore, can model the interaction between two or more predictors (x). The optimal model was selected using the backward procedure. An overfitted model was generated with all possible predictor variables, and the model was pruned by removing terms one by one, deleting the least effective term at each step until the best submodel was found. Model subsets were compared using the generalized crossed validation (GCV) criterion. GCV is an approximation to 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 will be that yielding the lowest GCV. GCV pMq “ řn i“1pyi´fMpxiqq2 n`1´pM n˘2(4) where yi are the observed values of the dependent variable, nis the number of observations, fMpxiq is the MARS model with Mbasis functions, xi are the observed values of the predictors included in the MARS model and pMis the number of parameters of the MARS model. To analyze the importance of the independent variables that contribute most to the final model, three different criteria were used: (i) nsubsets,i.e., the number of times that each variable is included in a subset (in the final model); (ii) sqr_rss,i.e., the decrease in the sum of square errors (RSS) for each Remote Sens. 2016,8, 369 7 of 19 subset relative to the previous subset; and (iii) sqr_gcv,i.e., the decrease in the GCV value for each subset relative to the previous subset. These criteria are normalized on a scale of 100, to facilitate the interpretation of the contribution of each predictive variable in the model. With the aim of evaluating the performance of the fitted model for each of the radiometric correction algorithms, the following statistics were determined to compare the AGB reference data for the monitoring plots with the estimates based on the images subjected to different corrections: root mean squared error (RMSE), the coefficient of determination (R 2 ) and the generalized coefficient of determination (GR2): RMSE “ g f f f e n ř i“1 pyi´ˆ yiq2 n´p(5) R2“1´ n ř i“1 pyi´ˆ yiq2 n ř i“1`yi´yi˘2(6) GR2“1´GCV GCVnull (7) where y i , ˆ y i and y are, respectively, the observed, estimated and mean values of the dependent variable; nis the total number of observations used to fit the model; pis the number of parameters to estimate; and GCVnull is the GCV of a model with a single independent term (β0). Because the quality of the fit does not necessarily reflect the quality of the prediction, assessment of the validity of the model with an independent dataset is desirable [ 62 ]. Therefore, cross-validation was carried out by splitting the dataset into different numbers of folds (2–10), and the goodness-of-fit statistics (Equations (5)–(7)) were calculated for each fold (test data) by using the models fitted to the other folds (training data). 2.8. Generation of Thematic Maps The final stage of the overall process illustrated was the calculation of AGB maps. These maps were derived by implementing the basic functions resulting from the MARS model fit for each of the atmospheric correction algorithms by using the raster calculator in ArcGIS®software [63]. An overview of the workflow described in the previous sections is presented in Figure 2. Remote Sens. 2016, 8, 369 7 of 18 subset relative to the previous subset. These criteria are normalized on a scale of 100, to facilitate the interpretation of the contribution of each predictive variable in the model. With the aim of evaluating the performance of the fitted model for each of the radiometric correction algorithms, the following statistics were determined to compare the AGB reference data for the monitoring plots with the estimates based on the images subjected to different corrections: root mean squared error (RMSE), the coefficient of determination (R 2 ) and the generalized coefficient of determination (GR 2 ): () pn yy n i ii − − =  =1 2 ˆ RMSE (5) () ()   = = − − −= n i ii n i ii yy yy R 1 2 1 2 2 ˆ 1 (6) =1−  (7) where y i , ŷ i and  are, respectively, the observed, estimated and mean values of the dependent variable; n is the total number of observations used to fit the model; p is the number of parameters to estimate; and GCV null is the GCV of a model with a single independent term (β 0 ). Because the quality of the fit does not necessarily reflect the quality of the prediction, assessment of the validity of the model with an independent dataset is desirable [62]. Therefore, cross-validation was carried out by splitting the dataset into different numbers of folds (2–10), and the goodness-of-fit statistics (Equations (5)–(7)) were calculated for each fold (test data) by using the models fitted to the other folds (training data). 2.8. Generation of Thematic Maps The final stage of the overall process illustrated was the calculation of AGB maps. These maps were derived by implementing the basic functions resulting from the MARS model fit for each of the atmospheric correction algorithms by using the raster calculator in ArcGIS ® software [63]. An overview of the workflow described in the previous sections is presented in Figure 2. Figure 2. Flow diagram of the processes involved in estimating AGB. Figure 2. Flow diagram of the processes involved in estimating AGB. Remote Sens. 2016,8, 369 8 of 19 3. Results In the first step, a visual comparison of the application of the different radiometric correction algorithms to the Landsat 5 TM image of the study area was done. All atmospheric correction algorithms were parameterized using the same type of “rural” aerosol, for zones not affected by urban or industrial sources and a visibility larger than 40 km. In comparison with the original image in DNs, the brightness of the other images differs due to the parameterization used in each algorithm to correct the radiometric and atmospheric effects (Figure 3). Remote Sens. 2016, 8, 369 8 of 18 3. Results In the first step, a visual comparison of the application of the different radiometric correction algorithms to the Landsat 5 TM image of the study area was done. All atmospheric correction algorithms were parameterized using the same type of “rural” aerosol, for zones not affected by urban or industrial sources and a visibility larger than 40 km. In comparison with the original image in DNs, the brightness of the other images differs due to the parameterization used in each algorithm to correct the radiometric and atmospheric effects (Figure 3). Figure 3. Comparison of the application of the different radiometric correction algorithms: (a) original image in DNs (multiband composition 5, 4 and 3); (b) ATCOR2; (c) COST; (d) FLAASH; (e) 6S; and (f) ToA. The mean reflectances of the 99 SPIFyS plots obtained for each band of the Landsat 5 TM sensor and radiometric correction algorithm were graphically compared (Figure 4). The differences vary according to the atmospheric effects corrected by each particular algorithm. In general, the spectral signature of the vegetation showed in all cases typical behavior: low reflectance in the bands of the visible spectrum (450–690 nm) and high reflectance in the near- and mid-infrared bands (760–2350 nm). However, a clear difference appears depending on whether the radiometric correction algorithms were applied to variables with physical meaning (e.g., ground reflectance) rather than to DNs. In this respect, the ATCOR2 and FLAASH algorithms overcorrected the atmospheric effects (underestimation of transmissivity) in the short wavelength regions and underestimated them in the near- and mid-infrared regions, while COST and ToA produced the inverse spectral pattern. Finally, the performance of 6S was intermediate between that of the other algorithms evaluated. Figure 3. Comparison of the application of the different radiometric correction algorithms: ( a ) original image in DNs (multiband composition 5, 4 and 3); ( b ) ATCOR2; ( c ) COST; ( d ) FLAASH; ( e ) 6S; and (f) ToA. The mean reflectances of the 99 SPIFyS plots obtained for each band of the Landsat 5 TM sensor and radiometric correction algorithm were graphically compared (Figure 4). The differences vary according to the atmospheric effects corrected by each particular algorithm. In general, the spectral signature of the vegetation showed in all cases typical behavior: low reflectance in the bands of the visible spectrum (450–690 nm) and high reflectance in the near- and mid-infrared bands (760–2350 nm). However, a clear difference appears depending on whether the radiometric correction algorithms were applied to variables with physical meaning (e.g., ground reflectance) rather than to DNs. In this respect, the ATCOR2 and FLAASH algorithms overcorrected the atmospheric effects (underestimation of transmissivity) in the short wavelength regions and underestimated them in the near- and mid-infrared regions, while COST and ToA produced the inverse spectral pattern. Finally, the performance of 6S was intermediate between that of the other algorithms evaluated. Remote Sens. 2016,8, 369 9 of 19 Remote Sens. 2016, 8, 369 9 of 18 Figure 4. Spectral performance of each of the radiometric correction algorithms considered. As the Shapiro–Wilk’s test indicated that the corrected values of the spectral bands of the sensor were not normally distributed, therefore a log transformation was used. The ANOVA applied to the data from the 99 SPIFyS plots detected significant differences for all algorithms in the visible region (Figure 5), except for the COST and FLAASH algorithms in Band 3. Significant differences were observed in the near-infrared region for all algorithms tested, except ATCOR2 and 6S, which were included in the same population group (Tukey’s test); they also shared characteristics with FLAASH and ToA. Finally, in the bands corresponding to mid-infrared regions, the algorithms ATCOR2, FLAASH and 6S did not provide significantly different results (p > 0.05). The importance of the predictor variables (spectral bands and indices and topographic variables) for each of the MARS models fitted to the different radiometric correction algorithms was calculated (Figure 6). Results pointed out the near-infrared as the most relevant spectral region for AGB prediction. In fact, among the predictor variables, Band 7 made the greatest contribution to the capacity for predicting AGB for the different algorithms analyzed in almost all the cases except in algorithm 6S, where the variable Band 5 appeared as the most important in accordance with the number of times that each variable is included in a subset. Table 3 shows the results of applying the MARS technique for estimating the AGB on the basis of different correction algorithms. The algorithm that yielded the best predictions was ToA (R 2 = 0.89, RMSE = 29.82 Mg·ha −1 ; 33.49% of the average stand biomass). However, after the application of cross-validation, the generalized coefficient of determination (GCV) decreased to 0.68. By contrast, the algorithm with the lowest predictive power was ATCOR2 (R 2 = 0.75 and GR 2 = 0.59). The algorithm with the highest generalized coefficient of determination was 6S (GR 2 = 0.71). Figure 4. Spectral performance of each of the radiometric correction algorithms considered. As the Shapiro–Wilk’s test indicated that the corrected values of the spectral bands of the sensor were not normally distributed, therefore a log transformation was used. The ANOVA applied to the data from the 99 SPIFyS plots detected significant differences for all algorithms in the visible region (Figure 5), except for the COST and FLAASH algorithms in Band 3. Significant differences were observed in the near-infrared region for all algorithms tested, except ATCOR2 and 6S, which were included in the same population group (Tukey’s test); they also shared characteristics with FLAASH and ToA. Finally, in the bands corresponding to mid-infrared regions, the algorithms ATCOR2, FLAASH and 6S did not provide significantly different results (p> 0.05). The importance of the predictor variables (spectral bands and indices and topographic variables) for each of the MARS models fitted to the different radiometric correction algorithms was calculated (Figure 6). Results pointed out the near-infrared as the most relevant spectral region for AGB prediction. In fact, among the predictor variables, Band 7 made the greatest contribution to the capacity for predicting AGB for the different algorithms analyzed in almost all the cases except in algorithm 6S, where the variable Band 5 appeared as the most important in accordance with the number of times that each variable is included in a subset. Table 3shows the results of applying the MARS technique for estimating the AGB on the basis of different correction algorithms. The algorithm that yielded the best predictions was ToA (R 2 = 0.89, RMSE = 29.82 Mg ¨ ha ´1 ; 33.49% of the average stand biomass). However, after the application of cross-validation, the generalized coefficient of determination (GCV) decreased to 0.68. By contrast, the algorithm with the lowest predictive power was ATCOR2 (R 2 = 0.75 and GR 2 = 0.59). The algorithm with the highest generalized coefficient of determination was 6S (GR2= 0.71). Remote Sens. 2016,8, 369 16 of 19 13. Ediriweera, S.; Pathirana, S.; Danaher, T.; Nichols, D.; Moffiet, T. Evaluation of different topographic corrections for Landsat TM data by prediction of foliage projective cover (FPC) in topographically complex landscapes. Remote Sens. 2013,5, 6767–6789. [CrossRef] 14. Gagliasso, D.; Hummel, S.; Temesgen, H. A comparison of selected parametric and non-parametric imputation methods for estimating forest biomass and basal area. Open J. For. 2014,4, 42–48. [CrossRef] 15. Hantson, S.; Chuvieco, E. Evaluation of different topographic correction methods for Landsat imagery. Int. J. Appl. Earth Obs. Geoinform. 2011,13, 691–700. [CrossRef] 16. Balthazar, V.; Veerle, V.; Eric, F.L. Evaluation and parameterization of ATCOR3 topographic correction method for forest cover mapping in mountain areas. Int. J. Appl. Earth Obs. Geoinform. 2012 ,18, 436–450. [CrossRef] 17. Cortés, L.; Hernández, J.; Valencia, D.; Corvalán, P. Estimation of above-ground forest biomass using Landsat ETM+, Aster GDEM and Lidar. Forest Res. 2014. [CrossRef] 18. Kelsey, K.C.; Neff, J.C. Estimates of aboveground biomass from texture analysis of Landsat imagery. Remote Sens. 2014,6, 6407–6422. [CrossRef] 19. Tian, X.; Li, Z.; Su, Z.; Chen, E.; Van der Tol, C.; Li, X.; Guo, Y.; Li, L.; Ling, F. Estimating montane forest above-ground biomass in the upper reaches of the Heihe River Basin using Landsat-TM data. Int. J. Remote Sens. 2014,35, 7339–7362. [CrossRef] 20. Vincent, R.K. An ERTS multispectral scanner experiment for mapping iron compounds. In Proceedings of the 8th International Symposium on Remote Sensing of Environment, Ann Arbor, MI, USA, 2–6 October 1972; pp. 1239–1247. 21. Chavez, P.S. Image-based atmospheric corrections-revisited and improved. Photogramm. Eng. Remote Sens. 1996,62, 1025–1036. 22. Richter, R. A spatially adaptive fast atmospheric correction algorithm. Int. J. Remote Sens. 1996 ,17, 1201–1214. [CrossRef] 23. Adler-Golden, S.; Berk, A.; Bernstein, L.S.; Richtsmeier, S.; Acharya, P.K.; Matthew, M.W.; Anderson, G.P.; Allred, C.L.; Jeong, L.S.; Chetwynd, J.H. FLAASH, a MODTRAN4 atmospheric correction package for hyperspectral data retrievals and simulations. In Procedings of the 7th Annual JPL Airborne Earth Science Workshop, Pasadena, CA, USA, 12–16 January 1998; Green, R.O., Ed.; pp. 9–14. 24. Anderson, G.P.; Felde, G.W.; Hoke, M.L.; Ratkowski, A.J.; Cooley, T.W.; Chetwynd, J.H., Jr.; Gardner, J.A.; Adler-Golden, S.M.; Matthew, M.W.; Berk, A.; et al. MODTRAN4-based atmospheric correction algorithm: FLAASH (fast line-of-sight atmospheric analysis of spectral hypercubes). In Algorithms and Technologies for Multispectral, Hyperspectral, and Ultraspectral Imagery VIII (Proceedings of SPIE); Shen, S.S., Lewis, P.E., Eds.; Society of Photo Optics: Orlando, FL, USA, 2002; pp. 65–71. 25. Vermote, E.F.; Tanre, D.; Deuze, J.L.; Herman, M.; Morcette, J.J. Second simulation of the satellite signal in the solar spectrum, 6S: An overview. IEEE Trans. Geosci. Remote 1997,35, 675–686. [CrossRef] 26. Saatchi, S.S.; Harris, N.L.; Brown, S.; Lefsky, M.; Mitchard, E.T.A.; Salas, W.; Zutta, B.R.; Buermann, W.; Lewis, S.L.; Hagen, S.; et al. Benchmark map of forest carbon stocks in tropical regions across three continents. Proc. Natl. Acad. Sci. USA 2011,108, 9899–9904. [CrossRef] [PubMed] 27. Moisen, G.G.; Frescino, T.S. Comparing five modelling techniques for predict-ing forest characteristics. Ecol. Model. 2002,157, 209–225. [CrossRef] 28. Güneralp, I.; Filippi, A.M.; Randall, J. Estimation of floodplain aboveground biomass using multispectralremote sensing and nonparametric modeling. Int. J. Appl. Earth Obs. Geoinform. 2014 , 33, 119–126. [CrossRef] 29. Filippi, A.M.; Güneralp, I.; Randall, J. Hyperspectral remote sensing of above-ground biomass on a river meander bend using multivariate adaptive regressionsplines and stochastic gradient boosting. Remote Sens. Lett. 2014,5, 432–441. [CrossRef] 30. Instituto Nacional de Estadística Geográfica e Informática (INEGI). Uso del Suelo y Vegetación Escala 1:250 000 Serie V, Información Vectorial; Instituto Nacional de Estadística Geográfica e Informática: Ciudad de Mexico, Mexico, 2012. 31. Corral-Rivas, J.J.; Vargas, B.; Wehenkel, C.; Aguirre, O.; Álvarez, J.G.; Rojo, A. Guía para el Establecimiento de Sitios de Inventario Periódico Forestal y de Suelos del Estado de Durango; Facultad de Ciencias Forestales: Universidad Juárez del Estado de Durango, Durango, México, 2009. Remote Sens. 2016,8, 369 17 of 19 32. Vargas-Larreta, B. Estimación del Potencial de los Bosques de Durango Para la Mitigación del Cambio Climático. Modelización de la Biomasa Forestal. Proyecto FOMIX-DGO-2011-C01-165681; Comisión Nacional de Ciencia y Tecnología (CONACYT): Durango, Mexico, 2013. 33. United States Geological Survey. Available online: http://glovis.usgs.gov (accessed on 11 January 2011). 34. Rouse, J.W.; Haas, R.H.; Schell, J.A.; Deering, D.W.; Harlan, J.C. Monitoring the Vernal Advancements and Retrogradation (Greeen Wave Effect) of Natural Vegetation; NASA/GFSC: Greenbelt, MD, USA, 1974. 35. Baret, F.; Guyot, G. Potentials and limits of vegetation indices for LAI and APAR assessment. Remote Sens. Environ. 1991,35, 161–173. [CrossRef] 36. Lavreau, J. De-hazing landsat thematic mapper images. Photogramm. Eng. Remote Ssens. 1991 ,57, 1297–1302. 37. Guyot, G.; Gu, X.F. Effect of radiometric corrections on NDVI-determined from SPOT-HRV and Landsat-TM data. Remote Sens. Environ. 1994,49, 169–180. [CrossRef] 38. Richter, R. Atmospheric correction of satellite data with haze removal including a haze/clear transition region. Comput. Geosci. 1996,22, 675–681. [CrossRef] 39. Carlotto, M.J. Reducing the effects of space-varying, wavelength-dependent scattering in multispectral imagery. Int. J. Remote Sens. 1999,20, 3333–3344. [CrossRef] 40. Song, C.; Woodcock, C.E.; Seto, K.C.; Lenney, M.P.; Macomber, S.A. Classification and change detection using Landsat TM data: When and how to correct atmospheric effects? Remote Sens. Environ. 2001 ,75, 230–244. [CrossRef] 41. Liang, S.; Fang, H.; Morisette, J.T.; Chen, M.; Shuey, C.J.; Walthall, C.L.; Daughtry, C.S.T. Atmospheric correction of Landsat ETM+ land surface imagery. II. Validation and applications. IEEE Trans. Geosci. Remote Sens. 2002,40, 2736–2746. [CrossRef] 42. Chuvieco, E. Teledetección Ambiental: La Observación de la Tierra, 1st ed.; Ariel Ciencia: Madrid, Spain, 2010. 43. Geosystems. Haze Reduction, Atmospheric and Topographic Correction. User Manual ATCOR2 and ATCOR3; Geosystems GmbH: Geneva, Switzerland, 2013. 44. ERDAS Inc. Erdas Imagine. Available online: http://www.hexagongeospatial.com/products/ERDAS-IMAGINE/details.aspx (accessed on 6 October 2014). 45. Chavez, P.S. An improved dark-object subtraction technique for atmospheric scattering correction of multispectral data. Remote Sens. Environ. 1988,24, 459–479. [CrossRef] 46. Moran, M.S.; Jackson, R.D.; Slater, P.N.; Teillet, P.M. Evaluation of simplified procedures for retrieval of land surface reflectance factors from satellite sensor output. Remote Sens. Environ. 1992,41, 169–184. [CrossRef] 47. EXELIS Inc. ENVI ® 5.1 Fast Line-of-Sight Atmospheric Analysis of Hypercubes (FLAASH); EXELIS Inc.: Boulder, CO, USA, 2013. 48. Huang, C.; Yang, L.; Homer, C.; Wylie, B.; Vogelman, J.; DeFelice, T. At Satellite Reflectance: A First Order Normalization of LANDSAT 7 ETM+ Images. Available online: http://digitalcommons.unl.edu/usgspubs/109/ (accessed on 6 October 2014). 49. Goslee, S. Radiometric and Topographic Correction of Satellite Imagery. Package Landsat. Version 1.0.8. Available online: http://cran.r-project.org/web/packages/landsat/landsat.pdf (accessed on 18 May 2015). 50. R Core Team. R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing. Available online: http://www.R-project.org/ (accessed on 18 March 2014). 51. Instituto Nacional de Estadística Geográfica e Informática (INEGI). México Continuo de Elevaciones Mexicano 3.0 (CEM 3.0). Available online: http://www.inegi.org.mx/geo/contenidos/datosrelieve/continental/descarga.aspx (accessed on 6 May 2014). 52. McNab, W.H. Terrain shape index: quantifying effect of minor landforms on tree height. Forest Sci. 1989 ,35, 91–104. 53. Roberts, D.W.; Cooper, S.V. Concepts and techniques of vegetation mapping. In Land Classifications Based on Vegetation: Applications for Resource Management; Ferguson, D., Morgan, P., Johnson, F.D., Eds.; General technical report INT-US Department of Agriculture, Forest Service, Intermountain Research Station (USA): Ogden, UT, USA, 1989; pp. 90–96. 54. Moore, I.D.; Nieber, J.L. Landscape assessment of soil erosion and non-point source pollution. J. Minn. Acad. Sci. 1989,55, 18–24. Remote Sens. 2016,8, 369 18 of 19 55. Wilson, J.P.; Gallant, J.C. Digital terrain analysis. In Terrain Analysis: Principles and Applications; Wilson, J.P., Gallant, J.C., Eds.; John Wiley & Sons, Inc.: New York, NY, USA, 2000; pp. 1–28. 56. Zevenbergen, L.W.; Thorne, C.R. Quantitative analysis of land surface topography. Earth Surf. Proc. Landf. 1987,12, 47–56. [CrossRef] 57. Hijmans, R.J. Raster: Geographic Data Analysis and Modeling. R Package Version 2.4-20. Available online: http://CRAN.R-project.org/package=raster (accessed on 6 October 2015). 58. Shapiro, S.S.; Wilk, M.B. An analysis of variance test for normality (complete samples). Biometrika 1965 ,52, 591–611. [CrossRef] 59. Levene, H. Robust tests for equality of variances. In Contributions to Probability and Statistics: Essays in Honor of Harold Hotelling; Olkin, I., Ed.; Stanford University Press: Stanford, CA, USA, 1960; pp. 278–292. 60. Friedman, J.H. Multivariate adaptive regression spline. Ann. Stat. 1991,19, 1–141. [CrossRef] 61. Milborrow, S. Earth: Multivariate Adaptive Regression Splines. R Package Version 4.2.0. Available online: http://CRAN.R-project.org/package=earth (accessed on 6 October 2015). 62. Myers, R.H. Classical and Modern Regression with Applications; PWS-Kent Publishing Company: Belmont, CA, USA, 1990. 63. ESRI Inc. ArcGIS Desktop: Release 10; Environmental Systems Research Institute: Redlands, CA, USA, 2011. 64. Fernández, S.C.; Rodenas, A.Q. Teledetección: Avances y Aplicaciones; Diputación de Albacete: Albacete, Spain, 1999. 65. Buho, H.; Kaneko, M.; Ogawa, K. Correction of NDVI calculated from ASTER L1B and ASTER (AST07) data based on ground measurement. In Advances in Geoscience and Remote Sensing; Jedlovec, G., Ed.; In-Tech Press: Vienna, Austria, 2009. 66. Richter, R.; Schläpfer, D. Atmospheric/Topographic Correction for Satellite Imagery. ATCOR 2/3 User Guide, Version 9.0.0.; ReSe Applications Schläpfer: Langeggweg, Switzerland, 2015. 67. Broszeit, A.; Ashraf, S. Using different atmospheric correction methods to classify remotely sensed data to detect liquefaction of the February 2011 earthquake in Christchurch. In Proceedings of the SIRC NZ 2013—GIS and Remote Sensing Research Conference, Dunedin, New Zealand, 29–30 August 2013. 68. Nazeer, M.; Nichol, J.E.; Yung, Y.K. Evaluation of atmospheric correction models and Landsat surface reflectance product in an urban coastal environment. Int. J. Remote Sens. 2014,35, 6271–6291. [CrossRef] 69. Mahiny, A.S.; Turner, B.J. A comparison of four common atmospheric correction methods. Photogramm. Eng. Remote Sens. 2007,73, 361–368. [CrossRef] 70. Hamdan, O.; Hasmadi, I.M.; Aziz, H.K. Combination of SPOT-5 and ALOS PALSAR images in estimating aboveground biomass of lowland Dipterocarp forest. In Proceedings of the IOP Conference Series: Earth and Environmental Science, Tomst, Russia, 7–11 April 2014. 71. Nguyen, H.C.; Jung, J.; Lee, J.; Choi, S.U.; Hong, S.Y.; Heo, J. Optimal atmospheric correction for above-ground forest biomass estimation with the ETM+ remote sensor. Sensors 2015 ,15, 18865–18886. [CrossRef] [PubMed] 72. Bilgili, A.V.; Van Es, H.M.; Akbas, F.; Durak, A.; Hively, W.D. Visible-near infrared, reflectance spectroscopy for assessment of soil properties in a semi-arid area of Turkey. J. Arid Environ. 2010 ,74, 229–238. [CrossRef] 73. Ghasemi, J.B.; Zolfonoun, E. Application of principal component analysis-multivariate adaptive regression splines for the simultaneous spectrofluorimetric determination of dialkyltins in micellar media. Spectrochim. Spectrochimi. Acta Part A Mol. Biomol. Spectrosc. 2013,115, 357–363. [CrossRef] [PubMed] 74. García, M.A.; Pérez, C.F.; De la Riva, J.; Fernández, J.; Pascual, P.E.; Herranz, A. Estimación de la biomasa residual forestal mediante técnicas de teledetección y SIG en masas puras de Pinus halepensis yP. sylvestris. In Proceedings of the IV Congreso Forestal Español, Zaragoza, Spain, 26–30 September 2005. 75. López-Serrano, P.M.; López-Sánchez, C.A.; Díaz-Varela, R.A.; Corral-Rivas, J.J.; Solís-Moreno, R.; Vargas-Larreta, B.; Álvarez-González, J.G. Estimating biomass of mixed and uneven-aged forests using spectral data and a hybrid model combining regression trees and linear models. iForest 2015. [CrossRef] 76. Lam, C.S.Y. Comparison of Flow Routing Algorithms Used in Geographic Information Systems. Ph.D. Thesis, University of Southern California, Los Angeles, CA, USA, 2004. 77. Hoechstetter, S.; Walz, U.; Dang, L.H.; Thinh, N.X. Effects of topography and surface roughness in analyses of landscape structure–A proposal to modify the existing set of landscape metrics. Landsc. Online 2008 ,3, 1–14. [CrossRef] Remote Sens. 2016,8, 369 19 of 19 78. Sesnie, S.E.; Dickson, B.D.; Rosenstock, S.S.; Rundall, J.M. A comparison of Landsat TM and MODIS vegetation indices for estimating forage phenology in desert bighorn sheep (Ovis canadensis nelsoni) habitat in the Sonoran Desert, USA. Int. J. Remote Sens. 2012,33, 276–286. [CrossRef] 79. Ali, M.; Montzka, C.; Stadler, A.; Menz, G.; Thonfeld, F.; Vereecken, H. Estimation and validation of rapideye-based time-series of leaf area index for winter wheat in the rur catchment (Germany). Remote Sens. 2015,7, 2808–2831. [CrossRef] 80. Yao, Z.; Liu, J.; Zhao, X.; Long, D.; Wang, L. Spatial dynamics of aboveground carbon stock in urban green space: A case study of Xi’an, China. J. Arid Land. 2015,7, 350–360. [CrossRef] © 2016 by the authors; licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC-BY) license (http://creativecommons.org/licenses/by/4.0/).