Estimation of forest structure metrics relevant to hydrologic modelling using coordinate transformation of airborne laser scanning data

An accurate characterisation of the complex and heterogeneous forest architecture is necessary to parameterise physically-based hydrologic models that simulate precipitation interception, energy fluxes and water dynamics. While hemispherical photography has become a popular method to obtain a number of forest canopy structure metrics relevant to these processes, image acquisition is fieldintensive and, therefore, difficult to apply across the landscape. In contrast, airborne laser scanning (ALS) is a remotesensing technique increasingly used to acquire detailed information on the spatial structure of forest canopies over large, continuous areas. This study presents a novel methodology to calibrate ALS data with in situ optical hemispherical camera images to obtain traditional forest structure and solar radiation metrics. The approach minimises geometrical differences between these two techniques by transforming the Cartesian coordinates of ALS data to generate synthetic images with a polar projection directly comparable to optical photography. We demonstrate how these new coordinatetransformed ALS metrics, along with additional standard ALS variables, can be used as predictors in multiple linear regression approaches to estimate forest structure and solar radiation indices at any individual location within the extent of an ALS transect. We expect this approach to substantially reduce fieldwork costs, broaden sampling design possibilities, and improve the spatial representation of forest structure metrics directly relevant to parameterising fully-distributed hydrologic models.


Introduction
Forested environments create unique microclimatic conditions that modulate a wide array of biophysical processes tightly linked to components of the hydrologic cycle.Structural and physiological characteristics of forests and their relationship to evapotranspiration, interception, soil moisture and energy fluxes have, therefore, been intensively studied to develop physically-based models capable of simulating water dynamics in diverse hydroclimate regimes (e.g., Wigmosta et al., 1994;Pomeroy et al., 2007;Kuraś et al., 2012).In snow-dominated regions, forests generally reduce the amount of snow present on the ground prior to the onset of spring due to snowfall interception and sublimation in the canopies.The attenuation of solar radiation as it passes through forest structural elements is also of particular importance as it controls the rate and timing of snow melt, and hence strongly determines flooding risk levels and seasonal water availability (Varhola et al., 2010a).
Modelling snow interception, radiation attenuation and other biophysical processes requires a detailed characterisation of vegetation structure.While the capacity of forests to intercept snow is primarily affected by snow density, stand architecture and branch flexibility (Parviainen and Pomeroy, 2000), spatiotemporal patterns of light transmission through the canopies are created by the interaction between local solar paths, the anisotropy of diffuse sky brightness, cloud cover and the three-dimensional distribution of all canopy elements (i.e., foliage, branches, boles and gap space) (Hardy et al., 2004).Variations of these factors can create an unlimited array of micro-environments within a forest, each with a distinctive gap distribution that ultimately determines how much of the falling snow and incoming radiation actually reaches the ground (Essery et al., 2007;Hardy et al., 2004).Characterising sub-canopy snow dynamics and radiation regimes within a specific spatial unit thus requires quantification of this structural complexity into numerical parameters readily available as inputs for hydrologic models.
Three of the metrics most commonly used to describe forest structure and its relationship to hydrologic processes are leaf area index (LAI), gap fraction (GF) and sky-view factor (SVF).Although LAI has several definitions (Bréda, 2003), it is generally described as the ratio of one-half of the total leaf area per unit of ground surface area (Chen et al., 1997).Despite known difficulties with the accurate estimation of LAI, one of its versions known as effective LAI or plant area index (PAI) has become a key input parameter in physicallybased hydrologic models because it directly affects rain and snow interception, wind speed reduction and radiation attenuation in forested environments (Ellis et al., 2010).GF is the fraction of view that is unobstructed by canopy elements in any particular angular direction (Welles and Cohen, 1996), equivalent to the probability of a light beam passing through the forest to reach a point near the ground (Danson et al., 2007).SVF, used to model absorption of longwave radiation by snow (Wigmosta et al., 2002), is usually defined in hydrologic models as the fraction of celestial (sky) hemisphere visible from a point near the forest floor (Sicart et al., 2004) and is calculated as a cosine-weighted 180 • integration of GF (Frazer et al., 1999).
Several methods have been developed to directly or indirectly estimate LAI, GF, and SVF in the field.One instrument frequently used is the LI-COR ® LAI-2000 Plant Canopy Analyser (LAI-2000), which can provide LAI and GF by simultaneously comparing incoming diffuse radiation above and below the canopy (Welles and Norman, 1991).Hemispherical photography (HP) is another popular alternative that uses skyward-looking images taken from beneath the forest to estimate various attributes of canopy structure and to model light penetration over periods of time (i.e., growing season).Both the LAI-2000 and HP are based on a hemispherical projection geometry usually comprising a wide field of view (AOV) (∼ 180 • for HP and 148 • for LAI-2000), which is fundamental to provide multi-angular estimates of GF and, in HP, to account for local solar paths and the angular variation in diffuse sky brightness.Advantages of HP over the LAI-2000 are that HP does not require above-canopy measurements of diffuse sky radiation to compute GFs and it provides a permanent image of the forest that can be processed with software tools to automatically obtain a variety of structural and site-specific radiation parameters (e.g., Gap Light Analyser (GLA) (Frazer et al., 1999) or Hemiview; Rich et al., 1999).Although HP is not free of bias in the presence of heterogeneous lighting conditions and is subject to certain subjectivity when manually binarising the images to separate canopy and sky pixels, it has been validated as a tool to accurately model radiation regimes beneath forest canopies, provided that a few basic local parameters are known (Coops et al., 2004).Hardy et al. (2004), for example, compared above-and below-canopy incoming global solar radiation measurements from pyranometers with radiation transmission estimates obtained from HP, and concluded that both agreed well enough to be interchangeably used in snow models.
One disadvantage of HP, however, is that image acquisition and processing are time-consuming and, therefore, cannot be easily applied to vast, remote areas (Essery et al., 2007).Light Detection And Ranging (LiDAR), on the other hand, is a remote-sensing technology capable of providing three-dimensional representations of canopy structure over large, continuous regions.LiDAR sensors actively emit laser pulses and record the distance between sensor and target, providing point cloud-type representations of the scanned objects.LiDAR systems are generally classified as either Terrestrial Laser Scanning (TLS) or Airborne Laser Scanning (ALS) according to platform, or as discrete or full-waveform according to the type of digitisation filter (Lim et al., 2003).Discrete ALS sensors mounted on helicopters or airplanes at low flying altitudes (500-1000 m) are currently the most widely used LiDAR systems in forestry (Lee et al., 2009) due to their extensive spatial coverage and sampling densities of one to several laser returns per m 2 (Wulder et al., 2008).
There is a significant body of literature investigating the application of ALS to predict traditional stand attributes such as tree density, diameter, height, timber volume, biomass and forest cover (e.g., Lim et al., 2003;Lovell et al., 2003;Naesset, 2002;Wulder et al., 2008), while only a few articles have directly compared ALS metrics with HP-derived stand parameters.Solberg et al. (2006), for example, parameterised models to estimate LAI from discrete ALS by fitting simple ALS return penetration ratios to LAI data obtained from HP and LAI-2000 measurements, with the aim of detecting and mapping defoliation caused by an insect outbreak in Norway.They found a strong linear relationship between the log-transformed inverse of vertical GF obtained from repeated ALS acquisitions and LAI estimated in the field.To improve the relationship, follow-up studies reapplied similar methodologies varying image pre-processing procedures (Hanssen and Solberg, 2007) and testing different ranges of ALS plot radii, tree species and ALS return configurations (Solberg et al., 2009;Solberg, 2010).These and comparable articles published by Riaño et al. (2004), Morsdorf et al. (2006), Jensen et al. (2008Jensen et al. ( , 2011) ) and Korhonen et al. (2011), all rely on regression-based estimates of LAI or GF using simple vertical ALS return ratios obtained from cylindrical plots as predictors.The studies recognise and conclude that the different perspectives and projection Hydrol.Earth Syst.Sci., 16, 3749-3766, 2012 www.hydrol-earth-syst-sci.net/16/3749/2012/ geometries associated with HP (upward-looking, angular) versus ALS (downward-looking, near-vertical) sensors make it difficult to establish an exact match between the two techniques.
The objective of this article is to develop a methodology to obtain HP-equivalent forest canopy GF, LAI, SVF and solar radiation transmission metrics at any location within a discrete ALS cloud of points.Our approach transforms the Cartesian coordinates of the ALS return cloud into a polar coordinate system to produce synthetic, upward-looking hemispherical images suitable for processing with specialised software (GLA).Metrics obtained from these images are then calibrated directly with real optical HP counterparts collected within a network of ground-reference sites.This novel approach has the following advantages compared to previous studies that have attempted to link ALS and HP metrics: (1) it is based on the same geometrical projection and, therefore, minimises calibration errors; (2) it takes advantage of the entire functionality of GLA or Hemiview, including the calculation of forest structure parameters and a variety of light indices for user-defined requirements; (3) it is less restricted to any particular spatial resolution associated with ALS cylinder size (Zhao and Popescu, 2009); (4) it does not require direct radiation measurements for validation due to the proven ability of HP to predict radiation regimes (Hardy et al., 2004); (5) it is based on a paired one-on-one comparison of hemispherical images rather than plot averages, allowing a more detailed exploration of the ideal physical representation of canopies by ALS and the direct input of point-level forest structure metrics into hydrological models to analyse relevant processes at the finest possible scale (e.g., Pomeroy et al., 2007); and, finally, (6) it only requires raw ALS data and HP without relying on manual ground measurements, complementary spectral remote sensing tools or sophisticated tree-reconstruction or stem mapping techniques (e.g., Roberts et al., 2005).
The analyses of this study are focused on obtaining the forest structure metrics that are currently used by most hydrologic models at any point within an ALS point cloud.The direct input of these remotely sensed variables into the models is not tested because the main benefit of this methodology is the better characterisation of canopy structure in space rather than the simulation of hydrologic processes at the point level, which can be achieved with traditional optical HP.Future work will take advantage of the opportunity to generate thousands of synthetic hemispherical images derived from ALS to assess the spatial distribution of forest structure metrics relevant to hydrologic modelling at the watershed level, and later fulfill the ultimate goal of allocating fully distributed, spatially explicit versions of these metrics to the models.

Study area
The study took place in central British Columbia (BC), Canada, near the cities of Quesnel and Vanderhoof (Fig. 1).For a decade, this area has been affected by an outbreak of mountain pine beetle (MPB) (Dendroctonus ponderosae) that significantly changed the landscape by defoliating and killing large continuous forests of lodgepole pine (Pinus contorta), the dominant species.The Interior Plateau of BC is characterised by cold, dry winters with snow cover for up to seven months every year; snow melt constitutes a main source of water during spring and is also associated with annual peak streamflow (Bewley et al., 2010).Since the physical processes that govern snow accumulation and ablation are highly sensitive to changes in forest cover, the impacts of forest disturbance on hydrologic regimes has recently been under intensive research in BC (Bewley et al., 2010;Boon, 2009;Coops et al., 2009;Teti, 2008Teti, , 2009;;Uunila et al., 2006;Varhola et al., 2010b).
Seven forested plots established by Teti (2008) provided the model calibration data for this study, and are described in detail in Table 1.The first four plots are located in the Baker Creek watershed, near Quesnel, while three are southwest of Vanderhoof.All are characterised by their low-gradient, relatively flat terrain and each consists of 36 sampling points separated by 10 m to create a squared 50 × 50 m grid, as shown on Fig. 1.The plots are representative of the main stand types and their relative predominance in the area at the time of their installation, namely: mature stands (height > 15 m) where most of the trees have been severely defoliated by MPB (BOD1, BOD3, VOD1 and VOD2); intermediate (height ∼ 10 m) stands affected by MPB, but with their trees still holding dehydrated, red foliage (BRC2); and young healthy-looking stands (height ∼ 3 m) (BRC1).An additional plot was located in a dense stand resulting from postfire regeneration (VYN), with high stem densities and 25 % mortality caused by within-stand competition and ice/snowrelated breakage rather than MPB.More information about these plots and the methodology for capturing inventory metrics and MPB defoliation are available on Teti (2008).Eleven additional plots of a different size and configuration, which included additional stand types such as healthy spruce, were used to validate our methodology at the plot-level.

Data acquisition for modelling
Hemispherical photographs were taken within 1 m of each of the 36 sampling points during the summer of 2008 using a Nikon Coolpix 4500 digital camera with a Nikkor FC-E8 auxiliary fisheye lens (view angle = 183 • ) mounted 110 cm above the ground as specified by Teti (2008).Although it is recommended to capture HP during overcast skies to favour maximum contrast with the canopy elements (Frazer et al.,  a Following codes by Teti (2008); first letter in code corresponds to the study area.b Defoliation percentage calculated as proportion of basal area falling into each health category, which are described by Varhola et al. (2010b).(Teti, 2008).
2001), these ideal conditions are rarely present in central BC.To prevent sunlight from directly hitting the lens, Teti (2008) used a small shading paddle which was later eliminated from the images by careful retouching.Although the sampling points were registered by a GPS with differential correction, the exact location of the optical camera when acquiring the HP still contained some error and thus the maximum estimated deviation between each point and the actual corresponding camera position was approximately 1.5 m.
ALS data were acquired in February 2008 by Terra Remote Sensing (Sidney, BC) with a TRSI Mark II discrete return sensor mounted on a Bell 206 Jet Ranger helicopter at a flying altitude of ∼ 800 m above ground level.The sensor's wavelength was 1064 nm with a pulse repetition frequency of 50 kHz, maximum off-nadir scan angle of 15 • , and a fixed beam divergence angle of 0.5 mrad.A 200 km × 400 m ALS transect was acquired over the ground plots in four separate sections, as shown on Fig. 1, with a resulting average foot-print size of 0.35 m and an average effective density of 4.8 returns m −2 .To estimate plot center elevations, a groundlevel digital elevation model (DEM) with 5 m pixels (25 m 2 ) was created by applying the ground filter algorithm used in FUSION ® software (McGaughey, 2010) to the ALS data as proposed by Kraus and Pfeifer (1998).
A selection of basic ALS metrics was obtained for each 50 × 50 m plot to explore the overall variability of forest structure and data configuration (Table 2).The total number of ALS returns per plot was separated into sub-canopy (below 0.5 m) and canopy (above 0.5 m) classes.ALS return density was calculated by dividing the total number of laser returns by plot area (2500 m 2 ), while ALS metadata provided mean absolute scan angle directly.Vertical GF was calculated for each plot as the ratio of sub-canopy returns (> 0.5 m) to total ALS returns.A mean canopy height proxy was the average height above ground of all canopy returns (different from mean tree height).In all cases, no distinction was made between first and other return types.

Synthetic ALS hemispherical image generation
ALS data was extracted for 75 m radius cylinders centered at each of the sampling points, a size chosen to ensure that enough ALS returns were included closer to the horizon to mimic the infinite viewing distance of optical HP, yet small enough so that the entire cylinders fitted into the 400 mwide data transect.Only the three northernmost rows of plot VOD1 were excluded as they were too close to the ALS boundary, thus reducing the sample to a total of 234 cylinders.All the ALS returns (first, intermediate, last) were included in the cylinders to maximise density as proposed by Todd et al. (2003) and Lee et al. (2009).
Each Cartesian (XYZ) reference position was set to 60 cm above the ALS DEM to account for the fact that ALS was collected during winter when ground returns recorded a snow layer averaging 50 cm of depth (Coops et al., 2009)   of view were eliminated from the 75 m cylinders to increase data processing efficiency.The XYZ positions of all the remaining returns were transformed with simple trigonometry to polar coordinates composed of angles of azimuth ( • ) and zenith ( • ), and distance (m) with respect to each HP reference position.Finally, angles of azimuth were flipped in an eastwest direction to reflect an upward-looking field of view as in HP.
The fine-scale representation of physical vegetation structure by individual ALS returns is not well understood, and several assumptions are therefore required to convert laser points into geometrically simplified, 3-D plant structures.We undertook a sensitivity analysis on a subsample of calibration plots to explore the impact of three main parameter settings related to projected canopy element size and shape: (1) projecting returns with a fixed or variable size (inversely proportional to distance), or a combination of the two; (2) minimum ALS return circle size for fixed projections; and (3) ALS return sphere size for variable projections.More details about these parameters and their implications are clarified below as the methodology for generating the synthetic images is explained.The optimal parameter settings were chosen by evaluating scatter plots and correlation coefficients of observed GFs in real versus synthetic images, while keeping the other parameters constant.Figure 2 shows an example testing three different minimum fixed projected sizes.
Based on the results of the sensitivity analysis, each ALS return was represented as an opaque sphere with a 15 cm diameter centered on the original ALS XYZ return location.These spheres were projected as black circles on a twodimensional plane to create one synthetic hemispherical image for each sampling cylinder.Calculating the diameter of each projected ALS return first required the selection of an arbitrary radius to the circular images (r) and a theoretical focal length (f ), both 10 cm.The ratio between the projected (r p ) and full-scale ALS return radii (r r ) is then assumed equivalent to the ratio between the focal length (f ) and the absolute distance between the return's centroid and the camera position (d): r p /r r = f/d.This is illustrated in Fig. 3a.
Optical distortions typical of hemispherical lenses were accounted for when generating the synthetic images.When viewed from a distance d, a sphere subtends an angle equal to the arctangent of the ratio between the sphere's diameter and d.A sphere located along the optical axis of a fisheye lens (i.e., at zenith = 0 • in this case) appears as a circle when projected on the image and gradually flattens into an ellipse as the azimuthal diameter is stretched in proportion to the zenith angle (Fig. 3b).This can be illustrated by considering the polar coordinate representation of three points in the celestial hemisphere -one at the zenith, one at the east horizon, and one at the north horizon.In polar coordinates, the line connecting the zenith with one of the points on the horizon has a length equal to the image radius (r), whereas the line connecting the two points on the horizon has a length of π ×r/2, even though the angular separation between all three points is 90 • .Thus, the apparent area of a feature in a hemispherical image increases from zenith to horizon according to Z/90 × (π/2 − 1), where Z ( • ) is the zenith angle of the feature (e.g., a sphere projected in the horizon (90 • ) appears with its azimuthal radii π/2 (57 %) larger than a sphere at the zenith).To simplify plotting synthetic images, we used circles to represent the modelled ALS returns with their areas increased to account for the stretching.A subsample of images with and without this radial geometric correction was generated to quantify the effect of optical distortions on LAI in different stand types.
Beam fraction Ratio of direct to total spectral radiation incident on a horizontal surface at the ground over a specified period, which is a function of cloud cover for supra-daily periods.Values calculated from cloudiness index as explained by Frazer et al. (1999).

Sky region brightness
Method for describing the intensity of the solar disk and diffuse sky.
The selected Standard Overcast Sky (SOC) assumes that the zenith is three times as bright as the horizon.

SOC SOC
Clear sky transmission coefficient Factor that describes the regional clarity of the atmosphere with respect to the instantaneous transmission of direct (beam) radiation.Value used recommended for the area by Frazer et al. (1999).
The radial location of ALS returns in the synthetic images followed the same equiangular projection produced by the optical lens system (FC-E8 fisheye converter) used to collect the real HP images (Inoue et al., 2004), where the radial distance from the centre of the image is directly proportional to the zenith angle (Rich et al., 1999).A projection based solely on variable diameters resulted in unrealistic-looking images with distant returns appearing too small to be detected as dark pixels.To avoid this, a minimum base constant return diameter (2.15 mm within a 20 cm diameter image, in this case; Fig. 2) was assigned to all returns in the images so that only those returns close enough to the HP reference to exceed this diameter were projected at variable, larger sizes.Finally, ALS returns closer than 75 cm to the reference were eliminated in accordance with the HP field procedures, which specified a minimum distance between foliage and camera lens for protection purposes and to eliminate the possibility of having a large proportion of the optical field of view obscured by a single foliage unit.
Each synthetic image was generated as a 750 × 750 pixel bitmap image file (BMP).All 234 files were automatically created in 5 min and 15 s using MATLAB ® (1.3 s/image) on a 2.66 GHz, 4.0 GB RAM, 64-bit computer.

HP analysis
Both optical HP and ALS synthetic images were analysed with GLA using the parameters listed on Table 3.The binarisation thresholds required for optical HP were provided by Teti (2008), while synthetic ALS photos were logically generated in black and white.GLA processing was done manually for each of the 468 images (optical and synthetic) at a rate of around 17 s per file (excluding HP binarisation).

LSA
Mean absolute ALS scan angle a All definitions of GLA-derived variables have been taken from Frazer et al. (1999), where they are described with more detail.b FLAI θ and LLAI θ are obtained from the appended output of GLA, where LAI 4 ring corresponds to θ = 60 • and LAI 5 ring corresponds to θ = 75 • .For additional information on LAI see Welles and Norman (1991) and Stenberg et al. (1994).
Table 4 describes the output variables obtained from the GLA analysis for each HP and synthetic image for the 234 sample locations.All of the variables except LAI were numerically integrated across the following zenith angle AOV (θ): 0-30, 0-45, 0-60, 0-75 and 0-90 • .Lower zenith angle limits were included here to allow later explorations of their relationships with hydrologic processes, based on previous findings.Teti (2003), for example, concluded that a 0-30 • zenith AOV was most effective at explaining differences in snow storage in the presence of different sized gaps.In addition, sky regions including the 0-45 • and 0-60 • zenith ranges contain most of the solar paths directly responsible for spring melt in our study area.FLAI θ and LLAI θ (defined on Table 4) automatically integrate LAI for zenith angles of 60 and 75 • (Welles and Norman, 1991;Stenberg et al., 1994).
In light of intensive discussions and debate about true LAI derivation for decades (Bréda, 2003), we highlight that this study is only focused on the version of LAI that has been widely used in hydrologic models: that obtained from LAI-2000 or HP, commonly known as effective LAI (eLAI) or plant area index (PAI) (e.g., Pomeroy and Dion, 1996;Bewley et al., 2010).It is not our goal to correct for clumping, isolate foliage from stems or formulate hypotheses about the angular distribution of leaves because more nuanced (or simply different) versions of LAI would also require the reparameterisation of existing hydrologic models.In this respect, LAI obtained from optical HP is considered as our groundtruth variable.

Standard ALS metrics
The traditional ALS metrics (without coordinate transformation) listed at the bottom of Table 4 were calculated from variable-size cylinders whose diameter hypothetically intercepted the canopies at each zenith AOV (θ) based on plotlevel maximum tree height (Table 1) (for θ = 90 • , cylinder diameters calculated for θ = 75 • were used).From these point clouds, vertically projected gap fraction (LVGF), mean laser canopy height (LMH), return density (LD) and scan angles (LSA) were estimated for each of the field plot locations (Table 2).LVGF and LMH are important descriptors of stand structure, while LD and LSA are principally related to ALS sensor configuration and data acquisition conditions, all potentially important sources of variation as shown in Table 2.More specifically: -LVGF represents a downward-looking ALS ground-tocanopy return penetration ratio known to be well correlated with upward-looking HP sky/canopy GFs (Solberg, 2010); -LMH is an indicator of stand height and provides substantial differentiation between stands; -LD accounts for the varying ALS return density among plots (see Table 2) produced by changing flight altitude and speed (Bater et al., 2011;Goodwin et al., 2006), which can alter the probability of ALS returns being intercepted by the canopy.Varying patterns of flight line overlap may also contribute to markedly different laser return densities throughout the spatial coverage, and must be explicitly accounted for when neighbouring or intersecting datasets result in higher numbers of returns (our ALS data was restricted to a single flight line and not subject to changes in density due to overlap); -LSA increases from the centre of the flight line (nadir) towards the swath edge and changes the probabilities of ALS hitting different canopy sections.For example, if scan angles are too large, laser pulses are less likely to penetrate the canopy because of a higher probability of being intercepted, resulting in a different spatial representation of the forest (Korhonen et al., 2011).

Regression modelling
GF, LAI, SVF and total direct and diffuse (global) radiation transmittance are GLA outputs directly applicable in hydrologic simulators (e.g., Wigmosta et al., 1994), and were, therefore, selected as the main response variables in this study (Table 4).To simplify the regression analyses, given the large number of variables and zenith angle of view (AOV) combinations, a specific modelling strategy was designed to ensure that all models: (1) shared a unique, stable structure and (2) were applicable to the full range of sampling points, avoiding the need to pre-identify different forest populations or use indicator (dummy) variables.A preliminary analysis showed no evidence suggesting the existence of nonlinear relationships between the variables, so all models were maintained linear and variable transformations were not necessary.
Both simple and multiple linear regression analyses were applied to predict the four HP-derived metrics: where θ is the maximum zenith AOV for metric aggregation (30, 45, 60, 75 and 90 • ); Y S and Y M can be either FGF θ , FRT θ , FSVF θ or FLAI θ ; X 1 is the ALS-derived counterpart of Y S or Y M (LGF θ , LRT θ , LSVF θ , or LLAI θ , respectively); X 2 corresponds to LVGF, X 3 is LMH; X 4 is LD; and X 5 is LSA.Please refer to Table 4 for a comprehensive definition of these variables.
Correlations between the dependent variables and each of the independent variables were examined through scatterplots, and a correlation matrix between all predictor variables (plus FGF θ ) was produced to ensure that variables with high inter-correlations were not added prior to fitting the models.
Equations ( 1) and (2) were fitted using ordinary leastsquares regression for each zenith AOV (θ) metric.Also, in order to justify the need and benefit of performing ALS coordinate transformations, a simple linear model was fitted to predict HP gap fraction (FGF θ ) with the vertical ALS gap fraction (LVGF) only: where Y S(θ) is the same as for Eq. ( 1) and X 1 corresponds to either simple LVGF or the transformation used by Solberg et al. (2006): ln(LVGF −1 ).
Choosing the best common multiple regression model structure followed a manual backward stepwise approach for variable selection.The intercept and all five predictor variables (Eq.2) were initially included in the regression to obtain a matrix of coefficient p-values that included the entire response variable -zenith AOV model combinations.Model coefficients of the supporting predictor variables (X 2 . . .X 5 ) that were significant less than twice in the matrix were removed until only those variables consistently showing statistical significance across all models were identified.The equations were then validated using the tests described below.
Goodness of fit was evaluated based on the models' adjusted coefficients of determination (r 2 and R 2 ) and three versions of root mean squared error: absolute (RMSE), leaveone-out cross-validation average derived by iteratively refitting the model with all but 1 of 9 randomly generated data groups (RMSE S ), and normalised by the range of observed values to enable comparisons between different variables (RMSE N ).All models were validated by observing the significance probability (p) of the regression coefficients and by performing Anderson-Darling and Shapiro-Wilk normality tests on the residuals, which confirmed the required normality if p > α (in all tests, α = 0.05).Linearity and variance stability was assessed through visual inspection of predicted vs. observed and predicted vs. residual scatterplots.
95 % prediction and confidence intervals were calculated and illustrated on the predicted vs. observed figures.For multiple regression, the confidence intervals were estimated with a quadratic function relating individual predicted values to their corresponding lower and upper limits as generated by the statistical software (I = c 0 + c 1 × P + c 2 × P 2 , where I = upper or lower interval limit, P = predicted value and c 0 , c 1 and c 2 are model coefficients).Variance inflation factors were also estimated to check for multicollinearity (Field, 2005).

Plot-level ground-truth model validation
Using the methodologies presented here, new and different grids of synthetic hemispherical images were generated in the same calibration plots and an additional set of 11 plots (other than those in Table 1) as part of a follow-up study.The resulting average GFs were compared with those estimated from pre-existing optical HP available in these plots.Although the number and distribution of synthetic and optical HP differed substantially within each plot, this comparison was very useful to validate our methodology with an independent dataset and justify the need for ALS coordinate transformation even if plot-level averages were required.This was done by contrasting the relationship between optical HP-derived GFs and both a simple ratio of raw ALS ground/canopy returns (calculated as LVGH, but for the extent of individual square plots) and GFs derived from synthetic hemispherical images.

Preliminary sensitivity analysis
The sensitivity analysis applied to a sub-sample of images indicated that the projection of ALS spheres using inverse distance-weighted diameters generally had little effect for the majority of synthetic images generated in mature stands, due to the large distances between the returns and the projection reference.Increasing the theoretical diameter of the spheres beyond 15 cm made returns close to the reference appear too big and subsequently blocked significant portions of the image, while substantially deteriorating the relationship between synthetic and optical estimates of GFs.However, variable-diameter projections were still necessary to produce an adequate representation of canopy structure in the young regeneration stand, where 15 cm spheres appeared optimal in all cases.On the other hand, varying the minimum size of projected returns is important for the calibration relationships because it affects the r 2 via changes in the regression slope, but will not significantly influence RMSE due to a lack of discernible scatter reduction (Fig. 2).

ALS-derived synthetic hemispherical photos
Figure 4 shows examples of ALS return clouds with corresponding synthetic and actual HP for six common stand structure types found in the study site.In general, there is good visual agreement between the HP and synthetic images across the stands, although the return density (∼ 5 returns m −2 ) of this ALS dataset makes identification of individual trees difficult.In the taller stands (BOD1-C5, BOD3-C4), larger and continuous gaps in the forest canopy are apparent at smaller zenith angles of both the actual and synthetic images.Denser and more homogeneous canopies (BRC2-C1, VYN-B4) show a more even distribution of canopy elements (i.e., ALS returns) across all the zenith angles.Markedly different to the other stands is the young regeneration (BRC1-A3), where images are dominated by sky and shorter, clumped vegetation leads to interception of ALS returns much closer to the projection reference.

Relationships between variables
The correlation matrix between GF as derived from the actual HP (FGF θ ) and estimated from the synthetic images (LGF θ ) for the five zenith AOV is shown in Table 5.The relationship between the ALS and HP-derived GFs is significant at all zenith angles and becomes stronger as zenith AOV increases (r = 0.75 for θ = 30 • and r = 0.93 for θ = 90 • ).
In addition, Table 5 summarizes correlations between the two angular GFs (FGF θ and LGF θ ) and the ALS simple metrics (LVGF, LMH, LD, LSA), which are generally poor.The best correlation occurred between mean LMH and LGF θ , suggesting that taller stands have a smaller GF across all zenith angles.The correlations between the predictor variables to be input in the multiple regression model (Eq.2) are generally weak, so redundancy is not likely introduced.
Scatterplots between FGF θ and LGF θ are shown for all zenith AOV on Fig. 5a-c.Differences in structure across the stands result in distinct population clusters clearly visible in the figures, which prevent fitting a single model to the data.
The relationship between LGF θ and FGF θ shows that, with the current parameter setup, the synthetic images underestimate GF when compared to HP, with the exception of stand BOD3.The young regeneration stand (BRC1) also deviates from the LGF θ -FGF θ general linear pattern and shows a considerably higher variability.

Simple linear regression
The simple regression model predicting FGF θ from ALS untransformed data (Eq.3) was weak across all zenith AOV (r 2 ranging from 0.31 to 0.41).Adjusted r 2 values increased (0.59-0.87) and RMSE decreased when the ALS transformed variable (LGF θ ) was used; however, nearly all the simple linear regression models failed residual normality tests.This is illustrated in Fig. 5d-f, where scatterplots of observed vs. predicted gap fractions indicate that a single regression line does not account for different populations of stand structures, particularly within the short regeneration stand (BRC1).

Multiple linear regression
After applying multiple regression to Eq. ( 2) with all the predictor variables included, LMH proved to be non-significant for θ values of 30 • and 45 • in all cases and for FLAI 60 and FLAI 75 , while LSA was consistently non-significant across all model specifications.The intercept (β 0 ) and LVGF were not significant in two cases only and were, therefore, maintained for a second run.After eliminating LMH and LSA from the regression, all of the remaining parameters were statistically significant with no exception.However, models for θ = 30 • did not pass the two residual normality tests, and the FLAI 60 model barely passed the Anderson-Darling test only.As suggested by Kutner et al. (2005)  non-normality of residuals in these cases while maintaining our modelling strategy.None of the transformations tested (i.e., inverse, square, square-root, log) solved the problem for FGF 30 , FRT 30 and FSVF 30 .However, using the square root of FLAI θ (FLAI 0.5 θ ) allowed the models to pass the Anderson-Darling test without losing variable significance while improving R 2 and RMSE.
Thus, the model that complied with the all the conditions was: Regression results for Eq. ( 4) are provided in Table 6, which shows that model behaviour was very similar for all dependent variables.Adjusted R 2 improves as zenith AOV increases due to more pixel aggregation that reduces the probability of canopy returns being assigned to the wrong sky region.RMSE N are very similar for FGF θ , FSVF θ and FRT θ with the same θ, while the prediction accuracy of all models is validated by the consistent similarities between RMSE and RMSE S (Kutner et al., 2005).The fitted parameter estimates of β 0 , β 1 , β 2 and β 4 shown in Table 6 can be readily used to predict FGF, FSVF, FRT and FLAI at location within the current ALS data.All models produced variance inflation factors ranging between 1.3 and 1.7, eliminating multicollinearity concerns (Myers, 1990).Scatterplots of predicted vs. observed values of FGF θ are shown in Fig. 5g-i for θ = 30, 60 and 90 • , respectively.When comparing Fig. 5d-f to Fig. 5g-i, it is evident that multiple linear regression was a successful tool to achieve more accurate predictions, especially by accounting for the distinct stand structure of the young regeneration stand (BRC1).Figure 5d and g also show that outliers might be preventing model validations for θ = 30 • .Scatterplots of predicted GFs and the model residuals appear visually satisfactory for all zenith AOV, with data points evenly distributed at both sides of the horizontal reference line of Fig. 6.  Figure 7b indicates that plot-level GF averages from optical HP are closely related to averages from new synthetic ALS images, both obtained for all the plots where ALS was available in the study area.The comparison between Fig. 7a and b constitutes strong evidence to justify coordinate transformation to accurately predict GF.

Discussion
The discussion regarding the methodology presented in this study is centred around the following questions: (1) how do synthetic hemispherical ALS images visually compare to their real HP counterparts and what are the main sources of error?; (2) how suitable is discrete ALS to represent the finescale canopy elements responsible for radiation transmission?; (3) how effective was the proposed modelling strategy aiming to predict HP-derived metrics from coordinatetransformed ALS?; (4) what are the perceived benefits of the methodology?; and ( 5) what lines of action are needed to improve and apply this approach in future research?

ALS synthetic hemispherical images
Our results indicate that coordinate transformation of ALS data produced synthetic HP images which were qualitatively and quantitatively similar to real optical HP.The use of a 75 m diameter ALS cylinder was deemed appropriate for this dataset, given that enough spheres appeared in the synthetic images at higher zenith angles to reproduce the saturation that occurs in real HP under closed canopy conditions.If needed, smaller cylinders could be tested for narrower ALS transects.The effects of the radial distortion correction were negligible in mature and medium stands, which showed a difference in resulting LAI of less than 1 % due to the increased overlap and saturation that occurs at higher view angles and larger distances between returns and reference.However, LAI was 20 % higher in images with geometric correction in young regeneration stands due to the abundance of returns closer to the reference.A visual inspection of the synthetic HP dataset showed that individual trees were difficult to identify in most stands, and that ALS returns occasionally appeared where canopy elements were absent in HP.There are three possible explanations for these differences: (1) density of the ALS point cloud was too low to capture basic crown-level structural details apparent in the optical HP; (2) HP was acquired six months after ALS and changes in stand structure (e.g., crown damage, tree fall, etc.) could have occurred in these stands affected by MPB, and (3) there were GPS positional errors in HP plot locations, camera orientation and image registration errors, as well as uncertain snowpack depths when ALS was collected.

Physical representation of canopy elements with ALS
Successful transformation of the ALS point cloud into realistic synthetic HP images depends on a number of factors.First, the density of laser returns needs to be high enough to capture the basic geometry of individual tree crowns and branches.Second, the size and shape of projected synthetic canopy elements (in this case spheres) must be closely related to some basic unit of light-intercepting foliage or branch structure found in real forests.Third, the density and distribution of laser returns must be relatively uniform throughout the entire survey area to avoid bias.It was shown here that a minimum constant projected size was necessary for all returns to resemble HP; however, inverse-distance-weighted variable projections were still necessary for returns closer to the reference, particularly for the short regeneration stand.ALS returns portrayed as opaque spheres represent a crude approximation of canopy structure as seen by a camera.Real canopy elements are far more complex, but because it is impossible for discrete ALS returns to accurately characterise fine-scale details of plant canopies, some arbitrariness is inevitable when assigning theoretical shapes and sizes to ALS returns.It must be highlighted that our methodology was not designed to reproduce the scale of detail found in real HP images as possible with higher density TLS (Côté et al., 2009), but to capture the basic patterns of canopy structure responsible for light interception and penetration that may in turn influence snow accumulation and melt.
The detection of canopy elements by ALS and HP are both dependent on optical properties of the canopy; however, the former technique is based on reflectivity while the latter on opacity.Another disparity between HP and ALS is that the downward, near-nadir view angle of ALS provides a biased vertical profile of forest canopies in which upper elements have a higher probability of being detected, leading to an under-representation of lower branches and stems (Hilker et al., 2010).This may be compensated in synthetic ALS hemispherical projections because image saturation increases towards the horizon mainly due to the corresponding exponential increase in the number of ALS returns, while in HP it is common to observe tree stems closer to the camera occluding the farther views as the main source of saturation.This will have an effect on the amount of unexplained variance in the regression models, but the strength of model fits (adj.R 2 = 0.8 to 0.92) for FGF, FSVF, FRT, and FLAI suggest that any bias in height distribution has little effect on the results.While some of these shortcomings might be minimised by more sophisticated individual tree-reconstruction routines, volumetric rendering or ray-tracing methods, the corresponding uncertainties here are partly masked and absorbed by the calibration models.
In this study, we made several assumptions about the size and shape of a specific combination of predictors and parameter estimates were chosen so that the empirical relationship between canopy metrics derived from synthetic and real HP and their visual similarity was maximised.This methodology, empirical in nature, is admittedly susceptible to interactions between parameters.For instance, a larger minimum projected circle size (Fig. 2) might be needed if ALS return density is lower, or the maximum sphere size could be reduced if returns are too close to the reference.Applying this methodology to a wider combination of forest stands and ALS datasets is required to evaluate parameter stability and optimisation.

Modelling strategy
Regression models using simple vertical gap fractions (LVGF or ln(LVGF −1 )) to estimate HP metrics generally had low r 2 , high RMSE N and produced model residuals that failed normality tests.These results contrast with those of Solberg et al. (2006Solberg et al. ( , 2009)), Hanssen and Solberg (2007) and others, in part because their statistical comparisons were based on the average of multiple photo plots rather than individual photo points, masking within-plot variability.Simple linear regression directly estimating HP metrics from their ALS-derived counterparts (X 1 ) was also unsuccessful because it failed to include other key explanatory variables (namely LVGF and LD) and describe the relationships among all stands in one single model, especially due to the deviations shown by BRC1 and BOD3 (Fig. 5).However, the need to perform coordinate transformations of ALS data to better predict HP-derived metrics in individual sampling points was strongly justified (Fig. 7).
Multiple linear regression was suitable to calibrate ALS metrics with HP by accounting for both forest structure and data configuration properties.The models' R 2 values above 0.80 and RMSE N below 10 % across all zenith AOV higher than 45 • suggest that confident predictions can be made throughout the entire ALS transect.This idea is also supported by the wide structural diversity of stands included in the regression dataset and the successful validation performed at plot-level averages in additional stands which represented even more diverse conditions (Fig. 7).Better HP geographical registrations and simultaneous HP/ALS data collection plus a detailed outlier analysis are required to fully validate the models for θ = 30 • .
A strong component of this study involved the use of a large network of individual ground-reference samples representing a heterogeneous collection of forest structure conditions that appropriately represented both within-and between-stand variability.The latter was particularly important to understand the relationship between GF derived from ALS synthetic images and HP across a broad range of GF estimates (e.g., 0.3 to 0.9 for θ = 60 • , Fig. 5h).A more complete sample of stand structures would have included mature, non-defoliated pine stands; however, this stand type was absent from the study area at the time of data collection.The accuracy of predicting HP metrics directly with ALS synthetic counterparts should be independent of stand health status in light of both ALS and HP being able to detect defoliation (Solberg et al., 2006).
The developed models (Eq.4, Table 6) proved suitable for our range of sampled forest structures and ALS data.Consequently, they need to be tested and validated for different stand types (species, densities, heights, health, etc.) and other ALS data acquisitions (e.g., return density, scan angle, footprint size, overlapping transects, return classes, etc.).Of particular interest is the application of this method to full-waveform (FW) LiDAR data, which will be increasingly used in the future and provides a more detailed profile of canopy elements and additional radiometric information (Pirotti, 2011).FW LiDAR also has the potential to assist in the improved estimation of the return dimensions by the analysis of target backscatter cross sections (Wagner et al., 2006(Wagner et al., , 2008)).However, given that discrete ALS has been used extensively in many regions, our methodology is not likely to become obsolete in the near future.(Côté et al., 2009).

Methodological advantages and applicability
There are a number of advantages associated with transforming ALS coordinates to generate hemispherical synthetic images.First, geometrical discrepancies between ALS and HP are minimised, allowing a direct comparison of structural and radiation metrics at the individual point level.Second, the methodology is simple because it is based on raw ALS pointcloud data and avoids the need for elaborate canopy models while minimising the number and complexity of geometrical parameters.Third, GLA or other specialised programmes can be used to directly estimate GF, LAI, SVF and local transmission of direct, diffuse or total radiation through forest canopies.We have chosen to generate synthetic hemispherical images from ALS which are then readily available for processing with GLA because hemispherical photography has become one of the most popular methods to obtain GF and associated forest structure metrics, and its outcomes have been systematically used to parameterise hydrologic models (e.g., Woods et al., 2006;Ellis and Pomeroy, 2007;Ellis et al., 2011).Alternative approaches might produce versions of these metrics that can deviate substantially from those used to develop existing process-based equations, thus requiring their parameters to be revised.Fourth, generating synthetic hemispherical images from ALS introduces unlimited flexibility in terms of sample size and experimental design layouts: any number of images can be obtained at user-defined spacing options and sub-pixel specific locations.Fifth, single calibration models for each variable proved to be applicable to a wide range of stand types.Finally, variables directly applicable to hydrologic modelling can now be obtained at any point within ALS datasets -significantly reducing fieldwork requirements while improving the parameterisation of vegetation classes at the landscape-level.The normal distribution of the calibration models' residuals suggests that our method is unlikely to produce systematically biased estimates of forest structure variables, which will benefit fully-distributed hydrologic modelling exercises.

Future work
Further research should focus on (1) improving the accuracy of the methodology by better geographical registration methods and coordinated data collection; (2) validating or reformulating the current models using different datasets and study areas (e.g., other species, mountainous topography) to evaluate parameter stability; (3) exploring the relationships between structural variables obtained from synthetic ALS hemispherical images and satellite-derived spectral indices for watershed-or landscape-level extrapolations; (4) developing alternative methodologies to parameterise hydrologic models with metrics directly obtained from ALS and other remote sensing technologies; and (5) improving the functionality of HP processing algorithms to estimate radiation components at sub-daily time steps (e.g., Leach and Moore, 2010).The latter represents a difficult challenge given the inaccuracies in camera orientation, anisotropy of sky brightness and atmospheric attenuation, among others; however, if achieved, it would allow the direct input of radiation transmission into point-based process simulation of hydrologic models and better performances if above-canopy radiation is available.
While GLA can directly estimate GF and radiation transmission, most hydrologic models have used LAI as the forest structure parameter input to calculate hourly or daily radiation components (e.g., Wigmosta et al., 1994;Pomeroy et al., 2007).Nevertheless, the difficulties of accurately measuring true LAI in the field are well known, and optical methods only measure the effective plant area index unless corrections are made for foliage clumping and the surface area contributed by branches and boles (Bréda, 2003).LAI-2000 or HP processed with GLA are also impacted by this bias and yet remain a popular method to estimate LAI by integrating log-transformed gap fractions through cosine-weighted zenith rings (Welles and Norman, 1991), only to be used as an intermediate parameter to simulate radiation transmission and other processes in hydrologic models.However, since radiation transmission and all the light indices available from GLA are also directly obtained from the simple sky/canopy pixel ratio defined here as gap fraction (GF), this variable might constitute a conceptually simpler and more parsimonious average forest structure parameter than LAI when modelling below-canopy radiation regimes.Since producing alternative physically-based models requires detailed measurements of above-and sub-canopy shortwave and longwave radiation, precipitation interception and evaporation, snow accumulation and depletion, among others, new studies are required to re-parameterise hydrologic models to substitute LAI with GF or other metrics directly obtained from remote sensing and quantify the resulting benefits or losses.A new era of research in hydrologic modelling should use alternative metrics derived from ALS, TLS and even satellite-derived spectral indices to develop entirely new process-based equations for multi-scale modelling of radiation transfer, precipitation interception, evapotranspiration and water routing, among others.
To our knowledge, this is the first study that has attempted to re-project discrete ALS individual returns in a polar coordinate system to directly model forest structure and radiation regimes with currently available HP image processing tools.Our results suggest that reprojection of the ALS point cloud is necessary if accurate HP-derived estimates of canopy gap fraction, LAI, SVF and solar radiation transmission are required at the point level.The goal of this study was not to provide prediction models with universal application across all forest types and ALS datasets, but to reveal the importance of coordinate transformation for the estimation of GF and other bulk-canopy metrics, and to demonstrate that these variables can be predicted from discrete ALS calibrated with HP.Our main research objective was fulfilled with the current approach as the models developed can operationally predict canopy GF, LAI, SVF and light indices with reasonable accuracy in any location within this ALS dataset, regardless of forest type.
When evaluating the strategies to improve hydrologic models with remotely sensed forest structure metrics, two alternatives arise: (1) redevelop process equations at the point level based on several years of detailed meteorological data collection to be directly linked with the 3-D capabilities of ALS or TLS to portray canopy structure, and (2) improve the characterisation of forest structure at the landscape level so that spatially explicit versions of the relevant variables can be directly input in existing fully distributed hydrologic models.While the two approaches need to be completed and this can be achieved in parallel, our study focuses on the latter mainly because it has the potential to directly improve the modelling of snow and streamflow processes on a larger scale using the current tools available.Here we present the first step towards obtaining fully distributed versions of the forest structure variables needed to parameterise current physically-based hydrologic models.Following the line of approach (2) mentioned above, the next stage consists of correlating these metrics -now available at any location within the spatial extent of an ALS dataset -with satellite-derived spectral indices in order to obtain fully distributed structural variables to replace the bulk and discrete vegetation classes that are predominantly used at present.The ultimate goal is to input these detailed, spatially explicit and continuous versions of the variables into fully distributed models to evaluate the changes in model efficiencies when estimating snow accumulation and ablation as well as streamflow generation.
As ALS is becoming increasingly available worldwide, this article represents a major contribution to hydrologic studies by facilitating the estimation of forest structure metrics relevant to model parameterisation through reduced fieldwork requirements and unlimited, spatially-explicit sampling designs.

Fig. 1 .
Fig. 1.Study area location within British Columbia (left) including ALS transects (black straight lines in the close-up) and ground plot locations in the Vanderhoof (top left corner ellipse) and Baker Creek (bottom right corner circle) areas; the 2500 m 2 square ground plots (right) are constituted by 36 individual stakes (spaced 10 m) labelled as letter and number combinations (A1 to F6)(Teti, 2008).

Fig. 2 .
Fig. 2. Effect of projected ALS return size on the relationship between observed and predicted gap fractions.The projected circle size of synthetic image examples (a, b, c) and corresponding relationships (d, e, f) is expressed as the fraction between the diameter of each ALS return and the diameter of the image: 0.0137 (a, d), 0.0176 (b, e), and 0.0215 (c, f).

Fig. 3 .
Fig. 3. Example of the projection of an ALS spherical return into a 2-D focal plane (a), and azimuthal radial distortion correction at representative zenith angles (b).

Fig. 4 .
Fig. 4. Representative examples of ALS point clouds (left), ALS synthetic hemispherical images (centre) and real optical hemispherical photographs (HP) (right) for each stand; azimuths ( • ) are shown on the hemispherical illustrations.

Fig. 5 .
Fig. 5. Relationship between ALS-derived (LGF θ ) and HP-derived gap fraction (FGF θ ) (a, b, c) and between predicted and observed values of gap fraction obtained from simple (d, e, f) and multiple (g, h, i) linear regression models across three representative zenith angle AOV (30 • , top; 60 • , centre; and 90 • , bottom); legends in sub-figures (a) and (i) apply to all.

Fig. 7 .
Fig. 7. Comparison between the relationship of plot-level average gap fractions obtained from optical HP and (a) vertical gap fraction estimated from untransformed ALS, and (b) gap fractions from calibrated synthetic hemispherical images.

Table 1 .
Stand locations and physical characteristics as of 2008.

Table 2 .
ALS simple metric summary for each major plot.

Table 5 .
Correlation matrix showing the coefficient of correlation (r) between variables used in multiple regression (gap fraction only, for simplicity); non-significant (p > 0.05) values shown in italic.

Table 6 .
Multiple linear regression results (refer to Table4for RMSE and RMSE S units).M )/main independent (X 1 ) variables; supporting variables X 2 and X 4 are common for all models.b Significant with 0.01 <= p < 0.05; c Significant with p < 0.01.
a Dependent (Y

Hydrol. Earth Syst. Sci., 16, 3749-3766, 2012 www
Despite the supporting ALS vertical variables increasing the significance of the model if applied to alternative datasets, new HP/ALS calibrations are required every time this approach is applied in a different area.This includes reassessing variable significance and full model validation.It is especially important to account for changes in ALS density caused by systematic overlapping multiple transects, where duplicate sampling might not be properly captured by LD alone.A voxelization of the ALS point cloud might minimise the effect of varying return densities, whereby each volume element (voxel) is coded as one if occupied by one or more ALS returns, or as zero if empty .hydrol-earth-syst-sci.net/16/3749/2012/