Airborne laser scanning of natural forests in New Zealand reveals the influences of wind on forest carbon

Forests are a key component of the global carbon cycle, and research is needed into the effects of human-driven and natural processes on their carbon pools. Airborne laser scanning (ALS) produces detailed 3D maps of forest canopy structure from which aboveground carbon density can be estimated. Working with a ALS dataset collected over the 8049-km2 Wellington Region of New Zealand we create maps of indigenous forest carbon and evaluate the influence of wind by examining how carbon storage varies with aspect. Storms flowing from the west are a common cause of disturbance in this region, and we hypothesised that west-facing forests exposed to these winds would be shorter than those in sheltered east-facing sites. The aboveground carbon density of 31 forest inventory plots located within the ALS survey region were used to develop estimation models relating carbon density to ALS information. Power-law models using rasters of top-of-the-canopy height were compared with models using tree-level information extracted from the ALS dataset. A forest carbon map with spatial resolution of 25 m was generated from ALS maps of forest height and the estimation models. The map was used to evaluate the influences of wind on forests. Power-law models were slightly less accurate than tree-centric models (RMSE 35% vs 32%) but were selected for map generation for computational efficiency. The carbon map comprised 4.5 million natural forest pixels within which canopy height had been measured by ALS, providing an unprecedented dataset with which to examine drivers of carbon density. Forests facing in the direction of westerly storms stored less carbon, as hypothesised. They had much greater above-ground carbon density for a given height than any of 14 tropical forests previously analysed by the same approach, and had exceptionally high basal areas for their height. We speculate that strong winds have kept forests short without impeding basal area growth. Simple estimation models based on top-of-the canopy height are almost as accurate as state-of-the-art tree-centric approaches, which require more computing power. High-resolution carbon maps produced by ALS provide powerful datasets for evaluating the environmental drivers of forest structure, such as wind.


Background
Forests are a key component of the global carbon cycle (Pan et al. 2013), storing and sequestering more carbon than any other ecosystem (Gibbs et al. 2007). Forest degradation and deforestation cause substantial releases of greenhouse gases to the atmosphere, estimated at 1-2 billion tonnes of carbon per year, which equates to 10% of global emissions (Baccini et al. 2012). Even if nations de-carbonise their energy supply chains within agreed schedules, a rise of 2°C in mean annual temperature is unavoidable unless forests are managed to store more carbon, which will involve protecting mature forests, restoring degraded forests, and taking marginal agricultural land out of production and reforesting it (Houghton et al. 2015). Accurate monitoring of forest carbon stocks underpins these climate change mitigation programmes (Agrawal et al. 2011) and most developed nations have reporting systems as part of international treaty commitments (e.g. New Zealand, Coomes et al. (2002)). The reporting systems track anthropogenic changes, but also natural process, which may be indirectly affected by humans. For example, disturbance of forest by wind and fire has enduring influences on regional carbon stocks and fluxes (e.g. Bradford et al. (2008), Coomes et al. (2012), Holdaway et al. (2014)) and appear to be increasing in response to changing climate.
Conventional approaches to national forest inventory are increasingly supplemented by active remote sensing technologies that provide the means to map forest structure over large spatial scales. In particular, Airborne Laser Scanning (ALS, also known as light detection and ranging or LiDAR) is now routinely used for mapping forests (e.g., Asner et al. (2010), Lefsky et al. (1999), Nelson et al. (1988), Wulder et al. (2012), Popescu et al. (2011)). The principle of ALS is that high-frequency laser pulses are emitted downwards from an aircraft, and a sensor records the time it takes for individual beams to reflect off surfaces (e.g., leaves, branches or the ground) and return to the sensor, thereby measuring the distance between the object and the airborne platform with sub-metre accuracy. Divergence of the beam means it is wider than individual leaves when it reaches the canopy, allowing some energy to penetrate through the upper canopy and reveal the foliage and the ground below, resulting in a 3D point cloud that captures the vertical structure of the forest (Lefsky et al. 1999). Conventional approaches to calculating national forest carbon stocks rely on measuring trees in inventory plots that encompass some tiny fraction of land area, and then calculating the mean carbon stock of this sample. ALS approaches also rely on permanent plot networks, this time to calibrate regression models that predict above-ground carbon density (ACD, as defined by Asner and Mascaro (2014)) from LiDAR information (Pan et al. 2011). The main advantage of using ALS is that it produces detailed maps of forest carbon rather than simply a national mean (Asner et al. 2010), and produces a vast number of ACD estimates, which supply statistical power to understand the drivers of carbon variation across landscapes (Getzin et al. 2017).
Area-based approaches are currently used to map carbon from ALS data over large scales. The principle is that estimates of ACD obtained from inventory plots are related to summary statistics derived from the ALS point cloud, such as the mean height of returns recorded within the plot areas, using regression analyses (Zolkos et al. 2013;Asner et al. 2010;Jubanski et al. 2013;Longo et al. 2016;Réjou-Méchain et al. 2015). Many summary statistics are used to construct area-based regression relationships, but one has proven particularly effective for modeling carbon from ALS: the top-of-the-canopy height, TCH, defined as the average height of the uppermost returns measured in the pixels comprising a plot (Asner and Mascaro 2014). TCH is power-law related to aboveground carbon of forests in many tropical forests (i.e. ACD = a · TCH b ). The values of a and b vary with forest type. In an attempt to generate a general model that applies across all forest types, Asner and Mascaro 2014 argued that carbon stocks should be modelled as ACD = a · TCH b · BA c · ρ d where BA is the basal area of the plot, ρ is the mean wood density of species in the plot (basal-area weighted mean), and TCH is a LiDAR-derived metric for height. This equation draws its inspiration from allometric equations used to estimate the biomass of individual trees from their basal area, wood density and height (Vincent et al. 2012). Once this equation has been fitted and parameters estimated, the next step is to fit submodel relating field BA and ρ to TCH. This two-stage modelling process results in a formula that relates ACD to TCH using seven parameters rather than two, defying the conventional wisdom that statistical models should be as simple as possible given the data. The reward of added complexity comes in the form of enlightenment about process. Asner and Mascaro show that contrasting forests, from wet and dry regions of the tropics, all have similar regression curves once adjustments are made for regional differences in the relationships between basal area, wood density and TCH (Asner and Mascaro 2014). It is also clear from this formulation that power-law models are only effective if basal area is closely related to top-of-canopy height, an assumption that is often supported in tropical studies, but not always in temperate regions Spriggs 2015). We have recently shown that model accuracy can be improved by including gap fraction alongside TCH in regression models developed for tropical regions Jucker et al. 2017), because it improves the prediction of basal area, but we have not explored whether this is also true in temperate regions.
There is current interest in developing individual-treebased approaches to make greater use of the 3D information contained in ALS data (Eysn et al. 2015;Ferraz et al. 2016;Vauhkonen et al. 2012), but these are computationally demanding. Technological advances have spawned many new algorithms for detecting individual trees in ALS point clouds, some working with rasterized surface models (e.g. Hyyppä et al. (2001), Chen et al. (2006), Yu et al. (2011)), and others with the complete point cloud (Morsdorf et al. 2004;Reitberger et al. 2009;Duncanson et al. 2014;Ferraz et al. 2016). Once the height and crown width of trees are detected, they can be used to estimate aboveground biomass of individual trees using allometric equations (Jucker et al. 2016) and these biomasses summed to give stand biomass. We recently showed that tree-centric approaches to carbon mapping perform well in Alpine coniferous forests (r 2 = 0.98 when field and ALS estimates of carbon stocks are compared) and that a correction factor can be applied to account for the small obscured trees that were invisible from the air . However, no advantage of individualbased modeling was found in a study of Malaysian tropical forests . The potential advantages of tree-based mapping include that: (i) it has a strong fundamental basis, being conceptually similar to allometric approaches used in field-based inventories; (ii) uncertainty in the estimation model is much less dependent on plot size, allowing calibration using individual trees and small plots ; (iii) narrow patches of forest with high conservation value can be mapped; (iv) growth and death of individual trees can be tracked, providing abundance data to parameterise simulation models of forest dynamics. However, these must be weighed against disadvantages, recognised as (i) high computational demand, (ii) delineation methods can only distinguish trees in the upper canopy leading to biomass underestimation, and (iii) over-or under-segmentation of large trees can result in bias Coomes et al. (2017). If these issues can be resolved, individual-based modeling could mark a fundamental shift in the way forests are monitored remotely Shugart et al. (2015). This paper explores the accuracy of area-based vs tree-centric approaches for mapping the carbon density of temperate forests in New Zealand and uses the maps to explore drivers of carbon stock variation. The forests are dominated by southern beeches, other evergreen broadleaf species and southern hemisphere conifers (podocarps). They differ in their structure from tropical forests, where TCH-based power-law models have been widely used, so there are questions about the appropriateness of this approach. The forests of concern were disturbed by cyclone Ita, one of the worst storms to hit New Zealand in decades, that brought strong winds and rain to the islands on 17 April 2014 and caused insured losses of US $42.9 million (https://hwe.niwa.co.nz/event/April_ 2014_New_Zealand_Storm). Forest damage was extensive, (e.g. Platt et al. (2014)). An ALS survey of Wellington Region was conducted partly to document the extent of damage and assess the feasibility of salvage logging. This ALS survey forms the basis of the study.

Site description and field plots
The study area included mountainous regions in the Wellington Region of New Zealand (forest map Fig. 1; altitude map Fig. 2; 175.3 • E 41.1 • S). Most mountains in the region retain a cover of natural forest, with natural sub-alpine shrublands growing at high altitude and successional forests and shrublands growing at lower altitudes following human or natural disturbance (Wiser et al. 2011); the lower lying regions are mostly occupied by agricultural land and urban areas. This study uses information from 35 permanent forest plots to generate equations for predicting ACD from ALS data, 32 of which were established as part of New Zealand's Land Use and Carbon Analysis System (LUCAS), while 3 others were established following the same protocol for other purposes. LUCAS is a national network of 20 m × 20 m plots set up for carbon accounting in indigenous forests and shrublands (Coomes et al. 2002;Holdaway et al. 2014;Holdaway et al. 2016). The plots systematically sample the country, following a 8 m × 8 km grid superimposed onto a map of indigenous forests and shrublands (see Fig. 1). Plot locations were then determined in the field using hand-held GPS devices (Garmin 62). One corner of the plot was geo-located, and N-S and E-W bearings were taken to create the plot edges. All trees > 2.5 cm stem diameter at 1.35 m height were tagged in the field, and their stem diameters (D) measured (Payton et al. 2004). The heights (H) of a subset of trees were measured and the heights of all other trees estimated from D using published allometries (Beets et al. 2012). Wood density estimates for each species were obtained from databases (Holdaway et al. 2014). The aboveground biomass of each tree was estimated from D, H and wood density of the species using functions developed specifically for New Zealand forests (Holdaway et al. 2014). Field estimates of above-ground carbon density (ACD plot , in Mg C·ha −1 ) were calculated by summing these tree biomasses, multiplying by the carbon content of wood (= 0.48 g C g −1 (Mason et al. 2012) and dividing by plot area A. The plots had dimensions of 20 m × 20 m on the ground's surface, while A ≤ 400 m 2 because it is calculated on a horizontal plane (i.e. area as seen on a map). Plots had a mean slope of 29.8 • , resulting in A being on average 14% less than field-measured areas. To calculate A, a GIS polygon was created from the four distances and bearing Fig. 1 The Wellington Region of New Zealand, which was surveyed by airborne LiDAR in 2013. Land-cover types, extracted from LCDB-IV, include natural old-growth forests (dark green), secondary forest (light brown), while plantations, urban and agricultural are left uncoloured. The location of LUCAS carbon monitoring plots are shown measured in the field, applying a Bowditch rule correction (Jones 1972) when those bearings and distances did not result in a closed polygon, and calculating the area of the polygons in the horizontal plane. Preliminary analyses of these plots show they are comprised primarily of Weinmannia racemosa (25% by basal area), southern beech species (Lophozonia menziesii 19%, Fuscospora fusca 17%, Fuscospora solandri 4%), conifers in the podocarpaceae (primarily Prumnopitys ferruginea 4%), and tree ferns (e.g. Cyathea smithii 3%). Successional woodlands are a mix of native species, such as Kunzea ericoides (2%), invasive species such as gorse (Ulex europaeus, 2%) and several broadleaf indigenous hardwoods. A further 84 species of tree and shrub comprise the remaining 21% of basal area.

Airborne laser scanning (LiDAR) survey
The ALS (LiDAR) data used in this study were acquired over the entire Wellington Region in 2013, surveying an area slightly larger than 8500 km 2 (Müller et al. 2015). The majority of the ALS survey was performed in early 2013 with some additional aircraft flights later in 2013 and 2014, depending on weather and data quality considerations. The LiDAR scanner was an Optech ALTM 3100EA flown at a nominal height of 1000 m. The target survey point density was 1.73 ppsm with 50% swath overlap to ensure the minimum raw data specification of 1.3 ppsm and a vertical accuracy of ± 0.15 m (i.e. ± 1 σ ). The 1261 LiDAR flight lines of raw point cloud data were merged, tiled, and automatically classified into ground and non-ground returns using the open-source LiDAR processing software, Sorted Pulse Data Software Library (Bunting et al. 2011). The tiles of ground classified points were then interpolated and mosaicked into a 1-m resolution digital terrain model (DTM), again using SPD (see Fig. 2). Similarly, the tiles of non-ground classified points were interpolated and mosaicked at 1-m resolution from which the DTM is subtracted to form a Canopy Height Model (CHM). The term CHM is nominal as it is primarily derived from canopy returns but will contain other identified non-ground features such as buildings in urban areas. A 5-metre median Gaussian smoothing filter was applied to the non-ground returns to discard outliers in the model. From the CHMs we calculated two metrics for each of the permanent field plots: top-of-canopy height (TCH, in m) and canopy cover at 10 m aboveground (Cover 10 ). TCH is the mean height of the pixels which make up the surface of the CHM within a particular region of interest (i.e. 20 m × 20 m plots in this paper). Canopy cover is defined as the proportion of area occupied by crowns at a given height aboveground (i.e., 1 -gap fraction). Cover 10 was calculated by creating a plane horizontal to the ground in the CHM at a height of 10 m aboveground, counting the number of pixels for which the CHM lies above the plane, and then dividing this number by the total number of pixels in the plot. A height of 10 m above ground was chosen because previous work has shown that the optimal height is about 1/2 of the mean TCH ). R packages raster (Hijmans 2015), rgdal ) and rgeos (Bivand and Rundel 2016) were used to extract LiDAR data from specific plots and calculate metrics from these data. Slope and aspect were calculated from the highresolution DTM (Fig. 2) using the terrain function within the raster package of R.

Area-based regression analysis to estimate carbon density
Models were fitted using maximum likelihood estimation. The simplest model we fitted is given by Eq. 1: where the residual error is ∼ N 0, σ 1 + σ 2 2 · TCH . The second variance term σ 2 2 in the error distribution allows uncertainty in ACD plot to increase with increasing forest height, as commonly observed. The parameters were estimated by maximum likelihood estimation (MLE), using the mle2 function in the bbmle library of R. MLE is a statistical curve fitting approach that has greater flexibility than least-squares regression techniques. While this particular model could have been fitted by regression with loglog transformed data and size-invariant uncertainty, more complex models (e.g. 3) cannot be fitted by regression, and so MLE was used for all model fitting for the sake of consistency. The error structure ∼ N 0, σ 1 + σ 2 2 · TCH was used consistently in all models.
Model accuracy was assessed by leave-one-out crossvalidation (Arlot and Celisse 2010). Plots were omitted one at a time, and their ACD values predicted by MLE using all other data. Goodness of fit was estimated using the normalised root mean square error RMSE% = (pred − obs) 2 /N 0.5 · 100/obs, where obs and pred are observed and predicted values, and obs is the mean of the observations. Normalised bias was calculated as bias% = (pred − obs) · 100/obs. The coefficient of determination was calculated as . Cross validation also provided uncertainty estimates for model coefficients (i.e. jackknife estimates of variance).
Having fitted the TCH-only model, Asner and Mascaro's modelling approach was followed (Asner and Mascaro 2014). First, the following model was fitted: The next stage is to fit the following submodels: Asner and Mascaro report coefficients from a model fitted to 14 contrasting forest types after correcting for regional difference in relationships between BA, ρ and TCH their "general model", ACD plot = 3.8358· TCH 0.2807 BA 0.9721 plot ρ 1.376 plot but recommend fitting local relationships to achieve the more accurate predictions (regional model). We explored both approaches. Finally, a model was fitted that included canopy cover (i.e. Cover 10 ) as well as TCH (see Jucker et al. (2016) for details).
A previous analysis of random measurement and modeling uncertainties associated with estimating forest carbon from LUCAS plots reported that plot-to-plot variation was the greatest source of uncertainty (SEM = 9.1% of mean ACD) while the propagation of all other errors resulted in only a 1% increase in overall uncertainty (i.e. SEM = 0.1%) (Holdaway et al. 2014). Inaccuracies in the geolocation of plot corners can introduce significant uncertainty into the relationship between estimated ACD and LiDAR metrics, particularly when field plots are small (Gobakken and Naesset 2009). This could be a particularly major issue with New Zealand data, because the field plots are small (0.04 hectares) and geopositioning was not conducted with a survey-grade (i.e. differential) GPS. To explore the extent of this problem, every plot was shifted by a random distance in both north-south and east-west directions (i.e. N(0, σ = 5) in both directions) from the recorded location, the LiDAR TCH calculated from canopy-height-model pixels that lay within the boundary of the new plot, and a regression line fitted to these data. A standard deviation of 5 m for the random shift corresponds approximately to accuracy of the Garmin 62 GPS used. This random shifting of plots was repeated 100 times to give a distribution of regression lines from which the uncertainty created by geolocation inaccuracies could be assessed, as well as estimates of uncertainty in plot-level TCH predictions.

Tree-centric mapping of carbon density
Delineation of individual tree crowns was performed using the itcSegment R package described in Dalponte and Coomes (2016), Coomes et al. (2017) and available through CRAN, which is an adaptation of the approach originally developed by Hyyppä et al. (2001). In essence, itcSegment finds local maxima in the canopy height model, using a moving window approach, and then finds the crown associated with each local maximum by searching locally for pixels of similar height, using a "region growing" subroutine based on various rules about tree geometry. The size of the window varies with forest height, in recognition that larger trees have wider crowns so a wider search window is needed. To assess the accuracy of delineation, histograms were produced of number of stems observed in the field versus those delineated in different height classes within the forest plots.
The crown width and height of delineated trees were used to estimate tree biomasses, and, by summation, the aboveground biomass of the 32 field plots. A global compilation of studies has reported that the following allometric function gives an unbiased estimates of individual tree biomass: AGB = α(CD · H) β , where CD is crown diameter and H is tree height of harvested trees, and α and β are parameters (Jucker et al. 2016). For angiosperm trees, which dominate the New Zealand forests in the Wellington Region, β was estimated at 2.013 (Jucker et al. 2016). Here, we estimated the CD of all delineated trees by assuming a circle equal in area the delineation polygon (i.e. CD = 2 √ CA/π , where CA is the delineated crown area. The total above-ground biomass of a plot is the sum of biomass estimates of delineated trees, multiplied by a correction factor accounting for unobserved trees hidden beneath the overstorey . Maximum likelihood estimation to fit the following relationship plot biomass and individual tree dimensions: Here parameter h combines the scaling coefficients from the tree biomass models with the subcanopy tree correction factor. Finally, we explored whether the mean biomass of overstorey trees is closely related to ACD by modelling

Removing an outlier
An outlier was identified in the LUCAS plot dataset and removed from further analyses. The predicted ACD for plot CN95 was almost double the field-estimated value, irrespective of which area-based model was fitted (e.g. 251 vs 551 Mg C·ha −1 for the TCH-only model). Bootstrap values obtained by cross-validation revealed that coefficients estimated when plot CN95 was removed from the MLE model were quite distinct from coefficients estimated after removing other plots, underscoring that this plot was unusual. Furthermore, random shifts of the plot's location around the recorded position resulted in larger variation in TCH estimates for CN95 than found in any other plot; this variation is expected when a large tree sits just outside a plot boundary and may, or may not, be included in TCH estimates when plot corners are randomly moved. This edge tree is apparent in the LiDAR imagery.

Carbon density mapping and trend analyses
An aboveground carbon density map was created by using the 1 × 1-m resolution CHM raster to calculate TCH and forest cover within 19 m × 19 m plots, a size chosen because the 20 m × 20 m LUCAS plots have a mean horizontal area of 335 m 2 . The best-supported area-based estimation model was used to predict ACD for each grid cell. The New Zealand Landcover Database comprises maps of land-cover types (see Fig. 1). Working with version 4.1 of the database, released February 2017 (https://lris.scinfo.org.nz/layer/423-lcdb-v41-landcover-database-version-41-mainland-new-zealand/), the carbon map was restricted to old-growth forest ( i.e., LCDB class 69 = indigenous forest, comprised mostly of natural forests at various phases of regeneration following disturbance by wind, earthquakes, bark-beetle outbreaks etc.) and secondary forests recovering from human disturbance (i.e. combining kanuka and/or manuka woodlands = LCDB class 52 and broadleaf indigenous hardwoods = LCDB class 54). Variation in ACD with altitude and aspect was calculated and displayed graphically.

Area-based modeling approaches
The simplest model, containing only TCH as an explanatory variable, has an R 2 of 0.75, similar to that achieved by more complex models. The regional model of Asner and Mascaro was: With 1 standard deviations of the coefficients of ± 0.41, ± 0.039, ± 0.029, ± 0.22 respectively (R 2 =0.86, bias = − 4%). The sub-models were: demonstrating a close relationship between TCH and BA, but no discernible relationship between TCH and ρ. Submodels 5 and 6 were used to predict BA and ρ from TCH for each plot (i.e. BA, ρ), and these predictions entered into Eq. 4 to predict ACD. This two-stage approach gave an R 2 = 0.76, slightly higher than the simple power-law model (Model 2 in Table 1). Predictions made from Asner and Mascaro's general model, using published coefficients obtained from tropical forests and New Zealand's BA-TCH and ρ-TCH relationships, had a similar R 2 value to the regional model (0.85) but showed much greater bias (-13% vs -4%). Introducing canopy cover (Cover 10 ) into the model did not improve goodness of fit (model 3 in Table 1), because Cover 10 and TCH were closely correlated (r = 0.82) and including both led to redundancy. To evaluate errors arising from geolocation inaccuracies, every plot was shifted by random distances from its recorded location, a ACD-prediction model was calculated using TCH values estimated for these shifted plot locations, and ACD values predicted for the plots using that model. The mean and standard deviation of ACD estimated from 100 such randomisations is shown in Fig. 3. The figure shows that geolocation errors vary with plot ACD: values of σ/μ × 100 are about 20% at low ACD is low, falling to 2.5% at intermediate values, and rising to 8% when ACD is high.

Tree delineation approach
The delineation algorithm was able to detect overstorey trees (i.e. those >12 m in height), although it tended to over-segment trees (i.e. to subdivide single trees), finding 16% more large trees than were actually present in the plots (Fig. 4). It was far less successful at detecting shorter trees, many of which were in the understory and thus invisible to ITCsegment (Fig. 4).
The model based on summing estimated biomasses of delineated trees (model 4 in Table 1) performed less well than the area-based approaches, having the lowest R 2 and highest RMSE % of all the models. However, the model based on mean(CD · H) values (model 5) was the best performing of all models.

Carbon mapping
With such an large dataset, computational efficiency is essential. Thus we opted to map carbon using Model 1 in Table 1, even though delineation methods produce slightly more accurate predictions. There were 4.5 million measurements of top-of-the-canopy height within the natural old-growth and secondary forests of Wellington Region, each estimated from the ALS canopy height model averaged over 25 m × 25 m. Figure 5 shows aboveground carbon density for oldgrowth and secondary indigenous forests across the entire Wellington Region, with a mountainous region of Tararua Forest Park shown in greater detail. It is apparent from these maps that the low-altitude forests around the coast store less carbon than those in the mountains.
The ACD probability distribution in old-growth forests is bell-shaped with a mean ACD of 301 Mg C·ha −1 while secondary forests have a highly skewed distribution, with a mean ACD of 105 Mg C·ha −1 (Fig. 6). Combining these datasets, the regions natural forests are estimated to store 216 ± 169Mg · ha −1 of carbon (mean ± 1 SD). Note that this standard deviation is based on spatial variation in ACD without propagation of errors associated with modelling or considering the effects of spatial autocorrelation, and would not be an appropriate estimate of uncertainty for reporting regional carbon storage. Terms in the models include basal area (BA) and wood density (ρ) recorded in plots, Top-of-the-Canopy Height (TCH) and residual variation in canopy cover (Cover 10 ) estimated from ALS data, as well as Crown Diameter (CD) and Height (H) of individual overstorey trees obtained by segmentation of the ALS dataset. Mean ± 1 SD of coefficient values are given. Goodness of fit is assessed using R 2 , RMSE (%) and bias (%) Fig. 3 Effects of introducing random geolocation errors into plot locations on ACD estimates for those plots. The mean and standard deviation of 100 randomisation are shown Strong trends in carbon density were observed with altitude and aspect (Fig. 7). Natural old-growth forest is tallest (and therefore has greatest carbon density) at mid-elevation, but this pattern is not seen in secondary forests. For both forest types we see greater carbon storage on east-facing slopes, which are protected from westerly storm damage, than for forests on west-facing slopes. Wind is expected to be greatest in the mountains, and it putative effects do indeed appear to be stronger at mid-to high-elevations than in the lowlands. Note Fig. 4 Evaluation of the accuracy of tree delineation by itcsegment; the number of trees observed in the field in four height classes are compared with numbers delineated from the ALS point cloud the standard errors of the mean, calculated as SEM = SD/(samplesize) 0.5 , are all very small because of the large sample sizes, so conventional inferences statistics such as F-tests give highly significant differences for all comparisons (not shown).

Airborne laser scanning of aboveground carbon density
Airborne Laser Scanning (LiDAR) provides precise measurements of forest structure from which forest carbon can be monitored, providing wall-to-wall mapping at higher resolution than any other remote sensing product, and producing vast datasets with which to study drivers of forest structure and dynamics with strong statistical power. Using ALS maps, the regions natural forests are estimated to store 216 ± 169 Mg · ha −1 of carbon (mean ± 1 SD), which is lower than estimated directly from the LUCAS plots (243 ± 173 Mg C·ha −1 ). By far the greatest source of uncertainty in field estimates of regional carbon storage is "sampling error", i.e. inherent variation between plots, with "model uncertainty" adding only 1% to uncertainty (Holdaway et al. 2014). We would expect ALS-estimated carbon density of plots to have greater uncertainty than field estimates. The coefficients of variation indicate that uncertainty is greater in the ALS estimate (0.78 for ALS, 0.71 for plot-based estimates) and uncertainty for ALS-based estimates would have been even higher if we had propagated the considerable uncertainty arising from geolocation of small plots (Chen et al. (2015) and see Fig. 3). However, the advantage of ALS is that it provides 4.5 million ALS samples, compared with just 31 from the field surveys. This enormous sample size outweighs the disadvantage of plot-level increases in uncertainty when it comes to estimating uncertainty in the mean. The standard error of the mean is tiny for ALS-estimated carbon compared with plot estimates (0.08 vs 32 Mg C·ha −1 ). Thus ALS delivers estimates of regional carbon storage with very narrow confidence intervals. It is important to realise, though, that ALS-estimated ACD may still be highly biased (i.e. contain systematic error) if there are biases in the estimation equations used to predict plot-ACD from field measurements and ALS measurements. In particular, ALS measure the sizes of trees and draw inferences from those sizes, but wood density can vary systematically across regions and at present that variation is not measured from aircraft or satellite (Avitabile et al. 2016). Hyperspectral remote sensing is able to classify forests with greater accuracy than any previous remote sensing approach ) and by linking each forest type with a community-weighted mean wood density value, it should in future be possible to reduce uncertainties in ACD arising from wood density variation.

Area-based vs tree-centric approaches
A simple estimation model that predicts ACD as a powerlaw function of top-of-the-canopy height (TCH) was almost as accurate as elaborate tree-centric approaches that require much more computational power (Table 1). The explanation for its success lies in the close relationship between TCH and forest basal area (Eq. 5; see ). Using the two-stage approach adopted by Asner and Mascaro 2014 helps put New Zealand's forest in a global perspective. In Fig. 8, we plot TCH-ACD relationships for New Zealand's forests against the relationships observed in 14 tropical regions, including areas of rain forest and dry forest (Asner and Mascaro 2014). To our surprise, we find that New Zealand's forests have a much greater carbon density for a given height, than any of these tropical forests (Fig. 8). The reason why they attain such high carbon densities is not the result of high community-weighted wood densities in New Zealand (mean = 0.49 g·cm −3 ); these are at the lower end of the range of tropical wood densities (0.48 -0.69 g·cm −3 ; see Fig. 9). Nor does it arise from a fundamentally different scaling relationship (i.e. Eq. 2), because the regional model fitted to the New Zealand data had similar goodness of fit to the general model reported by Asner and Mascaro. In fact, the explanation lies in the fact that New Zealand forests pack in a much greater basal area for a given height than tropical forest; for the relationship is BA = 6.44 · TCH, whereas for tropical forests the slope ranges from 1.13 to 4.37, with a mean of 1.88 (Fig. 9). The fundamental explanation for this dense Fig. 6 Probability distribution of aboveground carbon density (ACD, Mg C · yr −1 ) in natural and secondary forests, based on 2.6 million and 1.9 million measurements, respectively, of top-of-the-canopy height by airborne laser scanning. a Indigenous forest. b Secondary forest Fig. 7 Influences of aspect and altitude on aboveground carbon density (ACD) of (a) old-growth and (b) secondary forests in the Wellington Region of New Zealand, based on 4.5 million measurements of top-of-the-canopy height made by airborne laser scanning. Mean (± 1 standard error of the mean, as black arrows) are shown, but standard error of the mean are narrower than symbol widths because of high sample size. Note that virtually no secondary forests occur at high altitude packing of New Zealand forests remains unresolved, but could be related to strong winds in the region that keeps forests short but do not diminish their basal area. It would be fascinating to follow up this work by generating the relationships observed in Fig. 9 and for various regions exposed to strong and weak winds.
The tree-centric function that used mean(CD · H) as its explanatory variable was the most accurate of all Fig. 8 Comparison of the relationship between TCH and ACD observed in New Zealand's natural forests with those observed in 14 tropical regions, including wet and dry forests (Asner et al. 2014) the models tested, but the function that summed the biomasses of individual delineated trees -and is based on the fundamentals of forest inventory -performed poorly. Our tree-centric approach failed to detect many small trees, and over-segmented some large ones. The problem of undetected understorey trees can be resolved by applying a correction factor (Dalponte and Coomes 2016) but over-segmentation can create large uncertainties ). In the future, the answer may lie in more sophisticated algorithms (e.g. Ferraz et al. (2016)). Others have found mean-size metrics extracted from treecentric models are strong predictors of ACD (Singh et al. 2016); it seems that taking averages reduces the influence of segmentation issues described above.

The influences of wind storms on forest carbon density
Disturbance by wind, fire, snow storms, earthquakes, and insect outbreaks can have very strong influences on carbon cycling (Allen et al. 2010;Harcombe et al. 1998;Wardle 2002). In old-growth forests, the death of old trees creates large gaps that can take many years to refill (Zeide 2005), and understanding the process of gap creation and filling is essential to predicting carbon cycling (Korner and Körner 2003;Seidl et al. 2011;Coomes et al. 2012). ALS provides fresh insights into the effects of disturbance on forests. The damage wreaked by storms such as cyclone Ita, which hit New Zealand in 2013, is immediately obvious in their aftermath, in the form of wind-thrown trees and ripped-out tree crowns. The long-term influences of living in exposed versus sheltered sites are less obvious. ALS provides an ideal tool for exploring influences of environment on forest structure and dynamics, because of the unprecedented sample sizes it provides. Wind has a major influence on forests in the lower North Island of New Zealand (Fig. 7) (Zotov et al. 1938;Elder 1965). Strong north-westerly winds, coming off the Tasman sea, reduce forest stature in parts of the Tararua Range and other coastal forests (PJ Bellingham, pers. comm., (Zotov et al. 1938)). Funneling of the prevailing north-westerly winds through the Cook Straits magnifies their strength. The fact that ALS picks up differences between sheltered east-facing sites and exposed west-facing sites in both old-growth and secondary forests is reassuring, suggesting that wind is the general driver of the observed patterns, although of course it is impossible to rule out other drivers in this correlative analysis. The ALS survey provides an insightful "snapshot" into the effects of wind in region, but we know the long-term effects of wind are complex. For instance, a descriptive paper dating back to the 1930s reports that an unusual south-easterly storm caused extensive damage to forests in the Tararua Range in 1936 (Zotov et al. 1938). The extent to which prevailing westerly storms vs unusual south-easterly storm affect forests can only be understood by repeated surveying of the region.
A hump-shaped effect of altitude of ACD is also apparent. The decline in ACD at high altitude is consistent with ecological understanding of resource limitation. Towards the top of temperate mountains, productivity is limited by low temperatures, N-limitation resulting from there being fewer degree days for mobilisation, and P-limitation, e.g. (Harcombe et al. 1998). The explanation for low ACD in the low-altitude forests is most probably related to humans, because many areas have been modified through clearing, fire, invasive species, and logging, leaving a patchy distribution of undisturbed forests. However, we cannot rule out the influence of natural process in driving this low biomass, because warmer more fertile soils in the lowlands give rise to faster growing trees that turnover more frequently (Ferry et al. 2010;Coomes et al. 2005). More elaborate analyses using multiple environmental layers is not possible here, but might shed more light on the processes driving the observed patterns.
To conclude, this paper has illustrated how highresolution carbon maps produced by ALS are powerful for evaluating the environmental drivers of forest structure, demonstrating that forests are shorter on exposed west-facing slopes in New Zealand. Simple estimation models based on top-of-the canopy height were found to be almost as accurate as state-of-the-art tree-centric approaches that require more computing power. By working with Asner and Mascaro's approach, we showed that New Zealand forests have much higher ACD than tropical forests of comparable size, because they have extremely high basal areas for their height. It would be fascinating to expand this analysis to other regions, to gain a better understanding of the drivers of forest structure globally.

Conclusion
Simple estimation models based on top-of-the canopy height are almost as accurate as state-of-the-art treecentric approaches, which require more computing power. High-resolution carbon maps produced by ALS provide powerful datasets for evaluating the environmental drivers of forest structure, such as wind.