UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Assessing indicators of forest sustainability using lidar remote sensing Bater, Christopher William 2008

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
24-ubc_2008_spring_bater_christopher.pdf [ 9.16MB ]
Metadata
JSON: 24-1.0066311.json
JSON-LD: 24-1.0066311-ld.json
RDF/XML (Pretty): 24-1.0066311-rdf.xml
RDF/JSON: 24-1.0066311-rdf.json
Turtle: 24-1.0066311-turtle.txt
N-Triples: 24-1.0066311-rdf-ntriples.txt
Original Record: 24-1.0066311-source.json
Full Text
24-1.0066311-fulltext.txt
Citation
24-1.0066311.ris

Full Text

ASSESSING INDICATORS OF FOREST SUSTAINABILITY USING LIDAR REMOTE SENSING by  Christopher William Bater  B.A. (Hons.) University of Winnipeg, 2002  A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF  MASTER OF SCIENCE  in  The Faculty of Graduate Studies  (Forestry)  THE UNIVERSITY OF BRITISH COLUMBIA  (Vancouver)  March 2008  © Christopher William Bater, 2008  ABSTRACT The Province of British Columbia is developing a suite of attributes to assess and monitor forest sustainability. Each attribute is in turn evaluated using a variety of indicators. Recently, digital remote sensing technologies have emerged as both alternative and supplement to traditional monitoring techniques, with light detection and ranging (lidar) in particular showing great promise for estimating a variety of indicators. The goal of this thesis was to review and assess the ability of lidar to estimate selected indicators of forest sustainability. Specifically, digital elevation model (DEM) interpolation (from which indicators are extracted both directly and indirectly) and wildlife tree class distributions were examined.  Digital elevation models are a key derivative of lidar data, and their generation is a critical step in the data processing stream. A validation exercise was undertaken to determine which combination of interpolation routine and spatial resolution was the most accurate. Ground returns were randomly subsetted into prediction and validation datasets. Linear, quintic, natural neighbour, spline with tension, regularized spline, inverse distance weighting, and ANUDEM interpolation routines were used to generate surfaces at spatial resolutions of 0.5, 1.0, and 1.5 m. The 0.5 m natural neighbour surface was found to be the most accurate (RMSE=0.17 m). Classification and regression tree analysis indicated that slope and ground return density were the best predictors of interpolation error.  The amount and variability of living and dead wood in a forest stand is an important indicator of forest biodiversity. In the second study, the capacity of lidar to estimate the  ii  distribution of living and dead trees within forests is investigated. Twenty-two field plots were established in which each stem (DBH>10cm) was assigned to a wildlife tree (WT) class. For each plot, a suite of lidar-derived predictor variables were extracted. Ordinal logistic regression was then employed to predict the cumulative proportions of stems within the WT classes. Results indicated that the coefficient of variation of the lidar height data was the best predictor variable (r = 0.85, p <0.000, RMSE = 4.9%). The derived relationships allowed for the prediction of the proportion of stems within WT classes across the landscape.  iii  TABLE OF CONTENTS ABSTRACT........................................................................................................................ ii TABLE OF CONTENTS................................................................................................... iv LIST OF TABLES............................................................................................................. vi LIST OF FIGURES .......................................................................................................... vii ACKNOWLEDGEMENTS............................................................................................... ix CO-AUTHORSHIP STATEMENT ................................................................................... x 1 INTRODUCTION ........................................................................................................... 1 Lidar Background ........................................................................................................... 3 Use of Lidar to Estimate Indicators of Forest Sustainability.......................................... 5 Objectives ..................................................................................................................... 10 References..................................................................................................................... 12 2 EVALUATING ERROR ASSOCIATED WITH LIDAR-DERIVED DEM INTERPOLATION........................................................................................................... 16 Introduction................................................................................................................... 16 Methods......................................................................................................................... 23 Study Area ................................................................................................................ 23 Data Acquisition ....................................................................................................... 24 Data Processing......................................................................................................... 25 Error Analysis ........................................................................................................... 30 Error Prediction......................................................................................................... 31 Results........................................................................................................................... 32 DEM Validation........................................................................................................ 32 Error Analysis ........................................................................................................... 35 Error Prediction......................................................................................................... 39 Discussion ..................................................................................................................... 41 DEM Validation........................................................................................................ 41 Error Analysis ........................................................................................................... 42 Error Prediction......................................................................................................... 44 References..................................................................................................................... 45 3 TOWARDS THE ESTIMATION OF TREE STRUCTURAL CLASS IN NORTHWEST COASTAL FORESTS USING LIDAR REMOTE SENSING............... 50 Introduction................................................................................................................... 50 Methods......................................................................................................................... 54 Study Area ................................................................................................................ 54 Data Sets ................................................................................................................... 55 Lidar Data Processing ............................................................................................... 59 Logistic Regression for Estimating Wildlife Tree Class Proportions....................... 60 Predicting Wildlife Tree Class Across the Landscape.............................................. 63 Results........................................................................................................................... 63 Vertical Foliage Profiles ........................................................................................... 64  iv  Logistic Regression for Estimating Wildlife Tree Class Proportions....................... 66 Predicting WT Class Across the Landscape ............................................................. 72 Discussion ..................................................................................................................... 74 References..................................................................................................................... 78 4 CONCLUSION.............................................................................................................. 83 Future Research ............................................................................................................ 84 References..................................................................................................................... 87  v  LIST OF TABLES Table 1.1 An example of the specifications for the Optech ALTM 3100AE discrete return sensor, the Portable Airborne Laser System (PALS) laser profiler, the Laser Vegetation Imaging Sensor (LVIS) airborne full waveform sensor, and the Geoscience Laser Altimeter System (GLAS) space-borne full waveform sensor…………………..…4 Table 1.2 British Columbia Forest and Range Practices Act resource values, selected indicators, and how lidar data may be used to estimate them. Note that several of the resource values are currently under development by the provincial government...............6 Table 1.3 Examples of timber inventory and stand-level biodiversity indicators, including references. Note that the FRPA timber resource value indicators were in development at the time of this writing; the indicators listed below are typical variables measured in the field............................................................................................................7 Table 1.4 Additional indicators of stand-level forest biodiversity, and references to previous studies using lidar technology.................. ………………………………………9 Table 2.1 Interpolation routines and their parameterizations tested in this paper. Digital elevation models with spatial resolutions of 0.5, 1.0, and 1.5 m were created for all combinations, resulting in a total of 48 surfaces..……………………………………......22 Table 2.2 Global statistics summarizing validation errors for selected DEMs (n = 131,852). Only the most accurate parameterizations are shown for spline with tension, regularized spline, ANDUEM, and inverse distance weighting…………………………33 Table 3.1 Examples of previous research where ecological variables were estimated using lidar remote sensing in northwest coastal forests………………………………….52 Table 3.2 Summary statistics for sample plots by age class for stems with a DBH > 10 cm..……………………………………………………………………………………….56 Table 3.3 Lidar sensor and survey parameters…………………………………………..58 Table 3.4 Fit statistics for the ordinal logistic regression model using the lidar-derived variable (log transformed coefficient of variation of the vegetation return heights) as a predictor of wildlife tree class cumulative proportions………………………………….67 Table 3.5 Parameter statistics for the ordinal logistic regression model using the lidarderived variable (log transformed coefficient of variation of the vegetation return heights) as a predictor of wildlife tree class cumulative proportions. The intercept estimates are the coefficients used to predict the cumulative proportions of stems within the wildlife tree classes using Equation 1…………………………………………………………….67  vi  LIST OF FIGURES Figure 2.1 Hillshades derived from 1 m DEMs showing a meander of Lostshoe Creek. Note that all surfaces contain interpolation artefacts, particularly along the stream banks.…………………………………………………………………………………….34 Figure 2.2 Effect of increasing slope on mean absolute error of interpolation. Vertical lines indicate standard errors of the means. Note that the data points have been offset horizontally in order to increase legibility……………………………………………….35 Figure 2.3 Effects of vegetation structural class on mean absolute error of interpolation. Vertical lines indicate standard errors of the means. Note that the data points have been offset horizontally in order to increase legibility…………………………………...……36 Figure 2.4 Effects of ground return density on mean absolute error of interpolation. Vertical lines indicate standard errors of the means…………………………………..…37 Figure 2.5 Normalized difference vegetation indices (NDVI) derived from Quickbird and Landsat7 ETM+ imagery, and lidar ground return density. An increase in ground return density occurs at the highest NDVI values. Vertical lines indicate standard errors of the means…………………………………………………………………………………….38 Figure 2.6 Quickbird image (left), and CART-derived prediction uncertainty map with underlying hillshade for the 0.5 m natural neighbour DEM (right). For this surface, slope and ground return density were the best predictors of interpolation error……………….40 Figure 2.7 Simple classification tree resulting from CART analysis of absolute errors for the 0.5 m DEM created using natural neighbour interpolation..……………………..41 Figure 3.1 Typical field form used to estimate the wildlife tree class of conifers………57 Figure 3.2 Histogram displaying the proportions of WT classes within three forest plots. Note the relative increase in the proportions of the higher WT classes as seral stage increases from pole/sapling to old forest.………………………………………………..65 Figure 3.3 Typical vertical foliage profiles of pole/sapling (n=3) and old forest (n = 3) plots. Foliage is concentrated in the upper canopy of the pole sapling plots, whilst the more structurally complex old growth plots have vegetation distributed evenly through the profile. The higher percentage of dead trees within the old forest plots results in lower canopy densities.…………………………………………………..…………………......66 Figure 3.4 Predicted vs. observed values for cumulative proportions of stems within each wildlife tree class when using the logarithm of CVLidar as a predictor variable (r = 0.85, p <0.000, RMSE = 4.9%).....................................................................................................69  vii  Figure 3.5 Relationship between CVLidar and the predicted individual proportions of stems in each wildlife tree class. As the lidar-derived variable increases, the predicted proportion of healthy trees (i.e. class 1) decreases, whilst the proportions of trees within the other classes increase…….…………………………………………………………..70 Figure 3.6 The lidar-derived height percentiles plotted against the best subsets likelihood scores and the Cox and Snell R2 values. It is the lowest and median height percentiles that account for most of the variance in wildlife tree class..………………………………….71 Figure 3.7 Means and ranges of (1) lidar-derived heights of the 50th percentile, and (2) the percentage of dead trees, grouped by terrestrial ecosystem mapping structural class………………………………………………………………………………………72 Figure 3.8 Map displaying (A) a true colour Quickbird image of the study area (0.80 m spatial resolution); (B) the coefficient of variation (log-transformed) derived from the lidar vegetation returns; and (C) the predicted proportions of trees that are dead within a given pixel (e.g. WT class 3+). Maps B and C have spatial resolutions of 25 m………..73 Figure 3.9 Lidar ground and vegetation returns in pole/sapling (n = 1,794) and old forest (n = 1,969) plots. Dense canopies typical of the pole/sapling structural class occlude the lower canopy and ground, and as a result the majority of lidar returns are clustered in the upper parts of the canopy. Gaps that characterize old forest canopies, however, allow pulse more returns to penetrate to lower parts of the canopy, resulting in lower mean heights and higher variances…………………………………………………………..…76  viii  ACKNOWLEDGEMENTS This research is partly support by the British Columbia Forest Science Program (Y071024), NSERC grants to myself, Dr. Nicholas Coops, and Dr. Sarah Gergel, and funds from the BC Ministry of Forests and Range. I am deeply indebted to Dr. Nicholas Coops for his time, guidance, ideas, insight, editorial skills and honesty. Dr. Denis Collins provided ongoing support and valuable advice on methods and results. Thanks also to Trevor Jones, Shanley Thompson, Kate Kirby, and Sam Coggins for assistance in the field. I would like to thank Dr. Sarah Gergel, Dr. Valerie LeMay, and Dr. Michael Wulder for their advice and suggestions. Finally, I would like to thank the IRSS boys for their friendship and the good times.  ix  CO-AUTHORSHIP STATEMENT This thesis is the combination of (1) two submitted scientific papers for which I am lead author, and (2) extracts from a book chapter and a review paper for which I am second author. The initial project overview was proposed by my supervisor, Dr. Nicholas Coops, from which we identified digital elevation modelling and wildlife tree class estimation as key areas requiring further research. For these scientific journal submissions, I performed all research, data analyses, and interpretation of results, and prepared the final manuscripts. Co-authors provided advice on methodology and made editorial comments as required. Chapter 1 consists of my contributions to the book chapter and review paper.  x  1 INTRODUCTION1 Forests ecosystems cover approximately four billion hectares or 30% of the Earth’s land area, and are critical for conserving biological diversity, water and soil resources, and for providing supplies of wood and non-wood forest products (FAO, 2005; Siry et al., 2005). The wide variety of goods and services (both ecosystem and social) provided by forests has compelled societies, managers, and stewards to develop multi-functional sustainable forestry goals (Cubbage et al., 2007), and dramatically increased the need for information about forest status and condition that is not only spatially explicit and accurate, but also cost-effective and scientifically defensible (Noss, 1999).  Forest sustainability is typically monitored using a suite of indicators related to a variety of resource values, such as timber, soils, floral and faunal biodiversity, riparian areas, water and fish, and so on. Noss (1999) described a variety of indicators that may be used to monitor trends in forests, and their scales range from direct stand-level measurements made in the field, to remote sensing-based landscape-level estimates. According to Breckenridge et al. (1995) and Stone et al. (2000), to be useful, an indicator should: •  be easy to interpret;  •  be correlated with changes in ecosystem processes;  1  The tables in this section are exact reproductions of those I created for a recently submitted review article and a published book chapter. Wulder, M.A., Bater, C.W., Coops, N.C., Hilker, T., and White, J.C., Submitted. The role of lidar in sustainable forest management. The Forestry Chronicle. Wulder, M.A., Bater, C.W., Coops, N.C., Hirata, Y., and Sweda, T., 2007. Advances in laser remote sensing of forests. Invited Expert Commentary / Chapter. In, B.A. Larson (Ed.), Sustainable Development Research Advances. Nova Science Publishers, Inc.  1  •  have regional applicability;  •  show low temporal and spatial variability;  •  be quantifiable using synoptic or automated monitoring;  •  be responsive to change;  •  be anticipatory and provide potential for early warning;  •  have results that can be easily summarized and understood by non-experts; and  •  be cost effective.  The establishment of multi-functional sustainable forestry goals and the associated development of resource values and indicators is occurring in a number of Canadian provinces, including British Columbia. British Columbia contains approximately one half of Canada’s softwood lumber inventory, and in 2005 the forestry industry was responsible for 45% of the province’s manufacturing shipments (BC Stats, 2005). While forestry’s economic benefits are significant, extraction must be performed in a sustainable manner for future generations. In response to this need, the Province of British Columbia, under the auspices of the Forest and Range Practices Act (FRPA), has developed a suite of resource values to measure forest health and sustainability. The Forest and Range Evaluation Program (FREP) has been put in place as a multi-agency program to assess the effectiveness of the FRPA in achieving stewardship of the 11 resource values identified. Each resource value is assessed by monitoring a number of indicators which are traditionally measured using field-based approaches and aerial photography. Field assessments, however, can be expensive, labour intensive, provide small sample sizes and intensities, and often cover only small geographic areas, while  2  aerial photography suffers from time and cost issues, is prone to operator bias and subjectivity, and is limited by a shortage of trained interpreters.  Digital remote sensing systems are an evolving set of tools making available products that meet many of the criteria required of indicators: they provide data covering large areas, often at high spatial resolutions; they are capable of change detection through repeated observations; and they are becoming increasingly cost-effective. Stone et al. (2000), for instance, discussed the development of an indicator integrating a suite of space-borne image products to measure eucalypt forest health. Importantly, most sensors are passive and limited to imaging the upper portions of the forest canopy. Some sensors, however, actively emit radiation, which can improve canopy penetration and more effectively map forest structure. As a key example, light detection and ranging (lidar) employs an airborne laser to map terrain and vegetation simultaneously.  Lidar Background Lidar systems estimate distances between a sensor and a target based on half the elapsed time between a laser pulse emission and the detection of a reflected return. Lidar systems can be separated into two basic types: discrete return and waveform-recording (Lefsky et al., 2002; Lim et al., 2003). Discrete return sensors record single or multiple returns from a given laser pulse. As the laser signal is reflected back to the sensor, large peaks are interpreted to represent discrete objects in the path of the beam (e.g. the forest canopy, understorey, and ground). The sensor then records these peaks as discrete points in threedimensional space. Alternatively, full waveform instruments have a higher sampling rate  3  and record the full height distribution of the surfaces illuminated by the laser. Thus, within a forest canopy, discrete return instruments produce clouds of points representing intercepted surfaces, while full waveform sensors record the entire reflected signal for analysis (Lefsky et al., 2002). Generally, discrete return sensors use a small footprint (e.g. the laser’s circle of illumination on the ground) several decimetres in diameter, while waveform recording sensors employ a large footprint typically greater than 10 m in diameter (Table 1.1). Regardless of the type of system employed, lidar is capable of simultaneously mapping both vegetation height and vertical structure, and the underlying terrain’s morphology, with high accuracy.  Table 1.1 An example of the specifications for the Optech ALTM 3100AE discrete return sensor, the Portable Airborne Laser System (PALS) laser profiler, the Laser Vegetation Imaging Sensor (LVIS) airborne full waveform sensor, and the Geoscience Laser Altimeter System (GLAS) space-borne full waveform sensor. Specification  Optech 3100EA Airborne Discrete Return Scanner  PALS Airborne Laser Profiler  LVIS Airborne Full Waveform Sensor  GLAS Space-borne Full Waveform Sensor  Operating altitude (m)  80-3,500  <300  Typically 10,000  ~600,0000  Wavelength (nm)  1,050  905  1,064  1,064  Number of pulse returns  4, including last  1, first or last  N/A  N/A  Typical footprint diameter  ~0.3 m at 1,000 m flying altitude and 0.3 mrad beam divergence  0.3 m at 150 m flying altitude  ~25 m  ~65 m  Waveform digitization rate  NA  NA  500 Msamp/s  1 GHz  Scan angle  Variable, 0 o to ±25o  0o  ~ 12 o  0o  Maximum laser pulse repetition rate (Hz)  100,000  2,000  500  40  Optech, 2007  Nelson et al., 2003  Blair et al., 1999; NASA, 2007  Schultz et al., 2005  Reference  4  Airborne lidar surveys are typically designed to have a dense and evenly distributed point spacing. In areas of particularly dense forest canopy, however, the ground may be occluded from the sensor, resulting in a dataset containing a large number of vegetation returns and a relative paucity of terrain information. Under such conditions, the selection of an appropriate algorithm for digital elevation model (DEM) interpolation becomes an important decision, especially in uneven terrain, as differences in terrain model heights may directly affect the estimates of vegetation metrics. Nonetheless, lidar has been shown to be an effective remote sensing technology for measuring forestry-related variables.  Use of Lidar to Estimate Indicators of Forest Sustainability Lidar remote sensing offers the ability to accurately assess many indicators at the landscape level. Table 1.2 provides examples of FRPA resource values, some of their indicators, and comments on the capacity to estimate them using lidar. Many of the variables collected during field-based forest inventories may be accurately determined using lidar remote sensing (Table 1.3). Depending on lidar point densities, analyses may be performed at the individual tree- or plot-levels. A variety of tree height metrics may be estimated with accuracies rivalling field-based estimates (Næsset and Økland, 2002; Anderson et al., 2006). For instance, Magnussen and Boudewn (1998) estimated both mean and Lorey’s height of forest plots, while Næsset and Økland (2002) examined tree height and height to the base of crown of individual trees, and Lorey’s mean height and average height to the crown of forest plots.  5  Table 1.2 British Columbia Forest and Range Practices Act resource values, selected indicators, and how lidar data may be used to estimate them. Note that several of the resource values are currently under development by the provincial government. FRPA Resource Values  Indicators  Capacity of lidar to estimate indicator  Notes  Timber  See Table 1.3.  Stand-level biodiversity  See Tables 1.3 and 1.4. Plant community descriptions (e.g. vegetation species, browsing intensity)  Forage and associated plant communities  Visual quality  Limited  Limited ability to estimate vegetation species and browsing intensity.  Stream riparian functions  Moderate  DEM evaluation provides information on channel morphology at the reach scale. Use standard lidar vegetation metrics to monitor vegetation structure. Limited ability to estimate flow regimes, biotic communities, and water quality.  Erosion/deposition  Moderate  DEM evaluation, but scale of enquiry limited by vegetation cover and survey parameterization. Rills and gullies may not be detected.  Visual quality objective  High  Determine tree heights; may be useful for estimating volume removed. When combined with optical data, digital surface model may provide height information for 3D visualizations.  Forest health, invasive plants  Moderate  Identify changes in canopy structure (e.g. defoliation) using vegetation returns. Not useful for invasive plant identification.  Natural hazards  High  DEM evaluation for landslide chutes, slope mapping, tsunami hazard zones, and so on.  Size of retention areas surrounding karst resources  High  Map size and height using vegetation returns. Possible to estimate amounts of windthrow and harvesting within retention areas.  Soil disturbances, logging slash amounts  None  Condition of plant community  Moderate  Livestock management practices  None  Monitoring culturally sensitive natural resources, resource gathering areas  Moderate  Identify culturally modified trees  None  Point indicators (e.g. % moss, # of insect types)  None  Channel morphology  Moderate  Recreation resources  Resource features  Use standard lidar vegetation metrics to monitor riparian vegetation structure.  Water  Cultural heritage  Fish/riparian  Use standard lidar vegetation metrics to monitor forest structure.  DEM evaluation, map planform and cross-sectional morphology. Not useful for analysing sediments, bank characteristics and other channel characteristics.  6  Table 1.3 Examples of timber inventory and stand-level biodiversity indicators, including references. Note that the FRPA timber resource value indicators were in development at the time of this writing; the indicators listed below are typical variables measured in the field.  Indicator  Lidar sensor a  Sampling density (pulses/m2)  Individual tree, plot- or polygonbased  Species  Reference  2  Not applicable (waveform instrument; 10 m footprint)  Plot  Douglas-fir, western hemlock  Lefsky et al., 1999  1  1.2  Plot  Norway spruce, Scots pine  Næsset, 2002  Diameter  1  0.2  Plot  Douglas-fir  Magnussen and Boudewyn, 1998; Magnussen et al., 1999  1  0.6-2.3  Individual tree, plot  Norway spruce, Scots pine  Næsset and Økland, 2002  2  Not applicable (waveform instrument; 9 m footprint)  Plot, polygon  Black spruce, jack pine, aspen  Wulder and Seemann, 2003  1  6  Individual tree  Douglas-fir, ponderosa pine  Anderson et al., 2006  1  4-5  Individual tree  Norway spruce, Scots pine, silver birch, downy birch  Maltamo et al., 2004  Height  Volume Pathological indicators (e.g. scars, broken tops)  Not recommended using lidar  Growth  1  10  Individual tree  Scots pine  Yu et al., 2006  Species  1  4.7  Individual tree  Norway spruce, Scots pine  Holmgren and Persson, 2004  Height to base of crown  1  0.6-2.3  Individual tree, plot  Norway spruce, Scots pine  Næsset and Økland, 2002  Popescu et al., 2003  Næsset, 2002  Crown diameter  Number of trees per hectare  1  1.4  Individual tree  White oak, chestnut oak, northern red oak, southern red oak, yellow poplar, red maple,/Virginia pine, loblolly pine, shortleaf pine, pignut hickory, scarlet oak, black oak, blackgum, American beech  1  1.2  Plot  Norway spruce, Scots pine  Individual Norway spruce, Scots Pine, silver tree birch a 1 = Small footprint, discrete return; 2 = large footprint, full waveform 1  4-5  Maltamo et al., 2004  7  Predicting diameter and volume can be performed by relating back to height, typically by using regression models (e.g. Næsset, 2002) or other empirical allometric equations. Height, diameter, volume at the plot level, and to some degree, stocking densities, may be considered ready for operational application, while others, such as species determination, are still the subject of academic research and may not be feasible without supplemental data such as high spatial resolution optical imagery.  The assessment of forest biodiversity includes the examination of many of the same variables collected during typical forest inventories, but may be augmented by the analysis of additional species and structural information (Tables 1.3 and 1.4). The horizontal and vertical organization of forest canopies can provide managers with information relating to the development of plant communities, canopy function, and related habitat conditions for wildlife (Lefsky et al., 2002). It has been established that lidar can provide information pertaining to crown closure and canopy volume (e.g. Lefsky et al., 1999; Coops et al., 2007), while the estimation of windthrow, coarse woody debris, and wildlife tree class distribution is a matter of ongoing research. The identification of invasive plants and ecological anchors is not likely possible using lidar, and will continue to require field-based surveys.  When assessing water quality, lidar data may provide important ancillary information for field surveys, but most of the associated indicators must be measured on the ground (Table 1.2). Digital elevation models, for instance, may be used to map indicators such as  8  larger landslide events, but smaller-scale processes such as gullying may be difficult to resolve unless ground return densities are sufficiently high.  Table 1.4 Additional indicators of stand-level forest biodiversity, and references to previous studies using lidar technology. Indicator  Lidar sensor a  Sampling density (pulses/m2)  Individual tree, plot- or polygon-based  Species  See Table 1.3  Diameter  See Table 1.3  Height  Reference  See Table 1.3 2  Not applicable (waveform instrument; 10 m footprint)  Plot  Douglas-fir, western hemlock  Lefsky et al., 1999  1  0.7  Plot  Douglas-fir, western hemlock, western red cedar, red alder  Coops et al., 2007  2  Not applicable (waveform instrument; 20 m footprint, 20 m posting distance)  Polygon  Black oak, red maple, northern red oak, hemlock, red pine, white pine, white spruce  Weishampel et al., 2007  Crown closure, canopy volume  Windthrow  Species  Ecological anchors (e.g. bear den, cavity nest, wildlife trails)  Not recommended using lidar remote sensing  Invasive plants  Not recommended using lidar remote sensing  Coarse woody debris  Wildlife tree class  1  1  1.4  0.7  Plot  Lodgepole pine, unspecified spruce/fir  Seielstad and Queen, 2003  Plot  Western hemlock, western red cedar, amabilis fir, yellow-cedar, Sitka spruce, Douglas-fir, red alder  Bater et al., In review.  a  1 = Small footprint, discrete return; 2 = large footprint, full waveform  The technology can provide valuable information on riparian vegetation and changes in channel morphology at the reach scale or greater, but is of little utility for assessing indicators such as soil compaction or the impacts of livestock grazing. For the same reasons, natural hazard mapping is an area where laser scanning could be of great benefit,  9  as it provides accurate models of landscape morphology, typically with sub-metre vertical accuracy. Conversely, as mentioned earier, mapping channel erosion and deposition at the finest of scales is most likely unrealistic and will continue to require detailed field surveys.  Management of visual quality under FRPA seeks to ensure forest operations are consistent with public expectations. Management of this resource is governed by Visual Quality Objectives (VQOs) which are management objectives reflecting the public’s desired level of visual quality for a landscape. Lidar data can be used to produce very accurate topographic information and can be used to determine mean tree heights, information central to the visual impact assessment process. Specifically, the laser data can be used to enhance the accuracy of the 3D visualisations, helping forest managers assess whether proposed forest operations will achieve VQOs with greater confidence.  Objectives A review of the academic literature suggests that while many indicators of forest sustainability have been previously estimated using lidar remote sensing, certain areas still require additional investigation. Two of these include DEM creation, and the structural stage or wildlife tree class distributions within forest stands.  Digital elevation models are key derivatives of lidar data. Not only do they provide information on terrain morphology, but they also act as the reference to mean sea level against which vegetation heights are determined. Thus, literally all lidar-derived  10  variables, be they related directly to terrain, or to estimates of vegetation heights, will be affected by the accuracy of the DEM. The goals of Chapter 2 were threefold: first, by examining seven interpolation routines, identify the most accurate combinations of algorithm and spatial resolution for DEM creation; second, examine the effects of lidar ground return density, slope, vegetation cover, and vegetation structural class on interpolation errors; and third, use classification and regression tree analysis (CART) to identify important independent variables for error prediction, and create a rule-based classification for the derivation of prediction uncertainty maps.  Determining the distribution of living and dead trees, and the state of decay of the latter, is important for assessing the status of forest ecosystems. In British Columbia, the structural stage of living and dead trees is referred to as the decay class or wildlife tree (WT) class of a stem. Within the area of interest of this research, dead and decaying trees are crucial structural components of canopy architecture, and provide habitat for a variety of biota. The goal of Chapter 3 was to estimate the distributions of living and dead tree forms within forests by developing statistical relationships between plot-level distributions of living and dead trees, and commonly derived lidar vegetation metrics.  Finally, in Chapter 4, the previous chapters will be related to each other, and the significance of this work in relation to the field of lidar remote sensing and the estimation of indicators of forest sustainability will be discussed; the chapter will close with suggestions for future research.  11  References Anderson, H.E., Reutebuch, S.E., McGaughey, R.J., 2006. A rigorous assessment of tree height measurements obtained using airborne lidar and conventional field methods. Canadian Journal of Remote Sensing 32(5), 355-366. Bater, C.W., Coops, N.C., Gergel, S.E., Collins, D., In review. Estimating the distribution of tree structural classes in northwest coastal forests using lidar remote sensing. Remote Sensing of Environment. BC Stats, 2005. Quick facts about British Columbia. Ministry of Labour and Citizen’s Services. http://www.bcstats.gov.bc.ca/data/qf.pdf Accessed January 2008. Blair, J. B. Rabine, D. L., Hofton, M. A., 1999. The laser vegetation imaging sensor: a medium-altitude, digitisation-only, airborne laser altimeter for mapping vegetation and topography. ISPRS Journal of Photogrammetry & Remote Sensing 54, 115122. Breckenridge, R.P., Kepner, W.G., Mouat, D.A., 1995. A process for selecting indicators for monitoring conditions of rangeland health. Environmental Monitoring and Assessment 36, 45-60. Coops, N.C., Hilker, T., Wulder, M.A., St-Onge, B., Newnham, G., Siggens, A., Trofymow J.A., 2007. Estimating canopy structure of Douglas-fir forest stands from discrete-return LiDAR. Trees: Structure and Function 21(3), 295-310. Cubbage, F., Harou, P., Sills, E., 2007. Policy instruments to enhance multi-functional forest management. Forest Policy and Economics 9, 833-851.  12  FAO, 2005. Global Forest Resources Assessment 2005: Progress Towards Sustainable Forest Management. FAO forestry paper 147. ftp://ftp.fao.org/docrep/fao/008/A0400E/A0400E00.pdf. Accessed January 2008. Holmgren, J, Persson, A., 2004. Identifying species of individual trees using airborne laser scanner. Remote Sensing of Environment 90, 415-423. Lefsky, M.A., Cohen, W.B., Acker, S.A., Parker, G.G., Spies, T.A., Harding D., 1999. Lidar remote sensing of the canopy structure and biophysical properties of Douglasfir western hemlock forests. Remote Sensing of Environment 70, 339-361. Lefsky, M.A., Cohen, W.B., Parker, G.G., and Harding, D.J., 2002 Lidar remote sensing for ecosystem studies. BioScience 52(1),19-30. Lim, K., Treitz, P., Wulder, M., St-Onge, B., Flood, M., 2003. Lidar remote sensing of forest structure. Progress in Physical Geography 27(1), 88-106. Magnussen S., Boudewyn, P., 1998. Derivations of stand heights from airborne laser scanner data with canopy-based quantile estimators. Canadian Journal of Forest Research 28, 1016-1031. Magnussen, S., Eggermont, P., LaRiccia, V.N., 1999. Recovering tree heights from airborne laser scanner data. Forest Science 45(3): 407-422. Maltamo, M., Eerikäinen, K., Pitkänen, J., Hyyppä, J., Vehmas M., 2004. Estimation of timber volume and stem density based on scanning laser altimetry and expected tree size distribution functions. Remote Sensing of Environment 90, 319-330. NASA, 2007. Goddard Space Flight Center, Laser Vegetation Imaging Sensor, About LVIS. https://lvis.gsfc.nasa.gov/index.php. Accessed January 2008.  13  Næsset, E., 2002. Predicting forest stand characteristics with airborne scanning laser using a practical two-stage procedure and field data. Remote Sensing of Environment 80, 88-99. Næsset, E., Økland, T., 2002. Estimating tree height and tree crown properties using airborne scanning laser in a boreal nature reserve. Remote Sensing of Environment 79, 105-115. Nelson, R., Parker, G., Hom, M., 2003. A portable airborne laser system for forest inventory. Photogrammetric Engineering and Remote Sensing 69(3), 267-263. Noss, R.F., 1999. Assessing and monitoring forest biodiversity: a suggested framework and indicators. Forest Ecology and Management 115, 135-146. Optech. ALTM 3100AE Brochure. Optech Incorporated. February 2007. http://www.optech.ca/pdf/Brochures/ALTM3100EAwspecsfnl.pdf. Accessed March 2008. Popescu, S.C., Wynne, R.H., Nelson, R.F., 2003. Measuring individual tree crown diameter with lidar and assessing its influence on estimating forest volume and biomass. Canadian Journal of Remote Sensing 29(5), 564-577. Schultz, B.E., Zwally, H.J., Shuman, C.A., Hancock, D., DiMarzio, J.P., 2005. Overview of the ICESat mission. Geophysical Research Letters 32, L21S01. DOI: 10.1029/2005GL024009. Seielstad, C.A., Queen, L.P., 2003. Using airborne laser altimetry to determine fuel models for estimating fire behavior. Journal of Forestry 101(4), 10-15.  14  Siry, J.P., Cubbage, F.Q., Ahmed, M.K., 2005. Forestry Policy and Economics 7, 551561. Stone, C., Coops, N., Culvenor, D., 2000. Conceptual development of a eucalypt canopy condition index using high resolution spatial and spectral remote sensing imagery. Journal of Sustainable Forestry 11(4), 23-45. Weishampel, J.F., Drake, J.B., Cooper, A., Blair, J.B., Hofton, M., 2007. Forest canopy recovery from the 1938 hurricane and subsequent salvage damage measured with airborne LiDAR. Remote Sensing of Environment 109, 142-153. Wulder, M.A., Seemann D., 2003. Forest inventory height update through the integration of lidar data with segmented Landsat imagery. Canadian Journal of Remote Sensing 29(5), 536-543. Yu, X., Hyyppä, J., Kukko, A., Maltamo, M., Kaartinen, H., 2006. Change detection techniques for canopy height growth measurements using airborne laser scanner data. Photogrammetric Engineering & Remote Sensing 72(12), 1339-1348.  15  2 EVALUATING ERROR ASSOCIATED WITH LIDARDERIVED DEM INTERPOLATION2 Introduction Light detection and ranging (lidar) technology is an active remote sensing technique capable of simultaneously mapping the Earth’s surface and overlying features, including vegetation and buildings, with sub-metre vertical accuracy. Small footprint, discrete return airborne lidar typically employs a laser scanner slaved to an inertial measurement unit (IMU) and a global positioning system (GPS) for accurate measurements of aircraft orientation and position. The laser scanner emits discrete pulses to the ground from which multiple reflections or returns can be detected, allowing for the simultaneous mapping of the ground, vegetation, and other features. The distance estimates are based on the time between pulse emission and return detection.  Lidar data has been applied across a variety of disciplines, including archaeology, structural geology, geomorphology, engineering, resource management, and disaster assessment and planning. Recently, lidar has been gaining recognition in forestry and ecology as an effective tool for estimating a variety of vegetation metrics, including tree heights, biomass, crown size, leaf area index (LAI), and vertical canopy structure (e.g. Næsset, 1997; Lefsky et al., 1999; Dubayah and Drake, 2000; Lim et al., 2003; Coops et al., 2004). Regardless of the variable under investigation, estimates are typically based on their height above a continuous digital elevation model (DEM) representing the Earth’s bare surface. Airborne lidar surveys are typically designed to have a dense and evenly distributed point spacing. In higher leaf area canopies, however, ground visibility is 2  A version of this chapter has been submitted for publication. Bater, C.W. and Coops, N.C., In Review. Evaluating error associated with lidar-derived DEM interpolation. Computers & Geosciences.  16  reduced, resulting in datasets containing a large number of vegetation returns and a relative paucity of terrain information. This will have implications for not only the quality of derived DEMs and their representation of terrain morphology, but also for the accurate estimation of vegetation metrics.  Previous research has shown that the accuracy of DEMs varies with changes in terrain and land cover type (e.g. Adams and Chandler, 2002; Hodgson and Bresnahan, 2004; Hodgson et al., 2005; Su and Bork, 2006). By surveying ground returns at their x and y locations, Hodgson and Bresnahan (2004) decomposed lidar error into four components. Specifically and in decreasing order of importance, these included: lidar system measurements, interpolation error, horizontal displacement error, and survey error.  Adams and Chandler (2002) assessed the accuracy of a lidar-derived DEM of the Black Ven mudslide in Dorset, United Kingdom. The DEM had a spatial resolution of 2 m, although no information on the density of ground returns was available. Validation data consisted of a DEM generated from survey-grade points. Adams and Chandler (2002) reported an overall root mean square (RMS) error of 0.26 m, and found that the lidar data tended to increasingly underestimate terrain elevation as slope increased.  Hodgson et al. (2005) examined the effects of land cover and slope on DEM accuracy. Located in a watershed in North Carolina, USA, their study area consisted of gently rolling terrain. Land cover classes included grass and scrub/shrub, and pine, deciduous, and mixed forests. Lidar data were collected during leaf-off conditions, with an average  17  ground return posting distance of one point every 31.1 m (corresponding to density of 0.03 points/m2). Slope was then modelled by linear interpolation of a triangulated irregular network (TIN). Reference data consisted of 1225 survey grade points collected along 23 transects, and reference slope was calculated as the average slope of adjacent segments along survey transects. Hodgson et al. (2005) reported RMS errors of 0.145 m to 0.361 m for the different land cover classes, with higher errors occurring in areas with tall canopy vegetation. The scrub/shrub class, however, exhibited the largest RMS error. Little evidence was found for increased elevation errors in areas with slopes from 0o to 10o, but lidar-derived slope was generally under-predicted as terrain slope increased.  The selection of an appropriate algorithm and spatial resolution for DEM interpolation may become an important decision, especially in uneven terrain, as differences in terrain model heights may directly affect the estimates of vegetation metrics. Interpolation methods can be broadly defined as being deterministic or probabilistic (Maune et al., 2001). Deterministic methods are based only on the values of surrounding points, with algorithms using mathematical formulae to determine the influence of immediate neighbouring values. Inverse distance weighting, for example, applies a weighting that is derived from the distance of surrounding points; as a result, points that are close to the prediction location have more influence than points that are more distant. Probabilistic geostatistical methods rely on spatial autocorrelation between proximal points, and can account for both distance and direction when determining the importance of surrounding values (Maune et al., 2001). Kriging, for example, relies on the calculation of a semivariogram, which graphs the distance between point pairs and their difference in z-  18  values (e.g. lidar return heights). A model is then fit through the data from which weights for interpolation are obtained.  Interpolation algorithms vary widely in their complexity, ease of use, and computational expense. Previous authors have investigated DEM interpolation routines with varying results. Using lidar data collected over smoothly varying hillslopes near Humberside, United Kingdom, Lloyd and Atkinson (2002) employed cross-validation and a jack-knife approach (e.g. subsets of the lidar data) to test inverse distance weighted (IDW) interpolation and two types of kriging. Their original dataset consisted of 139 694 returns within a 500 m by 500 m area (corresponding to a mean point density of 0.56 pulses/m2). To assess the ability of the interpolation algorithms to predict heights where point densities were low, Lloyd and Atkinson (2002) systematically decimated their data, validating surfaces created using 50% (0.28 pulses/m2), 25% (0.14 pulses/m2), and 5% (0.03 pulses/m2) of the original dataset. Lloyd and Atkinson (2002) found that kriging was the more accurate method when point densities were low, but concluded that no benefit was derived from using the more sophisticated geostatistical approach when large amounts of data were available.  Yanalak (2003) examined 1 m resolution DEMs created using linear, natural neighbour, nearest neighbour, weighted average polynomial, minimum curvature and multiquadratic routines. The gridding methods were applied to five 1 ha test datasets, each containing 150 scattered reference points digitized from five theoretical test surfaces; the true heights of the test profile data were calculated using surface equations. Yanalak (2003)  19  found that the minimum curvature and multiquadratic routines were the most accurate based on test profiles across the surfaces.  Abrammov and McEwan (2004) tested linear interpolation, splining, nearest neighbour, and natural neighbour routines using regularly-spaced Mars Orbiter Laser Altimeter (MOLA) data collected over the Martian Korolev crater, and a simulated MOLA dataset collected in Iceland. While the former analysis was strictly qualitative, the latter employed a DEM produced by the Icelandic Geodetic Survey, the heights of which were sampled to simulate a typical MOLA survey. The point data were interpolated at three spatial resolutions using the four algorithms, and then the new DEMs were compared to the original. Abrammov and McEwan (2004) concluded that natural neighbour interpolation was not only the most accurate algorithm, but also generated the fewest visual artifacts. The most accurate spatial resolution varied by algorithm, with the high resolution surfaces being most accurate for natural neighbour and nearest neighbour, the medium resolution DEM for linear, and the low resolution DEM for splining.  Su and Bork (2006) tested DEMs generated from densely arranged lidar returns (densities averaged 0.75 points/m2). Using a jackknife approach similar to that of Lloyd and Atkinson (2002), Su and Bork (2006) validated DEMs created using splining, IDW, and kriging algorithms. They found that IDW was the most accurate interpolator, with RMS errors 0.02 m lower than those for splining and kriging. Su and Bork (2006) also collected reference points using a total station and a differential global positioning system (DGPS) to assess the effects of slope, vegetation and laser scan angle. They found that  20  the IDW-derived DEM overestimated heights in forested areas, and underestimated them in lowland meadows. Accuracies diminished as terrain slope angle increased, with RMS errors doubling in areas where slopes exceeded 10o, compared to areas where slopes were less than 2o. Laser scan angle, however, had little effect on DEM accuracy below 15o off nadir.  Based upon the results of this previous research, it is apparent that no interpolation method is universally superior. Ground return spacing (which is a function of lidar survey parameterization and vegetation density), raster spatial resolution, the complexity of terrain morphology, and the assumptions underpinning a given interpolator’s mathematical design affect the ability of interpolation algorithms to generate accurate DEMs.  In order to gain a better understanding of the potential errors introduced by these and other factors, the capacity to generate prediction uncertainty maps can provide the ability to visually explore and make comparisons between different prediction methods (e.g. Hengl et al., 2004). Prediction uncertainty maps can be a valuable tool for examining the spatial distribution of potential errors and their magnitude, and for identifying areas where a DEM may be of low quality. Thus, a simple method is needed to generate maps of prediction uncertainty for deterministic interpolation routines.  Within the context of the study area’s vegetative and topographic variability, the goals of this research were threefold: first, by examining seven interpolation routines, including  21  linear, quintic, natural neighbour, spline with tension, regularized spline, IDW, and ANUDEM (Table 2.1), identify the most accurate combinations of interpolation algorithm and spatial resolution for DEM creation; second, examine the effects of lidar ground return density, slope, vegetation cover, and vegetation structural class on interpolation errors; and third, use classification and regression tree analysis (CART) to identify important independent variables for error prediction, and create a rule-based classification for the derivation of prediction uncertainty maps. The routines selected are commonly used for terrain modelling, and are available in most geographic information system and image analysis software packages.  Table 2.1 Interpolation routines and their parameterizations tested in this paper. Digital elevation models with spatial resolutions of 0.5, 1.0, and 1.5 m were created for all combinations, resulting in a total of 48 surfaces. Type  Parameterization  TIN  References  Linear  n/a  Yes  ESRI, 2005  Quintic  n/a  Yes  ESRI, 2005  Natural Neighbour  n/a  Yes  Sibson, 1981; Sambridge et al., 1995  Spline with tension and regularized spline  Weights of 0, 0.1, 1.0 and 2.0 were tested  Inverse Distance Weighting  ANUDEM  No  Mitášová and Hofierka, 1993; Johnston et al., 2001  Powers of 0.5, 2.0 and 3.0 were tested  No  Johnston et al., 2001  Surfaces were created with drainage enforcement both on and off  No  Hutchinson, 1989, 1996  22  Methods Study Area The study area is a 4 x 7 km area in the Kennedy Flats, Clayoquot Sound on Vancouver Island, British Columbia (49o 0’ 35” N, 125o 37’ 21” W), approximately 15 km south of the town of Tofino. The area is classified as Coastal Western Hemlock (CWH) zone, based on the biogeoclimatic ecosystem classification (BEC) system (Meidinger and Pojar, 1991). Although flanked by the Vancouver Island Range, the topography is subdued and dominated by Pleistocene glacial deposits. Streams are generally meandering and alluvial, but most have short sections that show evidence of constraint. From 1971 to 2000, Tofino received an average of 3306 mm of precipitation each year, while mean daily minimum, average and maximum temperatures were 5.4oC, 9.1oC, and 12.8oC, respectively (Environment Canada, 20063).  Clayoquot Sound is a multi-use area and includes both recently harvested Crown land, and mature first and second growth forest in Pacific Rim National Park. Western hemlock (Tsuga heterophylla) is a dominant or codominant tree species throughout; western redcedar (Thuja plicata), amabilis fir (Abies amabilis), yellow-cedar (Chamaecyparis nootkatensis), Sitka spruce (Picea sitchensis), Douglas-fir (Pseudotsuga menziesii var. mensiesii), and red alder (Alnus rubra) also occur under differing conditions. The area has been previously mapped using the Province of British Columbia’s Terrestrial Ecosystem Mapping (TEM) classification system, which is derived from 1:20 000 to 1:50 000 scale aerial photography (Mitchell et al., 1989; Demarchi et al., 1990). A hierarchical classification system, TEM mapping stratifies the 3  Environment Canada, 2006. Canadian Climate Normals 1971-2000, Tofino A, British Columbia. http://www.climate.weatheroffice.ec.gc.ca/climate_normals/index_e.html  23  landscape into map units based on climate, physiography, surficial material, geology, soil and vegetation. Based on the TEM classification system, the study area encompasses a range of vegetation structural stages, including herb (1% of total area) shrub/herb (14%), pole/sapling (32%), young forest (4%), and old forest (46%).  Data Acquisition Small footprint laser data were acquired in July 2005 by Terra Remote Sensing (Sidney, British Columbia, Canada) using a Mark II discrete return sensor. Ground and nonground returns were separated using Terrascan v 4.006 (Terrasolid, Helsinki, Finland). Pulse frequency, flight speed, and altitude were optimized to achieve a nominal ground return spacing of one laser pulse return every 1.5 m or 0.7 returns/m2; ground return posting distances, however, were up to 20 m in the dense forest cover classes such as pole/sapling and young forest.  In order to estimate the vegetation cover across the study area, two sets of optical satellite data were obtained. First, high spatial resolution Quickbird data were acquired in June 2005 by DigitalGlobe (Longmont, Colorado, USA). The imagery consisted of four multispectral bands with a spatial resolution of 2.8 m in the blue, green, red, and near-infrared portions of the electromagnetic spectrum, with wavelengths of 0.45-0.52, 0.52-0.60, 0.630.69, and 0.76-0.90 µm, respectively. In addition, a moderate spatial resolution Landsat7 Enhanced Thematic Mapper plus (ETM+) scene, acquired during September 2000, was downloaded from the University of Maryland’s Global Land Cover Facility4. The image had been previously processed to the level 1G standard, which included radiometric and 4  University of Maryland, Global Land Cover Facility. http://glcf.umiacs.umd.edu/data/  24  geometric correction. The multispectral data had a spatial resolution of 28.5 m, and consisted of three visible bands and three infrared bands with wavelengths of 0.45-0.515, 0.525-0.605, 0.63-0.69, 0.75-0.90, 1.55-1.75, and 2.09-2.35 µm. Both the Quickbird and Landsat7 ETM+ data were acquired during cloud-free conditions while the forest canopy was in leaf-on condition.  Data Processing DEM Validation Of the interpolation algorithms tested, the linear, quintic, and natural neighbour methods employ TINs. A TIN is a vector terrain model using Delaunay triangles to join points in three-dimensional space (e.g. lidar ground returns). Linear interpolation calculates a surface based only on the coordinates of the Delaunay triangle within which a point is located, while the quintic method relies on the geometry of adjacent triangles and a fifthdegree polynomial to fit a smooth surface (ESRI, 2005). Natural neighbours employ a TIN to identify adjacent points, and Voronoi (Thiessen) polygons to determine the influence or weighting of each point (Sibson, 1981; Sambridge et al., 1995; Maune et al., 2001). Unlike linear interpolation, which is based on an unknown point’s distance from surrounding triangle nodes, natural neighbours determines the influence of surrounding points based on the proportion of overlap between an unknown point’s Voronoi region and the Voronoi regions of its neighbours.  In contrast to the TIN-based routines, the IDW, spline and ANUDEM algorithms are applied directly to a set of measured points, and do not rely on the generation of a TIN. Inverse distance weighted interpolation employs a weight that is an inverse function of  25  distance, so that adjacent values decrease in significance with increasing distance. Splines are conceptually comparable to fitting a thin flexible plate through a set of measured points, relying on a mathematical function to interpolate a smooth surface through a series of observations (Mitášová and Hofierka, 1993; Maune et al., 2001). Although the spline routines tested here are exact, meaning that the interpolated surfaces pass through the data points, some spline methods are also capable of smoothing, so that the interpolated surface passes near to rather than through the measured points. The only method tested here that is specifically intended for terrain modelling, the ANUDEM algorithm can use both point, line and polygon data to generate DEMs with realistic drainage characteristics (Hutchinson, 1989, 1996). The method is iterative, employing a finite difference interpolation algorithm that maintains a realistic drainage network. ANUDEM includes an optional drainage enforcement routine which attempts to remove all unidentified sinks or depressions. Most DEMs were generated using ArcGIS 9.1 Desktop; the ANUDEM and quintic routines were applied using ArcInfo workstation 9.2 (Redlands, California, USA). The implementation of the interpolation routines are detailed in Johnston et al. (2001) and ESRI (2005).  Lidar ground returns from across the entire study area were randomly subsetted into prediction and validation datasets consisting of 4 500 000 (97%) and 130 000 (3%) points, respectively. This proportion of prediction to validation data was used in order to test the quality of the DEMs without compromising the integrity of the lidar data. The prediction data were used to generate DEMs, while the validation data were used as an independent dataset to estimate vertical interpolation errors.  26  Critically, the exercise was not intended to test the absolute geodetic accuracy of the predicted surfaces, as the validation data were subject to the same degree of positional error as the prediction data. Instead, the validation data were intended to assess the ability of various interpolation routines to predict heights at locations lacking measured points in the prediction dataset. Thus, the results must be tempered because of the limitations inherent in the validation data.  In total, 48 DEMs were produced using the described interpolation routines at three spatial resolutions. The spline and IDW interpolators were parameterized to use the 12 nearest points in the vicinity of the location being estimated. Regularized spline and spline with tension were tested with weights set to 0, 0.1, 1, and 2; for IDW, surfaces were generated using power parameters of 0.5, 2, and 3. ANUDEM was tested both with and without drainage enforcement enabled. Hengl (2006) describes rules for selecting appropriate grid resolutions for interpolating point data, including those based on GPS horizontal error, map scale, point pattern geometry, and spatial autocorrelation. The lidar ground returns in the prediction dataset were randomly distributed, with a mean density of 0.12 points/m2 (corresponding to a mean posting distance of 2.9 m). The minimum distance between points, however, was 0.3 m, while 5% of the study area had point densities ≥ 0.5 points/m2, and 1% of the study area had point densities ≥1 points/m2. Applying Hengl’s (2006) rules indicated that spatial resolutions of 0.2 - 1.7 m were appropriate for the lidar data. Thus, for each combination of `interpolator and  27  parameterization, DEMs were generated with spatial resolutions of 0.5 m, 1.0 m, and 1.5 m (corresponding to an average of 0.04, 0.15 and 0.33 points/cell, respectively).  Slope To investigate the effects of slope on differences in interpolation routines and spatial resolution, a TIN was generated using all lidar ground returns in the dataset (e.g. prior to subsetting the data). A vector slope map was then derived from the TIN and reclassified into 10 categories, specifically: 0-1o, 1-2 o, 2-5 o, 5-10 o, 10-15 o, 15-20 o, 20-25 o, 25-30 o, 30-35 o and 35-40o. Although slopes of up to 85 o existed, the majority of ground returns were found between 0 o and 40 o.  Ground Return Density Lidar ground return density was calculated by counting the number of ground returns within each 1 m2 in the study area, and generating a raster dataset. Average filters were then applied to smooth the uneven distribution of ground returns and better characterize point density in a given area. A number of neighbourhood sizes were tested, including 3x3, 5x5, 7x7, 9x9, 21x21, and 31x31 pixel moving windows, resulting in six additional raster datasets mapping local mean ground return density across the study area.  Satellite Data Prior to analysis, the 2.8 m Quickbird image required both radiometric calibration and geometric correction. First, the multispectral bands were converted to top of atmosphere radiance units using calibration data accompanying the image (RSI 2006). Second, the  28  image was orthorectified both for accurate registration to the lidar data, and to remove distortions caused by terrain and the optical sensor. To do this, the lidar’s laser pulse return intensity values were interpolated into a continuous grayscale image, somewhat similar in appearance to an aerial photograph. Ground control points (GCPs) were then collected from identical locations in the lidar intensity image and the Quickbird scene. Using the GCPs and the lidar-derived DEM, the Quickbird image was then orthocorrected to match the lidar data, resulting in a planimetric image and precise spatial alignment between the two datasets. The coarser resolution Landsat7 ETM+ scene was not orthorectified to the lidar data, as the alignment between the two datasets was sufficiently accurate when compared together.  In order to provide an indication of vegetation cover throughout the study area, normalized difference vegetation indices (NDVI) (Tucker, 1979) were derived from both the Quickbird and Landsat7 ETM+ images. This commonly applied index is a well established as an indicator of LAI and vegetation cover, but is prone to saturation at high LAI values (Huete et al., 2002; Liang, 2004; Jensen, 2007). NDVI is calculated as:  NDVI =  NIR − red NIR + red  (1)  where NIR and red are bands in the near infrared and red portions of the electromagnetic spectrum, respectively.  29  Error Analysis Once the 48 DEMs were developed using the prediction dataset, the validation data were used to compare the biases and accuracies of the surfaces. For each DEM, vertical errors were calculated for every point in the validation data using the following equation:  E ( si ) = P( si ) − M ( si )  (2)  where E(si) is the error at location si, P(si) is the predicted value of the DEM at location si, and M(si) is the measured value from the validation lidar data at location si. Note that the z value of each DEM cell represents the estimated height at the cell centre, while validation points are unlikely to occur there. As a result, some additional error may be introduced when extracting heights from the interpolated surfaces, as it must be incorrectly assumed that elevation is constant across the surface of each cell. Global statistics were calculated to assess the overall performance of the interpolation routines. Mean error is a measure of bias, with positive values indicating the interpolation algorithm is overestimating the actual data. Root mean square error and mean absolute error (e.g. unsigned error) were calculated to measure the global accuracies of the surfaces.  In order to determine the influence of terrain, vegetation, and point density on interpolation errors for the various routines, the study area was stratified using the independent predictor variables, including slope, TEM vegetation structural class, Quickbird and Landsat7 NDVI, and lidar ground return density. To account for potential  30  misregistration of the TEM polygons and Landsat7 data, all lidar returns within 15 m of TEM structural class boundaries were removed. Mean absolute errors of interpolation were then summarized across the range of values for each independent variable.  Error Prediction Classification and regression tree analysis (Breiman et al., 1984; Steinberg and Colla, 1995) was used to develop maps of prediction uncertainty. A non-parametric classification tool using binary recursive partitioning, CART employs decision tree algorithms to determine how continuous and/or categorical variables can be used to predict or classify a dependent variable. Classification and regression tree analysis was chosen because it makes no assumptions about input data or their statistical distributions, and it is well suited to handling collinear datasets, outliers, and potentially insignificant predictors (Melendez et al., 2006; Schwalm et al., 2006). The output from a CART analysis is a series of logical if-then conditions ending in terminal nodes predicting the value of the response variable. Classification and regression tree analysis has commonly been applied to remote sensing data as a rule-based classification design (e.g. Lawrence and Wright, 2001; Brown de Colstoun et al., 2003; Coops et al., 2006).  Terrain slope, the Quickbird-derived NDVI image, the six mean lidar ground return density estimates, and TEM structural class were used as independent variables in a CART analysis. First, CART analysis was used as an exploratory data technique to identify those variables having the most influence on DEM error. Second, the conditional rules produced by the analysis were implemented using map algebra to generate categorical maps of prediction uncertainty.  31  Results DEM Validation Global statistics for the DEM validation are presented in Table 2.2. For all interpolation algorithms and all spatial resolutions, mean errors were sub-centimetre, indicating that interpolation biases were negligible. Root mean square errors ranged from 0.17 to 0.25 m, decreasing for all interpolation routines as spatial resolution increased from 1.5 to 0.5 m. Global mean absolute errors were also very consistent between the algorithms and varied more between spatial resolutions. The overall ranges of error were greatest for the IDW and splining algorithms, with the more conservative linear and natural neighbour interpolators having the lowest ranges. ANUDEM generally performed well, however, enabling drainage enforcement compromised its performance. For instance, RMS errors for the 1.5 m and 0.5 m drainage-enforced DEMs were 0.23 m and 0.22 m, respectively.  32  Table 2.2 Global statistics summarizing validation errors for selected DEMs (n = 131,852). Only the most accurate parameterizations are shown for spline with tension, regularized spline, ANDUEM, and inverse distance weighting. Interpolation Method (Weight or Power)  Resolution (m)  Linear  0.5  Mean Error (m) 0.0014  Quintic  0.5  Natural Neighbour  Min Error (m)  Max Error (m)  Range (m)  RMSE (m)  Mean Absolute Error (m)  -3.26  4.86  8.12  0.18  0.11  0.0015  -4.49  4.99  9.48  0.17  0.10  0.5  0.00059  -3.27  4.70  7.96  0.17  0.11  Spline with Tension (2)  0.5  -0.00030  -4.80  4.76  9.55  0.17  0.10  Regularized Spline (0)  0.5  -0.00022  -5.67  5.04  10.72  0.18  0.10  IDW (3)  0.5  -0.0045  -4.50  6.34  10.85  0.20  0.11  ANUDEM  0.5  -0.0035  -4.23  3.91  8.14  0.18  0.11  Linear  1.0  -0.00088  -3.53  4.82  8.34  0.19  0.12  Quintic  1.0  -0.00022  -4.69  4.95  9.64  0.19  0.11  Natural Neighbour  1.0  -0.0015  -3.55  4.74  8.29  0.19  0.11  Spline with Tension (2)  1.0  0.00061  -4.80  4.70  9.50  0.19  0.11  Regularized Spline (0)  1.0  0.00075  -5.75  4.94  10.69  0.19  0.11  IDW (3)  1.0  -0.0041  -4.68  6.28  10.97  0.20  0.12  ANUDEM  1.0  -0.0035  -5.79  4.38  10.17  0.20  0.12  Linear  1.5  0.0016  -3.53  5.12  8.64  0.19  0.12  Quintic  1.5  0.0026  -4.69  5.21  9.90  0.19  0.12  Natural Neighbour  1.5  0.0010  -3.55  4.90  8.45  0.19  0.12  Spline with Tension (2)  1.5  0.0021  -4.87  4.23  9.10  0.21  0.13  Regularized Spline (0.1)  1.5  0.0051  -6.76  5.09  11.84  0.25  0.15  IDW (3)  1.5  -0.00375  -4.61  7.46  12.07  0.21  0.13  ANUDEM  1.5  -0.0037  -5.01  4.35  9.36  0.22  0.13  All DEMs contained visually identifiable interpolation artefacts (Figure 2.1). Linear, and to a lesser degree, natural neighbour interpolation contained faceted surfaces where ground return densities were low, a result of their derivation from a TIN. The IDW routine generated dimpled surfaces, but more importantly the DEMs contained abrupt  33  Figure 2.1 Hillshades derived from 1 m DEMs showing a meander of Lostshoe Creek. Note that all surfaces contain interpolation artefacts.  steps in areas where ground returns were sparse and heights changed rapidly. So too, regularized spline and spline with tension sometimes contained spurious values in areas of rapidly changing slope, especially when an excessive weighting parameter was chosen. The ANUDEM algorithm generated the most visually appealing surfaces, which were generally smooth and free of obvious artifacts. Though not shown in Figure 2.1, enabling  34  ANUDEM’s drainage enforcement option resulted in a large number of spurious channels being cut into the DEM surfaces.  Error Analysis Increasing slope from 0o to 40o caused decimetre-level increases in mean absolute errors for all interpolation algorithms (Figure 2.2). The finer spatial resolution DEMs were generally more accurate than the coarser resolution surfaces. The 0.5 m surface created using spline with tension (weight = 2) had the lowest mean absolute errors across the range of slope classes.  Figure 2.2 Effect of increasing slope on mean absolute error of interpolation. Vertical lines indicate standard errors of the means. Note that the data points have been offset horizontally in order to increase legibility.  35  Vegetation structural class, as derived from the TEM data, also affected interpolation accuracy over a range of decimetres (Figure 2.3). The highest mean absolute errors occurred in the young forest class, which is characterized by vigorous growth and the differentiation of the canopy into distinct layers, and were typically 0.18 m to 0.23 m. In contrast, mean absolute interpolation errors in the old forest category were slightly better, at 0.10 to 0.13 m, a result of the numerous gaps in the forest canopy allowing a larger proportion of lidar returns to reach the ground. The 1.5 m DEMs tended to be slightly less accurate than the 0.5 m surfaces in all vegetation structural classes.  Figure 2.3 Effects of vegetation structural class on mean absolute error of interpolation. Vertical lines indicate standard errors of the means. Note that the data points have been offset horizontally in order to increase legibility.  36  Correlation analysis between the different local point density surfaces and the absolute errors of the DEMs indicated that a 5x5 cell mean filter was appropriate for generalizing ground return densities in a given area. Figure 2.4 shows the effects of lidar ground return density on interpolation accuracy. The greatest mean absolute errors were less than 0.20 m and occurred at the lowest densities; however, interpolation accuracy stabilized as point densities increased beyond approximately 0.7 points/m2. The 0.5 m DEMs performed slightly better than the coarser resolution surfaces regardless of the interpolation routine used.  Figure 2.4 Effects of ground return density on mean absolute error of interpolation. Vertical lines indicate standard errors of the means.  Figure 2.5 shows the relationship between the Quickbird and Landsat7 ETM+ NDVI images, and lidar ground return density. As NDVI (an estimate of vegetation cover 37  metrics such as LAI) increased, lidar ground return density decreased. At the highest NDVI values for both image datasets, however, an unexpected increase in ground return density occurred. As noted previously, NDVI may saturate at higher ranges of LAI. Histograms of NDVI values for both images, however, showed no indications of oversaturation, such as a distribution with a sharp right hand shoulder (e.g. Huete et al., 2002), nor did either image contain values exceeding 0.85.  Figure 2.5 Normalized difference vegetation indices (NDVI) derived from Quickbird and Landsat7 ETM+ imagery, and lidar ground return density. An increase in ground return density occurs at the highest NDVI values. Vertical lines indicate standard errors of the means.  38  Error Prediction An error prediction map was created for the 0.5 m natural neighbour DEM using CART analysis (Figure 2.6, Figure 2.7). The CART analysis indicated that the most important predictor variables were point density and slope, although the Quickbird NDVI could provide an alternative to ground return density. Vegetation structural class mapped using the TEM data was the least useful predictor. An example of the rule-based classification is as follows (refer to Figure 2.7): if point density was greater than 0.175 points/m2, and terrain slope was less than 12.5o, prediction uncertainty was very low.  39  Figure 2.6 Quickbird image (left), and CART-derived prediction uncertainty map with underlying hillshade for the 0.5 m natural neighbour DEM (right). For this surface, slope and ground return density were the best predictors of interpolation error.  40  Figure 2.7 Simple classification tree resulting from CART analysis of absolute errors for the 0.5 m DEM created using natural neighbour interpolation.  Discussion DEM Validation Our research confirmed that the linear, natural neighbour, quintic, the two spline algorithms, IDW and ANUDEM were similar with respect to their global RMS and mean absolute errors (Table 2.2). The linear and natural neighbor interpolators are the most conservative of the algorithms and tended to have the lowest overall range of errors. The spline and quintic algorithms will exceed local minima and maxima where unconstrained by measured values, resulting in a larger overall range of prediction errors, though not necessarily an increase in their RMS or mean absolute errors. Though not reported here, accuracies for the splining and IDW routines were heavily influenced by their  41  parameterizations, with outliers often exceeding ± 6 m. The steps evident in the IDW surfaces (Figure 2.1) are artefacts of the interpolation algorithm, and would have significant impact on any geomorphic analysis.  Results indicated that when parameters are optimized for a given interpolation routine, selecting a spatial resolution may be an equally important choice for DEM creation. As Hengl (2006) points out, there is no best resolution for interpolating point data, and choices should be based on point density, distribution, horizontal accuracy, spatial autocorrelation, terrain complexity, or a combination thereof. Consideration should also be taken for computational expense and data storage.  Ultimately, natural neighbour interpolation was favoured for the overall parsimonious nature of its performance; that is, its ease of use, simple parameterization, generally smooth and visually appealing surfaces, consistent accuracy, and conservative predictions made it the algorithm of choice for interpolating the ground returns in our data. If a natural neighbour approach is not available and linear interpolation is judged to be too simplistic, ANUDEM, quintic or spline with tension interpolation may provide good alternatives. It is recommended, however, that a validation exercise similar to the one carried out here be undertaken in order to identify optimal parameterizations and spatial resolutions.  Error Analysis Analysis of the various independent variables, including slope, TEM vegetation structural class, and point density showed that all caused decimetre-magnitude increases in mean  42  absolute errors for all DEMs. Critically, all algorithms estimated heights more accurately when finer spatial resolutions were used, implying that the choice of spatial resolution may be more important than the choice of interpolation method, assuming that the algorithms are properly parameterized. Although the spline and IDW algorithms appear to be more sensitive to changes in slope than other methods, it is recognized that these routines do not involve the initial derivation of a TIN. As a result, comparing a slope layer derived using a TIN against all of the interpolation methods may favour the TINbased routines compared to the splining and IDW approaches, potentially resulting in increased apparent accuracies for the linear, quintic and natural neighbour routines.  The linear, natural neighbour, quintic spline with tension, and ANUDEM interpolation routines were similar in their accuracies, both in terms of their global statistics and when stratified by a given predictor variable such as slope or point density. The regularized spline algorithm generally created surfaces which were too flexible, resulting in largemagnitude outliers. Inverse distance weighted did not perform well compared to the other algorithms, and the surfaces were the least realistic representations of the terrain.  An interesting result was the increase in ground return density at the highest NDVI values for both the Quickbird- and Landsat7-derived data (Figure 2.5). This observation would seem counterintuitive, as one would expect an increase in vegetation cover and subsequently, an increase in NDVI, would be associated with a reduction in ground return density. We suggest that where vegetation densities are especially high, canopy closure may be so complete that it becomes an impenetrable barrier to the lidar sensor and causes pulse return patterns similar to that generated by the ground. This may result 43  in the misclassification of vegetation returns into the ground category by the automated classification routines typically used for labeling lidar returns as ground or non-ground.  Error Prediction Classification and regression tree analysis provided a simple means to identify areas where DEM surfaces may be of low quality. The analysis identified the most important predictors of interpolation error, specifically point density and slope for the natural neighbour algorithm, and provided a set of rules to classify the study area into categories of prediction uncertainty. The map allows users to identify problematic areas and determine if the DEM is suitable for their needs, if necessary by conducting more rigorous field investigations. For instance, there is a high degree of uncertainty along stream banks, implying that estimates of riparian vegetation heights may be of lower quality than might otherwise be expected from lidar data. As expected, areas with a combination of high slope and low point density were the most prone to interpolation error. The advantage of the method is highlighted in areas where the combined effects of the two predictor variables were less intuitive, for instance, where one variable was high and the other low, or both were moderate. Ultimately, using prediction uncertainty maps may not only provide insights into the selection of an optimal DEM for a given combination of vegetation and terrain characteristics, but also into potential problems with lidar-derived vegetation height estimates in those areas where the DEM surface is suspect.  44  References Abramov, O., McEwan, A., 2004. An evaluation of interpolation methods for Mars Orbiter Laser Altimeter (MOLA) data. International Journal of Remote Sensing 25(3), 669-676. Adams, J.C., Chandler, J.H., 2002. Evaluation of lidar and medium scale photogrammetry for detection soft-cliff coastal change. Photogrammetric Record 17(99), 405-418. Breiman, L., Friedman, J., Olshen, R., Stone, C., 1984: Classification and regression trees. Wadsworth, Pacific Grove, Ca, 358 pp. Brown de Colstoun, E. C., Story, M. H., Thompson, C., Commisso, K., Smith, T. G., Irons, J. R., 2003. National Park vegetation mapping using multitemporal Landsat 7 data and a decision tree classifier. Remote Sensing of Environment 85, 316-327. Coops, N.C., Wulder, M.A., Culvenor, D.S., St-Onge, B., 2004. Comparison of forest attributes extracted from fine spatial resolution multispectral and lidar data. Canadian Journal of Remote Sensing 30(6), 855-866. Coops, N.C., Wulder, M.A., White, J.C., 2006. Integrating remotely sensed and ancillary data sources to characterize a mountain pine beetle infestation. Remote Sensing of Environment 105, 83-97. Demarchi, D. A., Marsh, R. D., Harcombe, A. P., Lea, E. C., 1990. The environment (of British Columbia), In: Campbell, R.W., Dawe, N.K., McTaggart-Cowan, I., Cooper, J.M., Kaiser, G.W., McNall, M.C.E. (Eds.) The Birds of British Columbia, Royal British Columbia Museum, Environment Canada, Canadian Wildlife Service, Victoria, pp. 55-142.  45  Dubayah, R.O., Drake, J.B., 2000. Lidar remote sensing for forestry. Journal of Forestry 98, 44-46. ESRI, 2005. ArcDoc Version 9.1. Environmental Systems Research Institute, Redlands, USA. Hengle, T., 2006. Finding the right pixel size. Computers & Geosciences. 32, 1283-1298. Hengl, T., Heuvelink, G.B.M., Stein, A., 2004. A generic framework for spatial prediction of soil variables based on regression-kriging. Geoderma 120, 75-93. Hodgson, M.E., Bresnahan, P., 2004. Accuracy of airborne lidar-derived elevation: empirical assessment and error budget. Photogrammetric Engineering and Remote Sensing 70(3), 331-339. Hodgson, M.E., Jensen, J., Raber, G., Tullis, J., Davis, B.A., Thompson, G., Schuckman, K., 2005. An evaluation of lidar-derived elevation and terrain slope in leaf-off conditions. Photogrammetric Engineering & Remote Sensing 71(7), 817-823. Huete, A.R., Didan, K., Miura, T., Rodriguez, E.P., Gao, X., and Ferreira, G., 2002. Overview of the radiometric and biophysical performance of the MODIS vegetation indices. Remote Sensing of Environment 83, 195-213. Hutchinson, M.F., 1989. A new procedure for gridding elevation and stream line data with automatic removal of spurious pits. Journal of Hydrology 106, 211-232. Hutchinson, M.F., 1996. A locally adaptive approach to the interpolation of digital elevation models. In: Proceedings, Third International Conference/Workshop on Integrating GIS and Environmental Modeling, Santa Fe, NM, National Center for Geographic Information and Analysis,  46  http://www.ncgia.ucsb.edu/conf/SANTA_FE_CDROM/sf_papers/hutchinson_michael_dem/local.html Jensen, J.R., 2007. Remote Sensing of the Environment: an Earth Resource Perspective (2nd Edition). Upper Saddle River, Pearson Prentice Hall. Johnston, K., Ver Hoef, J. M., Krivoruchko, K., Lucas, N., 2001. Using ArcGIS geostatistical analyst. Environmental Systems Research Institute, Redlands, USA, 300 pp. Lawrence, R.L., Wright, A., 2001. Rule-based classification systems using classification and regression tree (CART) analysis. Photogrammetric Engineering and Remote Sensing 67(10), 1137-1142. Lefsky, M.A., Cohen, W.B., Acker, S.A., Parker, G.G., Spies, T.A., Harding, D., 1999. Lidar remote sensing of the canopy structure and biophysical properties of Douglas-fir western hemlock forests. Remote Sensing of Environment 70, 339361. Liang, S., 2004. Quantitative Remote Sensing of Land Surfaces. Hoboken, John Wiley & Sons. Lim, K., Treitz, P., Wulder, M., St-Onge, B., Flood, M., 2003. Lidar remote sensing of forest structure. Progress in Physical Geography 27(1), 88-106. Lloyd, C.D., Atkinson, P.M., 2002. Deriving DSMs from lidar data with kriging. International Journal of Remote Sensing 23, 2519-2524. Maune, D.F., Kopp, S.M., Crawford, C.A., Zervas, C.E., 2001. Introduction, In: Maune, D.F. (Ed.) Digital Elevation Model Technologies and Applications: the DEM  47  Users Manual, American Society for Photogrammetry and Remote Sensing, Bethesda, pp. 1-34. Meidinger, D., Pojar, J., 1991. Ecosystems of British Columbia. BC Special Report Series No. 6. British Columbia Ministry of Forests, Victoria, Canada, 330 pp. Melendez, K.V., Jones, D.L., Feng, A.S., 2006. Classification of communication signals of the little brown bat. Journal of the Acoustical Society of America 120, 10951102. Mitášová, H., Hofierka, J., 1993. Interpolation by regularized spline with tension: II. Application to terrain modeling and surface geometry analysis. Mathematical Geology 25(6), 657-669. Mitchell, W.R., Green, R.N., Hope, F.D., Klinka, K., 1989. Methods for biogeoclimatic ecosystem mapping. Research Report RR89002-KL. British Columbia Ministry of Forests, Victoria, Canada, 33 pp. Næsset, E., 1997. Determination of mean tree height of forest stands using airborne laser scanner data. Journal of Photogrammetry and Remote Sensing 52, 49-56. RSI, 2006. ENVI Version 4.3. Research Systems Inc., Boulder, USA. Sambridge, M., Braun, J., McQueen, H., 1995. Geophysical parameterization and interpolation of irregular data using natural neighbours. Geophysical Journal International 122, 837-857. Schwalm, C.R., Black, T.A., Amiro, B.D., Arain, M.A., Barr, A.G., Bourque, C., Dunn, A.L., Flanagan, L.B., Giasson, M.A. Lafleur, P., Margolis, H.A., McCaughey, J.H., Orchansky, A.L., Wofsy, S., 2006. Photosynthetic light use efficiency of  48  three biomes across an east-west continental-scale transect in Canada. Agricultural and Forest Meteorology 140, 269-286 Sibson, R., 1981. A brief description of natural neighbor interpolation, In: Barnet, V. (Ed.) Interpreting Multivariate Data, John Wiley & Sons, New York, 21-36. Steinberg, D., Colla, P., 1995. CART: tree-structured non-parametric data analysis, Salford Systems, San Diego, CA, 336 pp. Su, J., Bork, E., 2006. Influence of vegetation, slope and lidar sampling angle on DEM accuracy. Photogrammetric Engineering and Remote Sensing 72(11), 1265-1274. Tucker, C.J., 1979. Red and photographic infrared linear combinations for monitoring vegetation. Remote Sensing of Environment 8, 127-150. Yanalak, M., 2003. Effect of gridding method on digital terrain model profile data based on scattered data. Journal of Computing in Civil Engineering 17(1), 58-67.  49  3 TOWARDS THE ESTIMATION OF TREE STRUCTURAL CLASS IN NORTHWEST COASTAL FORESTS USING LIDAR REMOTE SENSING5 Introduction Globally, forests are renewable resources which must be managed in a sustainable fashion whilst balancing economic, cultural, and social needs. The Canadian Province of British Columbia, which contains approximately half of the country’s softwood lumber inventory, and where the forestry industry is responsible for 45% of the province’s manufacturing shipments (BC Stats 2005), is no exception. While forestry’s economic benefits are significant, the extraction of forest resources must be performed in a sustainable manner.  In response to this need, the Province of British Columbia has developed a suite of resource values to monitor key attributes of forest health and sustainability, including biodiversity, timber, and soil, amongst others. Each resource value is assessed by monitoring a suite of indicators. To be useful, an indicator must be easy to interpret, correlated with changes in ecosystem processes, show low temporal and spatial variability, be applicable across different regions, and be cost effective (Breckenridge et al.,1995). Common examples include tree height, diameter at breast height (DBH), species richness, and wildlife tree (WT) class (or decay class), which are traditionally assessed using field-based approaches in association with aerial photography. Field assessments for both inventories and monitoring, however, can be expensive, labour  5  A version of this chapter has been submitted for publication. Bater, C.W., Coops, N.C., Gergel, S.E., and Collins, D. (In Review) Towards the estimation of tree structural class in northwest coastal forests using lidar remote sensing. Remote Sensing of Environment.  50  intensive, provide small sample sizes and intensity, and often cover only limited geographic areas, while aerial photography interpretation suffers from time and cost issues, is prone to operator bias and subjectivity, and is limited by a shortage of trained interpreters. As a result, there has been increased interest in augmenting provincial-level ecosystem and timber inventory mapping initiatives with digital remote sensing technologies, including research into light detection and ranging (lidar). Currently, acquisition, processing, and storage costs are high for lidar data, but its use in operational forest and resource management is becoming more common.  Lidar is an active remote sensing technology employing an airborne laser capable of simultaneously mapping terrain and vegetation heights with sub-metre accuracy. Various measures of forest structure and biodiversity have previously been estimated within the context of coastal northwest forests using laser altimetry (e.g. Lefsky et al., 1999; Hudack et al., 2002; Anderson et al., 2005; Lefsky et al., 2005a; Lefsky et al., 2005b; Coops et al., 2007). Table 3.1 presents examples of previous research where lidar data have been used to estimate ecological variables in northwest coastal forests. While tree and stand heights are obvious candidates for measurement using lidar, other variables such as DBH, basal area, leaf area index (LAI), crown morphology, and stem densities have also been successfully assessed using the technology. Furthermore, a wide range of methodologies have been used to process lidar vegetation height information, including calculating summary statistics and height percentiles, analyzing canopy height models, and extracting distribution parameters (e.g. Weibull shape and scale parameters) (Table 3.1).  51  Table 3.1 Examples of previous research where ecological variables were estimated using lidar remote sensing in northwest coastal forests. Reference  Type of Lidar Sensor  Species  Method  Variables  Magnussen and Boudewn (1998)  Small footprint, discrete return  Douglas-fir (Pseudotsuga menziesii)  Percentiles  Stand heights, canopy leaf area  Lefsky et al., 1999  Large footprint, full waveform  Douglas-fir (Pseudotsuga menziesii), western hemlock (Tsuga heterophylla)  Volume estimates  Biomass and leaf area index  Anderson et al., 2005  Small footprint, discrete return  Douglas-fir (Pseudotsuga menziesii), western hemlock (Tsuga heterophylla), red alder (Alnus rubra), maple (Acer spp.)  Percentiles, summary statistics  Crown fuel weight, crown bulk density, canopy base height, canopy height  Lefsky et al., 2005a  Large footprint, full waveform  Douglas-fir (Pseudotsuga menziesii), western hemlock (Tsuga heterophylla), Sitka spruce (Picea sitchensi), true firs (Abies sp.), ponderosa pine (Pinus ponderosa)  Canopy height profiles, volume estimates  Aboveground biomass; maximum, standard deviation, and mean DBH; basal area; mean, maximum and Lorey’s height; leaf area index; density of stems; basal area; fractional foliage cover  Coops et al., 2007  Small footprint, discrete return  Douglas-fir (Pseudotsuga menziesii), western hemlock (Tsuga heterophylla)  Distribution fitting, volume estimates  Mean height, stand volume, crown volume, leaf area index, basal area  One important ecological variable that has not been examined with comparable detail is the decay class or structural life stage of the tree. Decay class quantifies an individual tree, ranging from healthy trees, to live trees losing foliage, to standing dead stems (snags), to broken stems in various stages of decay. Within British Columbia, the form classification is known as the WT class, which when accumulated over a stand provides an indication of the amount of living and dead trees, and their state of decay. Seielstad and Queen (2003) previously discussed the ability of lidar to characterise fuel bed  52  roughness in forests in the western United States, and noted that the direct estimation of coarse woody debris loads may be achievable using the technology. To our knowledge, however, no attempt has been made to use lidar remote sensing to estimate the distribution of wildlife tree classes, including both living and dead stems and their state of decay, within forest stands.  The amount and variability of dead wood is an important indicator of forest structural diversity and biodiversity (Morrison and Raphael, 1993; Noss, 1999; McElhinny et al., 2005). Standing dead trees, or snags, are a critical component of coastal forests, increasing vertical and horizontal structural heterogeneity, and providing nesting platforms and denning sites for forest biota. Snags are often the centre of small openings in the canopy, allowing sunlight to reach understorey vegetation. Thus, forest managers in British Columbia have identified WT class as an important indicator of stand-level biodiversity and old growth forest health (e.g. Clayoquot Sound Scientific Panel, 1995; Delong et al. 2004; Province of British Columbia, 2005). Conventional assessments of dead trees are typically conducted in the field, where a WT class is assigned to each stem in a plot. Their numbers are then typically summarized as the number of dead trees (WT class 3+) per hectare (e.g. Province of British Columbia, 2005).  The goal of this research was to estimate the distributions of living and dead tree forms within forests by developing statistical relationships between plot-level distributions of WT class and commonly derived lidar vegetation metrics.  53  Methods Study Area The study area is a 4 x 7 km area in Clayoquot Sound on Vancouver Island, British Columbia (49o 0’ 35” N, 125o 37’ 21” W), approximately 15 km south of the town of Tofino. The area is classified as Coastal Western Hemlock (CWH) zone, based on the biogeoclimatic ecosystem classification (BEC) system (Meidinger and Pojar, 1991). Although flanked by the Vancouver Island Range, the topography is subdued and dominated by Pleistocene glacial deposits. Streams are generally meandering and alluvial, but most have short sections that show evidence of constraint. From 1971 to 2000, Tofino received an average of 3,306 mm of precipitation each year, while mean daily minimum, average and maximum temperatures were 5.4oC, 9.1oC, and 12.8oC, respectively (Environment Canada, 20066).  Clayoquot Sound is a multi-use area and includes both recently harvested Crown land, and mature first and second growth forest in Pacific Rim National Park. Western hemlock (Tsuga heterophylla) is a dominant or codominant tree species throughout; western redcedar (Thuja plicata), amabilis fir (Abies amabilis), yellow-cedar (Chamaecyparis nootkatensis), Sitka spruce (Picea sitchensis), Douglas-fir (Pseudotsuga menziesii var. mensiesii), and red alder (Alnus rubra) also occur under differing conditions.  6  Environment Canada, 2006. Canadian Climate Normals 1971-2000, Tofino A, British Columbia. http://www.climate.weatheroffice.ec.gc.ca/climate_normals/index_e.html  54  Data Sets Terrestrial Ecosystem Mapping for Stratification The Clayoquot Sound area has been previously mapped using the Province of British Columbia’s Terrestrial Ecosystem Mapping (TEM) classification system, which is derived from 1:20,000 to 1:50,000 scale aerial photography (Mitchell et al., 1989; Demarchi et al., 1990). A hierarchical classification system, TEM mapping stratifies the landscape into map units based on climate, physiography, surficial material, geology, soil and vegetation. Based on the TEM classification system, the study area encompasses a range of vegetation structural stages, including herb (1% of total area) shrub/herb (14%), pole/sapling (32%), young forest (4%), and old forest (46%).  Field Data Field data were collected in 2005 and 2006 from 22 forest plots ranging from pole/sapling to old forest based on the TEM classification (Table 3.2). Five of the old forest plots were located in variable retention harvest blocks. Data were collected from 625 m2 or greater rectangular plots, with plot centres and corners mapped at a horizontal accuracy of approximately 1-5 m using a post-processed differentially corrected GPS (Trimble GeoXT). For each stem with a DBH > 10 cm, distance and bearing from plot centre, tree height, DBH, and species were recorded, with crown dimensions measured for every fifth tree. For conifers, the WT class was estimated visually using a field sheet showing growth and decay stages ranked 1 through 9: class 1 are healthy trees with all foliage intact; class 2 are living trees exhibiting some defoliation; 3-5 are dead trees with hard wood, often with broken tops; 6 represents dead trees with broken tops and spongy  55  wood; 7 and 8 are dead trees with broken tops and soft wood; and class 9 represents dead and fallen trees (Figure 3.1).  Table 3.2 Summary statistics for sample plots by age class for stems with a DBH > 10 cm. One outlier was excluded from this summary and all subsequent analyses. Variable  Pole/Sapling n=6 (mean/range)  Young Forest n=3 (mean/range)  Old Forest n = 12 (mean/range)  Stems ha-1  1491 / 1544  1147 / 816  957 / 1391  Basal Area (m2 ha-1)  144.9 / 127.3  84.1 / 36.8  142.3 / 372.6  Mean Height (m)  19.3 / 5.3  18.3 / 3.9  12.6 / 12.6  Standard Deviation of Height (m)  6.1 / 2.0  5.1 / 1.3  6.33 / 12.1  Maximum Height (m)  32.5 / 18.0  25.8 / 4.7  27.0 / 30.4  Mean DBH (cm)  27.8 / 12.8  25.6 / 5.5  31.3 / 37.2  Maximum DBH (cm)  107.6 / 106.8  125.7 / 98.9  170.5 / 343.2  Standard Deviation of DBH (cm)  17.2 / 13.7  15.8 / 5.4  29.4 / 63.4  Dead Trees (WT Class 3+) (%)  12.0 / 18.1  13.1 / 9.0  19.6 / 12.1  In order to characterize the structural characteristics of the plots, vertical foliage profiles were measured using a modified point quadrat camera method. First described by Warren Wilson (1960, 1963), the point quadrat method originally relied on a series of vertical transects defined by plumb lines passing through the canopy. By recording the position of  56  leaves touching the plumb lines and then applying equations, LAI and vertical foliage profiles could be computed. MacArthur and Horn (1969) modified the technique by using a 35 mm camera equipped with a telephoto lens calibrated for range measurement, and a grid superimposed on the focusing screen. By assuming a random horizontal distribution of all foliage, data on the lowest leaf heights could be transformed into estimates of LAI for the canopy (MacArthur and Horn, 1969). Aber (1979a, 1979b) later applied the method to a range of northern temperate forests.  Figure 3.1 Typical field form used to estimate the wildlife tree class of conifers.  At nine points within a given plot (plot centre, and the cardinal and off-cardinal directions), a 35 mm camera was set up on a tripod and levelled. The ranges to the lowest canopy element (including stems, branches and leaves) within each cell of the 3 x 5 grid were then measured by focusing on the target, and then reading the range off of the telephoto lens’ focusing ring. The resulting 135 estimates of foliage height were  57  transformed into vertical foliage height profiles using equations described by Aber (1979b).  Lidar Data Small footprint laser data were collected during July 2005 by Terra Remote Sensing (Sidney, British Columbia), using a TRSI Mark II two-return sensor onboard a fixedwing platform. Flying at a mean height of 800 m above ground level, the survey was optimized to achieve a nominal point spacing of one laser pulse return every 1.5 m2 (Table 3.3). Ground and non-ground returns were separated using Terrascan v 4.006 (Terrasolid, Helsinki, Finland).  Table 3.3 Lidar sensor and survey parameters Sensor and Survey Parameters  Value  Sensor Type  TRSI Mark II discrete return sensor  Number of Returns  Two, first and last  Beam Divergence Angle (mrad)  0.5  Wavelength (nm)  1064  Mean Flying Height Above Ground (m)  800  Pulse Frequency (kHz)  50  Mirror Scan Rate (Hz)  30  Scan Angle (degrees)  ±23  Mean Footprint Diameter (m)  0.4  58  Lidar Data Processing A digital elevation model (DEM) was generated using methods described in Bater and Coops (2007). A number of different algorithms exist to interpolate lidar ground hits into a terrain surface, thus a validation exercise was undertaken to identify the optimal combination of interpolation routine and spatial resolution. First, the lidar data were randomly subsetted into a prediction dataset and a validation dataset. A suite of DEMs were then generated using linear, quintic, natural neighbour, regularized spline, spline with tension, a finite difference approach (ANUDEM), and inverse distance weighted interpolation routines, at spatial resolutions of 0.5, 1.0 and 1.5 m. Based on comparisons between the surfaces and the validation data, a 0.5 m spatial resolution digital elevation model was generated by applying a natural neighbour (Sibson, 1981; Sambridge et al., 1995) interpolation algorithm to the ground returns.  The heights of the vegetation returns above the ground were computed by subtracting the DEM heights from the vegetation return heights. A large number of variables were extracted from the lidar vegetation data based on Gobakken and Næsset (2005) and Næsset (2002; 2004). These variables attempt to capture vertical structure by classifying hits into percentiles based on their height distribution through the forest canopy, and included the 5, 10, 15… 95 percentiles, in addition to the means, maximums, standard deviations, and coefficients of variation of vegetation return heights within each plot. The natural logarithms of the cases of each variable were also computed.  59  Logistic Regression for Estimating Wildlife Tree Class Proportions The WT classes, ranging from 1 to 9, can be considered ordinal or ranked data, as the first two classes represent living trees, and the following seven represent dead trees with increasing amounts of decay. Thus, ordinal logistic regression was employed to estimate the cumulative proportions of WT classes present in plots using lidar-derived predictor variables. Ordinal logistic regression estimates probabilities for ranked data with more than two outcomes, based on one or more categorical and/or continuous predictors. The technique employs a cumulative odds model, which predicts the likelihood of a dependent variable being at or below a particular ranked category (Hosmer and Lemeshow, 2000; O’Connell, 2006). In the context of this research, the logistic model was used to predict the cumulative proportions of WT class 1, WT classes 1 and 2 combined, WT classes 1, 2 and 3 combined, and so on to WT classes 1– 9 combined, for each forest plot. The probabilities were calculated using the following equation:  P(WTclass ) =  1 1+ e  −( b0 +b1 X1 )  (1)  Where P(WTclass) = the probability of a wildlife tree class occurring. b0 and b1 = coefficients of a linear equation X1 = the lidar-derived predictor variable e = the base of the natural logarithm.  As the probabilities range from 0 to 1, the cumulative proportion of the classes could be calculated by simply multiplying P(WTclass) by 100.  60  Prior to building the logistic model, a variety of lidar-derived metrics were assessed using best subset techniques, which estimate likelihood scores for candidate predictor variables. Likelihood scores evaluate the statistical significance of parameter estimates using maximum likelihood methods (Hill and Lewicki, 2006). Once the best lidar-derived predictor variable was identified, ordinal logistic regression was used to predict the cumulative proportions of WT classes within the forest plots.  Models were assessed first by evaluating several diagnostics, and secondly by examining residuals between observed and predicted WT class cumulative proportions. First, a null logistic model was fitted to the WT class data using no predictor variable. This null model provided a baseline with which to evaluate the contribution of predictor variables by comparing the log-likelihoods (Tabachnick and Fidell, 2001) of models with predictor variables against the null case. A likelihood ratio test was used to assess the improvement of a model after a predictor has been added using the following equation:  X 2 = 2( LL New − LL Null )  (2)  Where χ2 follows a chi-square distribution, with degrees of freedom equal to the number of parameters in the new model minus the number in the null model. LLNew = the log-likelihood of a model incorporating predictor variables. LLNull = the log-likelihood of a model with no predictor variables.  61  Cox and Snell R2 values were obtained to estimate the amount of variance explained by the lidar-derived predictor variables, and were calculated using the following equation: 2 CS  R  = 1− e  ⎡ 2 ⎤ ⎢ − n ( LLNew − LLNull ) ⎥ ⎣ ⎦  (3)  2 Where RCS = the Cox and Snell R2.  n = the sample size.  Cox and Snell R2 values, however, are not the same as the coefficient of variation used in linear regression as they do not assess goodness of fit; rather, they are more suited to comparing the quality of competing models (Hosmer and Lemeshow, 2000). Thus, in addition to the likelihood scores described above, Cox and Snell R2 values were calculated for each lidar-derived predictor variable in order to compare their ability to estimate WT class.  In addition, the Wald statistic was calculated to test the significance of the coefficients of the predictor variables estimated by the logistic model, and test the null hypothesis that they differed significantly from 0. The Wald statistics were calculated using the following equation: Wald =  b SEb  (4)  Where SEb = the standard error of the regression coefficients.  62  Finally, observed and predicted WT class proportions were plotted, and the Pearson correlation coefficients (r) and root mean square errors (RMSE) of the residuals were examined.  Predicting Wildlife Tree Class Across the Landscape In order to demonstrate how results may be upscaled from the plot level to the landscape  scale, WT class was mapped using lidar vegetation data and the logistic regression models. The cumulative proportions of WT class 1, WT class 1 and 2 combined, WT class 1, 2 and 3 combined, and so on within a given area were estimated across the landscape using the lidar-derived predictor variables. A map of the relevant lidar-derived predictor variable was created at a spatial resolution of 25 m in order to approximate the scale of the forest plot data. The coefficients calculated by the logistic regression model were then applied to the lidar-derived surface to estimate the WT class proportions for each 25 m x 25 m pixel using Equation 1.  Results Initial examination of the field data indicated that one plot was an outlier, as it was located in a stand which had experienced significant disturbance, possibly from insect infestation, resulting in a stand structure not replicated in the dataset. This plot was considered an outlier and was excluded from further analyses. Figure 3.2 shows the WT class distributions for three plots, specifically one pole/sapling plot, one young forest plot, and one old forest plot.  63  Vertical Foliage Profiles Figure 3.3 shows the vertical foliage profiles for pole/sapling and old forest plots  measured using the modified point quadrat camera method. It is apparent that in the younger plots, which contain relatively low percentages of dead trees (Table 3.2 and Figure 3.2), the majority of the vegetation is concentrated in the upper canopy, where densities can be so high that the canopies are almost entirely closed. The dense overstories limit the amount of sunlight available to vegetation at lower levels, resulting in relatively low amounts of subdominant vegetation and open understories. The older plots contain a much more even distribution of vegetation down through their vertical profiles, and exhibit more open canopies. The open architecture of these forests allows for increased sunlight penetration, resulting in multiple canopies and dense understories.  64  Figure 3.2 Histogram displaying the proportions of WT classes within three forest plots. Note the relative increase in the proportions of the higher WT classes as seral stage increases from pole/sapling to old forest.  65  Figure 3.3 Typical vertical foliage profiles of pole/sapling (n=3) and old forest (n = 3) plots. Foliage is concentrated in the upper canopy of the pole sapling plots, whilst the more structurally complex old growth plots have vegetation distributed evenly through the profile. The higher percentage of dead trees within the old forest plots results in lower canopy densities.  Logistic Regression for Estimating Wildlife Tree Class Proportions The most significant lidar-derived predictor of the WT class proportions was the  coefficient of variation of the lidar vegetation return heights (log-transformed) (Table 3.4, Table 3.5), which had a likelihood score of 90.2 and a Cox and Snell R2 of 0.044. The full model chi-square statistic (χ2 = 94.51, p <0.000) indicated that the model was a better fit  66  than the null model. Table 3.5 displays model coefficients, and the Wald statistic indicates that the regression coefficient for the coefficient of variation of the lidar vegetation return heights (CVLidar) differs significantly from 0 (Wald = 3.981, p = 0.046).  Table 3.4 Fit statistics for the ordinal logistic regression model using the lidarderived variable (log transformed coefficient of variation of the vegetation return heights) as a predictor of wildlife tree class cumulative proportions.  n  2100  LLNull  -1919.77  LLNew  Chi-square (p)  Cox and Snell R2  -1872.52  94.51 (0.000)  0.044  Table 3.5 Parameter statistics for the ordinal logistic regression model using the lidar-derived variable (log transformed coefficient of variation of the vegetation return heights) as a predictor of wildlife tree class cumulative proportions. The intercept estimates are the coefficients used to predict the cumulative proportions of stems within the wildlife tree classes using Equation 1. Parameter  Estimate  Standard Error  Wald  p  Intercept 1  0.673  0.380  3.132  0.077  Intercept 2  1.050  0.396  7.040  0.008  Intercept 3  1.230  0.406  9.155  0.002  Intercept 4  1.480  0.426  12.095  0.001  Intercept 5  1.787  0.456  15.328  0.000  Intercept 6  2.039  0.489  17.405  0.000  Intercept 7  2.617  0.591  19.625  0.000  Intercept 8  6.541  3.611  3.282  0.070  CVLidar  -0.899  0.450  3.981  0.046  67  Examples of the form of the final logistic model predicting the probability that a tree within a forest plot will be WT class 1 (Equation 5), and the probability that tree will be WT classes 1 or 2 (Equation 6), when the logarithm of CVLidar =-0.34, are as follows: P(WTclass = 1) =  1 1+ e  P(WTclass = 1,2) =  − ( 0.673 + ( −0.899 * −0.34 ))  = 0.73  1 1+ e  − (1.050 + ( −0.899* −0.34 ))  = 0.80  (5)  (6)  By multiplying the results by 100, the probabilities can be converted to proportions.  Figure 3.4 shows predicted vs. observed values of the cumulative proportions of stems within the range of WT classes. The cumulative proportions of WT class predicted by the final logistic model were strongly correlated with those observed in the field (r = 0.85, p<0.000, RMSE = 4.9 %). Figure 3.5 shows how the predicted proportions of WT class are related to changes in the coefficient of variation of the lidar vegetation returns. As the lidar-derived variable increases, the predicted proportion of healthy trees (i.e. class 1) decreases, whilst the proportions of the other classes increase. The CVLidar increased as canopies became more structurally complex and open, which is at least partly the result of the presence of standing dead trees.  68  Figure 3.4 Predicted vs. observed values for cumulative proportions of stems within each wildlife tree class when using the logarithm of CVLidar as a predictor variable (r = 0.85, p <0.000, RMSE = 4.9%).  69  Figure 3.5 Relationship between CVLidar and the predicted individual proportions of stems in each wildlife tree class. As the lidar-derived variable increases, the predicted proportion of healthy trees (i.e. class 1) decreases, whilst the proportions of trees within the other classes increase.  As forest stands increase in structural complexity, and the percentage of dead trees and the number of canopy gaps increase, lidar pulse returns penetrate deeper through the forest canopy. Figure 3.6 demonstrates that the lower and median height percentiles were relatively strong predictors of WT class, with likelihood scores greater than 70, and Cox and Snell R2 values greater than 0.030. Figure 3.7 shows how the height of the 50th percentile decreases as the percentage of dead trees increases within the TEM structural classes. The strength of the relationship between the lidar height percentiles and WT class diminishes at the highest levels of the vertical distribution of vegetation returns, as maximum heights were poorly correlated with seral stage and canopy complexity. 70  Figure 3.6 The lidar-derived height percentiles plotted against the best subsets likelihood scores and the Cox and Snell R2 values. It is the lowest and median height  90  0.045  80  0.040  70  0.035  60  0.030  50  0.025  40  0.020  30  0.015  20 10  Cox and Snell R-Square  Likelihood Score  percentiles that account for most of the variance in wildlife tree class.  0.010 Likelihood Score Cox and Snell R-Square  0  0.005 0.000  5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 Lidar Height Percentile  The trend of the mean vegetation return height varied closely with that of the 50th percentile (r = 0.98), but virtually no correlation existed between the mean and the standard deviation of the vegetation return heights (r = 0.01). The lack of correlation between the mean and standard deviation of lidar vegetation return heights meant that the coefficient of variation was sensitive to changes in vertical canopy structure, and increased from approximately 0.15 to 0.3 for pole sapling and young forest, to 0.4-1.2 for old forest.  71  Figure 3.7 Means and ranges of (1) lidar-derived heights of the 50th percentile, and (2) the percentage of dead trees, grouped by terrestrial ecosystem mapping structural class.  Predicting WT Class Across the Landscape  Once the relationships between CVLidar and the field data were established, the WT class proportions could be predicted across the landscape. Figure 3.8 shows a Quickbird image of part of the study area, the lidar-derived estimates of the coefficient of variation (logtransformed), and the proportion of dead trees (WT class 3+) predicted to be within a given 25 m by 25 m pixel.  72  Figure 3.8 Map displaying (A) a true colour Quickbird image of the study area (0.80 m spatial resolution); (B) the coefficient of variation (log-transformed) derived from the lidar vegetation return heights; and (C) the predicted proportions of trees that are dead within a given pixel (e.g. WT class 3+). Maps B and C have spatial resolutions of 25 m.  73  Discussion The distribution of wildlife tree classes within a plot is important for quantifying the current structure of a forest stand, as well as for managing a stand for wildlife and biodiversity values. Because the proportion of WT classes from 1 to 9 within a plot was highly variable, examining the cumulative proportions of the observed frequencies was a convenient way to summarize the data for analysis. Once an understanding of how the proportions of wildlife tree classes vary over the landscape as a function of stand form, lidar data is employed to extrapolate over large areas.  The results presented here indicate the capacity of lidar to estimate the cumulative proportions of wildlife trees within forest plots. The lidar-derived coefficient of variation of the vegetation return heights (log transformed) was the most significant predictor variable. Given that CVLidar has been found to be a good descriptor of canopy structural attributes (e.g. Næsset, 2002; Næsset and Økland 2002: Anderson et al., 2005), this is not a surprising finding. This study also demonstrates that the lowest and median height percentiles were also good candidates for predicting WT class. This is likely the result of the direct linkage, noted by the Clayoquot Sound Scientific Panel (1995) and Clark et al. (2004), between tree mortality and overall stand structure.  Late seral stage forests in the study area are characterized by heterogeneous canopies interspersed with gaps containing dead trees and downed stems, and young trees regenerating through gap-phase replacement (Clayoquot Sound Scientific Panel, 1995). Many stands in the study area, however, are regenerating following harvest, and do not conform to the multiple-cohort distribution typical of old forests in Clayoquot Sound.  74  Instead, these regenerating stands are dominated by a single cohort, or age class, contain a lower proportion of dead trees, and tend to have dense canopies with foliage concentrated in their upper layers. Thus, stands exhibit the range of stand development described by Oliver and Larson (1996), including stand initiation, stem exclusion, understory reinitiation and old forest stages. The result is a varied mosaic of vegetation structural types and habitats across the landscape.  The vertical distribution and density of foliage has a direct impact on the vertical distribution of lidar returns within the plots. Gaps containing snags (defoliated, often limbless stems with very different structures than live trees), increased the mean penetration depth of lidar returns into the forest canopy. Conversely, the denser canopies of younger, regenerating stands tended to occlude the ground from the lidar sensor, resulting in very different vertical distributions of vegetation returns. Figure 3.9 shows the raw lidar pulse returns within a pole-sapling and an old forest plot. The lack of penetration to the lower canopy in the younger plot results in greater heights of the mean and 20th percentile, and a smaller standard deviation of heights, even though the old forest stand is 10 m taller. Critically, non-ground returns were not removed below a given height threshold, and though many may have actually intercepted the understorey, coarse woody debris, boulders, or the ground, their inclusion was nonetheless an important contribution to the analyses.  75  Figure 3.9 Lidar ground and vegetation returns in pole/sapling (n = 1,794) and old forest (n = 1,969) plots. Dense canopies typical of the pole/sapling structural class occlude the lower canopy and ground, and as a result the majority of lidar returns are clustered in the upper parts of the canopy. Gaps that characterize old forest canopies, however, allow more pulse returns to penetrate to lower parts of the canopy, resulting in lower mean heights and higher variances.  76  Increasing the number of plots across the full range of stand ages or structural class distributions is a necessary next step to both adequately capture the heterogeneity within and between the structural classes (especially old forests) found in the study area. Additional field data will also enable the application of multivariate statistical techniques, where more than a single predictor variable can be employed. Furthermore, additional research is required to determine if these techniques can be extended to different forest environments. We believe that lidar metrics such as the coefficient of variation and height percentiles, among others, can be useful proxies to ecologists and forest managers interested in augmenting their current mapping initiatives using lidar remote sensing.  77  References Aber, J.D., 1979a. Foliage-height profiles and succession in northern hardwood forests. Ecology 60, 18–23. Aber J.D., 1979b. A method for measuring foliage-height profiles in broad-leaved forests. The Journal of Ecology 67, 35–40. Anderson, H.A., McGaughey, R.J., and Reutebuch, S.E., 2005. Estimating forest canopy fuel parameters using LIDAR data. Remote Sensing of Environment, 94, 441-449. BC Stats, 2005. Fact sheet released by the British Columbia Ministry of Labour and Citizen’s Services, “Quick facts about British Columbia,” Victoria, Canada. http://www.bcstats.gov.bc.ca/data/qf.pdf (accessed 12 October 2006). Breckenridge, R.P., Kepner, W.G., and Mouat, D.A., 1995. A process for selecting indicators for monitoring conditions of rangeland health. Environmental Monitoring and Assessment 36 45-60. Bater, C.W., and Coops, N.C., (In Review). Evaluating error associated with lidarderived DEM interpolation. Computers & Geosciences. Clark, D.B., Castro, C.S., Alvarado, L.D.A., and Read, J.M., 2004. Quantifying mortality of tropical rain forest trees using high-spatial-resolution satellite data. Ecology Letters 7, 52-59. DOI: 10.1046/j.1461-0248.2003.00547.x Clayoquot Sound Scientific Panel, 1995. Sustainable Ecosystem Management in Clayoquot Sound. Report 5, Halfmoon Bay, British Columbia, Canada.  78  Coops, N.C., Hilker, T., Wulder, M.A., St-Onge, B., Newnham, G., Siggins, A., and Troymow, J.A., 2007. Estimating canopy structure of Douglas-fir forest stands from discrete-return LiDAR. Trees. DOI:10.1007/s00468-006-0119-6 DeLong, C.S., Burton, P.J., and Harrison, M., 2004. Assessing the relative quality of oldgrowth forest: an example from the Robson Valley, British Columbia. BC Journal of Ecosystems and Management 4(2), 1-16. Demarchi, D. A., Marsh, R. D., Harcombe, A. P. and Lea, E. C., 1990. The Environment (of British Columbia). In: R.W. Campbell, N.K. Dawe, I. McTaggart-Cowan, J. M. Cooper, G. W. Kaiser, and M. C. E. McNall (Eds.). The Birds of British Columbia. Royal British Columbia Museum, Environment Canada, and Canadian Wildlife Service, Victoria, British Columbia, Canada Victoria, pp. 55-142. Gobakken, T., and Næsset, E., 2005. Weibull and percentile models for lidar-based estimation of basal area distribution. Scandinavian Journal of Forest Research 20, 490-502. Hill, T., and Lewicki, P., 2006. Statistics: Methods and Applications. Tulsa, StatSoft. Hosmer, D.W., and Lemeshow, S., 2000. Applied Logistic Regression (2nd Edition). Toronto: John Wiley and Sons. Hudack, A.T., Lefsky, M.A., Cohen, W.B., Berterretche, M., 2002. Integration of lidar and Landsat ETM+ data for estimating and mapping forest canopy height. Remote Sensing of Environment 82, 397-416. Lefsky, M.A., Cohen, W.B., Acker, S.A., Parker, G.G., Spies, T.A. and Harding, D, 1999. LiDAR remote sensing of the canopy structure and biophysical properties of Douglas-fir western hemlock forests. Remote Sensing of Environment 70, 339-361.  79  Lefsky, M.A., Hudack, A.T., Cohen, W., and Acker, S.A., 2005a. Geographic variability in lidar predictions of forest stand structure in the Pacific Northwest. Remote Sensing of Environment 95, 532-548. Lefsky, M.A., Turner, D.P., Guzy, M., and Cohen, W.B., 2005b. Combining lidar estimates of aboveground biomass and Landsat estimates of stand age for spatially extensive validation of modeled forest productivity. Remote Sensing of Environment 95, 549-558. MacArthur, R.H., and Horn, H.S., 1969. Foliage profile by vertical measurements. Ecology 50, 802-804. Magnussen, S. and Boudewn, P. Derivation of stand heights from airborne laser scanner data with canopy-based quantile estimators. Canadian Journal of Forest Research 28(7), 1016-1031. McElhinny, C., Gibbons, P., Brack, C., and Bauhas, J., 2005. Forest and woodland structural complexity: its definition and measurement. Forest Ecology and Management 218, 1-24. Meidinger, D. and Pojar, J., 1991. Ecosystems of British Columbia. BC Special Report Series No. 6. BC Ministry of Forests, Victoria, British Columbia, Canada. Mitchell, W.R., Green, R.N., Hope, F.D. and Klinka, K., 1989. Methods for Biogeoclimatic Ecosystem Mapping. Rep. RR89002-KL. British Columbia Ministry of Forest Resources, Victoria, British Columbia, Canada. Morrison, M.L., and Raphael, M.G., 1993. Modeling the dynamics of snags. Ecological Applications 3(2), 322- 330.  80  Næsset, E., 2002. Predicting forest stand characteristics with airborne scanning laser using a practical two-stage procedure and field data. Remote Sensing of Environment 80, 88-99. Næsset, E., 2004. Practical large-scale forest stand inventory using a small-footprint airborne scanning laser. Scandinavian Journal of Forest Research 19, 164-179. Næsset, E., and Økland, T., 2002. Estimating tree height and tree crown properties using airborne scanning laser in a boreal nature reserve. Remote Sensing of Environment 79, 105-115. Noss, R.F., 1999. Assessing and monitoring forest biodiversity: a suggested framework and indicators. Forest Ecology and Management 115, 135-146. O’Connell, A., 2006. Logistic Regression Models for Ordinal Response Variables. London, Sage Publications. Oliver, C.D., and Larson, B.C., 1996. Forest Stand Dynamics. Toronto, John Wiley & Sons. Province of British Columbia, 2005. Monitoring and Evaluation Strategy. FRPA Resource Evaluation Program. British Columbia Ministry of Forests and Ranges, British Columbia Ministry of Environment, and British Columbia Ministry of Agriculture and Lands, Victoria, BC. Sambridge, M., Braun, J., and McQueen, H., 1995. Geophysical parameterisation and interpolation of irregular data using natural neighbours. Geophysical Journal International 122, 837-857. Seielstad, C,A, and Queen, L.P, 2003. Using airborne laser altimetry to determine fuel models for estimating fire behaviour. Journal of Forestry 101(4), 10-15.  81  Sibson, R., 1981. A brief description of natural neighbour interpolation. In V. Barnett (Ed.), Interpreting Multivariate Data. Wiley, Chichester, pp. 21-36. Tabachnick, B.G., and Fidell, L.S., 2001. Using Multivariate Statistics (4th Edition). Toronto, Allyn and Bacon. Warren Wilson, J., 1960. Inclined point quadrats. New Phytologist 59(1), 1-8. Warren Wilson, J., 1963. Estimation of foliage denseness and foliage angle by inclined point quadrats. Australian Journal of Botany 11, 95-105.  82  4 CONCLUSION High resolution DEMs are critical for terrain mapping, and can be used for monitoring indicators such as large mass wasting events, channel planform and cross-sectional morphology, erosion and deposition, and natural hazards (e.g. Table 1.2). They also act as the reference above mean sea level against which lidar vegetation return heights are estimated; thus their accurate derivation is a critical step which may impact the quality of both terrain- and vegetation-based indicators of forest sustainability. Chapter 2 described a methodology for assessing DEM interpolation errors by subsetting lidar data into prediction and validation datasets. While the natural neighbour routine was generally favoured, the intent was not to suggest that this will universally be the case. Instead, the important message is that when generating DEMs, a number of routines, parameterizations, and spatial resolutions should be explored and validated. Furthermore, by stratifying the study area by slope, ground return density, NDVI, and vegetation structural class, a CART-based method to examine interpolation error and map prediction uncertainty was developed. This method may provide users with the ability to assess the quality of their DEMs as it varies across the landscape using very little ancillary data. By mapping prediction uncertainty, users can determine if their data are appropriate for their intended purposes. Once a suitable DEM has been created, lidar vegetation data can be processed and analyzed.  A variety of vegetation-related indicators have been established, including those related to timber, biodiversity, riparian areas, and range management resource values. A key component of both structural and biological diversity, the distribution of wildlife trees is  83  an important indicator of forest health and habitat quality, but little attention has been paid to their estimation using remote sensing technologies. In Chapter 3, common lidar metrics were employed to estimate the cumulative proportions of wildlife tree classes within northwest coastal forests. Using ordinal logistic regression, a model was developed to predict the proportion of stems belonging to the range of WT classes at a given location. The lidar-derived coefficient of variation of the vegetation return heights (log transformed) proved to be the most reliable predictor variable. Ultimately, this research will contribute to our ability to characterize forest structure and health across the landscape.  Future Research Future lidar-related research will likely focus on two key areas. The first is the refinement of existing lidar metrics and the estimation of additional ecological indicators, in particular incorporating small footprint, full waveform data, and possibly lidar intensity data. The second will be the parameterization of lidar systems and surveys to maximize efficiency in terms of cost, coverage, data density, and accuracy.  Recently, digitizers have been developed for small footprint sensors which allow the return intensity to be recorded as a function of time. Full waveform data may provide more refined information about forest canopy architecture and understorey, and improved ground detection. Discrete return sensors record one or more returns per laser pulse, typically first and/or last, and possibly intermediate second, third and fourth returns. These sensors are limited by a reset time, which results in a minimum of several vertical  84  metres (typically 2-5 m) between the time when a given return has been detected and another can be measured. This arguably limits the ability of discrete return lidar to characterize short vegetation with great detail. Full waveform recorders, however, measure the amount of energy backscattered from surfaces within submetre resolution vertical bins, resulting in the capture of much more detail related to canopy and subcanopy architecture. Thus, combining the vertical detail of full waveform data with the recent increases in pulse repetition rates and prevalence of discrete return sensors promises to increase our ability to characterize forest structure.  In addition to height information, discrete laser pulse return data are usually accompanied by an intensity value, which is an indication of the strength or peak amplitude of the return. Some researchers in the lidar community, however, has been weary of intensity data and very little has been done to use it in any quantitative fashion, although Song et al. (2002) and Donoghue et al. (2007) are notable exceptions. Presently, lidar sensors are not calibrated to estimate intensity reliably, nor is there a general understanding about what the information really represents. The data tends to be noisy, as the values vary with platform altitude, sensor calibration, target bidirectional reflectance properties, and target reflectivity. At best, intensity may only be suitable for the classification of broad land cover types. However, the calibration of, and investigation into, intensity may lead to the development of new products that can be extracted from lidar data.  No standard exists for the collection of airborne lidar data. Surveys must be parameterized to insure the lidar data is appropriate for the indicators to be measured (for  85  instance plot based heights vs. individual tree heights). Flying height, pulse repetition rate, maximum scan angle, amount of sidelap, and footprint size influence not only the density of pulse returns, but may also affect the accuracies of ground and vegetation height estimates (Lim et al., 2005; Goodwin et al., 2006; Chasmer et al., 2006). State of the art sensors with pulse repetition rates in excess of 160 KHz are now in production (e.g. Optech Gemini), which may allow higher flying heights and larger swath widths. It remains to be seen, however, if the reduced laser power at these higher rates will have an adverse affect on data quality. Ultimately, managers will have to agree on a suite of indicators, and then research undertaken to identify optimal survey parameters that ensure data densities are appropriate, whilst balancing areal coverage and cost.  This research contributes to our ability to estimate indicators of forest sustainability in an accurate and spatially explicit fashion, enabling managers to carry out monitoring and evaluation programs at the landscape level. In particular, it contributes to the growing body of knowledge relating to the assessment of the status of forest ecosystems using remote sensing technologies.  86  References Chasmer, L., Hopkinson, C., Smith, B., Treitz, P., 2006. Examining the influence of changing laser pulse repetition frequencies on conifer forest canopy returns. Photogrammetric Engineering & Remote Sensing 72(12), 1359-1367. Donoghue, D.N.M., Watt, P.J., Cox, N.J., Wilson, J., 2007. Remote sensing of species mixtures in conifer plantations using LiDAR height and intensity data. Remote Sensing of Environment 110, 509-522. Goodwin, N.R., Coops, N.C., Culvenor, D.S., 2006. Assessment of forest structure with airborne LiDAR and the effects of platform altitude. Remote Sensing of Environment 103,140-152. Lim, K., Treitz, P., Wulder, M., St-Onge, B., Flood, M., 2003. Lidar remote sensing of forest structure. Progress in Physical Geography 27(1), 88-106. Song, J., Han, S., Yu, K., Kim, Y., 2002. Assessing the possibility of landcover classification using LiDAR intensity data. ISPRS Commission III, Graz, Austria.  87  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0066311/manifest

Comment

Related Items