Characterizing Forest Succession Stages for Wildlife Habitat Assessment Using Multispectral Airborne Imagery

In this study, we demonstrate the potential of using high spatial resolution airborne imagery to characterize the structural development stages of forest canopies. Four forest succession stages were adopted: stand initiation, young multistory, understory reinitiation, and old growth. Remote sensing metrics describing the spatial patterns of forest structures were derived and a Random Forest learning algorithm was used to classify forest succession stages. These metrics included texture variables from Gray Level Co-occurrence Measures (GLCM), range and sill from the semi-variogram, and the fraction of shadow and its spatial distribution. Among all the derived variables, shadow fractions and the GLCM variables of contrast, mean, and dissimilarity were the most important for characterizing the forest succession stages (classification accuracy of 89%). In addition, a LiDAR (Light Detection and Ranging) derived forest structural index (predicted Lorey’s height) was employed to validate the classification result. The classification using imagery spatial variables was shown to be consistent with the LiDAR derived variable (R2 = 0.68 and Root Mean Square Error (RMSE) = 2.39). This study demonstrates that high spatial resolution imagery was able to characterize forest succession stages with promising accuracy and may be considered an alternative to LiDAR data for this kind of application. Also, the results of stand development stages build a framework for future wildlife habitat mapping.


Introduction
Conventionally, the characterization of the physical habitat of wildlife is conducted by field work and generally is used at the local scale.For the past decades, remote sensing techniques have played a key role in mapping wildlife habitats, where large-scale spatial information is beneficial [1].Land cover classification maps derived from remote sensed data are commonly used to assess wildlife-habitat relationships and support wildlife conservation strategies and actions [2].Moderate spatial resolution satellite imagery (such as Landsat Thematic Mapper (TM) and Enhanced Thematic Mapper (ETM+) are widely available and have been used extensively for a range of species and habitat mapping applications [3,4].Land cover classification maps generated from Landsat imagery provide a broad scale delineation of different land cover types often including deciduous and coniferous forests [5].Applications of moderate resolution products to delineate wildlife habitat have revealed both benefits

Materials and Methods
The study area is located in the Hearst Forest near the town of Hearst, Central Ontario, Canada (49.7 • N, 83.7 • W) (Figure 1).Hearst Forest falls within the boreal mixed wood region and covers approximately 1.23 million ha, 1.00 million ha of which is productive forest.The main species within the forest are black spruce (Picea mariana), white spruce (Picea glauca), balsam fir (Abies balsamea), trembling aspen (Populous tremuloides), and black poplar (Populous balsamifera).
The optical imagery acquired using a Leica ADS-40 during July and August 2007 has four spectral bands (blue, green, red, and near-infrared) with a spatial resolution of 0.4 m by 0.4 m.The solar angle at the time of the data collection was 36.8 degrees.LiDAR data were also collected in between July and September 2007 during leaf-on conditions using an Optech ALS50 sensor mounted in a Cessna 310 aircraft with an average flying altitude of 2400 m above mean sea level.Average sampling density was 1.1 returns per square meter.The LiDAR data were used in this study to access the structural properties of forest canopies from different succession stages, and to allow for a comparison with those derived from optical imagery.Figure 1 shows the location of the study site and gives an overview of the multispectral imagery.Four general forest succession stages were adopted in this study [23]: stand initiation, young multistory, understory reinitiation, and old growth.The stand initiation stage indicates the stage where new vegetation becomes established and fully occupies a disturbed site.The following young multistory stage is mainly characterized by competition among the dominant trees.After some time when younger trees start to regenerate, the stand enters the understory reinitiation stage.In the last stage, the old growth stage, autogenic and gap-replacing processes create patches large enough for the stand reinitiation to begin.The optical imagery acquired using a Leica ADS-40 during July and August 2007 has four spectral bands (blue, green, red, and near-infrared) with a spatial resolution of 0.4 m by 0.4 m.The solar angle at the time of the data collection was 36.8 degrees.LiDAR data were also collected in between July and September 2007 during leaf-on conditions using an Optech ALS50 sensor mounted in a Cessna 310 aircraft with an average flying altitude of 2400 m above mean sea level.Average sampling density was 1.1 returns per square meter.The LiDAR data were used in this study to access the structural properties of forest canopies from different succession stages, and to allow for a comparison with those derived from optical imagery.Figure 1 shows the location of the study site and gives an overview of the multispectral imagery.
Four general forest succession stages were adopted in this study [23]: stand initiation, young multistory, understory reinitiation, and old growth.The stand initiation stage indicates the stage where new vegetation becomes established and fully occupies a disturbed site.The following young multistory stage is mainly characterized by competition among the dominant trees.After some time when younger trees start to regenerate, the stand enters the understory reinitiation stage.In the last stage, the old growth stage, autogenic and gap-replacing processes create patches large enough for the stand reinitiation to begin.
From the entire study area, 24 test areas representing four different succession stages were selected based on visual interpretation and the aid of tree height information from LiDAR.Each test site with a size of 501 × 501 pixels (200 m × 200 m) was subset.Examples of the test sites representing the four forest succession stages are shown in Figure 1a-d.
For each test site, to capture the spatial variation of forest canopies within a stand, calculations were conducted with a block size of 50 pixels by 50 pixels (i.e., 400 m 2 ).Such a block size was chosen as an adequate size to capture the variation within forest stands in each of the successional stages of this forest.To examine the ability of high spatial resolution imagery to differentiate forest stages, semi-variogram, GLCM, and the fraction of shadow and its spatial distribution were calculated through a self-written C++ program.
Semi-variograms measure the squared difference in brightness between pairs of pixels with a given distance apart and thus indicated the spatial variability of the observed tree canopies.The sill is the upper limit of the variance, which is equivalent to the overall variance in the image.The range is the distance at which the variance reaches the sill.The nugget effect is the variance at lag zero.From the semi-variogram, the variance values among pixels can be observed to increase with the lag distance between these pixels and become flattened out when the distance is large.The range was interpreted as the maximum crown size [24].The height of the sill of the variogram was related to the density of trees [25].To quantitatively relate the range and sill of a semi-variogram to the spatial variation of its corresponding stand, the experimental semi-variograms were calculated and then modeled using exponential model fitting, and range, sill, and nugget were determined.
To characterize forest canopies, the GLCM was used to calculate second order mean, standard deviation, contrast, Angular Second Moment (ASM), entropy, homogeneity, dissimilarity, and correlation [26][27][28].The first principal component from a PCA (principal component analysis) was used for the calculation of GLCM.For each pixel, four GLCMs were calculated within a neighborhood of 50 pixels by 50 pixels at four directions, 0 • , 45 • , 90 • , and 135 • .The displacement parameter in the GLCM calculation defines the scale at which the texture was analyzed.For our purpose, details of the spatial pattern were desired, and displacement distance was set equal to one pixel.The reason for selecting this value was to capture the variation within stands.When displacement value increases, the GLCM values become diverse and make it more difficult to distinguish between different stages.
In addition, shadows observed in imagery as a proxy for spatial attributes of forest canopies were investigated.In general, the distribution of shadows varies with tree height and size of the overstory, sun angle, and topography.As shown in Figure 2, the presence of shadows was different from stage to stage.Therefore, the shadow fraction that represents the shaded background was used as one descriptor to measure the structure difference of the forest stand development.In this study, shadow fraction was calculated as the ratio of shaded area to total area of the window size 50 × 50 pixels (20 m × 20 m).It was derived by using the near infrared band because the contrast between the tree stands and the shaded background was the strongest.First, a maximum-likelihood classification was made based on the distribution of the histogram.The threshold value was defined to separate the pixels into shaded and non-shaded groups.Then, the ratio between the number of shaded pixels and the total number of pixels within the window were calculated as the shadow fraction.
To further describe the forest succession stages in the measure of shadow distribution, a quadrant count (QC) method [29] was used.For a given site, the shadow maps from the previous step were used to calculate the degree of clumping to characterize the clumps of trees.Those maps were partitioned into m square regions called quadrats.The side length l of these quadrants was set to be 50 pixels to represent stands level for all sites.In each quadrat, the number of pixels that were marked as shadow was counted and the spatial pattern of the shadow for that site was described by a variance-to-mean ratio (VTMR) defined as: where x i is the number of pixels in each quadrat, and m is the number of quadrats in the layer.VTMR is a dispersion index, where the values equal to zero represent a constant random (not dispersed) case, and larger values represent greater clustering.The variables derived from semi-variogram, GLCM, shadow fraction, and degree of tree clumping were used in a Random Forest (RF) classification to investigate forest succession stage differentiation [30].The RF classifier consists of an ensemble of tree classifiers generated randomly from the input vector and outputs the classes that are voted by individual trees.It has been demonstrated in recent applications that a RF is an excellent classifier for hyperspectral remote sensing data, and multisource geographic data [31][32][33].In this study, RF was chosen due to the following two advantages, compared with other classification methods.First, RF classification was a non-parametric algorithm that was easy to interpret and free from dependence on the distribution assumption to achieve an accurate result.Second, each of the predictor variables was interpreted by calculating the influence of the model accuracy, which was used to determine the importance of each variable, and to select the optimal ones to further reduce the data redundancy [34].In this study, 50 trees were used in the classifier, as the analysis showed no significant difference between the forests after reaching that threshold.Analysis was conducted using the RandomForest package in the R statistical program [35].
Classification accuracy was assessed using a confusion matrix, where each column represented the predicted stage, and each row represented the actual stage based on visual interpretation.The class accuracy of each succession stage was calculated as the number of correctly predicted samples divided by the total number of samples in that stage, and the overall accuracy was calculated as the total number of correctly classified stages divided by the total number all tested samples.Also, to further validate the classification result, a LiDAR derived forest succession stage indicator was generated.The empirical model was built to link the LiDAR statistical derived metrics and the forest succession stage.Due to the lack of field surveys, no sufficient information was available to establish a direct relationship between the LiDAR point clouds and forest attributes at different stages.Thus, the model (Equation (2)) reported in van Ewijk et al. [12] to predict Lorey's height (PLH) for forest succession stage was adopted.Lorey's height (the stand tree height weighted by basal area) was found to be a good index in their study to separate the young multistory stage from an understory reinitiation stage.The reported model had an R 2 = 0.86, compared to the predicted quadratic mean diameter at breast height model which only achieved an R 2 = 0.68.PLH = 17.72 + 6.47 H_Skew + 0.94 H_Kurtosis + 0.80 H_ P9-16.59H_D6, ( where H_Skew is the Skewness of Heights, H_Kurtosisis the Kurtosis of Heights, H_ P9 is the 90th Height Percentile, and H_D6 is the 6th Height Decile.All returns of discrete LiDAR points were utilized for mapping the forest stage.First, the LiDAR data were classified by the vendor using the standard American Society for Photogrammetry and The variables derived from semi-variogram, GLCM, shadow fraction, and degree of tree clumping were used in a Random Forest (RF) classification to investigate forest succession stage differentiation [30].The RF classifier consists of an ensemble of tree classifiers generated randomly from the input vector and outputs the classes that are voted by individual trees.It has been demonstrated in recent applications that a RF is an excellent classifier for hyperspectral remote sensing data, and multisource geographic data [31][32][33].In this study, RF was chosen due to the following two advantages, compared with other classification methods.First, RF classification was a non-parametric algorithm that was easy to interpret and free from dependence on the distribution assumption to achieve an accurate result.Second, each of the predictor variables was interpreted by calculating the influence of the model accuracy, which was used to determine the importance of each variable, and to select the optimal ones to further reduce the data redundancy [34].In this study, 50 trees were used in the classifier, as the analysis showed no significant difference between the forests after reaching that threshold.Analysis was conducted using the RandomForest package in the R statistical program [35].
Classification accuracy was assessed using a confusion matrix, where each column represented the predicted stage, and each row represented the actual stage based on visual interpretation.The class accuracy of each succession stage was calculated as the number of correctly predicted samples divided by the total number of samples in that stage, and the overall accuracy was calculated as the total number of correctly classified stages divided by the total number all tested samples.Also, to further validate the classification result, a LiDAR derived forest succession stage indicator was generated.The empirical model was built to link the LiDAR statistical derived metrics and the forest succession stage.Due to the lack of field surveys, no sufficient information was available to establish a direct relationship between the LiDAR point clouds and forest attributes at different stages.Thus, the model (Equation ( 2)) reported in van Ewijk et al. [12] to predict Lorey's height (PLH) for forest succession stage was adopted.Lorey's height (the stand tree height weighted by basal area) was found to be a good index in their study to separate the young multistory stage from an understory reinitiation stage.The reported model had an R 2 = 0.86, compared to the predicted quadratic mean diameter at breast height model which only achieved an R 2 = 0.68.PLH = 17.72 + 6.47 H_Skew + 0.94 H_Kurtosis + 0.80 H_ P9-16.59H_D6, where H_Skew is the Skewness of Heights, H_Kurtosis is the Kurtosis of Heights, H_ P9 is the 90th Height Percentile, and H_D6 is the 6th Height Decile.All returns of discrete LiDAR points were utilized for mapping the forest stage.First, the LiDAR data were classified by the vendor using the standard American Society for Photogrammetry and Remote Sensing classification and Terra Solid's TerraScan software.Classified ground returns were used to generate a 1-m bare earth digital terrain model (DTM) raster, against which LiDAR vegetation returns were normalized.For stand level analysis, the entire area of 4 km by 4 km was arbitrarily cut into height bins in size of 20 m by 20 m, and the PLH within each of the height bins was calculated.Linear regression was modeled to determine if used variables and PLH were correlated with each other.

Semi-Variogram
The semi-variograms for the four stages revealed that the range of the variance for forest stands at the initiation stage was the smallest, and the range of the old-growth forest stands was the greatest (Figure 3).These values were close to the sizes of the dominant trees on the sites.The sill of the variance was related to the variation in the brightness value of the pixels, which in turn were related to percent canopy cover and vertical layering in the stands [36].As expected, stand initiation and young multistory stands have single canopy layers, and the corresponding sill values are low.The understory reinitiation and old-growth stands have multiple canopy layers with higher contrast and higher resulting sill values.Also, multiple peaks exhibit variance of forest stands at stand initiation and old growth stages.The results indicate gaps or clumps of trees in the stands at stand initiation stage, and mixed clumping pattern of trees in the stands at old-growth stage.

Semi-Variogram
The semi-variograms for the four stages revealed that the range of the variance for forest stands at the initiation stage was the smallest, and the range of the old-growth forest stands was the greatest (Figure 3).These values were close to the sizes of the dominant trees on the sites.The sill of the variance was related to the variation in the brightness value of the pixels, which in turn were related to percent canopy cover and vertical layering in the stands [36].As expected, stand initiation and young multistory stands have single canopy layers, and the corresponding sill values are low.The understory reinitiation and old-growth stands have multiple canopy layers with higher contrast and higher resulting sill values.Also, multiple peaks exhibit variance of forest stands at stand initiation and old growth stages.The results indicate gaps or clumps of trees in the stands at stand initiation stage, and mixed clumping pattern of trees in the stands at old-growth stage.The range value decreases from stands at the first stage (stand initiation) to the second (young multistory) and then increases to the third and fourth stage (understory reinitiation and old growth, respectively) (Table 1).At the early stage, when most of the area was bare soil, its spatial variation is rough and can be recognized as objects with an irregular shape in various sizes.From the young multistory stage to the old growth stage, the tree crowns are dominant.Therefore, the range value increases as the size of the crowns increase from young to old stages.This is validated by comparing the standard deviation of range value.In the second and third stages, the tree crown surface is more homogeneous, resulting in a smaller standard deviation.In contrast, both stand initiation and old growth stages had higher standard deviations, which may relate to greater canopy openness and associated variability.The range value decreases from stands at the first stage (stand initiation) to the second (young multistory) and then increases to the third and fourth stage (understory reinitiation and old growth, respectively) (Table 1).At the early stage, when most of the area was bare soil, its spatial variation is rough and can be recognized as objects with an irregular shape in various sizes.From the young multistory stage to the old growth stage, the tree crowns are dominant.Therefore, the range value increases as the size of the crowns increase from young to old stages.This is validated by comparing the standard deviation of range value.In the second and third stages, the tree crown surface is more homogeneous, resulting in a smaller standard deviation.In contrast, both stand initiation and old growth stages had higher standard deviations, which may relate to greater canopy openness and associated variability.To further demonstrate the relationship of range and sill for different stages, a scatter plot was created, with each stage assigned a different color.As shown in Figure 4, the stands at the stand initiation stage have a significantly smaller sill value, which can be used to distinguish it.Also, though there is some overlapping, the majority of young multistory and understory reinitiation values were separate from each other, while the old growth was mixed with them.initiation stage have a significantly smaller sill value, which can be used to distinguish it.Also, though there is some overlapping, the majority of young multistory and understory reinitiation values were separate from each other, while the old growth was mixed with them.

GLCM
Findings from the GLCM texture variable analysis indicated that no single variable can be used to separate all four successional stages; however, contrast showed the largest differences (Figure 5).Higher contrast values represent greater spatial variation between neighboring pixels.At the young multistory stage, the tree crowns were relatively small, and the neighboring pixels of different features (such as crown and soil) occur more frequently, leading to greater contrast.On the other hand, old growth stands have bigger crowns, and the neighboring pixels were more likely to be within the same features (crown or understory), and thus the contrast appeared to be the smallest.

GLCM
Findings from the GLCM texture variable analysis indicated that no single variable can be used to separate all four successional stages; however, contrast showed the largest differences (Figure 5).Higher contrast values represent greater spatial variation between neighboring pixels.At the young multistory stage, the tree crowns were relatively small, and the neighboring pixels of different features (such as crown and soil) occur more frequently, leading to greater contrast.On the other hand, old growth stands have bigger crowns, and the neighboring pixels were more likely to be within the same features (crown or understory), and thus the contrast appeared to be the smallest.
to separate all four successional stages; however, contrast showed the largest differences (Figure 5).Higher contrast values represent greater spatial variation between neighboring pixels.At the young multistory stage, the tree crowns were relatively small, and the neighboring pixels of different features (such as crown and soil) occur more frequently, leading to greater contrast.On the other hand, old growth stands have bigger crowns, and the neighboring pixels were more likely to be within the same features (crown or understory), and thus the contrast appeared to be the smallest.

Shadow Characteristics
Shadow fractions varied among the forest stages but the differences were not significant between adjacent stages (Table 2).6).The results demonstrated that shadow distribution could be used as an indicator for classifying between these three stages.Old growth had the highest VTMR and thus most clustered overstory, and young multistory had the lowest VTMR and most homogenous tree distribution.At the young multistory stand stage, most trees in early growth had crowns that were isolated, resulting in dispersed shadow areas.At old growth stand stage, big crowns from the mature stands and the multistory created clustered shadows.Bare soils were dominated in stands at stand initiation and caused shadow to be randomly distributed.
clustered overstory, and young multistory had the lowest VTMR and most homogenous tree distribution.At the young multistory stand stage, most trees in early growth had crowns that were isolated, resulting in dispersed shadow areas.At old growth stand stage, big crowns from the mature stands and the multistory created clustered shadows.Bare soils were dominated in stands at stand initiation and caused shadow to be randomly distributed.

Classification of Forest Succession Stage
In total, 32 variables as eight texture features per each direction, together with shadow fraction were used to perform RF classification.In total 1424 out of 1600 samples were correctly classified to achieve an accuracy of 89% (Table 3).

Classification of Forest Succession Stage
In total, 32 variables as eight texture features per each direction, together with shadow fraction were used to perform RF classification.In total 1424 out of 1600 samples were correctly classified to achieve an accuracy of 89% (Table 3).The forest successional stages derived from the image spatial pattern analysis agreed with the PLH from LiDAR data (Figure 7).All the areas with high PLH (brown-shaded area in Figure 7b) were classified as stands at old growth stage (red color areas in Figure 8a).Most of the forest stands identified from images as reinitiation stage areas identified from images (shaded green in Figure 7b) corresponded well to low PLH.The same can be noted for the stands at the young multistory stage too.Only the stands identified from the images as stands at an understory reinitiation stage appeared to be under classified, potentially related to its transitional nature between multistory and old growth.Also, the image texture of the stands at this stage was very similar to the stands at an old-growth stage.As evident from the confusion matrix of the RF classification (Table 3), most of the commission error of stands at understory reinitiation stage was due to these stands being classified as old-growth.
too.Only the stands identified from the images as stands at an understory reinitiation stage appeared to be under classified, potentially related to its transitional nature between multistory and old growth.Also, the image texture of the stands at this stage was very similar to the stands at an oldgrowth stage.As evident from the confusion matrix of the RF classification (Table 3), most of the commission error of stands at understory reinitiation stage was due to these stands being classified as old-growth.

GLCM Parameter Sensitivity Analysis
The effects of window size and displacement on GLCM derived texture variables were assessed using classification accuracies.The size of the processing window was a trade-off between finding enough spatial information to characterize the land features and limiting overlapping textures of different land features.A large window size captured a better representation of spatial patterns of each stage, but also contained patterns from other stages and introduced systematic errors.For our purpose, the average size of the sampling unit that described the within stand variation was used to accommodate a window size of 50 pixels by 50 pixels.To test if this value was the optimal value, GLCM texture variables were calculated with different window sizes and tested in the RF classification process.As shown in Figure 8a, the classification accuracy increased as the processing window size increased.When window size was greater than 40 pixels, the improvement of the classification accuracy started to diminish.It is worth mentioning that a standard window size in this study provides promising results, further analysis should include a segmentation process to extract canopy objects.The texture measurements should be processed within each object rather than predefined windows.The canopy openness, size of the crown/stands were all different at different stand development stages of the forest.At the early stage of the forest dynamics, the trees are generally more isolated, while in their mature stage, when they grow into each other, the overlaps of the tree crown made multiple crowns clustered as one single object.Thus, one window size cannot fit all forest stages.Also, to test if the displacement equal to one pixel was the optimal value, texture variables were calculated with different displacement distances and in the RF classification process.As shown in Figure 8b, the classification accuracy decreased with the increase of the displacement distance.When displacement equaled 11 pixels, a slight increase in the classification accuracy was observed, indicating a strong spatial relationship between two pixels that were 10 pixels apart, in our case meaning the dominant tree/tree clump size.

Feature Selection for Classification Variables Selection
Although promising classification accuracy was achieved, the number of variables used in the classification was large.To reduce the redundancy of the texture features, the omni-direction approach was tested.Different from calculating texture features from GLCM at one direction, the co- From a visual comparison, forest successional stages generated from image texture provide a similar classification map as one created from PLH generated from LiDAR point clouds.A linear regression fit to the suite of variables and the LiDAR derived PLH resulted in an R 2 = 0.68 and RMSE = 2.39.It should be noted that the R 2 was 0.86 for the PLH and field data for the Petawawa research forest.Differences in forest type exist between our study area and Petawawa research forest, and also between the real field data and regression model.Overall, the texture features used for classification were correlated with the LiDAR derived PLH.

GLCM Parameter Sensitivity Analysis
The effects of window size and displacement on GLCM derived texture variables were assessed using classification accuracies.The size of the processing window was a trade-off between finding enough spatial information to characterize the land features and limiting overlapping textures of different land features.A large window size captured a better representation of spatial patterns of each stage, but also contained patterns from other stages and introduced systematic errors.For our purpose, the average size of the sampling unit that described the within stand variation was used to accommodate a window size of 50 pixels by 50 pixels.To test if this value was the optimal value, GLCM texture variables were calculated with different window sizes and tested in the RF classification process.As shown in Figure 8a, the classification accuracy increased as the processing window size increased.When window size was greater than 40 pixels, the improvement of the classification accuracy started to diminish.
It is worth mentioning that a standard window size in this study provides promising results, further analysis should include a segmentation process to extract canopy objects.The texture measurements should be processed within each object rather than predefined windows.The canopy openness, size of the crown/stands were all different at different stand development stages of the forest.At the early stage of the forest dynamics, the trees are generally more isolated, while in their mature stage, when they grow into each other, the overlaps of the tree crown made multiple crowns clustered as one single object.Thus, one window size cannot fit all forest stages.Also, to test if the displacement equal to one pixel was the optimal value, texture variables were calculated with different displacement distances and in the RF classification process.As shown in Figure 8b, the classification accuracy decreased with the increase of the displacement distance.When displacement equaled 11 pixels, a slight increase in the classification accuracy was observed, indicating a strong spatial relationship between two pixels that were 10 pixels apart, in our case meaning the dominant tree/tree clump size.

Feature Selection for Classification Variables Selection
Although promising classification accuracy was achieved, the number of variables used in the classification was large.To reduce the redundancy of the texture features, the omni-direction approach was tested.Different from calculating texture features from GLCM at one direction, the co-occurrence matrix was generated in each direction and then added together for processing.In this case, the texture variables used for RF classification were lowered from 32 to 8. As a result, an accuracy of 82% was achieved.A confusion matrix of the classification result of both the four direction and omni-direction is presented in Table 4.As demonstrated in Table 4, sampling in all directions can significantly reduce the number of variables used in the classification while maintaining similar classification accuracy.Classification accuracy for young multistory stage conditions was similar between the two classifications but was lower for the other stages, potentially due to the solar angle causing the texture variation associated with shadows.The amount of the shadow varied when the tree height, the presence of understory, and size of the canopies were different.Stands at the young multistory stage had the most homogenous texture.Thus, the shadows were uniform and insignificant from GLCM in different directions.On the other hand, the shadows in the stands at the stand initiation and understory reinitiation stages were less uniform and had a direction associated with the solar angle in their spatial distribution.Therefore, texture features calculated from GLCM at the direction close to the direction of the shadow were quite different from those of the others.However, putting all directions together will minimize this difference and lower the classification accuracy.
To further reduce the variables used for classification, the first four texture features from the RF variable importance plot were selected to perform the classification.For our study, the most important variables are contrast, variance, mean, and dissimilarity (Figure 9).Both contrast and dissimilarity measured the spreading degree of the gray levels or the average gray level difference between neighboring pixels.The contrast values were higher for large local variations within the region.The local window grey-level variance was higher when there was a significant grey-level standard deviation in the local region.When forests develop, and multistory structures grow, the layers of tree crowns and understory reflect as the gray level variance in the texture.To test if such information is sufficient enough to differentiate the forest from young to mature, a RF classification was performed with an accuracy of 82%.
with shadows.The amount of the shadow varied when the tree height, the presence of understory, and size of the canopies were different.Stands at the young multistory stage had the most homogenous texture.Thus, the shadows were uniform and insignificant from GLCM in different directions.On the other hand, the shadows in the stands at the stand initiation and understory reinitiation stages were less uniform and had a direction associated with the solar angle in their spatial distribution.Therefore, texture features calculated from GLCM at the direction close to the direction of the shadow were quite different from those of the others.However, putting all directions together will minimize this difference and lower the classification accuracy.
To further reduce the variables used for classification, the first four texture features from the RF variable importance plot were selected to perform the classification.For our study, the most important variables are contrast, variance, mean, and dissimilarity (Figure 9).Both contrast and dissimilarity measured the spreading degree of the gray levels or the average gray level difference between neighboring pixels.The contrast values were higher for large local variations within the region.The local window grey-level variance was higher when there was a significant grey-level standard deviation in the local region.When forests develop, and multistory structures grow, the layers of tree crowns and understory reflect as the gray level variance in the texture.To test if such information is sufficient enough to differentiate the forest from young to mature, a RF classification was performed with an accuracy of 82%.To analyze the effect of semi-variogram variables range and sill on the overall classification, all variables derived from the GLCM and semi-variogram variables range and sill were used in a RF classification.This resulted in a slight classification accuracy improvement of 0.69%.However, among all the variables, sill was ranked the sixth most important variable and range was the ninth most important variable in the classification (Figure 10).Overall, the semi-variogram variables did not increase the classification significantly.Calculating range and sill values takes more effort than calculating GLCM variables.While range and sill values were useful for separating the stands at stand initiation stage from the rest, they were not worth calculating for the overall forest succession stage classification.
Adding shadow fractions, together with the variables derived from the GLCM, resulted in improved classification accuracy in a RF classification (89%), and shadow fraction was ranked the most important variable (Figure 11).Specifically, adding the shadow fraction improved the differentiation of the stands at the old growth stage from stands at the understory reinitiation stage, and the stands at the stand initiation stage from stands at the young multistory stage.However, this did not improve the classification between stands at the young multistory stage and understory reinitiation stage.among all the variables, sill was ranked the sixth most important variable and range was the ninth most important variable in the classification (Figure 10).Overall, the semi-variogram variables did not increase the classification significantly.Calculating range and sill values takes more effort than calculating GLCM variables.While range and sill values were useful for separating the stands at stand initiation stage from the rest, they were not worth calculating for the overall forest succession stage classification.Adding shadow fractions, together with the variables derived from the GLCM, resulted in improved classification accuracy in a RF classification (89%), and shadow fraction was ranked the most important variable (Figure 11).Specifically, adding the shadow fraction improved the differentiation of the stands at the old growth stage from stands at the understory reinitiation stage, and the stands at the stand initiation stage from stands at the young multistory stage.However, this did not improve the classification between stands at the young multistory stage and understory reinitiation stage.Adding shadow fractions, together with the variables derived from the GLCM, resulted in improved classification accuracy in a RF classification (89%), and shadow fraction was ranked the most important variable (Figure 11).Specifically, adding the shadow fraction improved the differentiation of the stands at the old growth stage from stands at the understory reinitiation stage, and the stands at the stand initiation stage from stands at the young multistory stage.However, this did not improve the classification between stands at the young multistory stage and understory reinitiation stage.RF classification was performed with the variables derived from the GLCM and VTMR derived from tree clumping.To provide spatial consistency, all metrics were calculated within a 50 pixels by 50 pixels window.Overall, adding the VTMR decreased the classification to 78%.Specifically, adding the VTMR improved the correctly classified stands at the stand initiation stage and old growth stage, but greatly affected differentiating the stands at the understory reinitiation stage from the other stages (from 72% down to 59%).It is postulated that this result occurred because the shadow distribution metrics best describe the scene of a large area.When calculated within a sampling unit (400 m 2 ), it was a poor indicator to separate the different stages.Since the VTMR failed to increase the overall classification accuracy, the variable was not included in the classification.
Finally, the classification scalability was tested at the coarser resolutions.To compare the classification accuracy quantitatively, a set of lower resolution images were resampled from the original image.As shown in Figure 12, the classification accuracy decreased as the image resolution decreased.This can be explained by the fact that in very high-resolution imagery, detailed information is revealed from the texture used to aid in distinguishing different classes.With coarser resolution imagery Forests 2017, 8, 234 14 of 16 information content is smoothed and made less distinct which cannot provide sufficient information to separate classes.However, it should be noted that high-resolution image computation was computer processing intensive and expensive for large area data collection, while coarser resolution satellite imagery is more readily available.Based on our findings, a resolution of 2 m was sufficient for large area habitat mapping, as the image size was smaller while the classification accuracy remains above 75%.
the overall classification accuracy, the variable was not included in the classification.
Finally, the classification scalability was tested at the coarser resolutions.To compare the classification accuracy quantitatively, a set of lower resolution images were resampled from the original image.As shown in Figure 12, the classification accuracy decreased as the image resolution decreased.This can be explained by the fact that in very high-resolution imagery, detailed information is revealed from the texture used to aid in distinguishing different classes.With coarser resolution imagery information content is smoothed and made less distinct which cannot provide sufficient information to separate classes.However, it should be noted that high-resolution image computation was computer processing intensive and expensive for large area data collection, while coarser resolution satellite imagery is more readily available.Based on our findings, a resolution of 2 m was sufficient for large area habitat mapping, as the image size was smaller while the classification accuracy remains above 75%.This study demonstrates that high spatial resolution imagery was able to characterize forest succession stages with promising accuracy and may be considered an alternative to LiDAR data with R 2 = 0.68 and RMSE = 2.39.We were unable to directly compare our results with the related literature because we did not find any work investigating the utilization of imagery texture features for forest succession stage classification.Also, the relevant studies assessing semi-variogram, GLCM, and tree shadows used satellite data with different spectral and spatial resolutions, and a different statistical approach that correlates the texture features and structural parameters.However, it is worth mentioning that our study supports the findings of the research conducted by Kayitakire et al. [19].In that study, forest structural parameters including basal area, top height, circumference, stand density, and age were modeled using image texture features derived from 1 m resolution IKONOS-2 imagery.The R 2 values of the best models ranged from 0.35 to 0.82 for these variables.Also, This study demonstrates that high spatial resolution imagery was able to characterize forest succession stages with promising accuracy and may be considered an alternative to LiDAR data with R 2 = 0.68 and RMSE = 2.39.We were unable to directly compare our results with the related literature because we did not find any work investigating the utilization of imagery texture features for forest succession stage classification.Also, the relevant studies assessing semi-variogram, GLCM, and tree shadows used satellite data with different spectral and spatial resolutions, and a different statistical approach that correlates the texture features and structural parameters.However, it is worth mentioning that our study supports the findings of the research conducted by Kayitakire et al. [19].In that study, forest structural parameters including basal area, top height, circumference, stand density, and age were modeled using image texture features derived from 1 m resolution IKONOS-2 imagery.The R 2 values of the best models ranged from 0.35 to 0.82 for these variables.Also, important GLCM features found are similar to what was demonstrated in Ozdemir and Karnieli's study [20].Their work examined forest structural parameters including number of trees, basal area, stem volume, diameter differentiation index, contagion index, gini coefficient, and standard deviation of diameters at breast height using correlation analyses with image texture features extracted from 2 m resolution WorldView-2 image.The results of the correlation analysis indicate that the forest structural parameters are significantly correlated with the contrast of red band (r = 0.75, p < 0.01), the entropy of blue band (r = 0.73, p < 0.01) and the contrast of blue band (r = 0.71, p < 0.01).

Conclusions
The study provided a comprehensive investigation of variables derived from texture analysis of GLCM, semi-variogram, and spectral variables of shadow fraction and shadow distribution.The result of RF classification demonstrated that using GLCM texture measurements and shadow fraction can differentiate stages of forest stand development with a classification accuracy of 89%.The key variables for classification were contrast, variance, mean, and dissimilarity from GLCM and shadow fraction.The development of forest stands and the change of the canopy surface created variation between canopies that was measurable from the spatial pattern analysis.Our investigation found no linear relationship between the semi-variogram parameters and forest succession stages.However, the sill value can be used to distinguish the stands at the stand initiation stage, and the range value can be used to separate stands at young multistory, understory reinitiation, and old growth stages.The cross-validation of image texture classification using LiDAR derived indices of PLH showed that the classification from young to mature corresponded to the PLH, providing an indication of agreement for classified forest stage with the LiDAR derived indices.Therefore, it is concluded that spatial pattern analysis of high spatial resolution imagery can be used as an alternative of LiDAR to derive forest successional stages.
In future work, we will investigate an index for characterizing dynamic forest development.This research used discrete stages to characterize forest development, where the continuous descriptor was more favorable.Also, within-stand variability that may exist within a stand development class will have individual stands deviated from the identified stages.The value of this research is in its contribution toward building a stand development framework that recognizes specific stand structural conditions as benchmarks for future wildlife habitat mapping.

Figure 1 .
Figure 1.False-color imagery of the study site, with the near infrared band printed as red, red as green, and green as blue.Test sites are marked with gray squares.The location of study area in Hearst Forest, Ontario, Canada is at lower left panel.The imagery of forest stands at the four development stages are as follows: (a) stand initiation; (b) young multistory; (c) understory reinitiation; and (d) oldgrowth stage.The false color composite of the optical imagery covering the study area is displayed with the near-infrared band as red, red as green, and green as blue.

Figure 1 .
Figure 1.False-color imagery of the study site, with the near infrared band printed as red, red as green, and green as blue.Test sites are marked with gray squares.The location of study area in Hearst Forest, Ontario, Canada is at lower left panel.The imagery of forest stands at the four development stages are as follows: (a) stand initiation; (b) young multistory; (c) understory reinitiation; and (d) old-growth stage.The false color composite of the optical imagery covering the study area is displayed with the near-infrared band as red, red as green, and green as blue.

Figure 2 .
Figure 2. Shadow map for forest stands at (a) stand initiation stage; (b) young multistory stage; (c) understory reinitiation stage; and (d) old growth stage, where black color indicates shaded background and white indicate sunlit background.

Figure 2 .
Figure 2. Shadow map for forest stands at (a) stand initiation stage; (b) young multistory stage; (c) understory reinitiation stage; and (d) old growth stage, where black color indicates shaded background and white indicate sunlit background.

Figure 3 .
Figure 3. Semi-variogram of the four forest succession stages: forest stands at stand initiation, forest stands at young multistory, forest stands at understory reinitiation, and forest stands at the oldgrowth stage.

Figure 3 .
Figure 3. Semi-variogram of the four forest succession stages: forest stands at stand initiation, forest stands at young multistory, forest stands at understory reinitiation, and forest stands at the old-growth stage.

Figure 4 .
Figure 4. Scatter plot of semi-variogram parameters range and sill for four forest succession stages

Figure 4 .
Figure 4. Scatter plot of semi-variogram parameters range and sill for four forest succession stages

Figure 5 .
Figure 5. Bubble plot of GLCM variables: second order mean, variance, contrast, Angular Second Moment (ASM), entropy, homogeneity, dissimilarity, and correlation for forest stands at four stages, where the location of the bubbles indicates the mean value of that stage and the size of the bubble represents the standard deviations.

Figure 5 .
Figure 5. Bubble plot of GLCM variables: second order mean, variance, contrast, Angular Second Moment (ASM), entropy, homogeneity, dissimilarity, and correlation for forest stands at four stages, where the location of the bubbles indicates the mean value of that stage and the size of the bubble represents the standard deviations.

Figure 6 .
Figure 6.Histogram of Variance-To-Mean Ratio (VTMR) of each quadrat from all four forest succession stages.

Figure 6 .
Figure 6.Histogram of Variance-To-Mean Ratio (VTMR) of each quadrat from all four forest succession stages.

Figure 7 .
Figure 7. (a) Color shaded forest succession stages predicted from multispectral image spatial patterns with brown, blue, green, and red representing the forest stands at (1) stand initiation, (2) young multistory, (3) understory reinitiation, and (4) old growth stage, respectively.(b) color shaded map of Predicted Lorey's Height (PLH) derived from Light Detection and Ranging (LiDAR) data with color from blue to red indicating PLH from low to high.

Figure 7 .
Figure 7. (a) Color shaded forest succession stages predicted from multispectral image spatial patterns with brown, blue, green, and red representing the forest stands at (1) stand initiation, (2) young multistory, (3) understory reinitiation, and (4) old growth stage, respectively.(b) color shaded map of Predicted Lorey's Height (PLH) derived from Light Detection and Ranging (LiDAR) data with color from blue to red indicating PLH from low to high.

Figure 8 .
Figure 8. Overall RF classification accuracy results of changing Gray Level Co-occurrence Measures (GLCM) parameters.(a) Classification accuracy versus window sizes.(b) Classification accuracy versus displacements.

Figure 8 .
Figure 8. Overall RF classification accuracy results of changing Gray Level Co-occurrence Measures (GLCM) parameters.(a) Classification accuracy versus window sizes.(b) Classification accuracy versus displacements.

Figure 9 .Figure 9 .
Figure 9. Variable importance plots for the forest succession stage classification using GLCM variables.Units are the percentage reduction of resulting classification accuracy from removing one given texture feature.

Figure 10 .
Figure 10.Variable importance plots for the forest succession stage classification using both GLCM variables and semi-variogram variables range and sill values.Units are the percentage reduction of resulting classification accuracy from removing one given texture feature.

Figure 10 .
Figure 10.Variable importance plots for the forest succession stage classification using both GLCM variables and semi-variogram variables range and sill values.Units are the percentage reduction of resulting classification accuracy from removing one given texture feature.

Figure 10 .
Figure 10.Variable importance plots for the forest succession stage classification using both GLCM variables and semi-variogram variables range and sill values.Units are the percentage reduction of resulting classification accuracy from removing one given texture feature.

Figure 11 .
Figure 11.Variable importance plots for the forest succession stage classification using GLCM variables and shadow fractions.Units are the percentage reduction of resulting classification accuracy from removing one given texture feature.

Figure 12 .
Figure 12.Overall RF classification accuracy results with image resolution changing from fine to coarse.

Figure 12 .
Figure 12.Overall RF classification accuracy results with image resolution changing from fine to coarse.

Table 1 .
Average and standard deviation of semi-variogram parameters.

Table 2 .
Shadow fraction derived from Near-InfraRed (NIR) band from four forest succession stage.differed significantly among young multistory, stand initiation, and old growth stages, as the mean values of the VTMR are 44.63,95.93, and 175.84 with stand deviations of 15.16, 37.59, and 64.96, respectively (Figure VTMR

Table 3 .
Confusion matrix from RF classification.

Table 4 .
Confusion matrix from RF classification.