Discriminating disturbance from natural variation with LiDAR in semi-arid forests in the southwestern USA

Discriminating amongst spatial configurations and climax size of trees in forests along varying physical gradients from time since last disturbance is a significant component of applied forest management. Understanding what has led to the existing vegetation's structure has important implications for monitoring succession and eco-hydrological interactions within the critical zone: the near-surface environment where rock, soil, air, and biota interact and regulate ecosystem services. This research demonstrates the utilities of local indicators of spatial association (LISA) to (1) quantify natural variation in forest structure as derived from aerial Light Detection and Range (LiDAR) across topographically complex landscapes at ecologically relevant scale, i.e., individual trees; (2) map previously recorded but poorly defined forest disturbances; and (3) link scalable topographic indices to observed tree size distributions. We first selected a priori undisturbed and disturbed stands scanned by aerial LiDAR t...


INTRODUCTION
Foresters and ecologists have long been challenged by the inherent scaling limitation of field techniques which limit their ability to relate plot level measurements to landscape scale processes. Although conventional forestry field protocols (i.e., Bechtold and Patterson 2005) provide abundant locally detailed information about composition and structure, plot-based spatial averaging results in large uncertainty at landscape scale (Chave et al. 2004, Tomppo et al. 2008, Mascaro et al. 2011. Conversely, remote sensing provides information at landscape level extents, but typically sacrifices the spatial resolution and detailed measurement information of plot-based inventory. This difficulty in linking landscape patterns to plot scale measurements is part of what Levin (1992) called 'the problem of pattern and scale.' Prior to airborne laser scanning (ALS), also called Light Detection and Range (LiDAR), remote sensing technology could not precisely resolve the most fundamental unit of ecology: the individual. With LiDAR foresters can characterize individual trees from the plot to the landscape scale (Falkowski et al. 2006, van Leeuwen and Nieuwenhuis 2010, Kaartinen et al. 2012, Vauhkonen et al. 2012, Maltamo et al. 2014). Today's dilemma is not how to collect enough detailed information to characterize a population with statistical significance, but how to analyze large LiDAR data sets in ways that answer the most beneficial questions at landscape scale.
Foresters have used 'site indices' to predict tree growth and to estimate potential maximum tree height from environmental factors at the plotscale for over a century (Woolsey 1911, Minor 1964, Clutter et al. 1983. Differences in growth rate and potential maximum tree size with age, i.e., the site index, in semi-arid forests are variable and emergent by elevation, slope and aspect (Whittaker andNiering 1975, Niering andLowe 1984). For example, in the southwestern USA at low elevation excessive heat and dryness result in low productivity desert environments; as elevation increases average daily temperatures decrease and average annual precipitation increases in relation to adiabatic lapse rate until the mountainsides are capable of supporting woodland and forest species. Biomass in southwestern USA semi-arid forests peak in the so-called mixed-conifer forest type where the growing season is not limited by freezing air temperatures and precipitation is abundant. At higher elevations temperature continues to decrease resulting in a progressively shorter growing season. Southfacing aspects (in the northern hemisphere) receive more direct solar radiation than do north-facing aspects. As a consequence south aspects are hotter and drier, and often support more xeric and light-tolerant vegetation relative to north aspects at the same elevations. In later successional stages after vegetation has established and shaded the soil surface, local microclimate conditions can increase soil wetness which allow for more mesic species to establish (Breshears et al. 1998). In semi-arid ecosystems availability of soil water with increased soil depth is key to sustaining increases in plant metabolism and thus greater biomass over time. Maximum tree height is known to vary with topographic position even in temperate and tropical forests (McNab 1993, Detto et al. 2013. Topographic wetness also affects total soil moisture and pH, which in turn affects the species richness of vascular plants (Zinko et al. 2005, Sørensen et al. 2006, Kopecký andČížková 2010, Moeslund et al. 2013). In the present study we explore the relationship of vegetation size and community complexity both by topographic position and wetness in relation to known disturbance histories.
Semi-arid forests are particularly vulnerable to episodic drought, insect outbreak, and wildfire disturbances that alter forest structure (Allen et al. 2002). Such threats will likely intensify as global warming progresses (Westerling et al. 2006, Williams et al. 2012, Stephens et al. 2013. Identifying and tracking the legacies of past disturbances are critical to managing semi-arid landscapes effectively (Allen et al. 1998, Swetnam et al. 1999. Across the late 19th and 20th century southwestern USA forests were significantly affected by wildfire suppression and anthropogenic land use practices, particularly logging and livestock grazing (Touchan et al. 1996, Swetnam et al. 1999, Allen et al. 2002. As a result, many southwestern forests' tree size distributions and age structure are today considered to be outside their 'historic range of natural variation' (Morgan et al. 1994). With incomplete historical records of disturbance, analyses of current forest condition and disturbance legacies cannot be placed into a developmental continuum useful for forest managers. For example, a contemporary mixedconifer stand could represent a climax state, a post-disturbance sere, or be the result of fire exclusion which led to a shift in dominance from pine-oak to more shade tolerant species (O'Con-nor et al. 2014a).
An important step in characterizing land use history along a developmental continuum is to differentiate whether existing vegetation represents the climax state or is an early-or midsuccessional stage following an unrecorded disturbance (Black et al. 1998, Foster et al. 1998, 2003. Land-use and disturbance histories are not particularly well documented in the southwestern USA prior to the 1900s (Bahre 1991, Turner 2003. Often only paper maps and reports survive, with an unknown number of records having been lost or discarded. Even if records of an event (either human or natural) exist, it is unlikely that they include information with precision sufficient to allow for geo-referencing. Further complications involve historical records that do not separate planned management actions from actual outcomes. Historical records of fires or logging contracts typically represent only the perimeter with little additional quantitative information (Johnson andWittwer 2008, Farris et al. 2010). These problems complicate the reconstruction of past disturbances and their lasting legacies on existing vegetation.
In cases where spatial information about historical disturbance is lacking, evidence can be obtained from tree rings. These include: tree establishment dates Wu 2005, Speer 2010), ring width variations caused by insect outbreaks (Lynch 2012), and fire scars (Margolis and Balmat 2009, Farris et al. 2010. Proxies such as firescar and tree-ring reconstructions can delimit where and when fires occurred, but they are difficult and expensive to acquire and are applicable only within a sampled area. More spatially explicit and practically easier to generate techniques are needed to help foresters and ecologists understand the existing community complexity of their forest landscapes. The goal of this research is to demonstrate the utility of local spatial statistics, called local indicators of spatial association (LISA; Anselin 1995), to discriminate variation and pattern at an individual tree scale. In the present study we use ALS-derived individual tree inventories ) to calculate the LISAs and delineate variation as it relates to spatial pattern. We further consider whether LISA relate to local-scale topographic variations and to recorded disturbances. Our objectives were to: (1) quantify variation in forest structure derived from ALS at an ecologically relevant scale, i.e., the individual tree; (2) map previously recorded but poorly defined disturbance locations; and (3) link variation in individual size and community composition to scalable indices of topographic wetness and topographic position.

MATERIALS AND METHODS
A large body of literature has shown trees in semi-arid forests tend to cluster in space by size (Sánchez Meador et al. 2009, Larson andChurchill 2012). To examine these patterns of clustering we considered how climax tree size distributions arise in relation to topographic variation and whether these can be detected by LISAs. We also examined whether recorded disturbances were readily apparent by comparing LISAs to where historical maps indicated logging contracts were awarded and along boundaries such as fence lines. We compared two related LISA statistics: Anselin local Moran's I (Anselin 1995), and the Getis-Ord G i statistic (Getis and Franklin 1987, Getis and Ord 1992, Ord and Getis 1995, 2001, O'Sullivan and Unwin 2010. We chose to use LISA over other statistical measures of forest structure such as mean canopy height (Detto et al. 2013), canopy cover, or lacunarity (Plotnick et al. 1996, Frazer et al. 2005) because we were interested in defining variation of local neighborhoods at the scale of individual trees.

Study areas
Our exemplar sites are (1) the Monument Canyon Research Natural Area (RNA) in the Jemez Mountains, New Mexico, which has never been logged and did not incur catastrophic disturbances over the entire 20th century, and (2) the Pinaleño Mountains, Arizona, which was opportunistically logged in the mid-20th century (Fig. 1).
The Jemez Mountains are located west of Santa Fe, New Mexico at 35.88 N, 106.58 W (Fig. 1). Their most dominant feature is a collapsed volcanic caldera with a rim approximately 19 km wide. Monument Canyon RNA is located on the south rim of the caldera on the Santa Fe National Forest (N.F.) at 35.88 N, 106.68 W (Fig. 1) v www.esajournals.org and is comprised of two mesas: Cat Mesa and San Juan Mesa, bisected by a central canyon and several tributary canyons draining to the southwest. Elevations vary between 2,460 and 2,550 m above mean sea level (amsl) (Figs. 2 and 3).
The Pinaleño Mountains are located southwest of Safford, Arizona at 32.78 N, 109.98 W, on the Coronado N.F. (Fig. 1) and reach 3,267 m amsl. The Pinaleños are the tallest of the Sky Island Archipelago mountains, characterized by a complex of steeply sided canyons with relatively gentle uplands above 2,700 m amsl ( Fig. 1) (Whittaker and Niering 1975).

Records of disturbance
Both sites have well characterized land use histories (Falk et al. 2007, Jones 2009, Swetnam et al. 2009). The disturbance regimes are similar to v www.esajournals.org other semi-arid forests across the southwestern USA and northern Mexico (Foster et al. 1998, Allen et al. 1998, Swetnam et al. 2009, Farris et al. 2010. Fire history reconstructions for the Pinaleño (Grissino-Mayer et al. 1994, O'Connor et al. 2014a) and Jemez (Touchan et al. 1996, Falk 2004, Falk et al. 2007, Dewar 2011) indicate frequent (3-10 year return interval) fire regimes with primarily low-and mixed-severity effects in pine and dry mixed-conifer forest types. There is also evidence of some high-severity wind-and topography-driven crown fires burning smaller areas in these ranges over longer intervals at higher elevations (Margolis et al. 2007, 2011, Margolis and Balmat 2009, O'Connor et al. 2014a). Other less frequent events (relative to surface fires) include semi-to multi-decadal drought Betancourt 2010, Williams et al. 2012), and insect outbreaks (Swetnam and Lynch 1993, Negrón et al. 2008, Lynch 2009, O'Connor et al. 2014b).
Monument Canyon RNA was preserved from all logging and wildfire since its establishment in 1935 (Falk et al. 2007). It is one of only a few stands of ponderosa pine that have not been logged in the southwestern USA. Overstory trees in Monument Canyon exceed 300 years in age (Falk et al. 2007). Part of Monument Canyon recently underwent a mastication thinning treatment in which trees ,20 cm DBH were removed along the San Juan Mesa. Sites with .308 slope were not thinned. Prior to thinning Monument Canyon was considered one of the worst examples of P. ponderosa 'dog-hair thicket' in the western USA (Allen et al. 2002, Falk 2004. Historically the whole RNA burned either at once or episodically oscillated between the two mesas across years with a mean fire interval of ,4 years (Falk et al. 2007. The Pinaleños logging campaign was prevalent   (Persson et al. 2002, Falkowski et al. 2006 in MATLAB 2012b (Mathworks 2012) (Figs. 5-7). The precision and accuracy of the VLM technique are reported in detail by  and , who found that the VLM had 68% error for overstory trees, with most commission errors due to multiple-leader canopy trees, and omission errors due to trees being located directly beneath a taller overstory tree.

Local indicators of spatial association (LISA)
Global spatial autocorrelations developed from weighted matrices assume homogeneity of the region of interest yielding a single statistic for a study area or population, and implicitly assume no variation across space (Moran 1950). In nature this condition is rarely satisfied, suggesting the need for statistics that vary over space (Legendre 1993). When data are heterogeneous, local spatial autocorrelation varies over space at the scale of the data (Anselin 1995, O'Sullivan andUnwin 2010). A global statistic is thus likely to be significantly different from autocorrelation measured at the local scale (O'Sullivan and Unwin 2010). In general, almost any standard statistic can be decomposed from a global to local statistic (Cressie 1993, O'Sullivan andUnwin 2010). Local statistics require a search radius or a number of near neighbors from which to draw a statistical distribution; and these groups of individuals Statistical tests of pattern analysis typically begin with a null hypothesis of complete spatial randomness: Rejection requires the existence of a dispersed or a clustered distribution based on a z-score (the standard deviation [SD] about the mean) where the p-value of statistical significance is typically set at 0.05 or 0.01 (a 95% or 99% confidence interval [CI]). Positive z-scores indicate clustering of large values and negative zscores indicate clustering of small values (e.g., a z-score of 2.58 indicates at the 99% CI a tree with tall neighbors, and À2.58 indicates a tree with short neighbors). Statistics that evaluate the existence of clusters in spatial arrangements are indicators of spatial association and include Moran's I (Moran 1950) and the Getis-Ord global G statistic (Getis and Ord 1992) (see Appendix A for equations).
For our study, Moran's I represents a summation of matrix cross products of z-scores of neighbor trees and the reference tree for all sample points, and thus can be illustrated by LISA where Anselin's local Moran's I (Anselin 1995) is calculated for statistical associations of each i, j pair at multiple spatial scales (Anselin 1995, Zhang andZhang 2007) (Appendix A). Similarly, Ord and Getis (1995) produced a LISA statistic of G, called G i * (Appendix A), to compare local concentrations of high to low zscores in an attribute where each is the center of a three-dimensional kernel. We applied these LISA statistics to the ALS derived individual tree inventory based on tree height and neighbor distance.
For the Getis-Ord G i * we set inverse Euclidean distance as the weight of the default neighborhood with a search threshold to the maximum observed tree height (49 m) with the 'Hot-Spot Analysis' Tool in ArcGIS 10.1 Toolbox (Mitchell 2005, ESRI 2012 (Mitchell 2005, ESRI 2012). The I i equation was also weighted by the height of trees using the inverse distance method with a maximum distance of 49 m. Maximum neighborhood area was based on the assumption that radial roots extend maximally out to the length of the tree's height (Woolsey 1911

Topographic indices
We also explored whether the topographic wetness index (TWI, Kirkby 1979, Quinn et al. 1991) and topographic position index (TPI, Guisan et al. 1999, Weiss 2001, Wilson et al. 2007) correlate with maximum tree size and the two LISA statistics. The TPI is calculated as: where r scale is the size of the moving window, Elevation is the value of the cell above mean sea level, l focal is the average of all elevations within the window, a is the annulus, r inner is the radius of the inner circle, and r outer is the radius of the outer circle (Weiss 2001). We calculated TPI in Quantum GIS 2.2 (QGIS 2013) with the GDAL-DEM tool set to the inverse distance weighting method (Wilson et al. 2007, GDAL 2013. Because we were interested in ecologically relevant scales including differences in soil depth and soil moisture, we set the r scale distance to 250 m, which corresponded to the average length scale of ridge-to-valley topography of the two featured landscapes (Figs. 2-4). Index values relate to topographic position (i.e., ridge, slope, valleybottom, etc.) with the most negative values in valley-bottoms and most positive values on ridge tops (Figs. 2-4). The TWI for estimating soil moisture was first demonstrated in the ''TOPMODEL'' (Beven and Kirkby, 1979) and is defined as where a is the local upslope area draining v www.esajournals.org through a digital elevation model (DEM) pixel or contour length unit and tanb is the tangent of the local slope b. All TWI calculations were performed in Quantum GIS 2.2 (QGIS 2013) with the extension for System for Automated Geoscientific Analysis (SAGA, v. 2.0.8). Kopecký anď Cížková (2010) tested 11 variations of TWI flow routing algorithms versus soil moisture and vegetation species diversity, they recommend a multi-flow routing algorithm (Quinn et al. 1991(Quinn et al. , 1995 for characterizing variation in vegetation. A multi-flow routing algorithm using the deterministic infinity (D'; Tarboton 1997) method was executed in the SAGA 'Catchment Area (Recursive)' and 'Topographic Wetness Index' module (Bö hner and Selige 2006). A convergence value equal to 1.1 was used. All DEMs were 10 m resolution (Figs. 2 and 4), which helped account for the lateral redistribution of water to the roots around each point located tree.

Individual tree isolations from the airborne LiDAR
For the Monument Canyon RNA (289 ha), the variable area local maxima algorithm identified 199,219 trees from the 0.333 pixel canopy height models. In the Pinaleños (166 ha, ALS Tile ID: F8405), the variable area local maxima algorithm identified 116,006 trees from the 0.333 pixel canopy height models. Tree heights range from a minimum 2-m cut-off to 49.3 m.

LISA of tree and stand level characteristics
Both LISA statistics discriminated various Tree heights are colored blue to red (in meters); black boundary circles are based on a wavelet fit to each colored blob and represents estimated canopy diameter. Anselin local Moran's I i (bottom two panels) is derived from the local maxima and weighted by tree height. Tall trees with tall neighbors are colored red, short trees with short neighbors are colored blue.
v www.esajournals.org structural differences and found strong evidence of clustering throughout Monument Canyon and the Pinaleño (Figs. 5-7). The negatively valued zscores tended to characterize clustering within the disturbed areas (e.g., within logged and thinned treatments). Recently established standing cohorts had either negative Anselin local Moran's I i z-scores (colored blue) or non-significant z-scores (gray), most evident in the Pinaleño logging deck where an aspect reversal takes place (Fig. 6).
For Getis-Ord G i * trees with short neighbors had significantly negative z-scores (z , À2.58, p , 0.01) and trees with tall neighbors had significantly positive scores (z . 2.58, p , 0.01) ( Table 2). Trees with z-scores between À2.58 and À1.65 and between 1.65 and 2.58 were also significantly positive at the 90-95% CI. Z-scores were non-significant (À1.65 , z , 1.65) in neighborhoods with a wide variation in heights, suggesting those trees were randomly distributed. For example, within the Pinaleño logging deck the average z-score changed rapidly from negative (green) to positive (red), with a brief non-significant transition (yellow) between the short and tall trees (Fig. 7, right panel). For the Anselin local Moran's I i z-scores the results are proportionately equivalent (Figs. 5 and 6), with negative z values for small neighbor trees and positive z values for tall neighbor trees ( Table 2). The added cross-product score of Anselin local Moran's I i , whether the object tree is also tall or short, resulted in a different set of classification values. The additional information further revealed more complex variations within the clusters of trees (Figs. 5 and 6, Table 2).

LISA patterns of logging and disturbance legacies
Both LISA statistics distinguished a variety of coherent patterns that aligned readily with known disturbance histories: timber contracts, fenced boundaries, recent fuels reduction treatments, historical fire perimeters; and changes in topography: aspect reversals, valley/canyon bottoms vs. ridge/mesa tops (Figs. 2-7).
Most of Monument Canyon's central drainage and the San Juan Mesa (southeast corner) are dominated by old-growth ponderosa pine (Falk et al. 2007); these stands are clearly defined by positive Anselin local Moran's I i z-scores (Fig. 6, red disks). Variations are defined by (1) the fence line along the southern boundary and (2) excluded areas from a recent fuels treatment operation. On San Juan Mesa (southeast corner of Fig. 5) south of the boundary fence there are mainly tall trees with short neighbors (yellow and green) that represent unthinned areas, while north of the fence tall trees with tall neighbors (red) dominate.
In the Pinaleños stand patterns are likely the result of logging within the known logging contract boundaries, e.g., clear cutting on the west side of the drainage, and patch cutting on the east side of the drainage (Figs. 6 and 7). Fig. 7. Orthophotograph (NAIP 2007) for Pinaleño sample area with historical logging records originally hand drawn on a 1960s orthophoto, later converted into a polyline (top left panel). A closer inspection of a patch-cut in the NAIP image with VLM local maxima represented (top right panel). Getis-Ord G i * (bottom panels) spatial statistics of the same area: negative values (green) represents clusters of small trees, positive values (red) represent clusters of tall trees. The west slope (left side) was clear cut, while the east slope (right side) was patch cut.
v www.esajournals.org Currently, forest structure within the logged area varies from unharvested old-growth mixed-conifer along the main riparian drainage (Big Creek) to recent reproduction in the logged areas. These differences are apparent in the ALS data, but are not evident in the 2007 National Agricultural Imagery Program (NAIP) orthophotography (Fig. 7).
Logged areas in the Pinaleños are in close proximity to skidding roads (Figs. 6 and 7). Within the clear cut sections, Anselin local Moran's I i reveals older 'seed-trees' amongst the young trees (Fig. 6, left side of the catchment). This is an important difference versus the Getis-Ord G i * (Fig. 7) which does not discriminate amongst the individuals because G i * does not look at the cross products of the neighbor trees.
Maximum tree size and LISA vs. topographic indices TPI and TWI were calculated with the same set of parameters using 10 m DEM models for both the Monument Canyon RNA (Figs. 2 and 3) and the Pinaleños (Fig. 4). The TPI for each site ranged from À5.5 to 3.6 ( Figs. 2 and 4). Most of the trees in this analysis were located on TPI between À3 and 3.5. The TWI surface with D' reported a range of values from 3.5 to 14.5, with most values for trees between 4 and 14 (Figs. 2  and 4).
Tree height increases asymptomatically with increasing wetness index over part of the range of TWI (!4 to approximately 7), above these values maximum tree height remains constant (Fig. 8, right panel). The tallest trees are in the valley/canyon bottoms (;49 m height), and the shortest are on ridge tops (;20 m height), with a nearly linear decrease in height as TPI increases in value from TPI ! À1.8 to TPI 3 (e.g., Fig. 8, left panel). In the Monument Canyon RNA oldgrowth stands are continuous across canyon bottoms up to the mesa tops, and are of similar ages (Falk et al. 2007) (Fig. 5). These stands end abruptly along the southern border fence, corroborating records that old growth trees south of the fence were logged in the 20th century. Stands of the tallest neighbors, given by the largest positive z-scores, versus TPI and TWI in Monument Canyon RNA (Fig. 9, black line) and the Pinaleños (Fig. 10, black line) skewed toward lower topographic positions (more negative TPI) and wetter sites (TWI . 6) than the population average (gray line).

DISCUSSION
Our results clearly demonstrate topographic position and wetness correlate to a majority of the variation in the climax tree height/size distribution. The tallest trees on the lowest topographic position and highest wetness indices was inversely correlated to the shorter trees on Table 2. Average z-scores and the comparative types of stand conditions for both Getis-Ord G i * (column 1) and Anselin local Moran's I i (columns two and three). In the far right column the first letter identifies the ith target tree, the second letter is the type of neighborhood it is surrounded by, e.g., H-L is a tall tree with short neighbors, L-L is a short tree with short neighbors, NS is not significant (the neighborhood trees are heterogeneous in structure), etc. with other forest studies in temperate and tropical forests (McNab 1993, Detto et al. 2013 suggesting a globally common forest feature driven by productivity gradients. Distinct variations in forest structure were also shown to be related to known disturbances and are clearly recognized by both LISA statistics. Logged areas were clearly defined versus a priori undisturbed sites. When topographic factors are accounted for, any departures from an expected climax state can be readily diagnosed. Our data are in agreement with Sánchez Meador et al. (2009), andLarson andChurchill (2012) having significant clustering in tree size across space. For Monument Canyon RNA, an essentially natural forest with a well-studied age distribution and fire regime (Falk et al. 2007, significantly positive LISA z-scores are found in both the canyon bottoms and on top of the mesas reflecting the nearly continuous distribution of old growth trees. In the Pinaleños, the most positive LISA z-scores are found in unharvested areas, as defined outside the recorded logging contract boundaries. Clusters of the tallest trees (z . 40) were also associated with negative TPI values and larger TWI values at both sites . The relationship of height and z-score with TPI and TWI appear to have minimum thresholds with respect to sites supporting taller trees .

Disturbance legacies characterized by LISA
Disturbances have left a lasting legacy on the tree size distribution and clustering patterns of these southwestern forests. LISA statistics clearly defined differences in forest structure more than fifty years after the logging disturbance (Figs. 6 and 7). In the Pinaleños, two different types of logging (clear cut vs. patch cut) were expressed by similar LISA scores but with different patch sizes on either side of the catchment.
Both LISA statistics revealed variation in forest structure not readily apparent in the orthophotography (Figs. 6 and 7). Anselin local Moran's I i provided additional information compared to the Getis-Ord G i * because Anselin local Moran's I i incorporates matrix cross products of z-scores of neighbor trees and the reference tree. In this regard Anselin local Moran's I i may be more informative than Getis-Ord G i * for characterizing different kinds of forest assemblages. For example, the z-score of a single large tree standing amongst a young cohort can be overwhelmed by the distribution of its short neighbors (and vice Fig. 9. Scatter plots and density distribution of zscores for TWI vs TPI in Monument Canyon. The distribution of high z-scores, i.e., the tallest clusters of trees, suggests tall trees aggregate in more negative TPI values, i.e., the drainage basins between ridge tops. The distribution of z-scores by TWI closely resembles the sampling distribution of TWI suggesting autocorrelation in those data. Fig. 10. Scatter plots and density distributions of zscores for the Pinaleño; same scheme as Fig. 9. v www.esajournals.org versa) for the G i *. In stands with mixed age structure, such as in the Monument Canyon RNA where dense dog-hair thickets exist, Anselin local Moran's I i was a key predictor for discerning where a tall tree is now surrounded by a new cohort (Fig. 5). Getis-Ord G i * clearly delineates where significant changes occur at the stand level, which may be more important when the objective is to delineate larger features, as in the logging perimeter (Fig. 7). In Appendix B we provide an extension to the Getis-Ord G i * technique, in which a Voronoi tessellation is used to create polygons of similar valued cohorts, equivalent to those used in USFS Common Stand Exams (http://www.fs.fed. us/nrm/fsveg/).
We tested these spatial relationships in two forest areas with known disturbance histories to evaluate fine scale variability in the ALS data. Cross validation to the known disturbance histories such as the logging record in the Pinaleños revealed other significant impacts not previously documented, mapped, or quantified explicitly. For both the Arizona and New Mexico examples of forest structure, Anselin local Moran's I i and Getis-Ord G i * showed coherent differences consistent with both known landuse history (Figs. 6 and 7), and interestingly to topographic variation (Figs. 8-10). Areas logged with either clear cuts or selective patch cuts were also easily differentiated (Fig. 7). Further, Anselin local Moran's I i and Getis-Ord G i * revealed variation not visible in the most recent orthophotography (Fig. 7). The presence of old logging roads which are most visible in the LiDAR bare earth model provided further confirmation logging was the disturbance agent in those locations.
Identifying structural changes in vegetation as a response to topographic variation versus disturbance is important for accurate reconstruction of the past disturbance regime. In Monument Canyon RNA, insignificant or negative zscores were shown to relate to young stands along steep slopes and the mesa tops. Based on reconstructed evidence from tree rings (Falk et al. 2007) frequent low-severity fires historically spread across the entirety of Monument Canyon area until the late 19th century. It is plausible these fires were intensified by slope alignments as they spread up out of the canyon bottoms which led to the mortality of trees along those locations preventing larger trees from establishing prior to fire exclusion. These regeneration trees, especially the dog-hair thickets, throughout Monument Canyon RNA are one of the consequences of P. ponderosa recruitment in the absence of episodic wildfire (Allen et al. 2002, Falk et al. 2007).

Conclusions and proposed applications
These examples illustrate how ALS data and LISA could be used in cases where disturbance history is less clearly documented, for example to distinguish residual old-growth forest from areas that are recovering from disturbance or some other unmeasured physical limitation. Given this potential, LISA statistics could be used to characterize forested stands with respect to: 1. Delineating existing stand boundaries, e.g., for common stand exams (see Appendix B for example); 2. Improving on historical maps where only a perimeter for a fire, insect disturbance, or timber management plan was recorded; 3. Quantifying site-indices based on scaled topographic indices, i.e., TPI and TWI, and using those to determine the level of departure from potential climax vegetation size; 4. Quantifying impacts of recent wildfires and insect outbreaks on stand structure versus those derived from other remote sensing products, e.g., NDVI mapping from satellite or aerial orthophotography, for which structure is only an implied variable; and 5. Comparing and contrasting forest-level responses to disturbance; e.g., using LISAs to evaluate the proportionality of mortality due to drought and/or an insect outbreak across different stands by topographic index variations.
Incorporating ALS data across landscapes and analyzing spatial properties of the resulting forest structure with LISA enables us to ask new questions about landscape patterns and processes. The forests in this study both have long histories of episodic disturbance, punctuated by logging and fire suppression. Where disturbance was known to have occurred LISAs showed skill at identifying spatial heterogeneity in the stand structure. These same methods can be applied to other forests, where an evaluation of current conditions relative to their historical range of variability is needed. LISAs can therefore be used to quantify the structural legacies of a range of natural and human disturbances (fire exclusion, high-severity wildfire, or logging). The methods demonstrated here indicate the potential to characterize semi-arid forest structure in the Southwestern USA and potentially other semi-arid forests globally.

APPENDIX B THIESSEN POLYGONS FOR COMMON STAND DELINEATION
Commonly shared characteristics of trees allow them to be aggregated into 'stands' of like-valued individuals based on (1) species, (2) local topographic features, (3) soil type, or (4) some combination of parameters. We were interested in delineating the pattern of the existing forest structure to create similarly classified stands.
We generated Thiessen polygons (Brassel and Reif 1979; also called Voronoi polygons [Boots 1986]) in ArcGIS 10.1 (ESRI 2012) about the individual points, which include the area of influence about a point relative to its neighbors. Thiessen polygons can be generated either before or after the Hot Spot Analysis (ESRI 2012), as the data are equivalently vector based (point or polygon).
We calculated the Getis-Ord G i * statistic to determine the break point in stand architecture and spatially connected the neighborhoods of Thiessen polygons. We merged the individual Thiessen polygons into a single polygon layer, which became the boundary of our common stand based on structural characteristics (Fig. B1).
In order to define common stand boundaries v www.esajournals.org we re-classified the continuous data using predefined sets of break-lines (e.g., manual, equal interval, quantile, standard deviation, etc., in ArcGIS 10.1 [ESRI 2012]). In our example here we used standard deviation. Once each break point was chosen we used the 'Select Tool' to create sets of polygons either greater than or less than the z-scores at a specified significance level. We chose to segment the example in Fig. B1 at the 95% level.