- Email: [email protected]
GLAS data: A case study in the New Forest National Park, UK
International Journal of Applied Earth Observation and Geoinformation 23 (2013) 109–118
Contents lists available at SciVerse ScienceDirect
International Journal of Applied Earth Observation and Geoinformation journal homepage: www.elsevier.com/locate/jag
A novel approach to estimate canopy height using ICESat/GLAS data: A case study in the New Forest National Park, UK Irfan Akhtar Iqbal a , Jadunandan Dash a , Saleem Ullah b,∗ , Ghayyas Ahmad c a b c
School of Geography, University of Southampton, Highfield, Southampton SO17 1BJ, UK Faculty of Geo-Information Science and Earth Observation (ITC), University of Twente, 7500 AA Enschede, The Netherlands Pakistan Forest Institute, Peshawar, Khyber Pakhtunkhwa, Pakistan
a r t i c l e
i n f o
Article history: Received 21 August 2012 Accepted 18 December 2012 Keywords: Canopy height LiDAR Remote sensing
a b s t r a c t The Geoscience Laser Altimeter System (GLAS) aboard Ice, Cloud and land Elevation Satellite (ICESat) is a spaceborne LiDAR sensor. It is the first LiDAR instrument which can digitize the backscattered waveform and offer near global coverage. Among others, scientific objectives of the mission include precise measurement of vegetation canopy heights. Existing approaches of waveform processing for canopy height estimation suggest Gaussian decomposition of the waveform which has the limitation to properly characterize significant peaks and results in discrepant information. Moreover, in most cases, Digital Terrain Models (DTMs) are required for canopy height estimation. This paper presents a new automated method of GLAS waveform processing for extracting vegetation canopy height in the absence of a DTM. Canopy heights retrieved from GLAS waveforms were validated with field measured heights. The newly proposed method was able to explain 79% of variation in canopy heights with an RMSE of 3.18 m, in the study area. The unexplained variation in canopy heights retrieved from GLAS data can be due to errors introduced by footprint eccentricity, decay of energy between emitted and received signals, uncertainty in the field measurements and limited number of sampled footprints. Results achieved with the newly proposed method were encouraging and demonstrated its potential of processing full-waveform LiDAR data for estimating forest canopy height. The study also had implications on future full-waveform spaceborne missions and their utility in vegetation studies. © 2012 Elsevier B.V. All rights reserved.
1. Introduction Geoscience Laser Altimeter System (GLAS), the first of a new generation spaceborne LiDAR onboard the Ice, Cloud and land Elevation Satellite (ICESat) mission, was launched on January 13, 2003. GLAS is the first spaceborne instrument which can fully digitize the backscattered waveform from incident surfaces (Wagner et al., 2006). GLAS offers near-global coverage by carrying three laser instruments (Laser 1, 2 and 3) with only one operational at a time. After the failure of Laser 1 in March 2003 (Schutz et al., 2005), GLAS followed a modified plan of 33 days sub-cycle with a 91 days repeat orbit (Abshire et al., 2005). Under the revised plan, Laser 3 acquired data during February–March, May–June, and October–November each year (Schutz et al., 2005), until its failure in October, 2008. Laser 2 also ceased emitting light on October 11, 2009 (NASA, 2009).
∗ Corresponding author at: Faculty of Geo-Information Science and Earth Observation (ITC), University of Twente, P.O. Box 217, 7500 AE Enschede, The Netherlands. Tel.: +31 53 4874 372. E-mail addresses: [email protected] (I.A. Iqbal), [email protected] (S. Ullah). 0303-2434/$ – see front matter © 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.jag.2012.12.009
After seven years of near-global coverage and 15 laser-operation campaigns, ICESat has stopped collecting science data. However, data from these historical acquisitions can be exploited to evaluate the potential of these data in other terrestrial applications and thus, can have useful implications for the next generation of sensor development, for example ICESat-2, DESDynI. Originally, ICESat was launched to monitor polar ice-sheet dynamics and its role in changing the global sea-level (Zwally et al., 2002). However, other scientific objectives of this mission also include surface reflectivity, precise measurement of land topography and vegetation canopy heights. Light Amplification by Stimulated Emission of Radiation (LASER) altimetry has supplemented the existing remote sensing techniques because of its high precision and capability to estimate canopy height. Light Detection and Ranging (LiDAR) directly measures the vertical component of vegetation, thus has the potential to measure the structural attributes of vegetation (Dubayah et al., 2000; Lefsky et al., 2002; Nelson, 2010). It has the capability to provide information at both canopy and subcanopy levels which is crucial for understanding forest health and regeneration capacity (Stone et al., 2000; Todd et al., 2003), forest fires (Andersen ˜ et al., 2003) and present and future carbon stock et al., 2005; Riano
110
I.A. Iqbal et al. / International Journal of Applied Earth Observation and Geoinformation 23 (2013) 109–118
estimation (Dean et al., 2004). Moreover, tree heights and vegetation density derived from LiDARs prove to be important inputs for aboveground biomass estimation, understanding global carbon cycle and ecosystem dynamics (Ranson et al., 2004). Forests, being important for their role in the global climate (UNFCCC, 1997, 2009a), need to be quantitatively as well as qualitatively analyzed. By estimating the structural attributes (e.g. aboveground biomass) of these ecosystems, the amount of carbon stored can be calculated as half the dry weight of biomass (Drake et al., 2002; Houghton, 2008; IPCC, 2003; UNFCCC, 2009b). Conventionally, aboveground biomass is estimated by using the DBH (Keller et al., 2001), though straightforward, estimation of aboveground biomass (using DBH) on ground is time consuming as well as expensive (Houghton, 2005) and mostly involves destructive sampling (Hiratsuka et al., 2003). Alternatively, due to its spatial and temporal coverage, remote sensing provides opportunity to derive vegetation biophysical variables and some of them are related to vegetation structure. Although DBH cannot be directly measured from remotely sensed data; other biophysical parameters like canopy height and canopy cover can serve as a proxy for estimating the amount of aboveground biomass. Several studies have attempted to use remotely sensed data (both passive and active (RADAR)) for estimating the biophysical parameters of vegetation and in turn, to estimate aboveground biomass e.g. (Dobson et al., 1992; Kaufman et al., 1990; Muukkonen and Heiskanen, 2007). However, the accuracy and sensitivity of these data types drop-off in areas of closed canopies and high-biomass (Waring et al., 1995; Mitchard et al., 2012).
2. Application of GLAS data in vegetation studies GLAS transmits short pulses (4 ns width equivalent to 60 cm in surface elevation) of infrared (1064 nm) and green (532 nm) light at a frequency of 40 Hz. The nominal diameter of the illuminated spot, the footprint, is ∼70 m, space at 172 m along the track; however, its size and ellipticity have varied through the time (Pang et al., 2008). The reflected photons are collected by a telescope of 1 m diameter, and their intensities are binned every nanosecond. Height (or range) can be worked out as the product of speed of light (3 × 108 m/s) and one-way travel time of the photons. Intensities of the reflected photons are stored in 544 bins over ice sheet and land, corresponding to a height of 81.6 m (Brenner et al., 2003). In steeply sloped areas and/or areas where feature heights exceed 81.6 m, GLAS waveform would truncate, making it inconvenient to derive range information. For this reason, in later operations height extent was increased to 150 m over land, using a ‘waveform compression scheme’ (Harding and Carabajal, 2005) according to which the first 152 bins correspond to 60 cm each. Fig. 1 shows the waveform components, used to estimate various physical attributes, which will be referred to in the succeeding text. The numbers 25, 50, 75 and 100 indicate height quartiles where respective percentage of waveform intensity is concentrated. GLAS products, typically GLA01 and GLA14, have been used for vegetation classification and biomass estimation (Boudreau et al., 2008; Lefsky et al., 2005; Ranson et al., 2004), and seasonal changes in vegetation (Duong et al., 2008), where most of the information is retrieved from GLA01 (the ‘raw’ waveform). Particular to vegetation studies, Ranson et al. (2004) found the front slope angle (Fig. 1) strongly correlated to canopy density and vertical variability of upper canopy. These relationships were further exploited in a classification of forests in central Siberia. Aboveground biomass estimation from GLAS data has been demonstrated to have agreement with field based or airborne measurements (Boudreau et al., 2008; Helmer and Lefsky, 2006; Lefsky
Fig. 1. Illustration of waveform components used to estimate canopy height. The waveform has been smoothed with a moving average of 5 to reduce noise and favour clarity. The numbers 25, 50, 75 and 100 indicate height quartiles where respective percentage of waveform intensity is concentrated. The distance between the two dotted blue lines is the waveform extent. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of the article.)
et al., 2005). Extraction of canopy height from raw waveform is complicated in areas of varied topography due to slope effects (Harding and Carabajal, 2005). GLAS returns from flat surfaces can be represented by a Gaussian curve; while flat areas with homogeneous vegetation can be approximated by two Gaussian curves. With increasing topographic variation and vegetation heterogeneity, the number of Gaussian peaks increases (Duncanson et al., 2010). To allow for such complexities, methods for canopy height estimation have been developed which are generally classified (Chen, 2010b) as (i) statistical methods and (ii) direct methods. In statistical methods, waveform processing is supplemented by additional information from the underlying terrain. So-called terrain indices are extracted either from the available Digital Terrain Models (DTMs) (Lefsky et al., 2005; Rosette et al., 2008) or derived from the waveforms (Lefsky et al., 2007; Pang et al., 2008). Direct methods concentrate on waveform processing for accurate identification of canopy top and the ground return with less, if any, assistance from auxiliary data, e.g. DTMs. Lefsky et al. (2005) derived canopy height as a linear function of waveform extent and terrain index. Waveform extent was taken as the difference between the first and last threshold crossings of the waveform; whereas terrain index was calculated, from Shuttle RADAR Topography Mission (SRTM) data, as the difference between minimum and maximum elevations within a footprint. Having both waveform extent and terrain index identified, canopy height was modelled as given in the equation: H = b0 (w − b1 g)
(1)
where, H is the measured canopy height, w is the waveform extent in m, g is the terrain index in m, b0 and b1 are coefficients applied to the waveform and terrain index, respectively. With their method, Lefsky et al. (2005) were able to explain a maximum variation of 68% with an RMSE of 4.85 m in canopy heights. They attributed the unexplained variance in canopy height estimation to the ‘measurement’ of waveform extent, more importantly, the identification of ground return (or peak) in the raw waveform. For efficient ground peak detection, Gaussian components are fitted to the raw waveform (Brenner et al., 2003) and iteratively reduced to six (Sun et al., 2008), using the Least Squares approach (Markwardt, 2009). Harding and Carabajal (2005) showed that out of the six peaks, the last peak with the highest amplitude is a representative of the ground. However, in densely vegetated areas, due to energy attenuation through the canopy (Parker et al., 2001),
I.A. Iqbal et al. / International Journal of Applied Earth Observation and Geoinformation 23 (2013) 109–118
the ground peak may have lower amplitude than adjacent peaks (Boudreau et al., 2008). In such cases, Boudreau et al. (2008) suggested to consider the peak with highest amplitude as the ground peak, even if it is not the last one. However, the method will lead to inaccuracies in estimating ground and therefore in canopy height. Rosette et al. (2008) reported an underestimation of 1.20 m in canopy height by considering the first or second highest amplitude peak as ground return. A study by Sun et al. (2008) demonstrated the utility of height quartiles (Fig. 1) to predict canopy heights. Quartiles at 25, 50, 75 and 100 were calculated on the basis of percentage of reflected energy. Quartile at 25, for example, represented the aboveground level at which 25% of the reflected energy is accumulated. Similarly, quartile at 100 is the distance between signal start and the ground return, which is essentially the waveform extent. Precise detection of the ground peak is important because it is used as a reference for calculating quartiles. Sun et al. (2008) used an automated approach for ground peak detection in which bin values (or amplitudes), from signal end to signal start, were compared to the two neighbouring values. The first peak, from signal end, with a value more than the half width of the transmitted laser pulse was taken as the ground peak. Their method explained a maximum variability of 83% in canopy height. However, signal strength was not considered which may label the first ‘significant peak’ as ground even if it has a lower energy. Ground or flat surfaces tend to register stronger responses. Another automated approach for peak identification with a search window of 5 ns was proposed by Duong et al. (2008). The search window operated on the normalized (amplitude at each instant divided by the sum of all amplitudes) and smoothed waveform, moving from beginning through end at a step of 1 ns. A peak was identified if the central bin had an incrementally higher amplitude value than the two bins on either side. Further, six Gaussian peaks were fitted to the resulting waveform, the last being representative of the ground return. Some irregular waveforms could not be adequately decomposed by Gaussian distributions and often missed a mode (Duncanson et al., 2010; Duong et al., 2009). The issues (of missed peaks) significantly limits the accuracy of estimated terrain and canopy surfaces. The last mode is critical because, principally, canopy height is the distance defined by the time difference between the signal start after noise and the centroid of the last peak (Sun et al., 2008). Although Gaussian fitting has shown considerable promise in various studies, its efficiency is ‘less satisfactory’ for pulses with high amplitudes (Wagner et al., 2006). Furthermore, in cases of low amplitude pulses, the extracted parameters are less accurate which demands stringent consideration of the waveform processing algorithms (Mallet and Bretar, 2009). There is a need to develop a robust waveform fitting method which circumvents the stated problems, allows a more accurate retrieval of canopy height and reduces uncertainties in derived variables such as aboveground biomass. Addressing the need, in this study, first a new method involving Fourier analysis of waveform processing for quantification of canopy height from GLAS data was developed and then the canopy heights derived from GLAS data were validated using field data.
3. Study area The study was selected within the New Forest National Park in Southeast England (50◦ 51 59 latitude and 01◦ 40 50 longitude) covering an area of 29262.36 ha. This site holds the status of “Site of Community Importance” since December 2004 (JNCC, 2006). Also, it was designated as a Special Area of Conservation in April 2005. The area is covered with deciduous broadleaved (oak/beach 29%) and coniferous woodland (17%). The remaining area is a mix of bogs,
111
marshes, scrub and heath. It is under the active management of the Forestry Commission of Great Britain. The study site is subdivided in compartments; a compartment is a managerial unit which has similar characteristics within the area; such as, soil types, vegetation, productivity and topography (Hubbard et al., 1998). Most of the compartments have pure and even-aged stands; however, in some compartments age-diversity does exist. Fig. 2 shows the GLAS data projected over the New Forest National Park. 4. Data 4.1. GLAS/ICESat data There are 15 GLAS data products available (Zwally et al., 2006). Data products GLA01 and GLA14 of release 28 and 29 were used in this study. Data were procured from the National Snow and Ice Data centre (NSIDC) through ICESat/GLAS data subsetter (URL2, 2009). GLA01 is a level 1A global altimetry product which contains intensities of transmitted and received waveforms. GLAS digitizes these intensities as counts which are then converted to volts using calibration tables in GLA01 header records. The transmitted pulses were arranged in time order, which means that the value of the first sample is for the sample closest to the spacecraft in time. Conversely, the received echo is in time reversed order, which means the value of the first sample is for the sample farthest from the spacecraft in time (NSIDC, 2008). GLA14 is a level 2 region-specific land surface altimetry data which is derived from level 1 products (GLA05 and GLA06). Elevations are given in m with respect to a reference ellipsoid accompanied by geoid elevations which allow users to calculate orthometric heights. Similar to GLA01, GLA14 also contains geolocation of footprint centroid, unique record and shot numbers, and a time stamp with Universal Coordinated Time (UTC) in J2000 (referenced from noon on January 1, 2000). As the transmitted and received pulses are assumed to resemble several Gaussian curves (Brenner et al., 2003), a set of Gaussian pulses (up to six) is fitted to the waveform during processing. For this purpose, GLA14 variables also include amplitudes, areas and standard deviations of 6 Gaussian distributions if the user is interested in Gaussian fitting. 4.2. Ancillary data In addition to GLAS data, CORINE land cover (CLC) classification 2000 (EEA, 2009) was used to identify forest patches in the study area and extract GLAS footprints over them. CLC 2000 is a geographic database which covers European countries and gives information on land cover classes at a spatial resolution of 100 m. Cloud free Landsat 7 (ETM+) imagery (URL1, 2009), dated April 2003, was used to visualize GLAS transects. A compartment and sub-compartment database of the study area was procured from the Forestry Commission of Great Britain (FCGB, 2009). For each compartment and sub-compartment in the study area the data-base provided information on: planting year, area of the compartment, type of species and yield class. The Forestry Commission has also designed yield tables to aid forest management in British conditions (Hamilton and Christie, 1971). From the species-specific yield tables and general yield class curves, crop height can be estimated if yield class and planting year of a compartment are known. 5. Methodology In the current study, a new method based on Fourier transforms, was used to process the GLAS waveform. Fourier transform converts a waveform in a time domain to a waveform in the
112
I.A. Iqbal et al. / International Journal of Applied Earth Observation and Geoinformation 23 (2013) 109–118
Fig. 2. Representation of ICESat tracts projected over the study area. Eastern tracks (in red box) were used for canopy height estimation, comprising of 145 footprints. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of the article.)
frequency domain with a real and an imaginary part (cosine and sine) (Bloomfield, 2000).
5.1. Data pre-processing The datasets used in this study were registered to different reference systems. To have the datasets coherent, they were all re-projected to a common reference system: the World Geodetic System 1984 (WGS1984). Information contained in GLA01 and GLA14 were extracted, using the IDL codes provided by the National Snow and Ice Data Centre (NSIDC, 2009). GLA01 waveforms were then linked to GLA14 on the basis of record and shot numbers. All GLAS footprints over the study area were identified and waveforms with high level of echo saturation were filtered out using a flag called “i gval rcv” (NSIDC, 2009a). This flag ranges from 0 to 250, with some values labelled as −999 indicating distorted waveforms. All waveforms having an i gval rcv of −999 were discarded. Further, CORINE land cover classification was used to identify forest patches in the study area and a total 145 GLAS footprints over these patches were extracted for analysis.
5.2. Field data collection Out of the total 145 footprints, 50 were randomly selected for field data collection to avoid bias and obtain representative samples (Lefsky et al., 2005). From these 50 identified footprints, 23 fell in the management boundaries of the New Forest National Park for which the Forestry Commission yield data was readily available. For the remaining 27 footprints, field work was conducted, in October 2009, to obtain height information for 19 points; remaining footprints could not be sampled due to logistic difficulties (e.g., limited time, weather conditions and refused access). A handheld Garmin eTrex GPS receiver was used, with an accuracy of about 7 m, to locate the footprint centres. Within each footprint, angle of elevation (top) and depression (base) of trees greater than 10 cm DBH were recorded with a clinometer. Also, horizontal distance from trees, ground elevations, species and description of footprint were recorded. Data collection was flexible in terms of the sample size; however, care was taken to characterize height variability. In compartments with homogeneous cover, fewer trees were measured than in uneven-aged plots. In uneven-aged areas, within footprints,
I.A. Iqbal et al. / International Journal of Applied Earth Observation and Geoinformation 23 (2013) 109–118
113
Fig. 3. Conceptual layout of sample plot for measuring canopy heights in the field.
data was collected in 4 annular plots, one in centre and three on the footprint margin at approximate azimuths of 0◦ , 120◦ and 240◦ . Radius of central plot was 10 m while the radii of smaller marginal plots were 5 m each. Radius of the central plot was greater than the marginal ones, to account for the reflected LiDAR energy which has a central maximum and fades away radially outwards (Duong et al., 2006b). Fig. 3 gives the conceptual sampling design for a footprint.
Fig. 4. Using a defined threshold level to discriminate between noise, signal-start and signal-end (as in Eq. (2)).
In this study, the first 200 bins were used to calculate the threshold parameters as shown below: Mean noise = E¯ =
5.3. Waveform processing 5.3.1. Waveform normalization Since different waveforms had different amplitude intensities due to changing atmospheric conditions and/or decay in the laser device (Duong et al., 2008), it would be improper to apply a single threshold value to all waveforms. Variation in amplitude intensities complicates comparison of waveforms from different campaigns and hence requires normalization, which was achieved by dividing the received energy at each instance (Ei ) by the total received 544 E (Duong et al., 2006a). energy (ET), where total energy is i=1 i
544
Mathematically, a normalized waveform is expressed as Ei /
E. i=1 i
200 i=1
Ei /200
Noise standard deviation =
(3)
200 (E − E) ¯ 2 i i=1
200
(4)
where, Ei is the received energy at time i. 5.3.4. Detecting ground, canopy-top and canopy height In order to estimate canopy height from the smoothed waveform, first, all peaks in the smoothed waveform were identified by a search window (1 ns) which moved from start to end of the waveform. The window calculated the difference between two succeeding amplitudes, and stopped until it found the point where the first derivative changed from positive to negative or vice versa.
5.3.2. Waveform smoothing and fitting The normalization step was followed by smoothing the normalized waveform and fitting a modelled waveform to it. For this purpose, Fast Fourier Transform was used, which is an efficient algorithm (Zadiraka and Igisinov, 1973) to compute discrete Fourier transform (analysis) and its inverse (synthesis). Power Spectral Density (PSD) was applied to find the frequency region of maximum signal strength. Based on the PSD, we used 30 harmonics for smoothing of and model-fitting to the normalized waveform.
5.3.3. Threshold and background noise The GLAS echo, telemetered to 544 bins, corresponds to a height of 81.6 m over land. Out of 544, the first 200 bins (30 m) were used to define the noise level (using Eq. (3)). The remaining 344 bins (51.6 m) were sufficient to represent the vegetation structure without truncation. A noise threshold level was used to single out the signal from noise and identify locations of signal-start and signal-end. Start and end of the signal are threshold crossings of the leading and trailing edges (Harding and Carabajal, 2005). We set the threshold level to the sum of mean noise and 4 times the noise standard deviation as in Lefsky et al. (2005) using Eq. (2). The resulting waveform is shown in Fig. 4. Threshold level = mean noise + 4 × noise standard deviation
(2)
Fig. 5. Identified peaks (in red circles) and troughs (in green circles) on the fitted waveform. Peak with maximum amplitude is also marked for determining ground, canopy top and canopy height. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of the article.)
114
I.A. Iqbal et al. / International Journal of Applied Earth Observation and Geoinformation 23 (2013) 109–118
Fig. 6. Final parameterization of GLAS waveform. Red and green circles show the waveform peaks and troughs, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of the article.)
Fig. 5 illustrates all peaks (in red circles), troughs (in green circles) and peak maximum amplitude. To detect canopy top, ground and canopy height, the peak with the maximum amplitude was used as a reference. From this reference point, the search window moved towards the signal-end and searched for a ground return peak. Peak with amplitude onefifth the maximum amplitude, above the defined threshold, was defined as the ground return. It was assumed, that due to reduced penetration and decay of energy through the canopy (Parker et al., 2001), the return energy can be as low as one-fifth of the maximum amplitude. If these conditions were not met, the maximum amplitude peak was designated as the ground return. To infer the position of the canopy, the window moved from the maximum amplitude peak towards the signal start to search for canopy point. A canopy point is defined here as the first peak after the signal start which is above the threshold level. In this study, we applied a threshold of 0.003 V, using Eq. (2). The leading edge of a canopy point indicates the canopy-top, because as soon as the signal interacts with the uppermost foliage and/or branches, it is excited. Canopy height was calculated as the distance between ground return and canopy top (Fig. 6). In case the ground return and the canopy top peaks coincided, the canopy height was assumed to be zero. 5.4. Field data analysis Angles of elevation and depression of tree top and bottom (respectively), and distance to tree from the point of observation were used to calculate tree heights (Eq. (5)). Tree height (field) = d × tan + d × tan ˚
(5)
where, d is the distance between the tree and the point of observation (in m), is the angle of elevation of the tree top in degrees, and ˚ is the angle of depression of the tree bottom in degrees. Because the waveform is sensitized by the uppermost foliage and/or branches, it represents the maximum canopy height inside the footprint. For this reason, only field heights of dominant trees were used for analysis. Moreover, there was a time lag of approximately 3–5 years between GLAS acquired data and field collected data, during which height increments were expected. Annual
Fig. 7. Correlation of canopy heights obtained from waveform analysis and those collected from field and the Forestry Commission data. Outliers are shown in red. Red dotted line shows 1:1 line. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of the article.)
height increments for each species were worked out from yield tables and subtracted from field calculated canopy heights, to get adjusted canopy heights. Canopy heights for the footprints falling in Forestry Commission sub-compartments were estimated from yield class curves in Hamilton and Christie (1971). Crop age was estimated as the difference between planting year and GLAS data acquisition year, and was used to estimate canopy height from yield class curves. 6. Results and discussion Two tailed Pearson’s correlation in Statistical Package for Social Scientists (SPSS v.16) was used to compare any two datasets resulting from the adopted methodology. R-value and RMSE were used as indicators of ‘goodness of agreement’ between corresponding variables. 6.1. Accuracy of GLAS derived canopy heights Canopy heights calculated from field and Forestry Commission data correlated with those obtained from GLAS waveforms showed a correlation of 0.59, with an RMSE of 5.18 m (Fig. 7). Heights derived from GLAS waveforms are hereinafter referred to as ‘GLAS estimated height’. The offset of about 6 m (in Fig. 7) represents positive bias in GLAS estimated canopy heights. Also, the data has a wider spread around the correlation line. The low correlation between the two datasets was traceable to six outliers, as shown in Fig. 7 (in red). In these six cases, the over and under estimation of canopy heights reached 11–13 m, respectively. Analysis revealed that over mixed and/or multiple canopies, the GLAS estimated heights were positively biased. In cases of thick vegetation cover, the heights were under-estimated. Biases in canopy height (both positive and negative) are also induced by the uncertainty in ground detection. For example, litter on forest floor is a possible reason of multiple scattering (Brenner et al., 2003; Mallet and Bretar, 2009). The extent of over and under estimations due to litter depends on its thickness, three-dimensional distribution and spectral properties (Kimes et al., 2006). Soil moisture may also be a factor to cause decay in the return pulse by absorbing some of the incoming energy (Brenner et al., 2003). However further studies are needed
I.A. Iqbal et al. / International Journal of Applied Earth Observation and Geoinformation 23 (2013) 109–118
Fig. 8. (a) Extreme case of over-estimation up to 11 m. Red and green circles show waveform peaks and troughs, respectively and (b) Footprint location of waveform (in (a)). The red dot shows the footprint centre, labelled with the record and shot number. Image taken from Google Earth (NOAA, 2009). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of the article.)
115
Fig. 9. (a) Extreme case of under-estimation up to 13 m. Red and green circles show waveform peaks and troughs, respectively and (b) location of waveform (in (a)). The footprint, centred at the red dot, covered young Douglas Fir crop. Image taken from Google Earth (NOAA, 2009). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of the article.)
to quantify the effect of soil moisture on LiDAR pulse absorption. After multiple scattering, a strong signal will under-estimate the ground; whereas a weaker signal will cause over-estimation due to increased path length (Brenner et al., 2003). Such uncertainties is ground detection are reflected in canopy height estimation.
dominant canopies have influenced the returned energy and resulted in slight rise at the signal start (Fig. 8a), which is introducing a potential source of error, thus failing to represent all vegetation within the footprint. Fig. 8b gives an aerial view of the location represented by waveform in Fig. 8a. The red circle approximates the GLAS footprint.
6.1.1. Extreme over-estimation Fig. 8a represents the waveform of a footprint which covered a sparse patch of oaks (196 years of age) dominated by low canopies. In the study area, oaks are not managed commercially; hence, height variations are expected. The LiDAR signal was excited by foliage of a few dominant trees, resulting in slight amplitude rise (red arrows in Fig. 8a); however, a strong response was registered from the trees in the under storey. For similar conditions, waveforms behaved the same way, resulting in over-estimated canopy heights. Vegetation structural heterogeneity is documented to complicate the relationship between LiDAR derived and field observed canopy heights (Clark et al., 2004). The incident LiDAR beam is directly affected by vertical (and horizontal) distribution of foliage and woody elements among other variables such as density and size of the canopy (Ni-Meister et al., 2001). LiDAR waveform structure, in particular, has been noted to be sensitive to canopy shape, reflective properties and associated photon interactions (North et al., 2010). A small number of
6.1.2. Extreme under-estimation Fig. 9a is a representative waveform where high underestimation (13 m) in canopy height determination was observed. This footprint covered a young crop (29 years of age) of Douglas Fir with dense canopy, because of which the signal was obstructed (Takeda, 2004). The signal, that penetrated the canopy, was weak and could hardly register a high amplitude peak. Another factor could be soil moisture, which is expected to be more under dense canopy cover, causing signal decay. Such conditions introduce uncertainties in precise ground detection, as discussed in Section 6.1, and ultimately result in biased canopy height estimates. Fig. 9b gives an aerial view of the location represented by waveform in Fig. 9a. The area is covered with young crop of Douglas Fir (29 years of age). The red circle approximates the footprint size. Another important factor, to which GLAS data is reportedly sensitive, is crop age (Helmer and Lefsky, 2006; Sun and Ranson, 2000) which might have introduced bias in GLAS estimated canopy heights. Note that Figs. 8a and 9a are representatives of footprints
116
I.A. Iqbal et al. / International Journal of Applied Earth Observation and Geoinformation 23 (2013) 109–118
method for retrieving ground elevation from Gaussian peaks (Chen, 2010a; Duong et al., 2008; Lefsky et al., 2005; Rosette et al., 2008). To date, Fourier method has not been applied to characterize LiDAR waveforms for canopy height estimation. Although, topographic influences cannot be overlooked, our method of waveform processing gives encouraging results in the study site. This implies the advantage of Fourier method of waveform analysis to automatically estimate canopy heights, in the absence of DTM, and less computational requirements. 7. Conclusion
Fig. 10. Correlation of field and GLAS estimated heights after removal of all six outliers.
covered with old oaks (196 years) and young Douglas Fir (29 years), respectively. It is interesting to note that four out of the six outliers fall within the Forestry Commission compartments, for which canopy heights were estimated from the yield class curves. Only by removing these four outliers from analysis, the correlation value improved from 0.59 to 0.72 and RMSE lowered from 5.18 m to 4.07 m. Yield tables are species-specific generalized models created for compartments and sub-compartments, covering larger spatial extent than covered by GLAS footprint. There is a possibility that variation in site qualities exist on (sub)compartment level, which affects tree growth. Furthermore, the Forestry Commission yield models are not dynamic and depict a mean trend on (sub)compartment level (Rosette et al., 2008). This seems to be a valid reason for biases in canopy heights extracted from yield class curves and those estimated from GLAS waveforms. Moreover, this analysis demonstrates robustness of the new method of waveform processing for canopy height estimation, showing ‘good’ agreement with field measured heights. Due to the uncertainties explained in Sections 6.1.1 and 6.1.2, these waveforms (six in total: 4 from Forestry Commission and 2 from field data) were removed to see the effectiveness of the method. Removal of these outliers further improved the coefficient of correlation from 0.59 to 0.79, and lowered the RMSE from 5.18 m to 3.18 m as shown in Fig. 10. With an automated method of ground identification based on iterative Gaussian fitting, Rosette et al. (2008) were able to explain 74% of the variance in canopy heights in Forest of Dean, UK. It is logical to compare results from Rosette et al. (2008) with this study, because the two study sites (Forest of Dean and the New Forest National Park) have similar environmental conditions and are managed by the Forestry Commission of Great Britain. It is important to mention that using DTM corrections, Rosette et al. (2008) were able to explain 89% of variability in canopy heights. Likewise, recently (Xing et al., 2010) came up with a modification to Lefsky et al. (2005) exploiting the non-linear relationship between waveform extent and maximum canopy height. Xing et al. (2010) were able to explain a maximum variation of 92% in canopy height in the complex terrain (0◦ –30◦ slope) of Wangqing forest, China. The method applied in this study, however, does not require any supplementary information from auxiliary data, i.e. the topographic correction. So far, Gaussian decomposition has been used for ground peak identification; even so, there is no widely accepted
Aim of the study was to evaluate the potential of ICESat/GLAS to estimate canopy height in the New Forest National Park, Southeast England. Canopy height from GLAS altimetry product (GLA01) was estimated and compared with field measured canopy height. Inaccuracies in ground peak detection have been reported to introduce uncertainties in canopy height estimation. This study demonstrated the ability of Fourier Transformation to tackle such issues. Moreover, an automated method for detecting major peaks in the ‘raw’ waveform was developed. Due to signal attenuation, it was assumed that the actual ground return can be the peak with energy 20% lower than the peak of maximum amplitude. Effects of forest litter and soil moisture could be possible reasons for signal weakening. Our presented method explained 79% of the variation in canopy height estimation compared to field data. Other factors which were not investigated in this study; like footprint eccentricity, representation of field measured height in GLAS data and energy decay between transmitted and received pulses may justify the unexplained variance. Future work will investigate the effect of footprint eccentricity, slope and energy decay between transmitted and receive pulse on retrieval of canopy height from GLAS data. The study succeeded to achieve encouraging results with minimum GLAS parameters. This implies that application of other ancillary variables and correction flags may enhance the results. Moreover, due to the possibility of losing some ‘significant’ peaks, the existing assumption of fitting six Gaussian modes to the raw waveform may be relaxed to include more Gaussian peaks. We strongly recommend application of the presented methodology in other similar studies to validate its repeatability and transferability. Acknowledgements The authors wish to thank the European Commission for awarding us with Erasmus Mundus – Scholarships. The authors also acknowledge the help of Dr. Zornitza Ivanova Yovcheva and Mr. M. Sohail Khan for editing the initial drafts and giving a critical feedback. References Abshire, J.B., Sun, X., Riris, H., Sirota, J.M., McGarry, J.F., Palm, S., Yi, D., Liiva, P., 2005. Geoscience Laser Altimeter System (GLAS) on the ICESat Mission: onorbit measurement performance. Geophysical Research Letters 32, L21S02, http://dx.doi.org/10.1029/2005GL024028. Andersen, H.-E., McGaughey, R.J., Reutebuch, S.E., 2005. Estimating forest canopy fuel parameters using LIDAR data. Remote Sensing of Environment 94, 441–449. Bloomfield, P., 2000. Fourier Analysis of Time Series: An Introduction (Wiley Series in Probability and Statistics). Wiley-Interscience, 78 pp. Boudreau, J., Nelson, R.F., Margolis, H.A., Beaudoin, A., Guindon, L., Kimes, D.S., 2008. Regional aboveground forest biomass using airborne and spaceborne LiDAR in Québec. Remote Sensing of Environment 112, 3876–3890. Brenner, A.C., Zwally, H.J., Bentley, C.R., Csath¢, B.M., Harding, D.J., Hofton, M.A., Minster, J.B., Roberts, L., Saba, J.L., Thomas, R.H., Yi, D., 2003. Derivation of range and range distributions from laser pulse waveform analysis for surface elevations, roughness, slope, and vegetation heights. In: Algorithm Theoretical Basis Document, Version 4.1. Goddard Space Flight Center, Greenbelt, MD, USA, p. 92, http://www.csr.utexas.edu/glas/atbd.html
I.A. Iqbal et al. / International Journal of Applied Earth Observation and Geoinformation 23 (2013) 109–118 Chen, Q., 2010a. Assessment of terrain elevation derived from satellite laser altimetry over mountainous forest areas using airborne lidar data. ISPRS Journal of Photogrammetry and Remote Sensing 65, 111–122. Chen, Q., 2010b. Retrieving vegetation height of forests and woodlands over mountainous areas in the Pacific Coast region using satellite laser altimetry. Remote Sensing of Environment 114, 1610–1627. Clark, M.L., Clark, D.B., Roberts, D.A., 2004. Small-footprint lidar estimation of subcanopy elevation and tree height in a tropical rain forest landscape. Remote Sensing of Environment 91, 68–89. Dean, C., Roxburgh, S., Mackey, B.G., 2004. Forecasting landscape-level carbon sequestration using gridded, spatially adjusted tree growth. Forest Ecology and Management 194, 109–129. Dobson, M.C., Ulaby, F.T., LeToan, T., Beaudoin, A., Kasischke, E.S., Christensen, N., 1992. Dependence of radar backscatter on coniferous forest biomass. IEEE Transactions on Geoscience and Remote Sensing 30, 412–415. Drake, J.B., Dubayah, R.O., Clark, D.B., Knox, R.G., Blair, J.B., Hofton, M.A., Chazdon, R.L., Weishampel, J.F., Prince, S., 2002. Estimation of tropical forest structural characteristics using large-footprint lidar. Remote Sensing of Environment 79, 305–319. Dubayah, R., Knox, R., Hofton, M., Blair, J.B., Drake, J., 2000. Land surface characterization using LiDAR remote sensing. In: Hill, M.J., Aspinall, R.J. (Eds.), Spatial Information for Land Use Management. Gordon and Breach Science Publishers, Singapore, pp. 25–38. Duncanson, L.I., Niemann, K.O., Wulder, M.A., 2010. Estimating forest canopy height and terrain relief from GLAS waveform metrics. Remote Sensing of Environment 114, 138–154. Duong, H., Lindenbergh, R., Pfeifer, N., Vosselman, G., 2009. ICESat full-waveform altimetry compared to airborne laser scanning altimetry over The Netherlands. IEEE Transactions on Geoscience and Remote Sensing 47, 3365–3378. Duong, H., Pfeifer, N., Lindenbergh, R., 2006a. Analysis of repeated ICESat full waveform data: methodology and leaf-on/leaf-off comparison. In: Workshop on 3D Remote Sensing in Forestry, Vienna. Duong, H., Pfeifer, N., Lindenbergh, R., 2006b. Full-waveform analysis: ICESat laser data for land cover classification. In: Remote Sensing: From Pixels to Processes, ISPRS Commission VII Mid-term Symposium, Enschede, the Netherlands, pp. 30–35. Duong, V.H., Lindenbergh, R., Pfeifer, N., Vosselman, G., 2008. Single and two epoch analysis of ICESat full waveform data over forested areas. International Journal of Remote Sensing 29, 1453–1473. EEA, 2009. Corine Land Cover 2000 (CLC2000) 100 m – Version 12/2009. Available: http://www.eea.europa.eu/data-and-maps/data/corine-landcover-2000clc2000-100-m-version-12-2009 (accessed 05.10.09). 2009. Forestry Commission in England. Available: FCGB, http://www.forestry.gov.uk/forestry/HCOU-5YHJS7 (accessed 22.09.09). Hamilton, G.J., Christie, J.M., 1971. Forestry Commission Booklet No. 34 Forest Management Tables (Metric). Her Majesty’s Stationery Office. Harding, D.J., Carabajal, C.C., 2005. ICESat waveform measurements of withinfootprint topographic relief and vegetation vertical structure. Geophysical Research Letters, 32. Helmer, E.H., Lefsky, M.A., 2006. Forest canopy heights in Amazon River Basin Forests as estimated with the Geoscience Laser Altimeter System (GLAS). In: AguirreBravo, C., Pellicane, P.J., Burns, D.P., Draggan, S. (Eds.), Monitoring Science and Technology Symposium: Unifying knowledge for Sustainability in the Western Hemisphere Proceedings RMRS-P-42CD. US Department of Agriculture, Forest Services, Rocky Mountain Research Station, Fort Collins, CO, pp. 802–808. Hiratsuka, M., Toma, T., Ramada, M., Heriansyah, I., Morikawa, Y., 2003. A general allometric equation for estimating biomass in Acacia mangium plantations. In: Proceedings of the International Conference on Tropical Forests and Climate Change: Carbon Sequestration and Clean Development Mechanism, Conference Committee, Traders Hotel, Roxas Blvd, p. xvii. Houghton, R.A., 2005. Aboveground forest biomass and the global carbon balance. Global Change Biology 11, 945–958. Houghton, R.A., 2008. Biomass. In: Sven Erik, J., Brian, F. (Eds.), Encyclopedia of Ecology. Academic Press, Oxford, pp. 448–453. Hubbard, W., Latt, C., Long, A., 1998. Forest Terminology for Multiple-use Management. University of Florida. IPCC, 2003. Good Practice Guidance for Land Use, Land Use Change and Forestry. http://www.ipccnggip.iges.or.jp/public/gpglulucf/gpglulucf.html Available: (accessed 25.08.09). JNCC, 2006. The New Forest SAC, Natura 2000 Standard Data Form. Available: http://www.jncc.gov.uk/ProtectedSites/SACselection/n2kforms/UK0012557.pdf (accessed 09.06.09). Kaufman, Y.J., Tucker, C.J., Fung, I., 1990. Remote sensing of biomass burning in the tropics. Journal of Geophysical Research 95, 9927–9939. Keller, M., Palace, M., Hurtt, G., 2001. Biomass estimation in the Tapajos National Forest Brazil: examination of sampling and allometric uncertainties. Forest Ecology and Management 154, 371–382. Kimes, D.S., Ranson, K.J., Sun, G., Blair, J.B., 2006. Predicting lidar measured forest vertical structure from multi-angle spectral data. Remote Sensing of Environment 100, 503–511. Lefsky, M.A., Cohen, W.B., Parker, G.G., Harding, D.J., 2002. LiDAR remote sensing for ecosystem studies. BioScience 52, 19–30. Lefsky, M.A., Harding, D.J., Keller, M., Cohen, W.B., Carabajal, C.C., Del Bom EspiritoSanto, F., Hunter, M.O., de Oliveira Jr., R., 2005. Estimates of forest canopy height and aboveground biomass using ICESat. Geophysical Research Letters 32, L22S02, http://dx.doi.org/10.1029/2005GL023971.
117
Lefsky, M.A., Keller, M., Pang, Y., De Camargo, P.B., Hunter, M.O., 2007. Revised method for forest canopy height estimation from Geoscience Laser Altimeter System waveforms. Journal of Applied Remote Sensing 1 (1), 013537, http://dx.doi.org/10.1117/1.2795724. Mallet, C., Bretar, F., 2009. Full-waveform topographic lidar: state-of-the-art. ISPRS Journal of Photogrammetry and Remote Sensing 64, 1–16. Markwardt, C.B., 2009. Non-linear least squares fitting in IDL with MPFIT. In: Pacific, A.S.o.t. (Ed.), Astronomical Data Analysis and Systems XVIII. San Francisco, p. 251. Mitchard, E.T.A., Saatchi, S.S., White, L.J.T., Abernethy, K.A., Jeffery, K.J., Lewis, S.L., Collins, M., Lefsky, M.A., Leal, M.E., Woodhouse, I.H., Meir, P., 2012. Mapping tropical forest biomass with radar and spaceborne LiDAR in Lopé National Park Gabon: overcoming problems of high biomass and persistent cloud. Biogeosciences 9, 179–191. Muukkonen, P., Heiskanen, J., 2007. Biomass estimation over a large area based on standwise forest inventory data and ASTER and MODIS satellite data: a possibility to verify carbon inventories. Remote Sensing of Environment 107, 617–624. NASA, 2009. ICESat Home Page. Available: http://icesat.gsfc.nasa.gov/ (accessed 23.01.09). Nelson, R., 2010. Model effects on GLAS-based regional estimates of forest biomass and carbon. International Journal of Remote Sensing 31, 1359–1372. Ni-Meister, W., Jupp, D.L.B., Dubayah, R., 2001. Modelling LiDAR waveforms in heterogeneous and discrete canopies. IEEE Transactions on Geoscience and Remote Sensing 39, 1943–1958. North, P.R.J., Rosette, J.A.B., Suárez, J.C., Los, S.O., 2010. A Monte Carlo radiative transfer model of satellite waveform LiDAR. International Journal of Remote Sensing 31, 1343–1358. NSIDC, 2008. GLAS Altimetry Product Usage Guidance. Available: http://nsidc. org/data/docs/daac/glas altimetry/pdf/NSIDC AltUserGuide Rel29.pdf (accessed 28.09.09). NSIDC, 2009. ICESat/GLAS Data: Tools. Available: http://nsidc.org/data/icesat/ tools.html (accessed 12.08.09). NSIDC, 2009a. GLAS Altimetry Data Dictionary [Online]. Available: http://nsidc. org/data/docs/daac/glas altimetry/data dictionary.html (accessed 29.09.09). Pang, Y., Lefsky, M., Sun, G., Miller, M.E., Li, Z., 2008. Temperate forest height estimation performance using ICESat GLAS data from different observation periods. In: Proc. ISPRS Congress, vol. 37, Beijing, pp. 1682–1750. Parker, G.G., Lefsky, M.A., Harding, D.J., 2001. Light transmittance in forest canopies determined using airborne laser altimetry and in-canopy quantum measurements. Remote Sensing of Environment 76, 298–309. Ranson, K.J., Sun, G., Kovacs, K., Kharuk, V.I., 2004. Landcover attributes from ICESat GLAS data in Central Siberia. In: Geoscience and Remote Sensing Symposium, 2004. IGARSS’04. Proceedings. 2004 IEEE International, pp. 753–756. ˜ D., Meier, E., Allgöwer, B., Chuvieco, E., Ustin, S.L., 2003. Modeling airborne Riano, laser scanning data for the spatial generation of critical forest parameters in fire behavior modeling. Remote Sensing of Environment 86, 177–186. Rosette, J.A.B., North, P.R.J., Suárez, J.C., 2008. Vegetation height estimates for a mixed temperate forest using satellite laser altimetry. International Journal of Remote Sensing 29, 1475–1493. Schutz, B.E., Zwally, H.J., Shuman, C.A., Hancock, D., DiMarzio, J.P., 2005. Overview of the ICESat mission. Geophysical Research Letters 32, L21S01. 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, 23–45. Sun, G., Ranson, K.J., 2000. Modeling lidar returns from forest canopies. IEEE Transactions on Geoscience and Remote Sensing 38, 2617–2626. Sun, G., Ranson, K.J., Kimes, D.S., Blair, J.B., Kovacs, K., 2008. Forest vertical structure from GLAS: an evaluation using LVIS and SRTM data. Remote Sensing of Environment 12, 107–117. Takeda, H., 2004. Ground surface estimation in dense forest. The International Archives of the Photogrammetry. Remote Sensing and Spatial Information Science 35, 1016–1023. Todd, K.W., Csillag, F., Atkinson, P.M., 2003. Three dimensional mapping of light transmittance and foliage distribution using LiDAR. Canadian Journal of Remote Sensing 29, 544–555. UNFCCC, 1997. Kyoto Protocol to the United Nations Framework Convention on Climate Change. Available: http://unfccc.int/resource/docs/convkp/kpeng.pdf (accessed 24.08.09). UNFCCC, 2009a. Copenhagen Accord. Available: http://unfccc.int/ files/meetings/cop 15/application/pdf/cop15 cph auv.pdf (accessed 08.01.10). UNFCCC, 2009b. Fact Sheet: Reducing Emissions from Deforestation in Developing Countries: Approaches to Stimulate Action. Available: http://unfccc.int/files/ press/backgrounders/application/pdf/fact sheet reducing emissions from deforestation.pdf (accessed 24.08.09). 2009. United States Geological Survey [Online]. Available: URL1, http://glovis.usgs.gov/ (accessed 14.08.09). 2009. ICESat/GLAS Data Subsetter [Online]. Available: URL2, http://nsidc.org/forms/glas subset form.html (accessed 11.08.09). Wagner, W., Ullrich, A., Ducic, V., Melzer, T., Studnicka, N., 2006. Gaussian decomposition and calibration of a novel small-footprint full-waveform digitising airborne laser scanner. ISPRS Journal of Photogrammetry and Remote Sensing 60, 100–112. Waring, R.H., Way, J., Hunt, E.R.J., Morrissey, L., Ranson, K.J., Weishampel, J.F., Oren, R., Franklin, S.E., 1995. Imaging radar for ecosystem studies. BioScience 45, 715–723.
118
I.A. Iqbal et al. / International Journal of Applied Earth Observation and Geoinformation 23 (2013) 109–118
Xing, Y., de Gier, A., Zhang, J., Wang, L., 2010. An improved method for estimating forest canopy height using ICESat-GLAS full waveform data over sloping terrain: a case study in Changbai mountains, China. International Journal of Applied Earth Observation and Geoinformation 12 (5), 385–392. Zadiraka, V.K., Igisinov, K., 1973. Analysis of the accuracy and efficiency of the Fast Fourier Transform algorithm and some of its applications. Cybernetics and Systems Analysis 9, 973–977.
Zwally, H.J., Schutz, B., Abdalati, W., Abshire, J., Bentley, C., Brenner, A., Bufton, J., Dezio, J., Hancock, D., Harding, D., Herring, T., Minster, B., Quinn, K., Palm, S., Spinhirne, J., Thomas, R., 2002. ICESat’s laser measurements of polar ice, atmosphere, ocean, and land. Journal of Geodynamics 34, 405–445. Zwally, H.J., Schutz, R., Bentley, C., Bufton, J., Herring, T., Minster, J., Spinhirne, J., Thomas, R., 2006. GLAS/ICESat L1&L2 Global Land Surface Altimetry Data, Global Land Elevation Data V028/V029, 2003-2006. National Snow and Ice Data Center. GLAS/ICESat Data Subsetter, Boulder, CO.
Contents lists available at SciVerse ScienceDirect
International Journal of Applied Earth Observation and Geoinformation journal homepage: www.elsevier.com/locate/jag
A novel approach to estimate canopy height using ICESat/GLAS data: A case study in the New Forest National Park, UK Irfan Akhtar Iqbal a , Jadunandan Dash a , Saleem Ullah b,∗ , Ghayyas Ahmad c a b c
School of Geography, University of Southampton, Highfield, Southampton SO17 1BJ, UK Faculty of Geo-Information Science and Earth Observation (ITC), University of Twente, 7500 AA Enschede, The Netherlands Pakistan Forest Institute, Peshawar, Khyber Pakhtunkhwa, Pakistan
a r t i c l e
i n f o
Article history: Received 21 August 2012 Accepted 18 December 2012 Keywords: Canopy height LiDAR Remote sensing
a b s t r a c t The Geoscience Laser Altimeter System (GLAS) aboard Ice, Cloud and land Elevation Satellite (ICESat) is a spaceborne LiDAR sensor. It is the first LiDAR instrument which can digitize the backscattered waveform and offer near global coverage. Among others, scientific objectives of the mission include precise measurement of vegetation canopy heights. Existing approaches of waveform processing for canopy height estimation suggest Gaussian decomposition of the waveform which has the limitation to properly characterize significant peaks and results in discrepant information. Moreover, in most cases, Digital Terrain Models (DTMs) are required for canopy height estimation. This paper presents a new automated method of GLAS waveform processing for extracting vegetation canopy height in the absence of a DTM. Canopy heights retrieved from GLAS waveforms were validated with field measured heights. The newly proposed method was able to explain 79% of variation in canopy heights with an RMSE of 3.18 m, in the study area. The unexplained variation in canopy heights retrieved from GLAS data can be due to errors introduced by footprint eccentricity, decay of energy between emitted and received signals, uncertainty in the field measurements and limited number of sampled footprints. Results achieved with the newly proposed method were encouraging and demonstrated its potential of processing full-waveform LiDAR data for estimating forest canopy height. The study also had implications on future full-waveform spaceborne missions and their utility in vegetation studies. © 2012 Elsevier B.V. All rights reserved.
1. Introduction Geoscience Laser Altimeter System (GLAS), the first of a new generation spaceborne LiDAR onboard the Ice, Cloud and land Elevation Satellite (ICESat) mission, was launched on January 13, 2003. GLAS is the first spaceborne instrument which can fully digitize the backscattered waveform from incident surfaces (Wagner et al., 2006). GLAS offers near-global coverage by carrying three laser instruments (Laser 1, 2 and 3) with only one operational at a time. After the failure of Laser 1 in March 2003 (Schutz et al., 2005), GLAS followed a modified plan of 33 days sub-cycle with a 91 days repeat orbit (Abshire et al., 2005). Under the revised plan, Laser 3 acquired data during February–March, May–June, and October–November each year (Schutz et al., 2005), until its failure in October, 2008. Laser 2 also ceased emitting light on October 11, 2009 (NASA, 2009).
∗ Corresponding author at: Faculty of Geo-Information Science and Earth Observation (ITC), University of Twente, P.O. Box 217, 7500 AE Enschede, The Netherlands. Tel.: +31 53 4874 372. E-mail addresses: [email protected] (I.A. Iqbal), [email protected] (S. Ullah). 0303-2434/$ – see front matter © 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.jag.2012.12.009
After seven years of near-global coverage and 15 laser-operation campaigns, ICESat has stopped collecting science data. However, data from these historical acquisitions can be exploited to evaluate the potential of these data in other terrestrial applications and thus, can have useful implications for the next generation of sensor development, for example ICESat-2, DESDynI. Originally, ICESat was launched to monitor polar ice-sheet dynamics and its role in changing the global sea-level (Zwally et al., 2002). However, other scientific objectives of this mission also include surface reflectivity, precise measurement of land topography and vegetation canopy heights. Light Amplification by Stimulated Emission of Radiation (LASER) altimetry has supplemented the existing remote sensing techniques because of its high precision and capability to estimate canopy height. Light Detection and Ranging (LiDAR) directly measures the vertical component of vegetation, thus has the potential to measure the structural attributes of vegetation (Dubayah et al., 2000; Lefsky et al., 2002; Nelson, 2010). It has the capability to provide information at both canopy and subcanopy levels which is crucial for understanding forest health and regeneration capacity (Stone et al., 2000; Todd et al., 2003), forest fires (Andersen ˜ et al., 2003) and present and future carbon stock et al., 2005; Riano
110
I.A. Iqbal et al. / International Journal of Applied Earth Observation and Geoinformation 23 (2013) 109–118
estimation (Dean et al., 2004). Moreover, tree heights and vegetation density derived from LiDARs prove to be important inputs for aboveground biomass estimation, understanding global carbon cycle and ecosystem dynamics (Ranson et al., 2004). Forests, being important for their role in the global climate (UNFCCC, 1997, 2009a), need to be quantitatively as well as qualitatively analyzed. By estimating the structural attributes (e.g. aboveground biomass) of these ecosystems, the amount of carbon stored can be calculated as half the dry weight of biomass (Drake et al., 2002; Houghton, 2008; IPCC, 2003; UNFCCC, 2009b). Conventionally, aboveground biomass is estimated by using the DBH (Keller et al., 2001), though straightforward, estimation of aboveground biomass (using DBH) on ground is time consuming as well as expensive (Houghton, 2005) and mostly involves destructive sampling (Hiratsuka et al., 2003). Alternatively, due to its spatial and temporal coverage, remote sensing provides opportunity to derive vegetation biophysical variables and some of them are related to vegetation structure. Although DBH cannot be directly measured from remotely sensed data; other biophysical parameters like canopy height and canopy cover can serve as a proxy for estimating the amount of aboveground biomass. Several studies have attempted to use remotely sensed data (both passive and active (RADAR)) for estimating the biophysical parameters of vegetation and in turn, to estimate aboveground biomass e.g. (Dobson et al., 1992; Kaufman et al., 1990; Muukkonen and Heiskanen, 2007). However, the accuracy and sensitivity of these data types drop-off in areas of closed canopies and high-biomass (Waring et al., 1995; Mitchard et al., 2012).
2. Application of GLAS data in vegetation studies GLAS transmits short pulses (4 ns width equivalent to 60 cm in surface elevation) of infrared (1064 nm) and green (532 nm) light at a frequency of 40 Hz. The nominal diameter of the illuminated spot, the footprint, is ∼70 m, space at 172 m along the track; however, its size and ellipticity have varied through the time (Pang et al., 2008). The reflected photons are collected by a telescope of 1 m diameter, and their intensities are binned every nanosecond. Height (or range) can be worked out as the product of speed of light (3 × 108 m/s) and one-way travel time of the photons. Intensities of the reflected photons are stored in 544 bins over ice sheet and land, corresponding to a height of 81.6 m (Brenner et al., 2003). In steeply sloped areas and/or areas where feature heights exceed 81.6 m, GLAS waveform would truncate, making it inconvenient to derive range information. For this reason, in later operations height extent was increased to 150 m over land, using a ‘waveform compression scheme’ (Harding and Carabajal, 2005) according to which the first 152 bins correspond to 60 cm each. Fig. 1 shows the waveform components, used to estimate various physical attributes, which will be referred to in the succeeding text. The numbers 25, 50, 75 and 100 indicate height quartiles where respective percentage of waveform intensity is concentrated. GLAS products, typically GLA01 and GLA14, have been used for vegetation classification and biomass estimation (Boudreau et al., 2008; Lefsky et al., 2005; Ranson et al., 2004), and seasonal changes in vegetation (Duong et al., 2008), where most of the information is retrieved from GLA01 (the ‘raw’ waveform). Particular to vegetation studies, Ranson et al. (2004) found the front slope angle (Fig. 1) strongly correlated to canopy density and vertical variability of upper canopy. These relationships were further exploited in a classification of forests in central Siberia. Aboveground biomass estimation from GLAS data has been demonstrated to have agreement with field based or airborne measurements (Boudreau et al., 2008; Helmer and Lefsky, 2006; Lefsky
Fig. 1. Illustration of waveform components used to estimate canopy height. The waveform has been smoothed with a moving average of 5 to reduce noise and favour clarity. The numbers 25, 50, 75 and 100 indicate height quartiles where respective percentage of waveform intensity is concentrated. The distance between the two dotted blue lines is the waveform extent. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of the article.)
et al., 2005). Extraction of canopy height from raw waveform is complicated in areas of varied topography due to slope effects (Harding and Carabajal, 2005). GLAS returns from flat surfaces can be represented by a Gaussian curve; while flat areas with homogeneous vegetation can be approximated by two Gaussian curves. With increasing topographic variation and vegetation heterogeneity, the number of Gaussian peaks increases (Duncanson et al., 2010). To allow for such complexities, methods for canopy height estimation have been developed which are generally classified (Chen, 2010b) as (i) statistical methods and (ii) direct methods. In statistical methods, waveform processing is supplemented by additional information from the underlying terrain. So-called terrain indices are extracted either from the available Digital Terrain Models (DTMs) (Lefsky et al., 2005; Rosette et al., 2008) or derived from the waveforms (Lefsky et al., 2007; Pang et al., 2008). Direct methods concentrate on waveform processing for accurate identification of canopy top and the ground return with less, if any, assistance from auxiliary data, e.g. DTMs. Lefsky et al. (2005) derived canopy height as a linear function of waveform extent and terrain index. Waveform extent was taken as the difference between the first and last threshold crossings of the waveform; whereas terrain index was calculated, from Shuttle RADAR Topography Mission (SRTM) data, as the difference between minimum and maximum elevations within a footprint. Having both waveform extent and terrain index identified, canopy height was modelled as given in the equation: H = b0 (w − b1 g)
(1)
where, H is the measured canopy height, w is the waveform extent in m, g is the terrain index in m, b0 and b1 are coefficients applied to the waveform and terrain index, respectively. With their method, Lefsky et al. (2005) were able to explain a maximum variation of 68% with an RMSE of 4.85 m in canopy heights. They attributed the unexplained variance in canopy height estimation to the ‘measurement’ of waveform extent, more importantly, the identification of ground return (or peak) in the raw waveform. For efficient ground peak detection, Gaussian components are fitted to the raw waveform (Brenner et al., 2003) and iteratively reduced to six (Sun et al., 2008), using the Least Squares approach (Markwardt, 2009). Harding and Carabajal (2005) showed that out of the six peaks, the last peak with the highest amplitude is a representative of the ground. However, in densely vegetated areas, due to energy attenuation through the canopy (Parker et al., 2001),
I.A. Iqbal et al. / International Journal of Applied Earth Observation and Geoinformation 23 (2013) 109–118
the ground peak may have lower amplitude than adjacent peaks (Boudreau et al., 2008). In such cases, Boudreau et al. (2008) suggested to consider the peak with highest amplitude as the ground peak, even if it is not the last one. However, the method will lead to inaccuracies in estimating ground and therefore in canopy height. Rosette et al. (2008) reported an underestimation of 1.20 m in canopy height by considering the first or second highest amplitude peak as ground return. A study by Sun et al. (2008) demonstrated the utility of height quartiles (Fig. 1) to predict canopy heights. Quartiles at 25, 50, 75 and 100 were calculated on the basis of percentage of reflected energy. Quartile at 25, for example, represented the aboveground level at which 25% of the reflected energy is accumulated. Similarly, quartile at 100 is the distance between signal start and the ground return, which is essentially the waveform extent. Precise detection of the ground peak is important because it is used as a reference for calculating quartiles. Sun et al. (2008) used an automated approach for ground peak detection in which bin values (or amplitudes), from signal end to signal start, were compared to the two neighbouring values. The first peak, from signal end, with a value more than the half width of the transmitted laser pulse was taken as the ground peak. Their method explained a maximum variability of 83% in canopy height. However, signal strength was not considered which may label the first ‘significant peak’ as ground even if it has a lower energy. Ground or flat surfaces tend to register stronger responses. Another automated approach for peak identification with a search window of 5 ns was proposed by Duong et al. (2008). The search window operated on the normalized (amplitude at each instant divided by the sum of all amplitudes) and smoothed waveform, moving from beginning through end at a step of 1 ns. A peak was identified if the central bin had an incrementally higher amplitude value than the two bins on either side. Further, six Gaussian peaks were fitted to the resulting waveform, the last being representative of the ground return. Some irregular waveforms could not be adequately decomposed by Gaussian distributions and often missed a mode (Duncanson et al., 2010; Duong et al., 2009). The issues (of missed peaks) significantly limits the accuracy of estimated terrain and canopy surfaces. The last mode is critical because, principally, canopy height is the distance defined by the time difference between the signal start after noise and the centroid of the last peak (Sun et al., 2008). Although Gaussian fitting has shown considerable promise in various studies, its efficiency is ‘less satisfactory’ for pulses with high amplitudes (Wagner et al., 2006). Furthermore, in cases of low amplitude pulses, the extracted parameters are less accurate which demands stringent consideration of the waveform processing algorithms (Mallet and Bretar, 2009). There is a need to develop a robust waveform fitting method which circumvents the stated problems, allows a more accurate retrieval of canopy height and reduces uncertainties in derived variables such as aboveground biomass. Addressing the need, in this study, first a new method involving Fourier analysis of waveform processing for quantification of canopy height from GLAS data was developed and then the canopy heights derived from GLAS data were validated using field data.
3. Study area The study was selected within the New Forest National Park in Southeast England (50◦ 51 59 latitude and 01◦ 40 50 longitude) covering an area of 29262.36 ha. This site holds the status of “Site of Community Importance” since December 2004 (JNCC, 2006). Also, it was designated as a Special Area of Conservation in April 2005. The area is covered with deciduous broadleaved (oak/beach 29%) and coniferous woodland (17%). The remaining area is a mix of bogs,
111
marshes, scrub and heath. It is under the active management of the Forestry Commission of Great Britain. The study site is subdivided in compartments; a compartment is a managerial unit which has similar characteristics within the area; such as, soil types, vegetation, productivity and topography (Hubbard et al., 1998). Most of the compartments have pure and even-aged stands; however, in some compartments age-diversity does exist. Fig. 2 shows the GLAS data projected over the New Forest National Park. 4. Data 4.1. GLAS/ICESat data There are 15 GLAS data products available (Zwally et al., 2006). Data products GLA01 and GLA14 of release 28 and 29 were used in this study. Data were procured from the National Snow and Ice Data centre (NSIDC) through ICESat/GLAS data subsetter (URL2, 2009). GLA01 is a level 1A global altimetry product which contains intensities of transmitted and received waveforms. GLAS digitizes these intensities as counts which are then converted to volts using calibration tables in GLA01 header records. The transmitted pulses were arranged in time order, which means that the value of the first sample is for the sample closest to the spacecraft in time. Conversely, the received echo is in time reversed order, which means the value of the first sample is for the sample farthest from the spacecraft in time (NSIDC, 2008). GLA14 is a level 2 region-specific land surface altimetry data which is derived from level 1 products (GLA05 and GLA06). Elevations are given in m with respect to a reference ellipsoid accompanied by geoid elevations which allow users to calculate orthometric heights. Similar to GLA01, GLA14 also contains geolocation of footprint centroid, unique record and shot numbers, and a time stamp with Universal Coordinated Time (UTC) in J2000 (referenced from noon on January 1, 2000). As the transmitted and received pulses are assumed to resemble several Gaussian curves (Brenner et al., 2003), a set of Gaussian pulses (up to six) is fitted to the waveform during processing. For this purpose, GLA14 variables also include amplitudes, areas and standard deviations of 6 Gaussian distributions if the user is interested in Gaussian fitting. 4.2. Ancillary data In addition to GLAS data, CORINE land cover (CLC) classification 2000 (EEA, 2009) was used to identify forest patches in the study area and extract GLAS footprints over them. CLC 2000 is a geographic database which covers European countries and gives information on land cover classes at a spatial resolution of 100 m. Cloud free Landsat 7 (ETM+) imagery (URL1, 2009), dated April 2003, was used to visualize GLAS transects. A compartment and sub-compartment database of the study area was procured from the Forestry Commission of Great Britain (FCGB, 2009). For each compartment and sub-compartment in the study area the data-base provided information on: planting year, area of the compartment, type of species and yield class. The Forestry Commission has also designed yield tables to aid forest management in British conditions (Hamilton and Christie, 1971). From the species-specific yield tables and general yield class curves, crop height can be estimated if yield class and planting year of a compartment are known. 5. Methodology In the current study, a new method based on Fourier transforms, was used to process the GLAS waveform. Fourier transform converts a waveform in a time domain to a waveform in the
112
I.A. Iqbal et al. / International Journal of Applied Earth Observation and Geoinformation 23 (2013) 109–118
Fig. 2. Representation of ICESat tracts projected over the study area. Eastern tracks (in red box) were used for canopy height estimation, comprising of 145 footprints. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of the article.)
frequency domain with a real and an imaginary part (cosine and sine) (Bloomfield, 2000).
5.1. Data pre-processing The datasets used in this study were registered to different reference systems. To have the datasets coherent, they were all re-projected to a common reference system: the World Geodetic System 1984 (WGS1984). Information contained in GLA01 and GLA14 were extracted, using the IDL codes provided by the National Snow and Ice Data Centre (NSIDC, 2009). GLA01 waveforms were then linked to GLA14 on the basis of record and shot numbers. All GLAS footprints over the study area were identified and waveforms with high level of echo saturation were filtered out using a flag called “i gval rcv” (NSIDC, 2009a). This flag ranges from 0 to 250, with some values labelled as −999 indicating distorted waveforms. All waveforms having an i gval rcv of −999 were discarded. Further, CORINE land cover classification was used to identify forest patches in the study area and a total 145 GLAS footprints over these patches were extracted for analysis.
5.2. Field data collection Out of the total 145 footprints, 50 were randomly selected for field data collection to avoid bias and obtain representative samples (Lefsky et al., 2005). From these 50 identified footprints, 23 fell in the management boundaries of the New Forest National Park for which the Forestry Commission yield data was readily available. For the remaining 27 footprints, field work was conducted, in October 2009, to obtain height information for 19 points; remaining footprints could not be sampled due to logistic difficulties (e.g., limited time, weather conditions and refused access). A handheld Garmin eTrex GPS receiver was used, with an accuracy of about 7 m, to locate the footprint centres. Within each footprint, angle of elevation (top) and depression (base) of trees greater than 10 cm DBH were recorded with a clinometer. Also, horizontal distance from trees, ground elevations, species and description of footprint were recorded. Data collection was flexible in terms of the sample size; however, care was taken to characterize height variability. In compartments with homogeneous cover, fewer trees were measured than in uneven-aged plots. In uneven-aged areas, within footprints,
I.A. Iqbal et al. / International Journal of Applied Earth Observation and Geoinformation 23 (2013) 109–118
113
Fig. 3. Conceptual layout of sample plot for measuring canopy heights in the field.
data was collected in 4 annular plots, one in centre and three on the footprint margin at approximate azimuths of 0◦ , 120◦ and 240◦ . Radius of central plot was 10 m while the radii of smaller marginal plots were 5 m each. Radius of the central plot was greater than the marginal ones, to account for the reflected LiDAR energy which has a central maximum and fades away radially outwards (Duong et al., 2006b). Fig. 3 gives the conceptual sampling design for a footprint.
Fig. 4. Using a defined threshold level to discriminate between noise, signal-start and signal-end (as in Eq. (2)).
In this study, the first 200 bins were used to calculate the threshold parameters as shown below: Mean noise = E¯ =
5.3. Waveform processing 5.3.1. Waveform normalization Since different waveforms had different amplitude intensities due to changing atmospheric conditions and/or decay in the laser device (Duong et al., 2008), it would be improper to apply a single threshold value to all waveforms. Variation in amplitude intensities complicates comparison of waveforms from different campaigns and hence requires normalization, which was achieved by dividing the received energy at each instance (Ei ) by the total received 544 E (Duong et al., 2006a). energy (ET), where total energy is i=1 i
544
Mathematically, a normalized waveform is expressed as Ei /
E. i=1 i
200 i=1
Ei /200
Noise standard deviation =
(3)
200 (E − E) ¯ 2 i i=1
200
(4)
where, Ei is the received energy at time i. 5.3.4. Detecting ground, canopy-top and canopy height In order to estimate canopy height from the smoothed waveform, first, all peaks in the smoothed waveform were identified by a search window (1 ns) which moved from start to end of the waveform. The window calculated the difference between two succeeding amplitudes, and stopped until it found the point where the first derivative changed from positive to negative or vice versa.
5.3.2. Waveform smoothing and fitting The normalization step was followed by smoothing the normalized waveform and fitting a modelled waveform to it. For this purpose, Fast Fourier Transform was used, which is an efficient algorithm (Zadiraka and Igisinov, 1973) to compute discrete Fourier transform (analysis) and its inverse (synthesis). Power Spectral Density (PSD) was applied to find the frequency region of maximum signal strength. Based on the PSD, we used 30 harmonics for smoothing of and model-fitting to the normalized waveform.
5.3.3. Threshold and background noise The GLAS echo, telemetered to 544 bins, corresponds to a height of 81.6 m over land. Out of 544, the first 200 bins (30 m) were used to define the noise level (using Eq. (3)). The remaining 344 bins (51.6 m) were sufficient to represent the vegetation structure without truncation. A noise threshold level was used to single out the signal from noise and identify locations of signal-start and signal-end. Start and end of the signal are threshold crossings of the leading and trailing edges (Harding and Carabajal, 2005). We set the threshold level to the sum of mean noise and 4 times the noise standard deviation as in Lefsky et al. (2005) using Eq. (2). The resulting waveform is shown in Fig. 4. Threshold level = mean noise + 4 × noise standard deviation
(2)
Fig. 5. Identified peaks (in red circles) and troughs (in green circles) on the fitted waveform. Peak with maximum amplitude is also marked for determining ground, canopy top and canopy height. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of the article.)
114
I.A. Iqbal et al. / International Journal of Applied Earth Observation and Geoinformation 23 (2013) 109–118
Fig. 6. Final parameterization of GLAS waveform. Red and green circles show the waveform peaks and troughs, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of the article.)
Fig. 5 illustrates all peaks (in red circles), troughs (in green circles) and peak maximum amplitude. To detect canopy top, ground and canopy height, the peak with the maximum amplitude was used as a reference. From this reference point, the search window moved towards the signal-end and searched for a ground return peak. Peak with amplitude onefifth the maximum amplitude, above the defined threshold, was defined as the ground return. It was assumed, that due to reduced penetration and decay of energy through the canopy (Parker et al., 2001), the return energy can be as low as one-fifth of the maximum amplitude. If these conditions were not met, the maximum amplitude peak was designated as the ground return. To infer the position of the canopy, the window moved from the maximum amplitude peak towards the signal start to search for canopy point. A canopy point is defined here as the first peak after the signal start which is above the threshold level. In this study, we applied a threshold of 0.003 V, using Eq. (2). The leading edge of a canopy point indicates the canopy-top, because as soon as the signal interacts with the uppermost foliage and/or branches, it is excited. Canopy height was calculated as the distance between ground return and canopy top (Fig. 6). In case the ground return and the canopy top peaks coincided, the canopy height was assumed to be zero. 5.4. Field data analysis Angles of elevation and depression of tree top and bottom (respectively), and distance to tree from the point of observation were used to calculate tree heights (Eq. (5)). Tree height (field) = d × tan + d × tan ˚
(5)
where, d is the distance between the tree and the point of observation (in m), is the angle of elevation of the tree top in degrees, and ˚ is the angle of depression of the tree bottom in degrees. Because the waveform is sensitized by the uppermost foliage and/or branches, it represents the maximum canopy height inside the footprint. For this reason, only field heights of dominant trees were used for analysis. Moreover, there was a time lag of approximately 3–5 years between GLAS acquired data and field collected data, during which height increments were expected. Annual
Fig. 7. Correlation of canopy heights obtained from waveform analysis and those collected from field and the Forestry Commission data. Outliers are shown in red. Red dotted line shows 1:1 line. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of the article.)
height increments for each species were worked out from yield tables and subtracted from field calculated canopy heights, to get adjusted canopy heights. Canopy heights for the footprints falling in Forestry Commission sub-compartments were estimated from yield class curves in Hamilton and Christie (1971). Crop age was estimated as the difference between planting year and GLAS data acquisition year, and was used to estimate canopy height from yield class curves. 6. Results and discussion Two tailed Pearson’s correlation in Statistical Package for Social Scientists (SPSS v.16) was used to compare any two datasets resulting from the adopted methodology. R-value and RMSE were used as indicators of ‘goodness of agreement’ between corresponding variables. 6.1. Accuracy of GLAS derived canopy heights Canopy heights calculated from field and Forestry Commission data correlated with those obtained from GLAS waveforms showed a correlation of 0.59, with an RMSE of 5.18 m (Fig. 7). Heights derived from GLAS waveforms are hereinafter referred to as ‘GLAS estimated height’. The offset of about 6 m (in Fig. 7) represents positive bias in GLAS estimated canopy heights. Also, the data has a wider spread around the correlation line. The low correlation between the two datasets was traceable to six outliers, as shown in Fig. 7 (in red). In these six cases, the over and under estimation of canopy heights reached 11–13 m, respectively. Analysis revealed that over mixed and/or multiple canopies, the GLAS estimated heights were positively biased. In cases of thick vegetation cover, the heights were under-estimated. Biases in canopy height (both positive and negative) are also induced by the uncertainty in ground detection. For example, litter on forest floor is a possible reason of multiple scattering (Brenner et al., 2003; Mallet and Bretar, 2009). The extent of over and under estimations due to litter depends on its thickness, three-dimensional distribution and spectral properties (Kimes et al., 2006). Soil moisture may also be a factor to cause decay in the return pulse by absorbing some of the incoming energy (Brenner et al., 2003). However further studies are needed
I.A. Iqbal et al. / International Journal of Applied Earth Observation and Geoinformation 23 (2013) 109–118
Fig. 8. (a) Extreme case of over-estimation up to 11 m. Red and green circles show waveform peaks and troughs, respectively and (b) Footprint location of waveform (in (a)). The red dot shows the footprint centre, labelled with the record and shot number. Image taken from Google Earth (NOAA, 2009). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of the article.)
115
Fig. 9. (a) Extreme case of under-estimation up to 13 m. Red and green circles show waveform peaks and troughs, respectively and (b) location of waveform (in (a)). The footprint, centred at the red dot, covered young Douglas Fir crop. Image taken from Google Earth (NOAA, 2009). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of the article.)
to quantify the effect of soil moisture on LiDAR pulse absorption. After multiple scattering, a strong signal will under-estimate the ground; whereas a weaker signal will cause over-estimation due to increased path length (Brenner et al., 2003). Such uncertainties is ground detection are reflected in canopy height estimation.
dominant canopies have influenced the returned energy and resulted in slight rise at the signal start (Fig. 8a), which is introducing a potential source of error, thus failing to represent all vegetation within the footprint. Fig. 8b gives an aerial view of the location represented by waveform in Fig. 8a. The red circle approximates the GLAS footprint.
6.1.1. Extreme over-estimation Fig. 8a represents the waveform of a footprint which covered a sparse patch of oaks (196 years of age) dominated by low canopies. In the study area, oaks are not managed commercially; hence, height variations are expected. The LiDAR signal was excited by foliage of a few dominant trees, resulting in slight amplitude rise (red arrows in Fig. 8a); however, a strong response was registered from the trees in the under storey. For similar conditions, waveforms behaved the same way, resulting in over-estimated canopy heights. Vegetation structural heterogeneity is documented to complicate the relationship between LiDAR derived and field observed canopy heights (Clark et al., 2004). The incident LiDAR beam is directly affected by vertical (and horizontal) distribution of foliage and woody elements among other variables such as density and size of the canopy (Ni-Meister et al., 2001). LiDAR waveform structure, in particular, has been noted to be sensitive to canopy shape, reflective properties and associated photon interactions (North et al., 2010). A small number of
6.1.2. Extreme under-estimation Fig. 9a is a representative waveform where high underestimation (13 m) in canopy height determination was observed. This footprint covered a young crop (29 years of age) of Douglas Fir with dense canopy, because of which the signal was obstructed (Takeda, 2004). The signal, that penetrated the canopy, was weak and could hardly register a high amplitude peak. Another factor could be soil moisture, which is expected to be more under dense canopy cover, causing signal decay. Such conditions introduce uncertainties in precise ground detection, as discussed in Section 6.1, and ultimately result in biased canopy height estimates. Fig. 9b gives an aerial view of the location represented by waveform in Fig. 9a. The area is covered with young crop of Douglas Fir (29 years of age). The red circle approximates the footprint size. Another important factor, to which GLAS data is reportedly sensitive, is crop age (Helmer and Lefsky, 2006; Sun and Ranson, 2000) which might have introduced bias in GLAS estimated canopy heights. Note that Figs. 8a and 9a are representatives of footprints
116
I.A. Iqbal et al. / International Journal of Applied Earth Observation and Geoinformation 23 (2013) 109–118
method for retrieving ground elevation from Gaussian peaks (Chen, 2010a; Duong et al., 2008; Lefsky et al., 2005; Rosette et al., 2008). To date, Fourier method has not been applied to characterize LiDAR waveforms for canopy height estimation. Although, topographic influences cannot be overlooked, our method of waveform processing gives encouraging results in the study site. This implies the advantage of Fourier method of waveform analysis to automatically estimate canopy heights, in the absence of DTM, and less computational requirements. 7. Conclusion
Fig. 10. Correlation of field and GLAS estimated heights after removal of all six outliers.
covered with old oaks (196 years) and young Douglas Fir (29 years), respectively. It is interesting to note that four out of the six outliers fall within the Forestry Commission compartments, for which canopy heights were estimated from the yield class curves. Only by removing these four outliers from analysis, the correlation value improved from 0.59 to 0.72 and RMSE lowered from 5.18 m to 4.07 m. Yield tables are species-specific generalized models created for compartments and sub-compartments, covering larger spatial extent than covered by GLAS footprint. There is a possibility that variation in site qualities exist on (sub)compartment level, which affects tree growth. Furthermore, the Forestry Commission yield models are not dynamic and depict a mean trend on (sub)compartment level (Rosette et al., 2008). This seems to be a valid reason for biases in canopy heights extracted from yield class curves and those estimated from GLAS waveforms. Moreover, this analysis demonstrates robustness of the new method of waveform processing for canopy height estimation, showing ‘good’ agreement with field measured heights. Due to the uncertainties explained in Sections 6.1.1 and 6.1.2, these waveforms (six in total: 4 from Forestry Commission and 2 from field data) were removed to see the effectiveness of the method. Removal of these outliers further improved the coefficient of correlation from 0.59 to 0.79, and lowered the RMSE from 5.18 m to 3.18 m as shown in Fig. 10. With an automated method of ground identification based on iterative Gaussian fitting, Rosette et al. (2008) were able to explain 74% of the variance in canopy heights in Forest of Dean, UK. It is logical to compare results from Rosette et al. (2008) with this study, because the two study sites (Forest of Dean and the New Forest National Park) have similar environmental conditions and are managed by the Forestry Commission of Great Britain. It is important to mention that using DTM corrections, Rosette et al. (2008) were able to explain 89% of variability in canopy heights. Likewise, recently (Xing et al., 2010) came up with a modification to Lefsky et al. (2005) exploiting the non-linear relationship between waveform extent and maximum canopy height. Xing et al. (2010) were able to explain a maximum variation of 92% in canopy height in the complex terrain (0◦ –30◦ slope) of Wangqing forest, China. The method applied in this study, however, does not require any supplementary information from auxiliary data, i.e. the topographic correction. So far, Gaussian decomposition has been used for ground peak identification; even so, there is no widely accepted
Aim of the study was to evaluate the potential of ICESat/GLAS to estimate canopy height in the New Forest National Park, Southeast England. Canopy height from GLAS altimetry product (GLA01) was estimated and compared with field measured canopy height. Inaccuracies in ground peak detection have been reported to introduce uncertainties in canopy height estimation. This study demonstrated the ability of Fourier Transformation to tackle such issues. Moreover, an automated method for detecting major peaks in the ‘raw’ waveform was developed. Due to signal attenuation, it was assumed that the actual ground return can be the peak with energy 20% lower than the peak of maximum amplitude. Effects of forest litter and soil moisture could be possible reasons for signal weakening. Our presented method explained 79% of the variation in canopy height estimation compared to field data. Other factors which were not investigated in this study; like footprint eccentricity, representation of field measured height in GLAS data and energy decay between transmitted and received pulses may justify the unexplained variance. Future work will investigate the effect of footprint eccentricity, slope and energy decay between transmitted and receive pulse on retrieval of canopy height from GLAS data. The study succeeded to achieve encouraging results with minimum GLAS parameters. This implies that application of other ancillary variables and correction flags may enhance the results. Moreover, due to the possibility of losing some ‘significant’ peaks, the existing assumption of fitting six Gaussian modes to the raw waveform may be relaxed to include more Gaussian peaks. We strongly recommend application of the presented methodology in other similar studies to validate its repeatability and transferability. Acknowledgements The authors wish to thank the European Commission for awarding us with Erasmus Mundus – Scholarships. The authors also acknowledge the help of Dr. Zornitza Ivanova Yovcheva and Mr. M. Sohail Khan for editing the initial drafts and giving a critical feedback. References Abshire, J.B., Sun, X., Riris, H., Sirota, J.M., McGarry, J.F., Palm, S., Yi, D., Liiva, P., 2005. Geoscience Laser Altimeter System (GLAS) on the ICESat Mission: onorbit measurement performance. Geophysical Research Letters 32, L21S02, http://dx.doi.org/10.1029/2005GL024028. Andersen, H.-E., McGaughey, R.J., Reutebuch, S.E., 2005. Estimating forest canopy fuel parameters using LIDAR data. Remote Sensing of Environment 94, 441–449. Bloomfield, P., 2000. Fourier Analysis of Time Series: An Introduction (Wiley Series in Probability and Statistics). Wiley-Interscience, 78 pp. Boudreau, J., Nelson, R.F., Margolis, H.A., Beaudoin, A., Guindon, L., Kimes, D.S., 2008. Regional aboveground forest biomass using airborne and spaceborne LiDAR in Québec. Remote Sensing of Environment 112, 3876–3890. Brenner, A.C., Zwally, H.J., Bentley, C.R., Csath¢, B.M., Harding, D.J., Hofton, M.A., Minster, J.B., Roberts, L., Saba, J.L., Thomas, R.H., Yi, D., 2003. Derivation of range and range distributions from laser pulse waveform analysis for surface elevations, roughness, slope, and vegetation heights. In: Algorithm Theoretical Basis Document, Version 4.1. Goddard Space Flight Center, Greenbelt, MD, USA, p. 92, http://www.csr.utexas.edu/glas/atbd.html
I.A. Iqbal et al. / International Journal of Applied Earth Observation and Geoinformation 23 (2013) 109–118 Chen, Q., 2010a. Assessment of terrain elevation derived from satellite laser altimetry over mountainous forest areas using airborne lidar data. ISPRS Journal of Photogrammetry and Remote Sensing 65, 111–122. Chen, Q., 2010b. Retrieving vegetation height of forests and woodlands over mountainous areas in the Pacific Coast region using satellite laser altimetry. Remote Sensing of Environment 114, 1610–1627. Clark, M.L., Clark, D.B., Roberts, D.A., 2004. Small-footprint lidar estimation of subcanopy elevation and tree height in a tropical rain forest landscape. Remote Sensing of Environment 91, 68–89. Dean, C., Roxburgh, S., Mackey, B.G., 2004. Forecasting landscape-level carbon sequestration using gridded, spatially adjusted tree growth. Forest Ecology and Management 194, 109–129. Dobson, M.C., Ulaby, F.T., LeToan, T., Beaudoin, A., Kasischke, E.S., Christensen, N., 1992. Dependence of radar backscatter on coniferous forest biomass. IEEE Transactions on Geoscience and Remote Sensing 30, 412–415. Drake, J.B., Dubayah, R.O., Clark, D.B., Knox, R.G., Blair, J.B., Hofton, M.A., Chazdon, R.L., Weishampel, J.F., Prince, S., 2002. Estimation of tropical forest structural characteristics using large-footprint lidar. Remote Sensing of Environment 79, 305–319. Dubayah, R., Knox, R., Hofton, M., Blair, J.B., Drake, J., 2000. Land surface characterization using LiDAR remote sensing. In: Hill, M.J., Aspinall, R.J. (Eds.), Spatial Information for Land Use Management. Gordon and Breach Science Publishers, Singapore, pp. 25–38. Duncanson, L.I., Niemann, K.O., Wulder, M.A., 2010. Estimating forest canopy height and terrain relief from GLAS waveform metrics. Remote Sensing of Environment 114, 138–154. Duong, H., Lindenbergh, R., Pfeifer, N., Vosselman, G., 2009. ICESat full-waveform altimetry compared to airborne laser scanning altimetry over The Netherlands. IEEE Transactions on Geoscience and Remote Sensing 47, 3365–3378. Duong, H., Pfeifer, N., Lindenbergh, R., 2006a. Analysis of repeated ICESat full waveform data: methodology and leaf-on/leaf-off comparison. In: Workshop on 3D Remote Sensing in Forestry, Vienna. Duong, H., Pfeifer, N., Lindenbergh, R., 2006b. Full-waveform analysis: ICESat laser data for land cover classification. In: Remote Sensing: From Pixels to Processes, ISPRS Commission VII Mid-term Symposium, Enschede, the Netherlands, pp. 30–35. Duong, V.H., Lindenbergh, R., Pfeifer, N., Vosselman, G., 2008. Single and two epoch analysis of ICESat full waveform data over forested areas. International Journal of Remote Sensing 29, 1453–1473. EEA, 2009. Corine Land Cover 2000 (CLC2000) 100 m – Version 12/2009. Available: http://www.eea.europa.eu/data-and-maps/data/corine-landcover-2000clc2000-100-m-version-12-2009 (accessed 05.10.09). 2009. Forestry Commission in England. Available: FCGB, http://www.forestry.gov.uk/forestry/HCOU-5YHJS7 (accessed 22.09.09). Hamilton, G.J., Christie, J.M., 1971. Forestry Commission Booklet No. 34 Forest Management Tables (Metric). Her Majesty’s Stationery Office. Harding, D.J., Carabajal, C.C., 2005. ICESat waveform measurements of withinfootprint topographic relief and vegetation vertical structure. Geophysical Research Letters, 32. Helmer, E.H., Lefsky, M.A., 2006. Forest canopy heights in Amazon River Basin Forests as estimated with the Geoscience Laser Altimeter System (GLAS). In: AguirreBravo, C., Pellicane, P.J., Burns, D.P., Draggan, S. (Eds.), Monitoring Science and Technology Symposium: Unifying knowledge for Sustainability in the Western Hemisphere Proceedings RMRS-P-42CD. US Department of Agriculture, Forest Services, Rocky Mountain Research Station, Fort Collins, CO, pp. 802–808. Hiratsuka, M., Toma, T., Ramada, M., Heriansyah, I., Morikawa, Y., 2003. A general allometric equation for estimating biomass in Acacia mangium plantations. In: Proceedings of the International Conference on Tropical Forests and Climate Change: Carbon Sequestration and Clean Development Mechanism, Conference Committee, Traders Hotel, Roxas Blvd, p. xvii. Houghton, R.A., 2005. Aboveground forest biomass and the global carbon balance. Global Change Biology 11, 945–958. Houghton, R.A., 2008. Biomass. In: Sven Erik, J., Brian, F. (Eds.), Encyclopedia of Ecology. Academic Press, Oxford, pp. 448–453. Hubbard, W., Latt, C., Long, A., 1998. Forest Terminology for Multiple-use Management. University of Florida. IPCC, 2003. Good Practice Guidance for Land Use, Land Use Change and Forestry. http://www.ipccnggip.iges.or.jp/public/gpglulucf/gpglulucf.html Available: (accessed 25.08.09). JNCC, 2006. The New Forest SAC, Natura 2000 Standard Data Form. Available: http://www.jncc.gov.uk/ProtectedSites/SACselection/n2kforms/UK0012557.pdf (accessed 09.06.09). Kaufman, Y.J., Tucker, C.J., Fung, I., 1990. Remote sensing of biomass burning in the tropics. Journal of Geophysical Research 95, 9927–9939. Keller, M., Palace, M., Hurtt, G., 2001. Biomass estimation in the Tapajos National Forest Brazil: examination of sampling and allometric uncertainties. Forest Ecology and Management 154, 371–382. Kimes, D.S., Ranson, K.J., Sun, G., Blair, J.B., 2006. Predicting lidar measured forest vertical structure from multi-angle spectral data. Remote Sensing of Environment 100, 503–511. Lefsky, M.A., Cohen, W.B., Parker, G.G., Harding, D.J., 2002. LiDAR remote sensing for ecosystem studies. BioScience 52, 19–30. Lefsky, M.A., Harding, D.J., Keller, M., Cohen, W.B., Carabajal, C.C., Del Bom EspiritoSanto, F., Hunter, M.O., de Oliveira Jr., R., 2005. Estimates of forest canopy height and aboveground biomass using ICESat. Geophysical Research Letters 32, L22S02, http://dx.doi.org/10.1029/2005GL023971.
117
Lefsky, M.A., Keller, M., Pang, Y., De Camargo, P.B., Hunter, M.O., 2007. Revised method for forest canopy height estimation from Geoscience Laser Altimeter System waveforms. Journal of Applied Remote Sensing 1 (1), 013537, http://dx.doi.org/10.1117/1.2795724. Mallet, C., Bretar, F., 2009. Full-waveform topographic lidar: state-of-the-art. ISPRS Journal of Photogrammetry and Remote Sensing 64, 1–16. Markwardt, C.B., 2009. Non-linear least squares fitting in IDL with MPFIT. In: Pacific, A.S.o.t. (Ed.), Astronomical Data Analysis and Systems XVIII. San Francisco, p. 251. Mitchard, E.T.A., Saatchi, S.S., White, L.J.T., Abernethy, K.A., Jeffery, K.J., Lewis, S.L., Collins, M., Lefsky, M.A., Leal, M.E., Woodhouse, I.H., Meir, P., 2012. Mapping tropical forest biomass with radar and spaceborne LiDAR in Lopé National Park Gabon: overcoming problems of high biomass and persistent cloud. Biogeosciences 9, 179–191. Muukkonen, P., Heiskanen, J., 2007. Biomass estimation over a large area based on standwise forest inventory data and ASTER and MODIS satellite data: a possibility to verify carbon inventories. Remote Sensing of Environment 107, 617–624. NASA, 2009. ICESat Home Page. Available: http://icesat.gsfc.nasa.gov/ (accessed 23.01.09). Nelson, R., 2010. Model effects on GLAS-based regional estimates of forest biomass and carbon. International Journal of Remote Sensing 31, 1359–1372. Ni-Meister, W., Jupp, D.L.B., Dubayah, R., 2001. Modelling LiDAR waveforms in heterogeneous and discrete canopies. IEEE Transactions on Geoscience and Remote Sensing 39, 1943–1958. North, P.R.J., Rosette, J.A.B., Suárez, J.C., Los, S.O., 2010. A Monte Carlo radiative transfer model of satellite waveform LiDAR. International Journal of Remote Sensing 31, 1343–1358. NSIDC, 2008. GLAS Altimetry Product Usage Guidance. Available: http://nsidc. org/data/docs/daac/glas altimetry/pdf/NSIDC AltUserGuide Rel29.pdf (accessed 28.09.09). NSIDC, 2009. ICESat/GLAS Data: Tools. Available: http://nsidc.org/data/icesat/ tools.html (accessed 12.08.09). NSIDC, 2009a. GLAS Altimetry Data Dictionary [Online]. Available: http://nsidc. org/data/docs/daac/glas altimetry/data dictionary.html (accessed 29.09.09). Pang, Y., Lefsky, M., Sun, G., Miller, M.E., Li, Z., 2008. Temperate forest height estimation performance using ICESat GLAS data from different observation periods. In: Proc. ISPRS Congress, vol. 37, Beijing, pp. 1682–1750. Parker, G.G., Lefsky, M.A., Harding, D.J., 2001. Light transmittance in forest canopies determined using airborne laser altimetry and in-canopy quantum measurements. Remote Sensing of Environment 76, 298–309. Ranson, K.J., Sun, G., Kovacs, K., Kharuk, V.I., 2004. Landcover attributes from ICESat GLAS data in Central Siberia. In: Geoscience and Remote Sensing Symposium, 2004. IGARSS’04. Proceedings. 2004 IEEE International, pp. 753–756. ˜ D., Meier, E., Allgöwer, B., Chuvieco, E., Ustin, S.L., 2003. Modeling airborne Riano, laser scanning data for the spatial generation of critical forest parameters in fire behavior modeling. Remote Sensing of Environment 86, 177–186. Rosette, J.A.B., North, P.R.J., Suárez, J.C., 2008. Vegetation height estimates for a mixed temperate forest using satellite laser altimetry. International Journal of Remote Sensing 29, 1475–1493. Schutz, B.E., Zwally, H.J., Shuman, C.A., Hancock, D., DiMarzio, J.P., 2005. Overview of the ICESat mission. Geophysical Research Letters 32, L21S01. 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, 23–45. Sun, G., Ranson, K.J., 2000. Modeling lidar returns from forest canopies. IEEE Transactions on Geoscience and Remote Sensing 38, 2617–2626. Sun, G., Ranson, K.J., Kimes, D.S., Blair, J.B., Kovacs, K., 2008. Forest vertical structure from GLAS: an evaluation using LVIS and SRTM data. Remote Sensing of Environment 12, 107–117. Takeda, H., 2004. Ground surface estimation in dense forest. The International Archives of the Photogrammetry. Remote Sensing and Spatial Information Science 35, 1016–1023. Todd, K.W., Csillag, F., Atkinson, P.M., 2003. Three dimensional mapping of light transmittance and foliage distribution using LiDAR. Canadian Journal of Remote Sensing 29, 544–555. UNFCCC, 1997. Kyoto Protocol to the United Nations Framework Convention on Climate Change. Available: http://unfccc.int/resource/docs/convkp/kpeng.pdf (accessed 24.08.09). UNFCCC, 2009a. Copenhagen Accord. Available: http://unfccc.int/ files/meetings/cop 15/application/pdf/cop15 cph auv.pdf (accessed 08.01.10). UNFCCC, 2009b. Fact Sheet: Reducing Emissions from Deforestation in Developing Countries: Approaches to Stimulate Action. Available: http://unfccc.int/files/ press/backgrounders/application/pdf/fact sheet reducing emissions from deforestation.pdf (accessed 24.08.09). 2009. United States Geological Survey [Online]. Available: URL1, http://glovis.usgs.gov/ (accessed 14.08.09). 2009. ICESat/GLAS Data Subsetter [Online]. Available: URL2, http://nsidc.org/forms/glas subset form.html (accessed 11.08.09). Wagner, W., Ullrich, A., Ducic, V., Melzer, T., Studnicka, N., 2006. Gaussian decomposition and calibration of a novel small-footprint full-waveform digitising airborne laser scanner. ISPRS Journal of Photogrammetry and Remote Sensing 60, 100–112. Waring, R.H., Way, J., Hunt, E.R.J., Morrissey, L., Ranson, K.J., Weishampel, J.F., Oren, R., Franklin, S.E., 1995. Imaging radar for ecosystem studies. BioScience 45, 715–723.
118
I.A. Iqbal et al. / International Journal of Applied Earth Observation and Geoinformation 23 (2013) 109–118
Xing, Y., de Gier, A., Zhang, J., Wang, L., 2010. An improved method for estimating forest canopy height using ICESat-GLAS full waveform data over sloping terrain: a case study in Changbai mountains, China. International Journal of Applied Earth Observation and Geoinformation 12 (5), 385–392. Zadiraka, V.K., Igisinov, K., 1973. Analysis of the accuracy and efficiency of the Fast Fourier Transform algorithm and some of its applications. Cybernetics and Systems Analysis 9, 973–977.
Zwally, H.J., Schutz, B., Abdalati, W., Abshire, J., Bentley, C., Brenner, A., Bufton, J., Dezio, J., Hancock, D., Harding, D., Herring, T., Minster, B., Quinn, K., Palm, S., Spinhirne, J., Thomas, R., 2002. ICESat’s laser measurements of polar ice, atmosphere, ocean, and land. Journal of Geodynamics 34, 405–445. Zwally, H.J., Schutz, R., Bentley, C., Bufton, J., Herring, T., Minster, J., Spinhirne, J., Thomas, R., 2006. GLAS/ICESat L1&L2 Global Land Surface Altimetry Data, Global Land Elevation Data V028/V029, 2003-2006. National Snow and Ice Data Center. GLAS/ICESat Data Subsetter, Boulder, CO.