Hierarchical Segmentation Framework for Identifying Natural Vegetation: A Case Study of the Tehachapi Mountains, California

Abstract

Two critical limitations of very high resolution imagery interpretations for time-series analysis are higher imagery variances and large data sizes. Although object-based analyses with a multi-scale framework for diverse object sizes are one potential solution, more data requirements and large amounts of testing at high costs are required. In this study, I applied a three-level hierarchical vegetation framework for reducing those costs, and a three-step procedure was used to evaluate its effects on a digital orthophoto quadrangles with 1 m spatial resolution. Step one and step two were for image segmentation optimized for delineation of tree density, which involved global Otsu’s method followed by the random walker algorithm. Step three was for detailed species delineations, which were derived from multiresolution segmentation, in two test areas. Step one and step two were able to delineating tree density segments and label species association robustly, compared to previous hierarchical frameworks. However, step three was limited by less image information to produce detailed, reasonable image objects with optimal scale parameters for species labeling. This hierarchical vegetation framework has potential to develop baseline data for evaluating climate change impacts on vegetation at lower cost using widely available data and a personal laptop.

1. Introduction

1.1. Motivation

Vegetation mapping is required for biological conservation and forest inventory; especially time series mapping is commonly used to detect transformations of species or suitable habitats (e.g., [1–3]) and evaluate impacts of climate change on species (e.g., [4,5]). Species distribution models (SDMs, also called habitat suitability models), which correlate environment variables with species sampling data to map species occurrences [6], have been used for vegetation mapping. However, the field-based approach is cost and labor-consuming, while pure remote sensing imagery interpretations are lacking in ability to detect individual species [7]. SDMs are viewed as static and equilibrium models capturing species-environment relations at large scales, but ignoring dynamic biological interactions, such as dispersal, migration, facilitation, competition, mutualism, and predation at local scales [8]. Mapping species directly using very high resolution imagery (1 to 2 m), and tracking vegetation transformations over time through repeated mapping provided the potential for building more dynamic vegetation mapping.

With advances of sensors, spatial resolutions of remotely sensed data have been improved to centimeter levels and provide more spatial information on species distributions. For example, WorldView-1, which was launched on 18 September 2007, collects panchromatic imagery at 0.5 m, and WorldView-2, which was launched on 8 October 2007, collects panchromatic imagery at 0.46 m and multispectral imagery at 1.84 m. GeoEye-1, which was launched on 6 September 2008, collects panchromatic imagery at 0.41 m and multispectral imagery at 1.65 m. As Nagendra and Rocchini [9] indicated, very high resolution imagery, which spatial resolution is fewer than 5 m, but only with four to eight bands or less bands, provides more accurate locations of tree canopies. The ability is better than hyperspectral imagery, which has a series of continuous spectral bands (over 100 bands), but coarser spatial resolution (20–30 at best). However, higher spatial resolutions may not have higher classification accuracies, due to higher variances within classes [10] (also called H-resolution problem, which occurs when land cover elements are larger than pixel sizes [11]). In addition, more pixels in very high resolution imagery result in large data sizes, and increase loadings of computer hardware requirements (i.e., segmentation procedures have very high memory and CPU requirements), so most studies only test in small areas, because of memory limitations of personal computers and long processing times. Object-based analysis can improve salt and pepper effects and increase classification accuracies over pixel-based image classifications, which ignore similarity of near pixels [12–18]. In addition, imagery can be stratified into smaller subsets [19] within the processing capabilities of personal computers.

Object-based image classification includes a two-step procedure, image segmentation and image classification. Image segmentation gathers several similar neighbor pixels together as objects, and image classification categorizes or labels objects as land cover types. In theory, image objects have equal internal variances at a common scale [20], and very high resolution imagery is helpful to derive correct locations. However, each class with varying sizes needs different scales to define appropriate objects [14,20,21], and multispectral bands are helpful to labeling or classification procedures.

The appropriate scale for identifying objects can be found by an iterative workflow [22] using hierarchical semantic models or knowledge [23]. In other words, there is no single optimal scale parameter, but a spatially-nested (multi-scale) structure [20–24] can be used to identify objects with different sizes and describe the object traits, especially in imagery without many spectral bands. For example, Yu et al., [18] pointed out vegetation alliances, which are more general than species types, are also critical to define tree species segmentation, and the similarity of near objects may decrease classification accuracies. Kim et al., [25] indicated that multiscale image classification involving both spectral and texture traits can increase classification accuracies. Nevertheless, those hierarchical-segmentation applications in very high resolution imagery are still limited to classify primary land cover categories [25] or estimate forest parameters [26–28] in small test areas, and there is lacking in a general framework for specific species in large areas.

Previous vegetation mapping frameworks have taken two approaches, the data-based orientation and vegetation-type orientation. In the data-based approach, integrating multiple remote sensing data layers can increase interpretations of vegetation traits [25,26,28,29]. In one example, Xie et al., [29] used this approach to identify exotic Australian pine with three-level segmentations: NDVI to distinguish vegetation from non-vegetation, tree heights, derived from LIDAR, to distinguish trees from short trees (lower than 3 m), shrubs and grasses, which did not belong to Australian Pine, and image traits using shape/color and smooth/compactness weighting parameters for target objects. However, higher data requirements raise costs, and may not be extensively adapted to other regions. In the vegetation-type approach, vegetation nested structures from tree, stands, forest types (e.g., pine, oak and red fir) and vegetation types (e.g., wetland, forest and grassland) are another solution [20]. This framework, incorporating object-based image segmentation [20] and labeling procedures [30], can be used to identify more general vegetation life form or land-cover types (e.g., conifer forest, hardwood forest, chaparral, soft chaparral) and within them, more specific vegetation types. As an example, the level of Associations (e.g., Jeffrey pine, black oak and coast live oak) is nested under the US National Vegetation Classification Systems ( http://usnvc.org). However, the labeling procedures for each individual species type, based on multispectral classification and image interpretation, may not all reach equally high accuracies (e.g., accuracies for conifer types ranged from 23% to 100% in Franklin et al., [30], but ecological studies may need to focus on certain species with lower accuracies).

1.2. Problem Statement and Literature Review

The goal of this study was to implement a hierarchical framework for segmenting specific tree species using three bands derived from 1 m digital orthophoto quadrangles (DOQ) and an evaluative framework to assess the segment results. For the purpose of time-series image processing, imagery, which governments started to conduct aerial surveys for, with only one single black/white band is the most widely available historical (long-term) sources. Thus, this study focused on the approach, which better used spatial information, rather than spectral analyses, as DOQ, which has very high spatial resolution, but the digital number may not be real. A further goal was to evaluate methods that can be implemented on personal computers (with conventional amounts of memory). Although identifying individual species can be demonstrated by machine learning methods (e.g., [31]), the hyperspectral and/or very high resolution imagery (less than 0.5 m) that these methods relied on is not always available. Thus, the hierarchical framework, which was built on widely available data and computing capabilities, and incorporated ecological knowledge of the study area, provided the ability to label segments within specific species association or species types, even if imagery data were not able to support species identifications.

Partitioning the image based on tree density was a key part of the natural vegetation framework used in this study (Figure 1). Previous studies emphasized intensive data requirements [25,26,28–30] and repetitive testing on the applications of vegetation mapping, but those approaches may not form robust procedures toward long-term image processing with uncertain data sources (i.e., validation points) and the lack of multiple attributes. Moreover, segmentation studies (reviewed by [32–37]) have not been extensively applied to natural vegetation using a robust approach. Among those segmentation methods, which were classified by Fu and Mui [32], edge-based algorithms were the most vulnerable to noise (i.e., heterogeneous pixels with higher variances, especially in sparse tree stands) [32], and threshold-based algorithms cannot deal with imagery complexities, even using local threshold approaches [38,39]. In contrast, region-based algorithms were widely incorporated with other algorithms, including the Woodcock and Harward [20] region growing algorithm, eCognition’s imagery merging and fractal net evolution approach [21,23,24], and multiscale object-specific segmentation (MOSS) using size constrained region merging [21,40]. Those algorithms could better deal with noise and avoid over-segmentation on vegetation mapping, compared to other algorithms, such as the watershed algorithm and the region growing algorithm [34,37,41]. However, two-step labeling procedures [30] or repetitive testing on parameters (e.g., scale, compact and shape parameters for multiresolution segmentation in eCognition) to construct hierarchical frameworks [21–24] were still time-consuming and data-intensive. Therefore, a robust alternative using image segmentation optimized for delineation of tree density, which distinguishes a specific threshold for separating meaningful species association and selecting segmentation algorithm or parameter, was necessary for simplifying those procedures and data requirements in time-series image analysis.

Furthermore, tree density was the main trait of the hierarchical vegetation framework for two reasons (Figure 1). Tree density distributions, which reflected specific ecological site factors, were more appropriate for occurrence of specific species association (e.g., shade-intolerant species and shade-tolerant species) [42–45]. In theory, shade-intolerant species tend to occur at sparse tree stands, while shade-tolerant species tend to occur at dense tree stands. Additionally, different tree density patterns were required to identify two kinds of image objects, tree species and vegetation stands, using two processing procedures. For identifying individual tree species, tree crown delineation algorithms and labeling classifiers using multispectral bands were widely used, especially for interpreting species counts and types, which can substitute for field-based surveys [46–50] to some extents. Those algorithms included valley-following, region growing and watershed segmentation [51–55]. Nevertheless, the fundamental assumptions were that tree crowns should be clearly separated in space (not overlap) and should be in regular shapes and similar sizes, according to algorithm functions [52,53]. Thus, no general rule in the literature to select algorithms is available [52,53]. Furthermore, this approach was less effective in heterogeneous and denser hardwood stands than conifer landscapes [51,56]. For vegetation stands, decomposing landscapes into smaller objects (vegetation stands, a group of trees) and labeling species by the majority classified pixels within segmentations [15], vegetation gradient model adding spectral mixture analysis [30] or non-parametric classifiers (e.g., [18,23]) were usually applied (e.g., [18,57–60]). Nevertheless, segmentations in sparse tree stands may not be partitioned well, due to effects of non-vegetation areas or shadows. In addition, precise vegetation boundaries were still challenging in denser tree stands, because region-based segmentation partitioned imagery, based on its variances, rather than true objects. To sum up, preferable algorithms in the literature for the two kinds of objects had their own problems, and no general solution was available. As a result, one of goals in this study was to test a flexible guideline for algorithm or parameter selections.

Specifically, the hypotheses tested in this study were:

• Environmental stratification, which environmental variables were separated by the global Otsu’s method for strata, will reduce the image variation of digital number (DN) and texture. It is assumed that reducing variance in image DN and texture indices simplified species composition based on studies showing that texture indices (local variances and second order textures) have high correlations with forest structural parameters (e.g., standard deviation of diameters and basal area) [61,62] and can clearly distinguish different forest parameters (e.g., stand ages) [63].

• Identifying tree-density patterns will be used to simplify species association in a given area where there was no recent disturbance. This hypothesis was based on assumptions that shade tolerance was related to variations in tree architectural parameters, such as stem and crown dimensions [42–45], and species life history traits (e.g., different seedling time between oaks and pines) and environmental conditions (e.g., topographic and climatic conditions), which constrained species regenerations, resulting in different species dominances (pine dominance and oak-pine codominance) in stands with different tree densities [64,65].

• Different tree density distributions will have different optimal scale parameters for image segmentation where the optimal scale parameter was that with the lowest segmentation variances and spatial autocorrelations [66–68]. The hypothesis was based on studies indicating that different image objects, such as sizes of tree crowns or clumps are a function of different spatial scales [10].

2. Methodology

Step one and step two, which involved the global Otsu’s method using environmental variables and the random walker algorithm on the imagery, was used to identify two tree density patterns. Step three was the multiresolution segmentation for extracting species objects (Figure 1). Environmental variables, based on elevation, slope or aspect, were used to partition the whole image into smaller subsets. Large images needed to be subset in order to reduce processing time and variation partitioning can help to statistically characterize certain components by increased stationarity or homogeneity within strata [69,70]. A non-parametric test was used to see which terrain variable was most effective at reducing variation in local DN and texture measures within strata. The reasons for using one environmental variable, rather than all three, were that one of goals in this study was to build a parsimonious procedure for reducing computations, and stratification based on a single variable was simple to implement. Then, the random walker algorithm was applied within the environmental strata to partition imagery into smaller subsets based on tree density patterns. The segmentation was evaluated by independent vegetation maps (the Timber survey and CALVEG) for segmentation accuracies. The goal was to evaluate whether the segmentation distinguishes image regions associated with particular species or species associates (Figure 1). Finally, the segmentation of tree densities was further partitioned by multiresolution segmentation with different scale parameters for species or stand level segmentation. The segmentation efficiency was evaluated by whether the detailed species segmentation related to specific guidelines of scale parameters on identifying species objects (Figure 1). So as to minimize costs, this study mainly used open source python libraries, scikit-image and other free software, and free DOQ imagery, and reduced processing time in commercial software to achieve the above goals.

2.1. Study Targets and Data Sources

Tejon Ranch, which belongs to Tejon Ranch Company, is located in the convergence of four eco-regions: the Mojave Desert, the Central Valley, the Sierra Nevada, and the Transverse Ranges [71]. The research area is located in the Tehachapi Mountains, elevation ranging from 368 to 2360 m [72] (Figure 2). Based on the only climate station within Tejon Ranch, 434.3 m elevation, the average yearly temperature is 59.61 °C (1895–2011) and yearly rainfall is 11.29 inches (1899–2011) [73]. This area has a typical Mediterranean climate. The main rainfall season concentrates in the winter from October to March of the next year, while the dry season is in the summer from July to September. The main vegetation type is oak woodlands, including canyon live oak, interior live oak, blue oak, California black oak, scrub oak, and others along with ponderosa pine and gray pine.

This study focused on two hardwood species, blue oak and California black oak, and two conifer species, ponderosa pine and grey pine. The four species are shade-intolerant [74], and tend to occur in sparse tree stands. Nevertheless, the shade tolerance of California black oak varies with age (higher shade-tolerant in saplings and less shade tolerance up maturity), so it can sustain in denser tree stands [74]. More importantly, owing to the long-term fire disturbance history, tree density does not always achieve its maximum possible value, and local adaptations to microenvironments and regeneration played more important roles in determining tree density. For example, ponderosa pines in southwestern United States were often found in dense stands changing from a range of 49–124 trees ha-1 at the time of Euro-American settlement to a range of 1235–2470 trees ha-1 now [75–77], due to high density of seedlings and saplings and fast growth near burned areas [65]. To sum up, species successions were not clearly identified in previous literatures and records and any rules about determinants of tree densities were not extensively applied to regions [74], so it is assumed that species have similar tree density patterns within a given area.

Data sources included a one-meter DOQ, SRTM 30 m Digital Elevation Data v.4 and two archived vegetation type maps, referred to in this study as Timber survey and the CALVEG. The DOQ was imagery with three 8-bit bands, which was rectified by digital terrain models and ground position points to remove terrain relief and camera tilt [78] so every land element was in corrected ground position. The horizon geometric accuracy, which was evaluate by a root-mean-square-error (RMSE), is 0.82 m [79]. The image was collected in 9 December 2009. This imagery is freely available from United States Department of Agriculture: Natural Resources Conservation Service Geospatial Data Gateway ( http://datagateway.nrcs.usda.gov/GDGOrder.aspx) or United States Geological Survey: The National Map Viewer and Download Platform ( http://viewer.nationalmap.gov/viewer/). However, because the DOQ was produced by aerial photos, the radiometric ranges cannot reflect the real spectral digital numbers. The SRTM 30m Digital Elevation Data, which were produced by National Aeronautics and Space Administration (NASA), were created by collecting elevation points from SRTM3 and a series of auxiliary digital terrain models for the purpose of interpolating voids to create seamless topography [72]. The data provided not only elevation information in 30 m spatial resolution but also slope and aspect gradients by calculating the rates of maximum change and their directions. The Timber survey was a field-surveyed map generated in 1980 for distributions of oak woodlands on Tejon Ranch [80]. The CALVEG (Classification and Assessment with Landsat of Visible Ecological Groupings) was produced by the USDA Forest Service using Landsat Thematic Mapper to construct a vegetation database across California in 1990s, and redeveloped in 2009 through 2010 using Landsat Thematic Mapper and Spot 5 panchromatic images [30]. Although the time of existing vegetation maps cannot correspond to the DOQ’s, this study assumed that the vegetation did not change much since 1980. The assumption was based on the fact of poor regeneration of California oak specie [81].

2.2. Environmental Stratification

To create environmental stratifications, a global Otsu’s method was used to separate elevation, slope and aspect images into smaller areas using scikit-image. The scikit-image ( http://scikit-image.org/), written by a community of volunteers, is one of python open-source peer-reviewed code libraries. The library is a collection of image processing algorithms in previous studies. The goal applying global Otsu’s method was to find an optimal global threshold in image histogram (one-dimension intensity measurements) to partition imagery into two groups. The function was to maximize the intra-group variances and minimize within-group variances [82] in order to categorize objects from backgrounds, such as tree crowns from barren lands. Preliminary analyses determined optimal thresholds; the three stratifications were partitioned by thresholds of 1169 m for elevation, 41° for slope and 149° for aspect (Figure 3).

The global Otsu’s method has been shown to be more effective than other binarization methods [35,39] with lower computational time. However, those binarization methods were not spatially separated [32,35]. The strata boundaries needed to be delineated by hand, depending on the spatial distributions of environmental variables beyond or below the thresholds. The manual boundaries of strata did show diverse patterns to reflect tree covers. Then, the random walker algorithm was implemented within each stratification for further partitions.

2.3. Image Segmentation Optimized for Delineation of Tree Density

The random walker algorithm was used for image segmentation. Segmentation was based on single first band. The goal of using a single band, instead of three bands, was to save computation time, if the single band can reach the goal to identify tree density. Additionally, DOQ may not represent the true spectral digital numbers. The random walker algorithm, which was proposed by Grady [83], was used to partition imagery into smaller subsets more robustly using scikit-image. The random walker algorithm originated from the graph theory to view the whole imagery as the combinations of vertexes (nodes) and edges (arcs), and the random walkers, which represented each individual pixel, were trying to formulate a path to their neighbors, according to certain probability distributions. The algorithm started from defining the markers, a group of seeds as the sampling of desired imagery objects. Then, this algorithm would assign each unseeded pixel a probability and a weight. The probability was assigned by absolute distances from seeds to each individual image pixels (unseeded pixels), while the weight was converted by the image intensity. Finally, the algorithm would assign a seed class to those unseeded pixels for cuts, mainly according to the probability, and the cuts might be adjusted by the weights to avoid crossing sharp image intensities. For example, if three out of four neighbors belong to one class, then the focal pixel is assigned to this class. Thus, the random walker algorithm can keep locally consistent boundaries, regardless of spatial extents and attribute ranges.

However, how to define the markers was not described in details in the original paper, and few studies have used the random walker algorithm to identify vegetation types, because the random walker algorithm was under the assumptions (also called supervised segmentation algorithms) that a series of pixels for desired objects and backgrounds were known and nearby pixels between desired objects and backgrounds can evolve to desired boundaries. To handle this issue, the grey-level image DN, which was assumed to be a function of tree densities, was used in the marker designs. The image DN was higher, and the tree densities were sparser. Using natural break optimizing, the map classification was to reduce within-group variances, in three classes to derive the thresholds for desired objects and backgrounds. As a result, barren lands or individual trees over the highest threshold were viewed as the desired objects, and denser tree stands below the lowest threshold were treated as the background (Figure 4). Then, the random walker algorithm could separate desired objects from backgrounds. In the second segmentation, the same procedures were applied to separate lower tree density stands from higher tree density stands (Figure 4). Finally, barren lands or individual trees and lower tree density stands were combined as lower tree density stands for species association labeling (Figure 4).

2.4. Species Segmentation

Detailed vegetation objects were further partitioned following the image segmentation optimized for delineation of tree density using multiresolution in eCognition 8, which has been extensively used in studies of object-base image classification for vegetation inventories (e.g., [18,58,59]). The algorithm of multiresolution segmentation was demonstrated by optimization procedures to minimize heterogeneity for each individual merging [24]. Although many algorithms were developed to delineate tree crowns and open for free, those algorithms are scene-dependent and time-consuming to test them all in various operational system and computer languages. Thus, eCognition can be the simplest way to see the effects of the hierarchical framework on species segmentation. Segmentation was produced by 17 attributes or features: RGB; first-order textures: mean, standard deviation, kurtosis, mean Euclidean distance, skewness, variance; second-order textures: mean, variance, homogeneity, contrast, dissimilar, entropy, second moment, and correlation [84,85]. For the purpose of understanding which attributes can form better segmentation, several attribute combinations were tested, especially for the second-order textures, which literature pointed out high correlations with vegetation characteristics [61,62]. Nevertheless, multiresolution segmentation counted on classifiers to label after segmentation procedures. Therefore, as a prior framework, the multiresolution segmentation cannot be implemented with appropriate parameters or ruleset beforehand, and applied the same setting to other images directly. Optimal scale parameters were needed to be identified. The tested scale parameters (without units, depending on image heterogeneity), which determined the object merges and object sizes for the purpose of reducing heterogeneity [24], included 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5, and 5.5.

2.5. Segmentation Evaluations

To test the hierarchical vegetation framework and the three hypotheses, a series of evaluations were used. First, image local DN and texture, mean and standard deviation in the first band, were used. The local mean and standard deviation were calculated based on window sizes, 3 by 3, 5 by 5, 7 by 7, 9 by 9, 15 by 15, 21 by 21 and 27 by 27. The goal was to see which terrain variable could be used to stratify the image most effectively, because image DN and texture correlates with forest structure parameters. A non-parametric test, Kruskal-Wallis test was applied using SPSS Statistics 20.0 to test whether variances (the local mean and standard deviation) were equal among groups (different environmental stratification).

Second, four focal species were chose to examine the effectiveness of the random walker algorithm by determining whether the two tree density patterns, which were combined from three tree density types (Figure 4), reflected different species association distributions as determined by the archival vegetation maps (reference data). In other words, the CALVEG and the timber survey were used as reference data for sampling four focal species types. The species samples were combined into species association and used to compare with tree density distributions within the study area. Indices of map agreement, including the kappa value, overall accuracies, user’s accuracies and producer’s accuracies, are based on the values in a confusion matrix [30,86]. Sampling the maps in order to calculate those indices were based on a stratified random approach accomplished using a “plug-in” sampling design tool in ArcGIS 10.0 to allocate 100 points proportionally by class areas. The strata were based on the three image segmentation optimized for delineation of tree density (Figure 4), where larger segmentation areas had more sampling points.

In order to evaluate the scale parameters used in the multiresolution segmentation, this study applied an objective (unsupervised) evaluation approach by using Moran’s I [67,68] as a measure of spatial autocorrelation to measure the similarity of segment-averaged attributes as a function of the distance between segments. There are several reasons to select the unsupervised evaluation approach. First, the approach could be applied in areas without abundant field data for supervised evaluations. The CALVEG was classified from SPOT (20 m) and Landsat TM (30 m), while the Timber survey was in coarser segmentation than the CALVEG. Also, the two reference maps were created at different scales [80], so it would not be possible to distinguish whether inaccuracies of vegetation segmentations were from boundary errors in maps made at other scales. In brief, it was questionable to sample points for validation of species distributions in 1-meter DOQ. Second, it is still challenging to evaluate results of object-based image analysis, compared to accuracy assessments in pixel-based image classifications, which applied sampling points to evaluate agreements [57,87,88]. Pixels within segmentations were not spectrally homogeneous [89], so the sampling points may not represent the whole objects. In addition, some of segmentation tools are proprietary software, such as Berkeley Image Segmentation (BIS) ( http://www.berkenviro.com/berkeleyimgseg/) and eCognition ( http://www.ecognition.com/). Although they have been described in the literature (e.g., [22–24,90]), there was no way to evaluate the segmentation accuracies by algorithms.

The basic idea of unsupervised evaluation was that attributes of optimal segmentations with clear boundaries (such as average spectral or texture values) should have low between-segment spatial autocorrelation, whereas over-segmentation (many small polygons), and under-segmentation (where the polygons are larger than the optimal segmentations), yielded segments whose attributes had higher between-segment spatial autocorrelation (Figure 5). Kim et al., [68] further pointed out that the segmentations in optimal scale parameters can contribute to better classification accuracies. Therefore, the segmentations based on different scale parameters were evaluated based on the Z values of Moran’s I in two test areas, sparse tree density stands and dense tree density stands (Figure 6). The goal was to see whether there were different optimal scale parameters for segmentation of two tree density patterns, because each has different sizes of individual trees and tree stands. This reason to use the Z values of Moran’s I, instead of Moran’s I in Kim et al. (2008) and Kim et al. (2009) [67,68] was that image segmentations with different scale parameters resulted in different weight matrices among images. For comparisons for an appropriate scale parameter, Moran’s I is required to standardized as the Z values of Moran’s I.

3. Results

3.1. Environmental Stratification

The aspect stratification had the best performance in terms of increasing stationarity within strata. Both DN and texture, mean and standard deviation of first-band digital number (DN) for different window sizes, showed significantly reduced variance within strata based on aspect, and the slope strata had the second best performance (Figure 7). Visually, the aspect stratification was most effective in separating tree crowns from barren lands, and the slope stratification again had the second performance (Figure 3). The slope stratification could be used to identify parts of barren lands, which had less steep slopes, but the elevation stratification was not useful for distinguishing woodland cover strata in this study area. Therefore, the aspect and slope variables produced strata that reduced variance in local measures of DN and texture, supporting the first hypothesis.

3.2. Image Segmentation on Tree Density

The two image segmentation optimized for delineation of tree density based on the random walker algorithm effectively separated the two species association. Ponderosa pines and California black oaks tended to occur in dense tree stands, while grey pines and blue oaks tended to occur in sparse tree stands within the study area. Agreement of species association with image segmentation optimized for delineation of tree density reached about 80% overall accuracies and 0.6 kappa values. In particular, using the aspect stratification for the thresholds of desired objects (firstly, for barren lands or individual tree stands and secondly for lower tree density stands) by natural break optimizing for the random walker algorithm did perform the best to produce image segmentation optimized for delineation of tree density (

Table 1. Agreements Between Species Association And Image Segmentation Optimized For Delineation Of Tree Density, measured by comparing two existing vegetation maps (Timber Survey, CALVEG) with the segmentations using Kappa, overall accuracy (percent correct classification), user’s accuracy (1-commission error), and producer’s accuracy (1-omission error).

(a) Timber Survey
	Elevation Stratification	Slope Stratification	Aspect Stratification
Kappa value	0.52	0.56	0.65
Overall accuracies	0.76	0.78	0.83
Users’ accuracies
Barren land	0.72	0.8	0.88
High density	0.81	0.76	0.76
Producers’ accuracies
Blue-oak and grey pine	0.82	0.8	0.81
Black oak and ponderosa pine	0.71	0.76	0.84
(b) CALVEG
	Elevation Stratification	Slope Stratification	Aspect Stratification
Kappa value	0.61	0.64	0.65
Overall accuracies	0.80	0.82	0.82
Users’ accuracies
Barren land	0.64	0.76	0.72
High density	0.96	0.88	0.92
Producers’ accuracies
Blue-oak and grey pine	0.94	0.86	0.9
Black oak and ponderosa pine	0.74	0.79	0.77

). Visually, step one and step two were effective to separate tree density patterns into three tree density categories, barren lands or individual trees, lower tree density stands and higher tree density stands, and the two selected regions of two tree densities had similar performances on separations of tree densities (Figure 8).

3.3. Species Segmentation

The evaluation of scale parameters suggested that multiresolution segmentation may not delineate appropriate detailed vegetation segmentation based on this imagery. Neither lower tree density stands nor higher tree density stands had an optimal scale parameter. The optimal scale parameter indicated that certain scale parameter would reach a minimum spatial autocorrelation, compared to its larger and smaller scale parameters. The Z values of Moran’s I continuously declined as the scale parameters increase, showing no local minimum (Figure 9). The difference between results for sparse tree density stands and dense tree density stands was the magnitude of Z values of Moran’s I in most attribute combinations. The sparse tree density stands had higher values, while the dense tree density stands had lower values (Figure 9). High variances of sparse tree stands, especially in large areas of barren lands, dominated the patterns of segmentations, while mixed species and overlapping tree crowns in dense tree stands resulted in similarity among segmentations. However, several segmentations, which were based on single second-order texture, especially homogeneity, dissimilarity and entropy, in both sparse tree density stands and dense tree density stands did not perform well (Figure 9).

Visually, similar neighbor segmentation was the main problem in either denser tree density stands or sparse tree density stands, regardless of scale parameters. In the sparse tree density stands, although the main tree crowns could be identified, the boundaries of tree crowns were not precise. Concretely, near segmentations of barren lands confused with parts of tree crowns using smaller scale parameter (over segmentation), while tree species segmentations included parts of barren lands using larger scale parameter (under segmentation). Under the tested scale parameters, there was no optimal points to balance the over segmentation and under segmentation. In the dense tree stands, the similarity among segmentation was still high. Segmentations using smaller scale parameters resulted in very small, complex segmentation. Through spectral difference segmentation, original segmentations with larger scale parameters merged into very large segmentations (Figure 10).

4. Discussion

4.1. Implications of Image Segmentation Optimized for Delineation of Tree Density

Under the limitations of long-term available data and computation abilities of personal computers, image segmentation optimized for delineation of tree density demonstrated in this study provided an alternative framework for vegetation mapping, instead of data intensive analyses. The approach involved both the global Otsu’s method and the random walker algorithm using terrain variables and 1 m very high resolution imagery with only three bands as limited inputs. Those procedures can be carried out in Window 7 64-bit operating system using Intel(R) Core(TM) i7-2860QM CPU with 12 GB memory. The concrete goals were to add environmental variables, which were correlated with tree covers, and reduce variances within environmental strata and image segmentations. Specifically, those segmentation results reflected the applications of ecological understandings, instead of intensive data requirements and repetitive tests on segmentation and labeling. For one thing, water availability, which west-side slopes, confronting the ocean, had more precipitations than east-side slopes, might be the reason for the best performances of the aspect stratification, although no high-density climate stations or climate mapping at high spatial resolution can be used to validate. For another thing, because species associations are correlated with forest density in the study area, image segmentation optimized for delineation of tree density could be labeled by specific species associations, even though the third-step procedure of detailed species segmentations did not work well at identifying individual species segmentation, due to coarse spatial resolutions and less attributes.

In this study area, although four target species are shade-intolerant, they do not all occur in early-successional or low density stands as expected for shade-intolerant species, because of their varying shade tolerance at different ages (i.e., California black oak) and the human management history (i.e., ponderosa pine). Ponderosa pines occur in dense tree stands as do California black oaks, while blue oaks and grey pines occur in sparse tree stands. The agreements between image segmentation optimized for delineation of tree density and species association reached about 80% (overall accuracy). Further testing of alternative species labeling strategies in different ecosystems is needed for general applications.

Indices summarizing the confusion matrix were used to examine hypothesis 2, species association labeling. However, those common indices, especially the kappa value, have confronted harsh challenges in recent decades. The main criticism was that the kappa value was constructed by comparisons between the reference map and classification map, but the two maps may not be meaningfully comparable. Overall, the kappa value requires that both maps followed some assumptions, such as fitting normality and not being affected by other covariates, so adjustments on each map were the key, such as the weighted kappa coefficient [91]. In remote sensing cases, the kappa value was only a summary statistics through randomness sampling as a baseline, not a meaningful index to indicate quantity disagreement and location disagreement [92]. This study applied the kappa values and the overall accuracy as indices, because the stratified random approach using environmental variables, instead of wholly random sampling was effective to emphasize species locations (e.g., [93]), and tree density focused on whether species association labels were in specific image segmentation optimized for delineation of tree density, rather than improving detailed land use mapping or species mapping. Thus, simple indices were adequate for comparing the effects of environmental stratifications and assessing the correctness of species association locations.

4.2. Implications of Species Segmentation

The implications of species segmentation were to test the applications of the vegetation hierarchical framework, and whether the results related to specific guidelines of algorithm or parameter selections. The results employed an unsupervised evaluation (the Z values of Moran’s I under 999 permutations to evaluate whether the attributes of segmentations (average within segmentation) were significantly different from the neighbor segmentation. Although all showing decreases of the Z values of Moran’s I with increases of the scale parameters using different attribute combinations for segmentation (Figure 9), and image objects in two test areas both violated the object assumption with equal internal variances [20], the two tree density patterns did have fairly different reactions in most of attribute combinations on defining detailed species objects.

The unsuccessful segmentation at the third level resulted from two aspects in object definitions. In dense tree density stands, individual tree crowns could not be delineated, since the 1 m spatial resolution imagery was still too coarse. The spatial resolution is larger than the tree crowns. The result was similar to the outcomes of Woodcock and Strahler [10], where the coarse spatial resolutions did not have an optimal peak to reflect tree crown sizes (30 m imagery presented an asymptote, which local variances decreased as spatial resolutions increased, while 0.75 m imagery showed a local peak). The Z values of Moran’s I also showed an asymptote, which the Z values of Moran’s I decreased as scale parameters increased. In addition, homogeneity throughout the whole image resulted in very complex segmentations, especially using smaller scale parameters. As a result, adding higher spatial resolution image data or ancillary data, especially for oak woodlands, which often have lower classification accuracies (e.g., [56]), and using larger scale parameters for homogeneous areas may be a better alternative.

In the sparse tree density stands, the segmentations using multiresolution segmentation also could not derive correct tree crown objects, although individual tree crowns were clearly separated visually. The reason came from not enough information to delineate tree crown objects from the neighbor segmentation. Even if involving multiple texture features, segmentations were still limited to characterizing individual tree crowns. The segmentation may improve performances by implementing other segmentation algorithms, such as the watershed algorithm, as the example in Baatz et al., [22] or the valley-following algorithms and adding attributes, which can reflect vegetation characteristics (e.g., near-infrared band).

Although the results showed that multiresolution segmentation cannot deal with detailed species segmentations based on the DOQ imagery, the third evaluation still provided guidelines. As Fu and Mui [32] pointed out, region-based algorithms may differ in segmentation results, according to the order of region-merging, even if region-based algorithms were widely used and had higher classification accuracies than pixel-based imagery classifications. As a result, applying global information as nested models within certain classes is necessary. In other words, a single algorithm may not apply to whole hierarchical frameworks, even with the multi-scale concept.

4.3. Applications of Hierarchical Vegetation Framework

The hierarchical vegetation framework, developed in this study, is more time- and cost-saving for species association labeling. For one thing, data can be collected and automatically processed at lower costs. Very high resolution imagery in the form of digital photography is widely available, and environmental variables help to improve imagery processing for broad-area investigations. For another thing, the approach can incorporate current imagery processing methods, and does not require repetitive testing, like procedures in eCognition, and was effective to reduce time on using eCognition.

Concretely, aerial photos with single imagery DN are more available data sources for vegetation mapping, and sometimes, aerial photos with very high spatial resolution are in archives, either national-level imagery datasets (e.g., USDA: Natural Resources Conservation Service Geospatial Data Gateway or ASO Taiwan Image supplier and services System) or international-level imagery datasets (e.g., EarthExplorer). For example, Corona Lanyard (1963) and Corona KH-4B (from 1967 to 1972) both provided 1.8 m imagery worldwide, and Corona imagery can be bought through EarthExplorer website ( http://earthexplorer.usgs.gov/) with very small costs. Meanwhile, SRTM, a digital terrain dataset has worldwide coverage in 90 m spatial resolution, and can be used to assist image processing.

In contrast, Landsat imagery with coarser spatial resolutions has been widely used in long-term land use/land cover mapping and monitoring [94], because of its long-term consistent platforms, since 1972 and multiple bands, which reflect vegetation distributions (near-infrared) and temperature (thermal-infrared). Nevertheless, except in the United States, Landsat satellites may not offer enough imagery with continuous or regular time intervals and good imagery qualities, such as low cloud covers in wetter tropical areas. Furthermore, other detailed ancillary data, which were helpful to assist imagery interpretations, may be available at national-level datasets, such as U.S. National Elevation Dataset 10 Meter ( http://datagateway.nrcs.usda.gov/GDGOrder.aspx), U.S. General Soil Map (STATSGO) ( http://websoilsurvey.nrcs.usda.gov/app/WebSoilSurvey.aspx), or LiDAR Topography Data ( http://opentopo.sdsc.edu/gridsphere/gridsphere?cid=datasets). However, those datasets cannot reflect temporal changes of species distributions and are limited in survey areas.

This research provided a robust way for time-series image processing. The approach can be achieved by a single-band, very high resolution image at lower time, labor and facility costs. However, this approach still suffers several limitations. First, the global Otsu’s method applying to the aspect factor did not consider areas where the terrains mostly face north, since the situation (when 360 deg. = 0 deg) was not included. As the aspect factor performed the best as the first step to reduce variances with ecological meaning (tree covers), the global Otsu’s method is also the best way of fast computation. Second, species association labeling may only be effective in certain areas, as the kappa values only reached 0.6. Both disturbance frequencies and intensities may violate the second hypothesis. Especially, the approach cannot apply to areas suffering large, infrequent disturbances, which may result in more heterogeneous species distributions [95]. The goal of this study was not to provide an ultimate solution for general ecosystem labeling, but to propose an examined framework for an efficient approach on time-series images. Third, compared to commercial software, open-source solution cannot always provide equal functions, regardless of personal expertise. As open-source tools may be developed for usage in other fields (e.g. medical image processing, instead of vegetation mapping), and code quality is not certain among different libraries, it takes much time and effort for researchers to write their own codes or test the workable functions for research questions. For these limitations, there are several reasons to use the open-source tools as the solution. For one thing, this research selected a peer-reviewed image processing library, scikit-image, which reduced code uncertainty. For another thing, this research was not trying to fill the gap between commercial software and open-source tool, but provided the best way of using current available data and tools, especially for long-term image processing. To be brief, regardless of limited image information, this research intended to offer a starting point of using open-source solutions on time-series image processing.

Future studies should focus on two aspects. One is finding a new method for species association labeling in different ecosystems. Ecological understanding and local knowledge should be applied to vegetation mapping. The other focus of future research should be the techniques themselves. Algorithms of tree crown delineations can be tested at sparse tree density stands and appropriate object definitions should be established for dense tree density stands. In particular, individual tree density patterns reflect different biological interactions over times. For example, at the species spreading front, reduced intra-specific competition in dense populations may increase population growth rates and migrations, but sparse populations may suffer decreased migrations [96]. However, highly overlapping tree crowns make it difficult to detect changes in density of individual species (e.g., [97]), and tree crowns in sparse tree stands cannot be identified using multiresolution segmentation. Therefore, based on the hierarchical framework, two separate algorithms with more precise object definition are required.

5. Conclusions

Object-based analysis with a multiple scale framework is one of solutions toward time-series image processing. However, the data-based orientation approach required more data sources and the vegetation-type orientation approach was limited to improve classification accuracies of specific species types for ecological modeling. In particular, long-term image sources, which have less image information and are lacking in detailed validation data are still challenging to deal with. Also, large amounts of tests on labeling procedures and obtaining the best multiresolution framework are still required. As a result, a robust segmentation approach mainly using open-source tools with less data requirements (i.e., SRTM, digital orthophoto quadrangles with 1 m spatial resolution and existing vegetation maps as validation data) was proposed.

The author applied a three-step procedure for the above issues. At the first level-segmentation, the global Otsu’s method was applied to characterize vegetation traits and reduce computation loads by gathering similar components together into strata on three environmental variables (i.e., slope, aspect and elevation). With the evaluation by the Kruskal-Wallis test, this step, especially the aspect factor was proved effective to reduce image variances in local DN and texture and more importantly, characterize similar components of tree covers within strata. At the second step, image segmentation optimized for delineation of tree density was used with the least computation requirements. Species association labeling for four target species, which were grouped into two species associations was tested for effectiveness, and indicated about 80% in overall accuracy and 0.6 in kappa values between image segmentation optimized for delineation of tree density and species associations, which were derived from two existing vegetation maps. In the third step, the image segmentation, optimized for delineation of tree density, was further tested for its applications on partitioning species segmentation using the Z values of Moran’s I. However, this step using multi-resolution segmentation did not bring an optimal scale parameter, due to limited information on a single-band imagery. More importantly, the results, which showed two reflections on the Z values of Moran’s I, indicated the suggestions for segmentation algorithm selections and developments.

Although the image information is limited to identifying individual species, and the species association labeling for the evaluation of image segmentation optimized for delineation of tree density is area-based, the study still provides the following contributions to vegetation mapping. First, the hierarchical vegetation framework for image segmentation provided the possibility of long-term change detections, at least to compare certain ecological indices over time. For example, worldview-3 or Geo-eye provide more opportunities to delineate species segmentation, while Landsat MSS, TM, ETM+ or Corona are unable to segment species, as compared to recent imagery. As more and more spatial-temporal modeling approaches are developed, long-term imagery change detections in object (entity) forms provide the potentials to extend current models to the recent past (over 70 years)[98]. Second, the hierarchical vegetation framework with ecological meanings is effective to reduce computation loads and indicate the processing procedures. In the former two steps, the whole image dataset (22393 × 19458) was processed, while the third step was limited by the commercial software, which required extra facility, such as the eCognition server ( http://www.ecognition.com/products/ecognition-server) for big data processing. As more historical images are available on the Internet, this approach can provide solutions to process large amounts of image sources with diverse quality. Third, based on the current object definition, the limited image information was proven unable to produce detailed, reasonable image segmentation with an optimal scale parameter. As the Z values of Moran’s I are more proper to be used as the evaluation index for appropriate segmentations, the index was also used to evaluate the optimal segmentation level. In this study, the second level segmentation is optimal, since the near segmentations at the third level are similar. Although multiresolution segmentation can produce segmentations, those segmentations may not be meaningful, as spectral difference segmentation further aggregated all segmentations together. At this point, the open-source solution using threshold-based segmentation is helpful to deal with those kinds of images.