Predictive Ecosystem Mapping of South-Eastern Australian Temperate Forests Using Lidar-Derived Structural Profiles and Species Distribution Models

Abstract

Modern approaches to predictive ecosystem mapping (PEM) have not thoroughly explored the use of ‘characteristic’ gradients, which describe vegetation structure (e.g., light detection and ranging (lidar)-derived structural profiles). In this study, we apply a PEM approach by classifying the dominant stand types within the Central Highlands region of south-eastern Australia using both lidar and species distribution models (SDMs). Similarity percentages analysis (SIMPER) was applied to comprehensive floristic surveys to identify five species which best separated stand types. The predicted distributions of these species, modelled using random forests with environmental (i.e., climate, topography) and optical characteristic gradients (Landsat-derived seasonal fractional cover), provided an ecological basis for refining stand type classifications based only on lidar-derived structural profiles. The resulting PEM model represents the first continuous distribution map of stand types across the study region that delineates ecotone stands, which are seral communities comprised of species typical of both rainforest and eucalypt forests. The spatial variability of vegetation structure incorporated into the PEM model suggests that many stand types are not as continuous in cover as represented by current ecological vegetation class distributions that describe the region. Improved PEM models can facilitate sustainable forest management, enhanced forest monitoring, and informed decision making at landscape scales.

1. Introduction

Understanding the composition and structure of forest ecosystems is critical for ecological understanding, and for developing effective forest management strategies. The increased availability of geographic information and remote sensing data for forested ecosystems, in combination with advances in computational modelling, has facilitated the development of predictive ecosystem mapping (PEM) models [1,2,3]. PEM models are defined as methods that identify ecological-landscape relationships from spatial environmental data and field observations (as available) to predict vegetation composition across a landscape [1,2,4,5,6,7]. Using similar approaches but with varying objectives and units, PEM models have been broadly defined as predictive vegetation mapping [1], which also encompass species distribution models (SDMs) [2,3,8,9], bioclimatic envelope models [10,11], habitat suitability or decision-support models [12,13,14], and ecological niche models [15,16]. Such models originated from ecological niche theory whereby vegetation distribution is predicted using variables that either correlate with or define tolerance ranges of species [1,8,9,10,17,18].

There are three typical approaches to modelling community distributions from species and environmental data: (1) assemble species into communities and then predict, (2) predict species individually then assemble into communities, and (3) assemble and predict species together [19,20,21,22,23]. The first approach uses some form of classification, ordination, or aggregation method to generate communities from individual species survey data. The presence-absence locations of each community are then compared with environmental predictors to describe their distributions. The second approach models each species distribution individually as a function of available environmental predictors. A community-level output is then generated from a classification, ordination, or aggregation of all the extrapolated individual species distribution models. The third approach characterizes each community and associated relationships with environmental data in a single process incorporating each species and all environmental data simultaneously. No single approach will be optimal in all circumstances, with the selection dependent on available data and study objectives [21,22]. Ohmann et al. [24] argued that any approach selected should aim to retain information of individual species and patterns of co-occurrence in the final predictions. Joint species distribution modeling (JSDM) is an approach that accounts for both presence-absence and abundance data jointly across species to classify communities [25,26]. Combining independent SDM model outputs may be just as effective as JSDM though where the species of interest do not interact with one another [25]. These techniques all aim to develop predictive models that utilize multidimensional descriptive variables to predict the distribution of vegetation across landscapes.

The use of multiple data types to describe environmental variables and forest structure can result in complex, often non-linear data, which challenges the accuracy of predictive models developed using traditional parametric approaches [27,28,29]. Machine learning techniques have emerged from synergies between computer science and the identification of ecological processes and patterns, known as ecological informatics, which directly address issues of data complexity and non-linearity to generate accurate predictive models [28,29,30]. Some examples of machine learning techniques that have been used for species distribution modelling include artificial neural networks, k-nearest neighbors, support vector machines, boosted regression trees, and random forests [6,8,27,31,32,33,34]. The random forest algorithm, developed by Breiman [35], has two key advantages for forest distribution mapping; (1) the novel variable importance measure, and (2) the proximity measures of similarity among data points [27].

Many predictor variables for PEM models can be measured in the field or collated from various geographic information system layers and remote sensing data [3]. Selection of predictor variables is critical to PEM accuracy [36,37,38,39], with each variable considered in terms of its importance to the distribution of an ecosystem and the ecological basis for its inclusion [18,37]. These variables can include a number of environmental gradients: direct gradients influence growth but are not consumed (i.e., temperature, pH); indirect gradients have a location-dependent correlation with forest distribution but have no functional effect on plant growth (i.e., latitude, altitude, slope); and resource gradients are consumed by vegetation (i.e., light, water) [2,5,18,36]. Interactions between structural and compositional gradients can influence stand development and are important to capture in ecosystem classification [40]. These gradients of structure and composition describe the state of vegetation at a specific point in time (i.e., height, cover, strata), and can be broadly defined as ‘characteristic gradients’. Characteristic gradients are increasingly being included in PEM models with inputs summarizing forest structural complexity into a single metric or multiple metrics when three-dimensional data are available [12,41,42].

Characteristic gradients can be quantified at broad spatial scales using remote sensing compared to spatially limited and cost-intensive field observations. Systematic field observations are often individual locations spread across large areas with limited sensitivity to the vertical assemblage of vegetation. Both horizontal and vertical structure can be obtained across a spatial continuum using light detection and ranging (lidar) technology. Lidar is an active remote sensing technique that uses laser ranging to determine the distance between the sensor and an object. Distance is calculated from half the time-lapse between when the laser pulse is emitted and the detection of the returned pulse after striking an object [14,43,44]. Key sensitivity and accuracy limitations of commonly used optical and radar imagery [45] can be overcome using lidar data [44]. Lidar pulses measure the presence or absence of foliage and topographic structure in three dimensions [46,47,48,49], where optical and radar data responses are generally canopy and mid-canopy driven, respectively [43]. Results from Kane et al. [50] highlight that lidar can accurately characterize forest successional stages in the absence of field measurements. The ability of lidar to exploit small gaps in the forest canopy and capture understorey structure allows for the inclusion of lower strata attributes in predictive vegetation models [51,52,53,54]. Predictive models that include characteristic variables from lidar often use discrete derivatives such as canopy height, foliage cover, leaf area index (LAI), and vegetation density [12,29,50,51,55,56,57,58,59,60]. Recent studies using the vertical continuum of lidar data have typically focused on structure only classifications using random forest [49,61], excluding the use of ecological data pertaining to species composition dominating landscape-scale PEM.

In south-eastern Australia, eucalypt forests have traditionally been delineated based on dominant Eucalyptus species without consideration for associated Acacia, rainforest and understory species [62]. The classification of rainforest focuses on the floristic composition at maturity and excludes changes in species composition over time due to changes in the composition that may be driven by disturbance (i.e., ingress of Acacia and Eucalyptus following fire) [62]. Cameron [62] argued that, while this traditional approach has advantages for interpreting community niche distributions, it is problematic at the ecotones that share dominant species, or where dominant canopy species are replaced by species characteristic of another stand type (e.g., patches of eucalypts within a rainforest mosaic after fire). The potential need to consider both climatic and structural gradients in PEM models for this region was further highlighted in a recent study that indicated both types of variables were important for explaining beta-diversity in the temperate forests of south-eastern Australia [63]. Accordingly, the cool temperature forest landscapes of the Central Highlands region in south-eastern Australia are an ideal case study for examining the potential to improve traditional PEM approaches using structural information such as characteristic gradients. The landscape contains wet sclerophyll forest, rainforest and an ecotone between these two stand types described as ‘cool temperate mixed forest’ [62], which has so far been excluded from the principal vegetation class system in the region (‘Ecological Vegetation Class’, EVC) [64,65].

This study seeks to examine the utility of combining recent developments in stand type predictions, using characteristic gradients in lidar structural attributes with complementary sources of ecologically meaningful information. The previously developed lidar-derived stand type map [61] demonstrates the assemble-and-predict approach, which evaluates characteristic gradients (i.e., the vertical profile of vegetation density) across the landscape to predict the distribution of rainforest, ecotone and eucalypt stands. The distribution of species which best characterize stand types in the region were predicted individually and subsequently assembled into community groups in order to provide important ecological context. The lidar-derived and SDM-derived stand type classifications, in addition to auxiliary criteria which further stratify these stand types (i.e., specific individual species, elevation, canopy height) were combined using a series of decision trees to generate a highly detailed PEM model of the study region. This approach to PEM demonstrates a hierarchical, stepwise approach to combining lidar-derived stand type classifications with auxiliary information which provide important ecological context for predicting the spatial distribution of ecological communities. The PEM model demonstrates the potential for defining the spatial distribution of ecological communities at a much higher resolution than existing EVC-based classifications in the study region.

2. Materials and Methods

2.1. Study Area

The Central Highlands region of south-eastern Australia is located in the state of Victoria approximately 100 km north-east of Melbourne (Figure 1). Three spatially contiguous stand types of varying structure dominate the region’s forests: cool temperate rainforest (herein referred to as ‘rainforest’); wet sclerophyll eucalypt forests (‘eucalypt’ forest); and the ‘ecotone’ where rainforest and eucalypt forest species mix [66]. Rainforest stands are dominated by 30–40 m tall Nothofagus cunninghamii, Acacia melanoxylon and Acacia dealbata with understorey species, including up to 20 m tall Atherosperma moschatum. Tree ferns Cyathea australis and Dicksonia antarctica dominate the lower 5 m strata with ground ferns at less than 1 m. Rainforest stands may occasionally include emergent 60 m tall Eucalyptus regnans. Eucalypt stands are dominated by up to 80 m tall E. regnans and occasionally Acacia-dominated understories reaching 40 m in height. Smaller trees and single-stemmed tall shrubs frequently occupy the lower 20 m of these forests. Ecotone forests occur in areas between distinct rainforest and eucalypt stands when eucalypts are co-dominant with rainforest trees in the overstorey above a rainforest understorey. Cameron [62] describes this ‘ecotone stand type’ as ‘cool temperate mixed forest’; a seral community where fire or other recurring disturbances prevent secondary succession to eucalypt-free rainforest.

Other stand types that may be dispersed within or surrounding rainforest, eucalypt and ecotone forests, include a range of eucalypt-dominated forests. Dry sclerophyll forests occur at lower elevations (up to approximately 900 m) and on warm and dry sites (elevated ridge lines and exposed northerly aspects). These dry eucalypt forest canopies are typically dominated by a mix of stringybark, peppermint, and gum-barked eucalypt species, excluding E. regnans and/or E. delegatensis. Similar in dominant overstorey species to wet forest stands, damp eucalypt forest occurs from 200 to 1000 m in elevation on warm, moist sites and contain local abundances of E. obliqua and E. regnans with both species replaced at higher elevations (approximately > 1000 m) by E. delegatensis. Areas of shrubby wet forest are similar to damp eucalypt forest but more consistently dominated by E. delegatensis and a dense understorey of tall shrubs, including Polyscias sambucifolia and Oleria phlogopappa. Montane eucalypt forests include both wet and damp forests, typically at elevations greater than 1000 m, often with a dominance of E. delegatensis that replaces other lower elevation eucalypt species. Sub-alpine woodlands occur at elevations above 1200 m, frequently receive snow, and are dominated by E. pauciflora. Riparian eucalypt forests occur as narrow buffers around inland waterways at lower elevations and occur as patches amongst rainforest at their ecotone. Riparian stands commonly contain E. regnans, E. viminalis, A. melanoxylon and a dense assemblage of understorey vegetation.

2.2. Species Composition of Stand Types

Ecological vegetation classes (EVCs) are used across Victoria to classify vegetation types in different bioregions, defined by climate, geomorphology, geology, soils and vegetation [64,65,67]. EVCs are characterized by floristics, lifeforms and ecological characteristics, with an inferred association to distinct environmental conditions [64,65,67]. Many of these characteristics were either from previously modelled data, derived from Landsat imagery, or from floristic information collected in 2004 and 2005 from an unspecified number of field plots across Victoria [65]. Multiple EVCs correspond to the stand types in this study, including cool temperate rainforest, wet forest, damp forest, shrubby wet forest, riparian forest, montane wet forest, montane damp forest, and sub-alpine woodlands. To facilitate analysis, the wet forest, damp forest, and shrubby wet forest were considered eucalypt stands, and the montane wet and montane damp forest classes were considered as montane eucalypt stands. Cool temperate rainforest, riparian forest, and sub-alpine woodlands retain their classification. There is currently no distinct EVC for ecotone stands. Benchmark lists of species were identified as ‘typical’ to each EVC [68], and were filtered to include only tree, shrub, and tree fern lifeforms.

2.3. Floristic Sampling

Floristic surveys were conducted in 71 field plots (after excluding burnt/logged/observably disturbed plots) that represented our three main stand types (rainforest, ecotone, and eucalypt) from four field campaigns in the study area between April 2011 and October 2014. The presence or absence of all native tree, shrub, and tree fern species were recorded across sites over a broad range of environmental conditions for each stand type (see Section 2.1 Study Area) in fixed-area plots of sizes 20 × 20, 30 × 30, or 10 × 40 m² (Table S1). Differences in species composition (dependent variable) among stand types (fixed factor) were tested using permutational ANOVA (PERMANOVA), using PERMANOVA+ Version 1.0.6 (PRIMER-E Ltd., Plymouth, UK) [69]. Species composition was based on the Bray-Curtis similarity matrix of the presence or absence of all tree, shrub, and tree fern species. Significance testing of the Bray-Curtis similarity measures and post hoc comparisons were made using 9999 permutations. Permuted residuals were calculated using the Type III (partial) partitioned sum of squares [69]. The similarity percentages (SIMPER) procedure of PRIMER (Version 6.1.15, PRIMER-E Ltd., Plymouth, UK) [70] was used to identify those species that contributed up to 90% similarity within stand types (Table S2) and also those species contributing most to differences among the three main stand types (Table S3). The results of the SIMPER analysis were used to differentiate stand types based on the presence or absence of key species. Species that had an average abundance value of <0.2 (occurs less than 20% of the time) in at least one stand type and also present on the EVC species lists were identified as species which best separate stand types. An exception was made to include E. delegatensis in future analyses as it exists in montane eucalypt forest and typically occurs at elevations between 500 m and 1200 m above sea level [71,72] but was under-represented in the field sampling.

2.4. Lidar-Derived Digital Elevation Model (DEM)

Airborne lidar data were commissioned by the state environment department and were acquired from late 2007 to early 2008 over a subset of the Central Highlands region (Figure 1) from a fixed wing aircraft (Data agreement number DQ201206071138; see

Table 1. Lidar data flight acquisition details and sensor configurations.

Lidar System Configurations
Acquisition date range	19 November 2007–0 January 2008
Sensor type	OptechALTM3100EA
Wavelength (nm)	1064
Scan rate (kHz)	71
Scan angle (°)	±25
Mean footprint size (m)	0.26
Pulses (m2)	0.90
Maximum returned signals	4
Horizontal accuracy (cm)	±35
Vertical accuracy (cm)	±50

for flight details and sensor settings). There should be limited impact of the date discrepancy between the lidar capture and floristic sampling due to the slow change in forest structure of these mature stand types (>100 years old) and the exclusion of disturbed field sites in the study [45,51,73,74]. This is further supported by a previous study where a time lag of six years between field-data collection and lidar acquisition only had a minimal influence on the predicted distribution of bird species modelled using structural variables in an undisturbed coniferous forest [74]. All mapping in this study was constrained within the lidar coverage. Data were provided in 2 km × 2 km tiles, including a ‘bare-earth’ digital elevation model (DEM) at 1 m resolution that was generated by the data provider using the ground-classified points. Cells were aggregated and the mean value was used to generate a 20 m resolution DEM for each tile. The reported DEM vertical (±50 cm) and horizontal accuracies (±35 cm) were the same as the lidar data (Table 1). Additional assessments of vertical accuracy were conducted using differential GPS coordinates of six field plots in each of the rainforest, ecotone, and eucalypt stand types of varying vegetation density (total of 18 point locations). Differential GPS (DGPS) coordinates were acquired using a Leica GNSS Viva GS10 with an AS10 antenna and post-processed with ground-based reference stations. The overall root mean square error (RMSE) of the DEM compared to the DGPS recorded elevation at the 25 locations was 0.49 m.

2.5. Environmental and Satellite Data

A series of topographic, climate and satellite-derived characteristic variables were considered for modelling individual species distributions (

Table 2. Environmental and satellite variables for the Central Highlands region (bold variables were included as potential predictors in species distribution models, non-bold variables were excluded based on strong correlations, r ≥│0.7│, with included variables). Maximum, minimum and mean values are provided for each variable within the lidar acquisition boundary.

Variable	Minimum	Maximum	Mean
Environmental—Topography
Elevation (m above sea level)	103.9	1567.9	652.0
SlopeA (°)	0	68.1	17.1
Topographic Position IndexA (TPI; 100 m radius)	−8.6	14.2	0
Heat load indexB (HLI)	0.2	1.0	0.8
Potential direct incident radiationB (PDIR; MJ cm−2 yr−1)	0.2	1.0	0.8
Proximity to waterwaysC (m)	0.0	6882.7	1817.5
Environmental - Climate
Annual mean temperature (BIO1; °C)	5.9	14.9	11.6
Mean diurnal range (BIO2; °C)	5.1	13.4	9.2
Isothermality (BIO3; %)	28.8	51.8	43.4
Temperature seasonality (BIO4; standard deviation; °C)	3.7	5.0	4.3
Maximum temperature of warmest month (BIO5; °C)	16.2	29.1	23.8
Minimum temperature of coldest month (BIO6; °C)	−2.4	5.7	2.6
Temperature annual range (BIO7; °C)	16.7	27.0	21.2
Mean temperature of wettest quarter (BIO8; °C)	1.6	13.5	6.8
Mean temperature of driest quarter (BIO9; °C)	11.1	19.8	16.7
Mean temperature of warmest quarter (BIO10; °C)	11.1	20.8	16.9
Mean temperature of coldest quarter (BIO11; °C)	0.7	9.5	6.3
Annual mean precipitation (BIO12; mm)	667.1	2089.5	1290.3
Precipitation of wettest month (BIO13; mm)	70.7	274.9	164.9
Precipitation of driest month (BIO14; mm)	31.8	84.2	53.0
Precipitation seasonality (BIO15; coefficient of variation; %)	14.8	47.2	31.9
Precipitation of wettest quarter (BIO16; mm)	202.3	745.4	454.3
Precipitation of driest quarter (BIO17; mm)	109.7	324.3	194.7
Precipitation of warmest quarter (BIO18; mm)	110.5	324.3	198.5
Precipitation of coldest quarter (BIO19; mm)	199.5	690.0	446.6
Annual mean vapour pressure deficit (hPa)	2.0	6.2	4.0
Summer mean vapour pressure deficit (hPa)	3.6	10.7	6.9
Annual heat moisture indexD (AHMI)	7.7	35.9	17.8
Satellite
Bare ground fraction (% + 100)	97.7	194.1	103.5
Green vegetation fraction (% + 100)	98.0	203.1	169.0
Non-green vegetation fraction (% + 100)	95.9	202.0	126.0

^A Calculated using SAGA GIS (Version 4.0.1) [79]. ^B Calculated using Hyperniche (Version 2.2) [80]. ^C Euclidean distance. ^D AHMI = (BIO1 + 10)/(BIO12/1000) [81].

). All topographic variables were generated using the 20 m lidar-derived DEM. Seasonal fractional cover products (bare ground, green vegetation, non-green vegetation) derived from Landsat imagery (30 m cell resolution) and acquired between December 2007 and February 2008 [75,76,77] were used as satellite-derived characteristic variables. The combination of the environmental and satellite data at varying resolutions reflects the scale of phenomena at which species distributions are driven or constrained by each variable [3,78]. Each of the satellite-derived characteristic variables and climate variables were resampled to 20 m resolution using bilinear and nearest neighbor interpolation, respectively, for all data analysis and spatial predictions.

Climate variables were developed using a combination of pre-existing and interpolated mean monthly datasets (250 m cell resolution) centered on the 1981–2010 reference period. A series of 19 bioclimatic indices commonly used in SDM analyses were computed using the dismo package [82] in R (Version 3.3.0) [83]. The annual and summer (from December to February) mean vapor pressure deficit (VPD), and an annual heat moisture index (AHMI) were also calculated using the available climate data. Maximum and minimum temperature layers were sourced from a previous study [84], which used a topographic index of relative elevation (minimum temperature only) and standardized moderate resolution imaging spectroradiometer (MODIS) land surface temperature (LST) data to improve interpolation performance. Precipitation and vapor pressure were interpolated using ANUSPLIN 4.4 [85,86] following Stewart et al. [87], with the addition of a regression procedure [84,88] to minimize bias in the monthly climate normals prior to model calibration. Pooled cross-validation statistics for mean monthly precipitation and vapor pressure are presented in Table S4. The mean monthly vapor pressure deficit was then computed by subtracting vapor pressure from the saturation vapor pressure. Saturation vapor pressure (SVP) was calculated using temperature data in conjunction with the August-Roche-Magnus equation [89,90]:

$$SVP=6.1094\times exp{ }^{\frac{17.625\times T}{243.04+T}}$$

(1)

where T is the mean monthly temperature.

Correlation matrices were used to detect collinearity and to select a reduced set of predictor variables for SDMs (Table S5). Variables were considered as highly correlated where r ≥ |0.7| [91]. Three climate metrics (annual heat-moisture index [81], temperature seasonality [92], temperature annual range) were selected after excluding highly correlated variables. The majority of climate variables were highly correlated with one another, including the mean annual temperature and the annual precipitation despite smoothing to account for the lack of short-range (5–10 km) correlation between precipitation and elevation [93,94]. As an alternative, the annual heat moisture index was used as it reflects moisture availability as a function of both temperature and precipitation [81]. The three climate variables selected were considered suitable for species distribution modelling, as they reflect both mean climates and seasonal variability, quantifying important components of the ecological niche. Correlation analysis of the optical satellite characteristic gradients showed a strong correlation between the green vegetation and non-green vegetation fractions (r = −0.98). Only the non-green vegetation fraction and bare ground fraction were retained as they were the only pairwise combination below the correlation threshold (r = 0.56). As each of these three cover components (approximately) sum to one, the green vegetation fraction can otherwise be inferred by the bare ground and non-green vegetation cover and therefore was unlikely to have deleterious effects on model performance. The final set of predictor variables included seven environmental variables and two satellite-derived fractional cover products (Table 2).

2.6. Species Distribution Modelling

SDMs for each of the five species which best differentiated stand types were fit using presence and absence records obtained across the Central Highlands region with the randomForest package [95,96] in R. Records were collated from a collaborative network of field plots among the listed authors, (including the 71 sites used for the SIMPER analysis) in addition to project-specific survey data obtained from the Victorian Biodiversity Atlas [97]. Each of the plots used for the analysis maintained a minimum separation distance of at least two pixels (at 20 m resolution), ensuring that no replicates were present in the complete dataset. A six-fold, spatially blocked cross-validation strategy was applied using the blockCV package in R [98] in order to minimize the effect of spatial autocorrelation on the model evaluation. Autocorrelation range in the prediction surfaces was evaluated by fitting variograms to subsampled points across the landscape. A median autocorrelation range of approximately 25 km was selected based on the variogram analysis and was used to stratify the plot data into six sets of training and testing data. Two random forest model variants were run for each species; one using environmental variables only and the other using environmental and satellite variables combined. The cross-validated predictive performance of each model was assessed using the area under the receiver operating characteristic (AUC), overall accuracy (%), sensitivity (%; true positive rate), specificity (%; true negative rate) and true skill statistic (TSS; sensitivity + specificity − 1). The optimal threshold for binary classification (presence or absence) of species distributions was selected by maximizing the sum of sensitivity and specificity of cross-validated predictions [99]. The final set of predictor variables used to map the distribution of each species was selected based on the overall accuracy of the classified cross-validated predictions. The spatial distributions of each species were calculated as the average value predicted by models generated during cross-validation and were classified into binary maps using the optimal threshold.

2.7. Predictive Ecosystem Mapping Model

The PEM model was developed using a combination of SDM-derived (landscape species model) and lidar-derived (landscape structure model) stand type classifications (rainforest, ecotone or eucalypt), individual SDMs, elevation, and a canopy height model (CHM) to identify non-forest areas. The landscape species model was developed using individual species distributions to minimize model complexity while maximizing the number of species presence and absence records used to support stand type classification. Structured survey records could therefore be sourced from ecological surveys with varying objectives and locations, which would not otherwise have been possible with a JSDM approach. We intentionally excluded the full vertical profile of vegetation from the modelling of individual species distributions as they reflect a complex assemblage of multiple species. Each modelled species had between 316 and 494 combined presence and absence records sourced from a selection of different field campaigns conducted in the region. The landscape species model was generated with the predicted distributions of species which best differentiated stand types, guided by a decision tree developed using the SIMPER analysis results (Figure 2). The presence or absence of E. regnans or E. delegatensis was used to differentiate between ecotone and rainforest stands in the landscape species model.

The landscape structure model refers to a specific landscape-scale map derived using structurally based characteristic gradients derived from the full vertical profile of lidar data [61]. The landscape structure model, previously developed by Fedrigo et al. [61] for the same study extent, provides a detailed delineation of forest classes into rainforest, ecotone and eucalypt stands, specifically recognizing the coalescing of stand types at the landscape scale. The model uses known stand types in conjunction with the first 10 principal components of 1 m strata plant area volume density (PAVD) profiles to predict rainforest probability across the landscape. The retained principal components were most sensitive to PAVD for strata between 1 and 10 m. The highest rainforest probabilities were associated with high PAVD in these strata. The classification was then performed by selecting probability thresholds, which were used to separate each stand type [61]. The landscape structure model achieved an overall accuracy of 83.8% and achieved a Cohen’s kappa coefficient of 0.62 during cross-validation, indicating substantial model agreement. Stand type rainforest probability thresholds were adjusted from Fedrigo et al. [61] to ranges that were predicted for field plots only. The rainforest probability range thresholds for eucalypt, ecotone and rainforest stands were <18%, 18% to 70%, and ≥70%, respectively. New thresholds were selected to eliminate potential bias from modelled EVC classifications, which were used to identify thresholds in Fedrigo et al. [61]. The new thresholds resulted in substantial model agreement with the same overall accuracy of 83.8% and a slightly increased Cohen’s kappa coefficient of 0.63. The species distribution modelling and hierarchical decision trees were designed to provide an ecological basis for constraining and further stratifying the lidar-derived stand type maps.

The landscape species model was combined with the landscape structure model to find areas of agreement based on a second decision tree (Figure 3). Congruence between the landscape species and landscape structure models was indicative of a higher degree of confidence in the modelled predictions. E. regnans, E. delegatensis, and elevation were used to separate montane eucalypt, sub-alpine woodland, and other forest types (inclusive of dry forest) from the broader eucalypt class. The CHM was used to filter out non-forest regions and was generated from the lidar data by identifying the maximum height of returns for each cell. Areas of recent logging, within the last 10 years before lidar acquisition, were used to identify a CHM threshold to identify areas of non-forest. A threshold of 6 m was defined as the first quartile CHM value in recently logged areas. All areas on the landscape with a CHM value of less than 6 m were classified broadly as non-forest.

A summary of the preceding steps, and the workflow for developing the final PEM model is illustrated in Figure 4.

The final PEM model was compared against the EVC distributions by summarizing the predicted spatial coverage (km²) of each stand type. The potential distribution of the ecotone stand type amongst the existing EVCs was also examined, as ecotones were not previously defined in the EVCs. The total area of rainforest, ecotone, and eucalypt stands predicted by the landscape species model and the landscape structure model were also compared to highlight how each of these preliminary classifications constrained one another. Total area and percentage of the total ecotone area was determined for each EVC. All comparisons were made separately for two 5 km × 5 km extents (Figure 1) and the full study area.

3. Results

3.1. Species Composition of Stand Types

A total of 19 species occurred in both the collated EVC list and the field-based list of species that were identified as contributing up to 90% of the similarity within each stand type (Table S2). The overall species composition differed significantly among stand types (P ≤ 0.001). Four species, including A. dealbata, A. moschatum, N. cunninghamii, and P. bicolor, were identified as those which were best able to separate stand types based on the SIMPER analysis (Table S3). E. regnans and E. delegatensis were also included to separate ecotone from rainforest stand types in the landscape species model. Rainforest stands were best separated from ecotone and eucalypt stands by the presence of A. moschatum and P. bicolor and the absence of A. dealbata (Table S3). The presence of N. cunninghamii also distinguished rainforest from eucalypt stands. Ecotone stands were best separated from the other stand types by the presence of N. cunninghamii and A. dealbata. While A. dealbata was identified as important in separating rainforest from ecotone and eucalypt stands, it was excluded from the categorization of stand types due to its ecological niche covering the majority of the study area.

3.2. Species Distribution Modelling

Cross-validated error statistics for the five SDMs used to generate the landscape species model and final PEM model are shown in

Table 3. Cross-validated (spatially blocked across six folds) random forest species distribution models for species used to delineate stand types in the Central Highlands region. Model variants selected for predictive ecosystem mapping are bolded for each species.

Species	Predictors	P	A	T	TPR	TNR	OA	TSS	AUC
Atherosperma moschatum	Environmental	38	423	0.22	0.84	0.80	0.80	0.64	0.87
	Environmental and Satellite	38	423	0.14	0.82	0.78	0.79	0.60	0.83
Eucalyptus delegatensis	Environmental	57	344	0.24	0.98	0.81	0.84	0.79	0.92
	Environmental and Satellite	57	344	0.41	0.96	0.88	0.89	0.84	0.92
Eucalyptus regnans	Environmental	92	273	0.25	0.86	0.73	0.76	0.58	0.81
	Environmental and Satellite	92	273	0.20	0.99	0.72	0.79	0.71	0.85
Nothofagus cunninghamii	Environmental	129	365	0.15	0.95	0.64	0.72	0.59	0.83
	Environmental and Satellite	129	365	0.21	0.90	0.68	0.74	0.58	0.84
Pittosporum bicolor	Environmental	52	264	0.16	0.88	0.66	0.70	0.54	0.81
	Environmental and Satellite	52	264	0.19	0.87	0.72	0.75	0.59	0.80

Environmental predictors = topography: slope, topographic position index, potential direct incident radiation, proximity to waterways, climate: temperature seasonality, temperature annual range, annual heat moisture index; Satellite predictors = bare ground fraction, non-green vegetation fraction (see Table 2); P = number of sites where species was present; A = number of sites where species was absent; T = mean threshold for binary classification, calculated by maximizing the sum of sensitivity and specificity for each cross-validation fold; TPR = true positive rate; TNR = true negative rate; OA = overall accuracy; TSS = true skill statistic, calculated as TPR + TNR − 1; AUC = mean area under the receiver operating characteristic, calculated for each cross-validation fold.

. Models using environmental and satellite variables achieved a 2% or greater improvement in overall accuracy during cross-validation relative to environmental variables only for E. delegatensis, E. regnans, N. cunninghamii and P. bicolor. Alternative metrics generally supported the selection of the environmental and satellite models based on the trade-off between the true positive and negative rates, and AUC scores (although AUC for P. bicolor was slightly higher for environmental variables). The cross-validation performance was higher across all metrics when using environmental variables for A. moschatum. The annual heat moisture index and fraction of non-green vegetation were consistently among the most important variables across cross-validation folds (see Figure S1 for variable importance plots). The predicted distributions of each species for a subset of the study area are mapped alongside canopy heights in Figure 5.

3.3. Predictive Ecosystem Mapping and Comparison With Ecological Vegetation Classes

The landscape species model and the landscape structure model both predicted much larger areas of rainforest and eucalypt stands in comparison to the PEM model (

Table 4. Area of each stand type (km²) as predicted by the landscape species model and the landscape structure model across the landscape and for two 5 km × 5 km subset extents (see Figure 1 for extents). Values in parenthesis are the percentage difference in area compared to the final predictive ecosystem mapping model.

Stand Type	Extent 1		Extent 2		Landscape
	Species Model	Structure Model	Species Model	Structure Model	Species Model	Structure Model
Rainforest	7.8 (107.6)	4.1 (10.2)	7.5 (267.5)	2.3 (15.3)	783.51 (509.2)	219.45 (70.6)
Ecotone	11.5 (36.0)	12.4 (45.6)	11.5 (−6.3)	15.7 (27.8)	477.30 (−26.9)	2921.27 (347.4)
Eucalypt	5.7 (−11.5)	8.5 (32.3)	6.0 (−20.2)	6.2 (−17.4)	3911.82 (234.5)	1712.16 (46.4)

). This overestimation was particularly excessive for the landscape species model, where the predicted area was 509.2% more for rainforest, and 234.5% more for eucalypt forest in comparison to the PEM. The landscape structure model also overestimated the total amount of rainforest and eucalypt forest relative to the PEM by 70.6% and 46.4%, respectively. The total ecotone area predicted by the landscape species model was somewhat lower than the PEM (−26.9%), but several times larger for the landscape structure model (347.4%).

These differences are illustrated by the stand type distributions predicted by a spatial subset of the landscape species model and landscape structure model illustrated in Figure 6 alongside corresponding PEM and EVC maps. The distribution of ecotone predicted by the landscape structure model was strongly constrained by the ecologically based landscape species model. Conversely, the structural information provided in the landscape structure model constrained the distribution of rainforest and eucalypt stands which would have otherwise been expected using the landscape species model alone. When combining these classifications, in conjunction with the auxiliary criteria listed in the decision trees (i.e., E. regnans and E. delegatensis distributions, elevation and canopy height), the final PEM model predicted a highly complex distribution of forest across the landscape (see Figure S2 for the full study region). The PEM model indicates a realistic mosaic of stand types with many examples of rainforest transitioning to ecotone, which is interspersed within eucalypt forest.

The EVCs (see Figure S3 for the full study region) predict much larger, more uniform clumps of stand types. The distribution of montane eucalypt and sub-alpine woodland predicted by the PEM model corresponds well to the EVC distributions of equivalent classes. Thick corridors of riparian forest indicated by the EVCs were much thinner, or not present, in the PEM model. Non-forest and regions of other stand types predicted by the PEM model indicate the absence of E. delegatensis and E. regnans, or a canopy height less than 6 m.

Differences in the area covered by each PEM model class and EVCs for the complete study area and both 5 km × 5 km extents were calculated (

Table 5. Area of each stand type (km²; percentage of total landscape in parenthesis) as predicted by the predictive ecosystem mapping (PEM) model compared with the ecological vegetation class (EVC) map. Values are for the entire study region (landscape) and two 5 km × 5 km extents (see Figure 1 for extents).

Stand Type		Extent 1		Extent 2		Landscape
PEM	EVC	PEM	EVC	PEM	EVC	PEM	EVC
Rainforest	Cool Temperate Rainforest	3.7	5.4	2.0	4.5	128.6 (2.5)	101.8 (2.0)
Ecotone	Not classified (NC)	8.5	NC	12.3	NC	652.9 (12.6)	NC
Eucalypt	Wet Forest	6.4	8.8	7.6	13.1	1169.4 (22.6)	717.5 (13.9)
	Damp Forest	-	0.0	-	0.0	-	970.9 (18.8)
	Shrubby Wet Forest	-	0.0	-	0.0	-	<0.1 (0.0)
Riparian	Riparian Forest	0.4	0.0	0.3	0.0	90.8 (1.8)	204.3 (3.9)
Montane Eucalypt	Montane Wet Forest	5.3	10.8	2.0	7.4	132.9 (2.6)	458.0 (8.8)
	Montane Damp Forest	-	0.0	-	0.0	-	211.5 (4.1)
Sub-Alpine Woodland	Sub-Alpine Woodland	0.6	0.0	0.0	0.0	90.6 (1.7)	81.3 (1.6)
Other Forest	All Other Vegetation Classes	0.0	0.0	0.1	0.0	2586.6 (50.0)	2135.0 (41.3)
Non-Forest	Non-Forest	0.0	0.0	0.7	0.0	323.1 (6.2)	294.5 (5.7)

). The PEM model predicted 56.2% of the landscape as covered by other forest or non-forest (2909.7 km²), in contrast to the EVCs which predicted 47.0% as covered by all other vegetation classes or non-forest (2429.5 km²). The wet forest, damp forest, and shrubby wet forest EVCs covered a larger area of the landscape (1688.5 km²) than the associated eucalypt class of the PEM model (1169.4 km²). The PEM model had higher area estimates for rainforest (128.6 km² vs. 101.8 km²) and the only estimate for ecotone (652.9 km²). The EVCs had higher area estimates for riparian forest (204.3 km² vs. 90.8 km²) and the combined montane wet forest (458.0 km²) and montane damp forest (211.5 km²) when compared against the associated montane eucalypt class of the PEM model (132.9 km²). The PEM model predicted a slightly larger area of sub-alpine woodland than the EVCs (90.6 km² vs. 81.3 km²). Ecotone stands identified by the PEM model contained multiple EVCs, but were most commonly found in the wet forest, damp forest and montane wet forest EVCs (35.3%, 13.9% and 26.4% of the total ecotone area, respectively;

Table 6. Area covered (km²; percentage cover in parenthesis) by ecological vegetation classes within the ecotone stand type as determined by the predictive ecosystem mapping model. Values are provided for the entire study region (landscape) and two 5 km × 5 km extents (see Figure 1).

Ecological Vegetation Class	Extent 1	Extent 2	Landscape
Cool Temperate Rainforest	2.2 (25.8)	2.2 (17.9)	42.8 (6.5)
Wet Forest	3.6 (43.0)	7.4 (60.6)	230.6 (35.3)
Damp Forest	0.0	0.0	90.5 (13.9)
Shrubby Wet forest	0.0	0.0	0.0
Riparian Forest	0.0	0.0	31.1 (4.8)
Montane Wet Forest	2.6 (31.2)	2.6 (21.4)	172.1 (26.4)
Montane Damp Forest	0.0	0.0	33.8 (5.2)
Sub-Alpine Woodland	0.0	0.0	9.0 (1.4)
All Other Vegetation Classes	0.0	0.0	41.9 (6.4)
Non-Forest	0.0	0.0	1.0 (0.1)
Total Ecotone Area	8.5	12.3	652.9

). The least common EVCs located within the ecotone region were riparian forest, montane damp forest, cool temperate rainforest and sub-alpine woodland (4.8%, 5.2%, 6.5%, and 1.4% of the total ecotone area, respectively).

A visual comparison of the PEM model and the EVC distributions in each 5 km × 5 km extent (see Figure 1 for locations) indicated good general agreement for rainforest and eucalypt forest areas (Figure 7). In general, the finer-scale resolution of the PEM model differentiated stand type patches in contrast to the more continuous distributions of the EVC map. A network of rainforest patches, evident in the bottom left of the first PEM subset, was absent from the EVC map. In contrast, rainforest in the second subset was patchy and covered less area than the cool temperate rainforest EVC.

4. Discussion

The fusion of statistical analyses, lidar-derived structural profiles, and an assemblage of SDMs demonstrates an integrated approach to PEM that was highly sensitive to spatial variations in vegetation structure. The field observations are complemented by the continuous coverage of remote sensing information that indirectly quantifies characteristics of the stand function [3]. The potential importance of capturing this spatial variability is highlighted when comparing the PEM model against the EVCs that are commonly used in the region. The PEM model shows a much greater level of detail in comparison to the EVCs, which are often spatially homogenized with large clusters of uniform stand types. In the little detailed information that has been published about the development of EVCs, it has been suggested that elements of structure were included in the classification process [65]. Structure (spatial arrangement of forest components) variables in modelling are often used interchangeably with attributes describing function (forest processes) and composition (species abundance and diversity), where details related to one can be used as a proxy for another [41]. The extent, level of detail, and number of variables by which forest structure can be defined by lidar differs greatly in magnitude from the in situ measurements and those derived from optical remote sensing data like Landsat data [54].

The high spatial variability of the PEM model suggests that many stand types do not have as continuous a cover as suggested by the EVC distributions, and this is consistent with our extensive field observations of these landscapes. Despite these differences in spatial resolution, there was generally a good level of agreement in the extent and locations of similar stand types predicted by the PEM model and the EVC maps. Our analysis highlights how complex structural and ecological information can complement one another through different modeling pathways. The lidar-derived landscape structure model was effective at constraining the much broader predictions of rainforest and eucalypt stands by the landscape species model. Conversely, the landscape species model was crucial in providing an ecological basis for constraining predictions of ecotone as indicated by the landscape structure model. Leveraging the congruence between predictions driven by characteristic and environmental gradients shows promise for mapping ecotone stands, as they are likely to have the most dynamic distributions.

Lidar-derived characteristic gradients provide a continuous representation of the full vertical vegetation profile, which is useful for identifying the structural gradient between rainforest, ecotone, and eucalypt stands. Lidar can resolve understorey structure with a level of detail that would not otherwise be possible using characteristic gradients derived from optical remote sensing [3,29,54]. This vertical complexity can be measured at very broad spatial scales using lidar data, and repeat acquisitions can provide valuable opportunities to refine model performance and identify important structural changes in vegetation between specific points in time across the landscape. While field observations were limited to undisturbed stands, the output PEM model may be subject to error in the landscape structure model due to the time difference between field observations and lidar capture. Future studies would benefit from coincident field observation and lidar capture. Furthermore, airborne lidar provides a consistent measurement of structure at a scale that is not possible from the ground, including in regions with limited ground access. Nonetheless, one of the main limitations of using purely structural information to predict stand types is that the tolerance of the individual species that comprise these ecosystems to environmental conditions is not considered [49,61]. For example, eucalypt forest and montane eucalypt forest share similar structural characteristics but are strongly delineated by the environmental conditions that are suitable for either E. regnans or E. delegatensis. By incorporating SDMs into the PEM model, an ecological basis for stand type classification can be determined independently of isolated lidar-derived classifications [96].

Lidar-derived structural profiles were intentionally excluded as predictors for SDMs, as the PAVD profiles represent assemblages of vegetation and their inclusion would have resulted in a large increase in model complexity [49,61]. The use of lidar in SDMs has typically included several lidar-derived metrics [12,56,57,59] rather than some derivative that leverages the entire stream of data collected along the vertical profile. Zimble et al. [12] used airborne lidar-derived tree heights to classify western United States forests into single- and multi-layered stands with > 90% accuracy to be used as the structural component to a variation of the United States Forest Service PEM called habitat decision-support systems. Recent studies that have utilized these profiles have expressed similar challenges in distinguishing stand types with similar vertical vegetations structure [49]. In this study, the use of maximum canopy height as a lidar metric was not useful for distinguishing stand types because E. regnans occurs across most of the landscape and grows to consistent maximum heights.

Recent studies [this study, 29,51,55] have identified how remotely sensed gradients constrain the compositional space of stand or community types based on observed changes in forest structure at different successional stages. Hakkenberg et al. [29] explores the use of an alternative approach to modeling community continua by incorporating lidar and hyperspectral remote sensing. Their study utilizes compositional ordination to summarize plot species information into axes of maximum floristic variation (similar to Simonson et al. [51]). Hakkenberg et al. [29] then classified communities using an unsupervised classification approach utilizing goodness-of-clustering evaluators and dissimilarity matrices. This approach is different to our study but achieves similar outcomes to the SIMPER analyses we used by identifying species that characterize forest communities based on underlying matrices of floristic dissimilarity. The use of dissimilarity matrices (e.g., Ferrier et al. [100]) in both species information [this study, 29] and lidar-only metrics [55] is increasingly being used to classify forest stand types. Compositional modelling in this study, and both the Hakkenberg et al. [29] and Moran et al. [55] studies utilized random forest for classification due to its generalizability and maximization of predictive accuracy with its ability to balance limited training data and high data dimensionality. Key differences in this study include the use of stand type classifications using the vertical continuum of lidar data through dimension-reduced principal components as predictor variables and the removal of highly correlated predictor variables (i.e., the landscape structure model). The ordination approach detailed in Simonson et al. [51] and Hakkenberg et al. [29] are ideally suited for landscapes with limited a priori knowledge about community/stand species composition while the approaches detailed in Moran et al. [55] and Fedrigo et al. [61] are suitable when considering classification using lidar-only metric and vertical continuum data, respectively. This study utilizes a priori knowledge of stand type composition, combined with the SIMPER approach to species classification at the landscape scale prior to the use of lidar data, for further stand type identification. Other studies that have performed similar classifications have been in less complex forests in North America [29] and Europe [51], while our study demonstrates the success of similar PEM approaches in structurally complex (multi-layered, dense, closed canopy temperate) forests of south-eastern Australia. All recent PEM studies highlight concerns with mapping stand or community types as discrete units, when their distribution across the landscape are ultimately continuous and can often overlap in areas of transition [29,55,61].

The PEM model developed as part of this research provides the first spatial prediction of dominant stand types in the study region that considers both the ecological niche and the vertical continuum of vegetation structure. Lidar can penetrate forest canopies and can resolve understorey structure with a level of detail that would not otherwise be possible using characteristic gradients derived from optical remote sensing. Lidar-derived characteristic gradients are therefore particularly useful for identifying ecotone stands that are characterized by eucalypt overstorey with rainforest associated understorey. Cameron [62] defined these ecotone stands as ‘cool temperate mixed forest’ but until now they have not been formally delineated. The ecotone may be the most likely to change in distribution due to the regeneration dynamics that occur in these stands with changes in microclimate, disturbance regimes and light conditions [62,63]. The succession of forest stands from one type to another suggests that these stand type maps reflect dynamic processes and will require regular assessment to evaluate their agreement with ground observations [50,58]. Repeat lidar acquisitions over time will provide opportunities to detect changes in the structure and distribution of these stand types in response to natural and anthropogenic disturbances at very high spatial resolutions.

5. Conclusions

This study demonstrates a novel approach to PEM that integrates a range of data types including field-based, optical remotely-sensed, and lidar data to model the continuous distribution of key stand types in the temperate forests of south-eastern Australia. By exploring the use of varying environmental and characteristic predictor variables and identifying those species which best separated stand types, we generated a landscape species model based on the combination of multiple SDMs. Characterizations of stand types based on significantly different species compositions and lidar structural profiles were critical to identifying and mapping the predicted distribution of ecotone between rainforest and eucalypt stands. Using a combination of SDMs and remote sensing data, the ecotone forest was clearly mapped for the first time, offering potential for finer-scale mapping of forests in landscapes beyond the temperate forests of south-eastern Australia. The PEM model allows processes influencing forest development patterns to be captured which should enhance classification accuracy and provide a baseline for evaluating changes in forest development over time. More detailed and high-resolution stand type mapping for this region and others will serve to improve our knowledge of the composition and structure of forested landscapes, and provide a valuable tool for supporting sustainable forest management decisions.