Rapid Characterisation of Vegetation Structure to Predict Refugia and Climate Change Impacts across a Global Biodiversity Hotspot

Identification of refugia is an increasingly important adaptation strategy in conservation planning under rapid anthropogenic climate change. Granite outcrops (GOs) provide extraordinary diversity, including a wide range of taxa, vegetation types and habitats in the Southwest Australian Floristic Region (SWAFR). However, poor characterization of GOs limits the capacity of conservation planning for refugia under climate change. A novel means for the rapid identification of potential refugia is presented, based on the assessment of local-scale environment and vegetation structure in a wider region. This approach was tested on GOs across the SWAFR. Airborne discrete return Light Detection And Ranging (LiDAR) data and Red Green and Blue (RGB) imagery were acquired. Vertical vegetation profiles were used to derive 54 structural classes. Structural vegetation types were described in three areas for supervised classification of a further 13 GOs across the region. Habitat descriptions based on 494 vegetation plots on and around these GOs were used to quantify relationships between environmental variables, ground cover and canopy height. The vegetation surrounding GOs is strongly related to structural vegetation types (Kappa = 0.8) and to its spatial context. Water gaining sites around GOs are characterized by taller and denser vegetation in all areas. The strong relationship between rainfall, soil-depth, and vegetation structure (R2 of 0.8–0.9) allowed comparisons of vegetation structure between current and future climate. Significant shifts in vegetation structural types were predicted and mapped for future climates. Water gaining areas below granite outcrops were identified as important putative refugia. A reduction in rainfall may be offset by the occurrence of deeper soil elsewhere on the outcrop. However, climate change interactions with fire and water table declines may render our conclusions conservative. The LiDAR-based mapping approach presented enables the integration of site-based biotic assessment with structural vegetation types for the rapid delineation and prioritization of key refugia.


Introduction
Considerable changes in the distribution and ecology of species and ecosystems are likely to be ongoing over the coming decades in response to anthropogenic climate change [1][2][3]. Identifying refugia (habitats that facilitate species persistence during largescale and long-term climatic change [4]), is increasingly important in conservation planning as a critical climate change adaptation strategy. Persisting in refugia may provide an important means of in-situ survival for many species [5][6][7].
Identifying the location of refugia requires a spatially explicit understanding of the relationships between biodiversity and the environment (including climate) at appropriate scales and through time. The current reliance on species distribution models (SDMs) is most often applied at coarse spatial scales, but refugia may occur at relatively fine spatial scales [8][9][10]. For example, the globally significant South-West Australian Floristic Region (sensu [11]; henceforth SWAFR) is predicted to experience a decrease in precipitation (e.g. [12,13]), and coarse SDMs predict large impacts on species distributions [14,15]. However, none of these models take into account fine scale environmental heterogeneity, and as a consequence are unable to identify refugia at finer scales -the scales likely to enable local persistence under predicted changes, though see [16,17].
Emerging technologies such as LiDAR (Light Detection and Ranging) and RADAR (Radio Detection and Ranging) systems are powerful tools for the spatially explicit modelling of environment and biodiversity. The increasing availability of these tools enables ready mapping of vegetation structure including overstorey and understorey characteristics [18][19][20]. Delineation of spatial patterns based on structural characteristics can be related to vegetation types on the ground [21,22]. Characterisation of vegetation structure further allows extraction of the key vegetation attributes of height and crown density [23], and can be used to quantify structural heterogeneity at local scales, and to identify environmental constraints and specific habitats. However, to identify refugia under projected climate change, vegetation structural characteristics (as measured by LiDAR or other remote sensing approaches) must be linked to predictive environmental variables. It may then be feasible to link structural vegetation mapping with local process-based measures to identify key vegetation types and places that are representative of refugia at finer spatial scales.
Ongoing changes and those projected under anthropogenic climate change are particularly important for mediterraneanclimate ecosystems [24,25]. These mediterranean-climate ecosystem regions occur on six continents [26], harbour a substantial proportion of the Earth's vascular plant flora [27], and are all recognized as global biodiversity hotspots [28]. The SWAFR is the least topographically complex of the five mediterranean-climate ecosystem regions with little opportunity for contraction to mountain refugia as the climate warms [26]. It is also predicted to be the most adversely affected by projected climate change, with consensus among global climate models that rainfall will continue to decline [13,25,29].
The SWAFR is characterised by the ancient granite-based landscapes of the Yilgarn Craton and Albany Fraser Orogen [30]. Granite inselbergs or outcrops (GOs) are topographically complex in comparison with the subdued surrounding landscape. Hence microclimatic variation, due to topographic and indirect effects of soil moisture variability [17], within GOs could buffer against regional climate change, and could continue to provide habitats for species occurring on or around them. Granite outcrops are generally rich in biodiversity, and are therefore of great conservation importance [31,32]. In south-western Australia, at least 1200 vascular plant taxa are found on GOs [31], and GOs harbour a considerable proportion of the region's invertebrate, reptile, bird and mammal faunas [33,34]. The importance of GOs is even greater in disturbed agricultural landscapes, where they constitute important habitat remnants for the biota [35,36], and gene pools for surrounding landscapes under restoration.
The elevated nature and geological constitution of GOs means that they channel water, nutrients and plant residues to the fringes of the rock [37,38], where growing conditions are more favourable for plants [39,40]. This capacity may be important in southwestern Australia, where moisture deficits, nutrient impoverishment, and acidity are typical features of local soils [41][42][43]. Weathering on exposed GOs provides nutrients and sediments to associated colluvial and alluvial fans surrounding the outcrops [40], reducing local constraints on plant growth. In addition, the slope and shallow soils of GOs reduce waterlogging, and basement rock beneath the fringe prevents water seeping away into deeper aquifers. Therefore, it is predicted that the fringes of GOs should have denser vegetation, greater biomass and higher productivity than the surrounding landscape [39].
Spatially explicit modelling of vegetation structure in conjunction with environmental variables will allow investigation of the interactions between vegetation characteristics and climate. A consistent characterization of local vegetation structure within a regional context should enable the quantification of habitats [44,45] and the identification of environmental constraints influencing growth. Such an approach is based on the strong relationships between environment, vegetation type, and density of vegetation [46][47][48]. The five mediterranean-climate regions are commonly cited examples of convergent evolution in vegetation structure and function [27,49]. Hence, some consistency in vegetation structure could be expected in similar environments across the SWAFR region despite remarkable floristic diversity (see for example [50]). Canopy height has been widely used to assess habitat condition and conservation status [21]. While correlations between canopy height/cover and environmental attributes are well known, there has been no previous attempt to provide finescale spatial realization of those relationships.
In this paper we establish relationships between vegetation structure and environmental variables to identify refugia using a case study of GOs across the SWAFR. We sought to generate detailed maps of local structural vegetation type as a means to relate growth to environmental variables indicative of local resource availability and growth constraints in topographically complex areas. This would need to be at a scale relevant to conservation management and the identification of refugia as safe havens for biodiversity. We have four specific aims and associated hypotheses.
1) Delineate and map local-scale vegetation structure across the region. We expect a wide variety of consistent vegetation structural forms applicable across the region reflecting variation in resource availability, despite local variation in species composition. 2) Compare the local spatial distribution of the types of vegetation structure near GOs across the rainfall gradient. We expect the tallest and densest vegetation to be confined to run-on areas at the base of GOs, because vegetation in these run-on areas have access to additional nutrients and water from the GO when compared to the surrounding areas. 3) Quantify the relationship between environmental variables and habitat attributes derived from vegetation structure. We expect a reduction in canopy height and cover of vegetation with reduced rainfall and soil depth. 4) Portray the vegetation structure predicted for sites under climate change projections for the region. We expect sites closest to GOs to retain proportionally denser and taller vegetation than sites further away from them. We also expect areas of tallest vegetation to be most affected by climate change.

Materials and Methods
Airborne LiDAR data and Red, Green and Blue (RGB) imagery were acquired by AAM Pty Ltd (Perth, Australia) from flights covering the areas around 28 targeted GOs across the SWAFR rainfall-gradient from mesic to low rainfall environments (Fig. 1). The selection of GOs was based on the inclusion of large and iconic GOs, while covering the full rainfall gradient with a limited number of flights, opportunistically including additional GOs in the flight path. The area bounded by the polygon connecting the surveyed sites was 296,361 km 2 . Within this region, LiDAR and RGB imagery were obtained over a total area of 95,485 ha.

LiDAR data processing
A description of the LIDAR processing is provided more fully elsewhere [22]. Overlapping areas of adjacent runs produced strips with a denser point cloud, these variations in point density being undesired [51]. Hence all layers were initially processed in a 464 m raster and, in each cell with LiDAR returns from overlapping runs only points were included from the run with the smallest mean scan angle. This produced grid cells with up to 32 returns per cell, although 11-15 returns were typically found in vegetated areas with trees.
A 262 m digital elevation model (DEM) of the terrain was derived from triangulation of all ground points, which was consequently used to determine the height above ground for all returns classified as non-ground (see [51] for a description of the procedure). Similarly, a 262 m canopy height layer was determined by subtracting the elevation of the ground from the elevation of each return. The DEM and selected layers of the GOs are downloadable (http://refugia.curtin.edu.au/). Other datalayers can be made available upon email request.

Volumetric pixels
Subtle variations in point density occurred throughout the flight. Therefore, presence and absence of vegetation was recorded in 3D volumetric pixels (voxels) as these were less sensitive to these point cloud density variations than calculated percentages of returns per defined vertical layer. For each vertical layer of 1 m height, percentages of filled voxels (PVF) were computed within a 363 m window, producing a vertical PVF profile for each grid cell. These PVF profiles were smoothed in the vertical direction to better define the top and bottom of canopy layers by applying a simple Gaussian filter (with weights of 0.27, 0.46 and 0.27). Following the approach of Reitberger [52], a threshold of 20% (i.e. at least 2 out of 9 voxels) was initially used to trigger the start and end of a vegetation layer within a PVF profile. For the first two of these vegetation layers, canopy height, layer thickness, mean coverage and mean intensity were recorded. Mean intensity was based on only first and single returns, as return number has a profound influence on the recorded intensity [53]. The smoothed PVF profiles were further sampled at 18 height intervals from 1 to 80 m with intervals ranging from 1 to 20 m, increasing with height. These were stored for further processing.

Identification of PVF profile classes
Local spatial heterogeneity in PVF profiles may identify patchiness that is significant for understorey and midstorey canopy layers. A hybrid classification procedure was used, aimed at using this spatial heterogeneity to identify vegetation-structural types. The sampled PVF profiles and canopy heights, cover, intensity and layer height of the top canopy layer (i.e. crown thickness) were Detailed plot-based floristic surveys were carried out at 16 of these sites (filled squares with abbreviated names, see also Table 1). Isohyets are also shown. Closed triangles are regional centres. doi:10.1371/journal.pone.0082778.g001 transformed into principle components using standard image processing software (ERDAS ER MapperH, Intergraph, Alabama, USA). Significant principal components, capturing 99% of the variation were used to identify unique PVF profiles for Mt Frankland National Park, and for Boyagin and Chiddarcooping Nature Reserves. These sites represented mesic, intermediate and dry environments, respectively. The iterative, self-organising (ISO, also known as the migrating means) clustering technique [54], was used to identify significant classes covering at least 1% of the area using ERDAS ER MapperH software. For each of these classes, based on all pixels within these classes, minimum, maximum and median values per vegetation layer were determined. Also, typical vertical profiles were constructed and used to derive a qualitative description for each class. Profiles of the three areas were combined and similar profiles were removed, resulting in 54 unique PVF profile classes.
LiDAR data were further classified according to the type and elevation of the returns [44], where four vegetation layers were identified: low (,1 m above the ground), medium (1-3 m), high (3-10 m) and top canopy (.10 m) vegetation layers. The PVF profile classes were used in a pixel-based supervised classification using a minimum distance classifier for all areas. LiDAR intensity was excluded from the classification as the LiDAR intensity corrections for scan angle and flying height did not fully correct intensity differences between scanned areas. The minimum distance classifier is non-parametric, simple and fast and does not require dispersion statistics that need to be derived from training data [54].

Conversion of PVF profile classes into structural vegetation classes
At selected locations in Mt Frankland, Boyagin and Chiddarcooping, vegetation types characterised by a similar floristic composition and vegetation structure were described and the presence of dominant overstorey and understorey species was recorded. Field based geocoded photographs were taken at the same locations. These locations were selected on transects covering a wide range of PVF profile classes to ensure that visited locations were representative of variation in structural vegetation types. Variations in the PVF profiles within an area with similar vegetation typically represented a range of understorey conditions and disturbance history. Typical combinations of 54 PVF profile classes were identified (based on similarity in vertical profiles and co-occurrence within a single vegetation type after examination of RGB imagery), and merged into 27 meaningful broader structural vegetation classes. Thus each of these structural vegetation classes were related to at least one of the vegetation types observed.
A segmentation step was used to derive so-called objects or polygons based on canopy height at local (i.e. a large tree) and slightly larger ''vegetation'' scale by adjusting the scale parameter (eCognition DeveloperH 8, Trimble Geospatial Imaging, München, Germany). These polygons were classified according to the structural vegetation class with the largest relative coverage. Bare areas without above-ground LiDAR returns within native vegetation areas were further subdivided into smaller polygons based on RGB brightness (i.e. the sum of the digital numbers in all three bands). Thresholds based on local brightness differences could be used to differentiate bare rock from moss-mats for those outcrops in higher rainfall areas where moss-mats are darker in colour. Thus, the brightness, derived from the RGB values, within these smaller scale polygons was compared to the brightness of vegetation scale polygons. They were assigned to the moss-mat class if between particular thresholds, with thresholds iteratively adjusted for each outcrop. For dryer areas, moss mats were less prominent and also less visible in summer when the aircraft was flown. There was also considerable variation in illumination across any individual outcrop. Therefore brightness of moss-mats either did not contrast strongly, or differences were small compared to variation in brightness across particular outcrops. Consequently, an approach based on local brightness differences could not be used, and moss-mats were not further classified for these drier areas.

Supervised classification
Delineations based on structural characteristics are most valuable when they also identify transitions between vegetation types. Structural vegetation classes can effectively be assigned to vegetation types when combined with local expert knowledge [22]. However, within a regional context, structural vegetation classes are not specific for a vegetation type. Therefore field records for validation were collected in transects at Boyagin (135 records), Chiddarcooping (25 records), Porongurups (106 records) and Mt Frankland (168 records) to determine whether identified structural vegetation classes were related to floristic vegetation types, or to transitions between them. The observations included locations with various periods of recovery after fire. Kappa coefficients, indicating classification accuracy [54], were determined for Boyagin, Porongurups and Mt Frankland.
Means were determined for canopy height and for ground cover for each structural vegetation class, based on all classified polygons within Boyagin, Mt Frankland and Chiddarcooping.

Floristic surveys and plots
A total of 16 GOs were selected for detailed study (Fig. 1). Between 25 and 36 plots (see below) were established for most GOs with the exception of Boyagin Rock where 71 plots were established. At each of these 16 GOs, plots were divided between three major habitat types: 1) sites either with shallow soils with water off-flow or seasonally inundated soil-filled rock-pools, also known as gnammas (OF, herbfield vegetation); 2) sites with moderately deep soils or with vegetation with access to cracks in the rock with on-and off-flow of water (INT, usually with shrubs or low trees); and 3) on-flow areas at the base of the outcrops (ON, typically with forest vegetation -particularly in higher rainfall areas).
Sizes of plots were based on accumulation curves for species and the size of taxa under investigation (i.e., 161 m for OF plots, 565 m for INT plots and 20620 m for ON plots). Plot locations were recorded with a hand-held GPS. Floristic composition (all vascular plants) and soil depth was recorded for each plot, using a soil depth probe hammered or pushed into the ground by hand until maximum depth (or at least 50 cm for deeper soils), at 5 random locations in the plot. These were used to calculate the probability of a soil deeper than 0.5 m (pDS). A total of 494 plot locations were available across the 16 sites with associated LiDAR data.

Statistics for plot locations
Statistics concerning canopy height and ground coverage were calculated based on a 4 m buffer around the geocoded plot locations, intersected with classified polygons. The buffer was used to account for the error made when recording the position using a hand-held GPS. Means of ground coverage and maximum canopy height were determined for each plot (including buffer). In addition, weighted values were determined from structural vegetation class means, based on all intersected polygons that cover at least a quarter of a plot.
The influence of environmental variables on canopy height and ground cover was tested with a simple linear regression model. Multivariate linear models (Y = aX+e) were fitted with interactions included. P-values were evaluated for each variable in the explanatory X-block and variables were excluded when not significant. Potential explanatory variables included annual precipitation (interpolated values from the WorldClim dataset based on means of the years 1950-2000 [55]), elevation range (difference in elevation between highest and lowest elevation plots for each GO) and pDS. Elevation range may be a proxy for potential runoff and may also influence rainfall, particularly in coastal zones. Robustness of relationships was evaluated with a leave-one-out (LOO) validation, using the Q 2 statistic [56].

Refugial capacity
For the SWAFR, the A1F1 scenario, including high CO 2 emissions predicts a greater than 40% chance of exceeding a 20% reduction in rainfall when compared to the 1961-1990 reference period for climate in 2070 [57]. We used this 20% reduction to illustrate potential changes in vegetation structure and the utility of the methodology derived.
The surveyed plots were located along the rainfall gradient from mesic in the High Rainfall Province (.600 mm rain p.a. -both Yilgarn Craton and Albany Fraser Orogen sites) to the inland side of the Transitional Rainfall (Yilgarn) and eastern edge of the Southeast Coastal (Albany Fraser) Provinces (both 300-600 mm rain p.a. [11]), with differences in vegetation structure reflecting these gradients. To enable use of these relationships for deriving projections around GOs, spatially explicit values for each polygon were needed for all terms included. For pDS, these values were estimated using regressions based on current vegetation structure and rainfall, using linear regressions based on mean values derived from the assigned structural vegetation class covering plot areas. Estimated pDS values, maximum supported canopy height and ground cover were determined for each polygon using the present mean vegetation height and ground coverage.
Under the assumption that current relationships between environmental variables and vegetation structure that vary along a spatial gradient can be used to predict an in-situ change over time, these relationships can be used to assess the impact of new climate regimes across the region. This means virtually ''relocating'' current outcrops to new climate regimes equivalent to current climates in lower rainfall areas. For each polygon, the current structural vegetation classification and means of environmental variables were exported (to a simple spreadsheet). The equations describing relationships between environmental variables and canopy height for plot locations, as described above, were used to predict future canopy height under an A1F1 scenario with a 20% rainfall reduction. This was repeated for ground cover. The distance between predicted canopy height and ground cover and mean values of each structural vegetation class were calculated and used to reclassify each polygon. Individual classes are heterogeneous and mean canopy height and ground coverage of a single polygon can be much higher than the mean of the class where the polygon was assigned to. This because many other features were also included in the supervised classification. Thus a small change in canopy height or ground cover may lead to a reclassification into a class with larger mean canopy height or ground coverage. Reclassification was therefore restricted and only vegetation structure classes with a lower or equal mean canopy height and ground cover than the class currently assigned to the polygon could be selected. From these, the vegetation structure class with the lowest minimal distance was selected and linked to the polygon for display.

Results
The canopy height of vegetation in plots on or near GOs (see Fig. 1. for the locations of the GOs) responded more strongly to precipitation than to ground coverage, with annual precipitation ranging between 314-1208 mm yr 21 for the GOs included in this study (Table 1). For most GOs, mean and maximum canopy height was taller and ground coverage higher in on-flow plots than in intermediate or of-flow plots.

Structural vegetation types
The vertical vegetation profiles generated showed a wide range in canopy height, distribution of canopy elements and ground cover across the region (Fig. 2), reflecting the range of structural vegetation types present. Although each individual PVF was distinct, groups of PVFs with similar vertical distribution profiles, but differences in ground cover, can be recognized. Some structural classes included a wider range of PVFs than others (e.g. the Open woodland class represents a wide range of vegetation types). However, these PVFs and structural classes can be linked to structural vegetation types when combined with topography and local expert knowledge [22]. In some areas, classes with a distinct vegetation structure can be directly linked to floristic vegetation types. For example, tall open-forest dominated by Eucalyptus diversicolor (Myrtaceae, karri) had a very distinct profile (Fig. 2), which in some areas can be further linked to the age of stands. Thus, in the area of Mt Frankland, substantial areas of karri regrowth show a distinct profile in comparison with oldgrowth stands dominated by the same species. However, in other situations, distinct profiles may be associated with multiple vegetation types, reflecting differences in landscape position and an array of species composition. For example, combinations of these profiles can be directly linked to vegetation types when combined with local descriptions for Mt Frankland and the Porongurups (Kappa coefficients of 0.86 and 0.78 respectively), and when combined with local expert knowledge for Boyagin, see [22].

Comparative spatial distribution of vegetation types
For all GOs, vegetation is taller and denser in on-flow areas at the base of the outcrops (Fig. 3). Tall karri trees (to 70 m in height) dominated on-flow plots with highest annual rainfall (i.e. Mt Chudalup and Mt Frankland), and at sites with relatively large topographic relief (i.e. the Porongurups which influences local climate). However, karri is restricted to the highest rainfall and least seasonal end of the High Rainfall Province. Hence vegetation height was much lower elsewhere, even in on-flow sites.
On GOs, low and open vegetation associated with shallow and rocky soils is dominant, usually covered with shrubs, herb fields or moss mats. Vegetation types with similar structure can be found across the rainfall gradient. Canopy height in off-flow areas for Mt Frankland and Mt Chudalup were taller than may have been expected. However, taller and denser vegetation occurs further from outcrops in higher rainfall parts of the Transitional Rainfall and Southeast Coastal Provinces but is confined to narrow fringes in on-flow areas near outcrops in lower rainfall areas.
In the High Rainfall Province, low and open vegetation in shallow soils occurs only on granite outcrops. However, this vegetation structure (within the study sites) is replaced by denser and taller vegetation where soil depth increases, in soil pockets on the GOs, and in areas surrounding them. In high rainfall and swampy areas, dense shrublands occur in lower landscape positions (e.g. surrounding Mount Chudalup and Mount Frankland). At the lower rainfall end of the Transitional Rainfall and Southeast Coastal Provinces, low, open vegetation with patchy or scattered shrubs is also common on deeper sandy soils with limited water holding capacity further from the outcrop.

Environmental variables and habitat attributes
The average height and ground cover in the polygons overlapping the plot locations were related to environmental factors (rainfall, elevation range and soil depth). The R 2 values of the simple regression models ranged from 0.48 to 0.91 for maximum canopy height, and from 0.39 to 0.84 for ground cover ( Table 2). For canopy height in on-flow plots the elevation range is a significant term, and in combination with rainfall, reflects the importance of moisture regimes in these sites. Separate models run only with the inclusion of off-flow and intermediate plots had low explanatory power (R 2 ,0.1, not shown).
When plots were grouped according to structural vegetation type, strong relations with environmental factors emerged (R 2 value of 0.92 and 0.84). The models were also robust, with a large positive value for Q 2 in the leave-one-out validation of the relationships. A combination of soil depth and annual rainfall resulted in relationships with R 2 values of 0.75 for ground cover and 0.93 for canopy height (Fig. 4), with each point in this relationship representing a different structural class in plots on or near outcrops.
The absolute changes in canopy height and ground cover are larger per mm of rainfall reduction for patches with a larger pDS. With an equally strong change in ground cover for shallow soils, the influence on structural vegetation types may be much larger for shallow rocky areas. For example, in an area with 500 mm rainfall, a 100 mm reduction in rainfall on soils with a pDS of 0.5 results in a change in canopy height of about 1 m. This reduction in available water can be compensated by moving to a site with a pDS of 0.7 (see Fig. 4).
We found that the pDS can be modelled using attributes of the structural vegetation maps combined with rainfall (R) values: pDS = 0.55+(10.896CH+8.936GC20.736R)/1000. Canopy height (CH), ground cover (GC) derived from LiDAR and rainfall explained 85% of the variation in pDS (N = 22, R 2 = 0.85, LOO Q 2 values of 0.76). This relationship allows estimation of the pDS of the wider surroundings of GOs, as the vegetation of those areas will also reflect their environment. Some structural vegetation types represent a transition (e.g. karri regrowth represents an intermediate phase in the life of the forest stand) in vegetation structure.

Refugial capacity
The means of canopy height and ground coverage within structural classes were directly related to environmental factors along the climate gradients in this study. Relationships were strong enough to meaningfully substitute space for time in assessing the impact of new climate regimes across the region.
There is a significant reduction in the area covered by vegetation with taller and denser canopies with a 20% reduction in rainfall (Fig. 5). Higher rainfall areas show a proportionally much greater change in both height and cover with this rainfall reduction (Fig. 6). For example, vegetation from Mt Frankland changes from Tall open-forest dominated by karri to scattered areas of karri forest. This structural vegetation type is projected to contract to water gaining areas, where currently the tallest trees are supported. For example, much of the current vegetation in the Open woodland 20-25 m class, typically including E. guilfoylei (yellow tingle) in association with Corymbia calophyla (Myrtaceae, marri) and E. marginata (jarrah), is projected to contract to narrow fringes surrounding the area of tall karri forest. Low, dense vegetation and shrubland found on Mt Frankland is projected to be replaced by a low and more open vegetation structure. However, dense shrub vegetation was not found on areas surrounding the rock in future climates, as taller vegetation was still supported on deeper soils. A strong reduction in proportional area for this habitat type is expected. However, in the wider landscape, the proportional reduction in dense shrub was much smaller (Fig. 6) than that of taller vegetation.
For Boyagin, the areas classified in the Trees 10-15 m class, with Allocasuarina huegeliana (Casuarinaceae, rock sheoak) and Eucalyptus accedens (Myrtaceae, powder-bark wandoo), are projected to be replaced by shrubs with scattered trees, a vegetation class now including kwongan vegetation [58]. The total proportion of this open shrubland vegetation is expected to decrease greatly (Fig. 6), and the shallow gravelly soils currently supporting kwongan vegetation will likely become low vegetation. There is a major projected reduction in patches with dense shrubland or scattered trees on Boyagin Rock (Fig. 6 and Fig. 7), with these structural vegetation types contracting to the base of the outcrop.
Areas with deeper soils along waterways below GOs are also important for relatively dense and tall vegetation. For example, in Chiddarcooping, Eucalyptus salubris (gimlet) or E. salmonophloia (salmon gum) occur in deeper soils, classified as Trees 10-15 m (Fig. 5). Areas further from streams will only support vegetation classified as Open tall shrubs or Scattered shrubs, reflecting a much more scattered and open vegetation, with likely changes in composition. The proportional changes in vegetation types in the wider landscape are expected to be more pronounced than those Table 2. Linear regression and leave-one-out (LOO) validation statistics of multiple linear regressions with parameter estimates of environmental factors determining maximum canopy height and ground coverage. directly surrounding the outcrop (Fig. 6). The dense vegetation in narrow fringes around the base of GOs, typically including Acacia lasiocalyx (Fabaceae, rock acacia), and rock sheoak, is projected to further contract and become more open. On the GOs, Patchy or Scattered shrubs on shallow soils are likely to be replaced by low vegetation or herbfields with annual plants (Fig. 7).

Discussion
We have provided a rapid characterisation of local-scale vegetation structure in relation to environment across the SWAFR. This has enabled the identification of fine-scale patterns of vegetation structure, and projections under anthropogenic climate change. A warmer, drier climate means that ecophysiological thresholds of some species may be reached locally in areas due to spatial variation in topography and radiation [9]. In this light, prudent management for conservation is likely to focus on areas such as refugia [59], where biodiversity may be able to persist for longest [4], although careful experimentation on ecophysiological thresholds is needed to test this hypothesis.
LiDAR-based mapping of vegetation structure can highlight specific areas for potential conservation and protection within a broader regional context. However, the assumption that spatial relationships between vegetation structure and environment can be used to predict in-situ temporal changes needs further critical testing. Our combination of spatially explicit mapping of structural vegetation types with environmental variables and site-based biotic assessment enabled the rapid delineation and prioritization of key potential refugia.

Structural vegetation types
We found that LiDAR can be used to quantify differences in vegetation structure at a local scale [44,53,60,61]. Canopy height in off-flow areas for GOs in very high rainfall areas was taller than expected due to the proximity of plots to tall vegetation, affecting the LiDAR derived height estimates. Tall vegetation strongly influences height in neighbouring raster cells due to the triangulation of top canopy returns, and some of these cells may have been included in the 4 m buffer that was used, affecting the height estimate of these herbfield plots.
The voxel-based characterization used here [22], allowed explicit characterization of vertical profiles for every pixel, despite low density point clouds. The identified classes reflected differences in canopy height, density and the vertical distribution of vegetation. The wide range of LiDAR instruments and techniques available for canopy characterization [20] provide the means to translate plot-based assessments to larger areas in a wide range of studies.  Typically, single vegetation types include several vertical profiles, mostly occurring in regular spatial sequences (e.g. a dense tree adjacent an open shrubland), and a single vertical profile may occur in several vegetation types [22]. Vegetation structure may be disturbed in fire-prone landscapes [62,63]. However, the structural classes recognised here were predominantly based on the overstorey, with the understorey being of secondary importance. Hence, the accuracy of the overall classification was not strongly influenced by the presence of various recovery periods after low intensity fire as observed in the visited locations, although high intensity or frequent fire may have very different consequences.
The importance of water gaining on-flow areas at the base of granite outcrops has long been recognised [38,64,65]. We found that these areas support denser and higher vegetation when compared to the immediate surroundings. However, as expected, they also share structural similarities to the vegetation in higher rainfall areas.

Environmental variables and habitat attributes
This study provides a methodology to link vegetation structure with environment at fine spatial scales over a broad geographic area. As expected, we found that rainfall and soil depth had a significant influence on vegetation height and ground cover in plots located on and around outcrops. Other environmental variables that were not considered in our study, e.g. the amount of water influx, and availability of cracks in the rock [66], may also be important.
The novelty of our approach is that it makes the relationships between environment and vegetation structure spatially explicit at a fine scale, and reveals potential associated patterns in relation to predicted climate change. For on-flow plots, significant interactions between elevation range and rainfall were found for canopy height, indicating the importance of runoff from the outcrop. However, elevation range was not significant when explaining differences in ground cover. The relationships based on values from individual plots were different between GOs in the High Rainfall Province and other areas. This difference was accounted for when averaged over structural class means, demonstrating that our approach can be applicable across the wider SWAFR.
The classification of vegetation structure at local scales enables a quantification of environmental drivers for important habitat characteristics such as vegetation height and cover. There was a direct relationship between rainfall, soil-depth and vegetation structure, suggesting that water availability is the major driver of vegetation structure in these environments. The strength of this relationship may be illustrated by the strong crown decline in response to the reduction in rainfall in recent years [67,68], with a higher incidence of crown dieback on soils with stony profiles and low water holding capacity during the 2010/2011 summer which was the hottest and driest on record [69].
The vertical distribution of canopy elements is strongly correlated with the diversity of vascular plant species [70] and with faunal diversity [60]. A change in vertical structure due to a reduction in rainfall therefore has direct implications for biodiversity. Monitoring of vegetation structure with LiDAR provides a means of assessing overall habitat condition, and ecophysiological response to changing climate [70]. These strong relationships between vegetation structure and climate indicate that the structural vegetation map may also be used to identify environmental constraints within the regional context for areas directly around GOs. It should, however, be noted that further away from outcrops, other constraints, such as waterlogging or salinity, may be of greater importance than the proximity of the GO in determining vegetation height. Therefore, care is required when extrapolating the relationship between rainfall, soil depth and vegetation structure where the response to climate change may be very different.

Comparative spatial distribution of vegetation types
The comparison of vegetation structure on plots across a rainfall gradient has provided a means of understanding spatial patterns within the landscape context that is an essential element of the identification of refugia [4]. Similar vegetation structure was found across the rainfall gradient, but their landscape position varied predictably in relation to water availability. Vegetation types occurring within the broader landscape in the mesic end of the gradient were confined to on-flow areas at the lower rainfall end. These on-flow areas have access to more water [38], and consequently may also have unique microclimates resulting from the topography and vegetation structure [10,17,71]. In more mesic areas, low and open vegetation is confined to outcrop areas, whereas these structural vegetation types were dominant on GOs in the lower rainfall areas of the region. Under projected climate change, reductions of up to 20% in rainfall in comparison with the base period 1960-1990 can be expected by 2070 [57], emphasising the importance of so called drought refugia [6]. This indicates that granite outcrops exhibiting a wider range of habitats and water gaining on-flow areas may facilitate species persistence -important characteristics of refugia under climate change.
A close connection between rainfall and catchment groundwater-storage has been documented in the area encompassing the High Rainfall Province within the Yilgarn Craton [72]. The significant decline (a reduction of 14% in May-July was already observed for the years 1975-2004 when compared to  in autumn and early winter rainfall in the area since the 1970s [12,13,73,74] has been accompanied by a shift from perennial to ephemeral streams, a regional decline in water-tables [72], a decline in the runoff coefficient, and the development of a new hydraulic regime [75]. Recent canopy death has been observed in the overstorey in this area [68,76], generally on shallow soils around granite outcrops. This may be contradictory to the suggestion that on-flow sites at the base of granite outcrops may serve as refugia for taller vegetation in lower rainfall areas. However, on-flow areas with deep soils below granite outcrops will have greatest access to moisture via rainfall and runoff as the regolith continues to dry. The open forest of this area currently accesses moisture from deep, highly weathered lateritic soil profiles that store a large proportion of winter rains [77,78]. As the water table further declines, forests on the shallowest soils of the region that have low water holding capacity, will be first affected [67,68]. The hydrology of the area may be completely transformed when the regolith dries and groundwater becomes disconnected from the stream zone [73] and water gaining on-flow sites become increasingly important for biodiversity.

The influence of disturbance history
In mediterranean-climate ecosystems fire regimes have a considerable influence on vegetation structure and also composition [63,79]. Thus, where fire is sufficiently frequent it can determine species composition [79] and may become a major driver of vegetation structural change under climate change [80,81]. The sites included in this study have an array of disturbance histories that may affect understorey vegetation structure. This should not be expected to greatly influence PVF structure in tall forest, where overstorey height and cover provides a distinct signature. However, in areas of high intensity fire (e.g. Mt Cooke, January 2003, 8 years prior to imagery being flown) or with low vegetation height and frequent fires e.g. Chiddarcooping [82,83], fire may have a significant influence on PVF that may be age-since-fire specific, or related to particular fire regime history. Recognition that disturbance in these fire-prone landscapes may temporarily (or permanently under regime shifts) change vegetation structure [62,63,79] must be accommodated in vegetation mapping [84].

Conclusion
Our study demonstrates the utility of enabling technologies such as LiDAR for identifying and mapping putative climate change refugia. Using granite outcrops in the SWAFR as a case study we found that the vegetation around outcrops included a wide range of structural classes, reflecting differences in local topography, soil depth and water influx associated with the large diversity of habitats found there [31,85]. Under a rainfall reduction scenario predicted for the region, we were able to identify areas where vegetation structure may be likely to persist for longest, therefore providing safe havens for the biota under climate change. However, we acknowledge that interactions such as fire and declining water tables also influence response to climate change [76,[79][80][81]86]. In addition, our projections are likely to be conservative since current vegetation structure may not yet reflect the major changes in rainfall reduction experienced in the years 2000-2010.