Comparing 3D Point Cloud Data from Laser Scanning and Digital Aerial Photogrammetry for Height Estimation of Small Trees and Other Vegetation in a Boreal–Alpine Ecotone

Changes in vegetation height in the boreal-alpine ecotone are expected over the coming decades due to climate change. Previous studies have shown that subtle changes in vegetation height (<0.2 m) can be estimated with great precision over short time periods (~5 yrs) for small spatial units (~1 ha) utilizing bi-temporal airborne laser scanning (ALS) data, which is promising for operation vegetation monitoring. However, ALS data may not always be available for multi-temporal analysis and other tree-dimensional (3D) data such as those produced by digital aerial photogrammetry (DAP) using imagery acquired from aircrafts and unmanned aerial systems (UAS) may add flexibility to an operational monitoring program. There is little existing evidence on the performance of DAP for height estimation of alpine pioneer trees and vegetation in the boreal-alpine ecotone. The current study assessed and compared the performance of 3D data extracted from ALS and from UAS DAP for prediction of tree height of small pioneer trees and evaluated how tree size and tree species affected the predictive ability of data from the two 3D data sources. Further, precision of vegetation height estimates (trees and other vegetation) across a 12 ha study area using 3D data from ALS and from UAS DAP were compared. Major findings showed smaller regression model residuals for vegetation height when using ALS data and that small and solitary trees tended to be smoothed out in DAP data. Surprisingly, the overall vegetation height estimates using ALS (0.64 m) and DAP data (0.76 m), respectively, differed significantly, despite the use of the same ground observations for model calibration. It was concluded that more in-depth understanding of the behavior of DAP algorithms for small scattered trees and low ground vegetation in the boreal-alpine ecotone is needed as even small systematic effects of a particular technology on height estimates may compromise the validity of a monitoring system since change processes encountered in the boreal-alpine ecotone often are subtle and slow.


Introduction
The world's climate will undergo distinct alterations over the coming decades leading to rapid changes in basic growth factors, such as temperature and precipitation. This will influence the alpine and tundra transition zones of the boreal forests. These areas are characterized by steep temperature-productivity gradients because the trees exist close to their tolerance limit in terms of temperature. Even a moderate increase in temperature may therefore lead to a rapid increase in growth of existing trees [1,2] as well as colonization of treeless areas and migration of the tree lines [3]. Other drivers of vegetation change also exist at the forest-tundra and forest-alpine interfaces, like herbivory [4][5][6][7]. Summer farms and grazing by domestic animals have been common in many montane areas. In areas where summer farming and grazing animals have previously limited tree growth, the tree line and the forest may expand towards higher elevations and latitudes when this activity is reduced. Migration of the alpine and northern tree lines will influence future carbon pools. A need therefore exists to monitor vegetation changes in these areas [8].
Airborne laser scanning (ALS) has been proposed as a technique to monitor subtle changes near the tree lines, such as colonization of treeless areas [9], to assist in prediction and estimation of tree height [9][10][11], and to estimate subtle changes in tree and vegetation height [12]. Numerous studies [9,11,[13][14][15][16][17][18] have shown that ALS data with point densities of 7-11 points m −2 may be applied to detect individual pioneer trees in the alpine tree line. With such pulse densities, about 90-100% of the trees with heights greater than 1 m are likely to be hit by laser pulses resulting in echoes with height values greater than zero, i.e., located above the terrain surface. Echoes with heights >1 m are mainly tree echoes. However, because laser beams of pulse lasers will tend to penetrate into the canopies before an echo is triggered, even the maximum of the recorded echo heights will typically underestimate the height of small trees. For example, Ref. [9] reported an underestimation of tree height by 0.43 m to 1.01 m, depending on tree species and tree height. The tree heights in their dataset ranged from 0.11 m to 5.20 m, which implies that the smaller trees will tend to have recorded laser echoes with height values equal to zero even when a tree is hit by a laser pulse. This will limit the sensitivity of ALS as a tool for early detection of recently established trees.
Three-dimensional (3D) photogrammetric point data produced in digital aerial photogrammetry (DAP) from imagery acquired by aircraft and unmanned aerial systems (UAS) have in recent years become a viable alternative to 3D data from ALS in forest and vegetation studies. DAP is an alternative to ALS in, for example, operational inventory of forest resources (e.g., [19,20]), although with somewhat lower accuracy than by use of ALS [21], but still economically competitive to ALS in forest inventory when both cost and utility of the data are taken into account [22]. 3D data from DAP based on UAS imagery have shown great promise in terms of precision of estimates of parameters such as volume and biomass in different forest types -from boreal forests ( [23]) to dry African savannah [24].
It is usually cheaper to acquire 3D data from DAP than from ALS and UAS offer greater flexibility for small areas than acquisition from larger airborne platforms. Use of 3D data extracted by DAP from imagery acquired by UAS may therefore be a viable option to assist in detection and estimation of height of small trees in smaller areas in the forest-alpine or forest-tundra ecotones. Further, because a passive sensing technique like optical imaging may better capture the properties on the outer surface of a tree canopy than lasers, which tend to penetrate the canopy, DAP may in fact offer additional advantages compared to ALS for tree detection and estimation of height [19].
Still, to our knowledge, there is no scientific evidence of the performance of 3D data from DAP for small tree detection and height estimation in the forest-alpine or foresttundra ecotones. A recent study by Ref. [25] may, however, give an interesting perspective on the potential usefulness of 3D data from UAS DAP for estimation of properties of small trees. Ref. [25] estimated mean tree height using DAP from imagery acquired with UAS over young forest stands under regeneration. The mean plot and stand height in their dataset ranged from 0.5 m to 13.0 m with an average value of 2.5 m. They found that the RMSE of the mean height estimate produced with assistance of 3D DAP data was substantially smaller than obtained with ALS. Although this study only addressed mean values of groups of trees (plots, stands) rather than individual trees, their findings may suggest that DAP also may be useful to assist in quantifying properties of individual small trees. In an assessment of UAS 3D point clouds from laser and DAP for individual tree height estimation over tree plantations with an average tree height of 2.6 m, Ref. [26] reported that DAP underestimated height to a greater extent than laser and that the accuracy was greater for laser than for DAP. However, their 3D data had an exceptionally high point density (443-939 points m −2 for DAP and 325-649 points m −2 for laser), which makes comparisons with previous findings based on data from airborne platforms and coarser resolution imagery difficult.
There are currently no commercial interests associated with small trees and other vegetation in the tree line ecotones, as there are in the productive forests where small trees represent young forest in an early stage of a rotation after clear-felling (cf. [25,26]). The primary purposes of quantifying and monitoring small trees and other vegetation in the tree line ecotones is therefore partly to keep track of how climate change affects this climate sensitive ecological environment and to enable consistent analysis of the net climate feedback of tree line migration caused by changes in biomass and soil carbon pools and changes in albedo. Ref. [12] proposed several statistical estimators to estimate changes in height of trees and other vegetation using bi-temporal data from ALS in a boreal-alpine ecotone. The estimators were applied to an observation period of six years and it was demonstrated that statistically significant increases in height could be found for relatively small monitoring units (1.5 ha primary monitoring units). The estimation framework was proposed as an operational methodology that could be applied to monitoring over vast tracts of land, and it was based on a so-called model-dependent approach to statistical inference. The method relies on bi-temporal 3D data from remote sensing and temporally consistent ground observations of heights of trees and other vegetation. Because the precision (confidence interval) of the height change estimates will determine the sensitivity of the method to detect subtle changes in height, it is important to quantify to what extent the source of the 3D information (ALS vs. DAP) influence the precision of the estimates.
The current study focused on the boreal-alpine ecotone in particular. The objectives were twofold. (1) We assessed and compared the performance of 3D data from ALS and DAP for prediction of tree height of small pioneer trees and evaluated how tree size and tree species affected the predictive ability of the two types of 3D data. (2) We compared the precision of vegetation height estimates (trees and other vegetation) across the chosen study area using 3D data from ALS and from DAP using the estimators proposed by Ref. [12]. As part of the latter objective, we also evaluated the different sources of uncertainty in the model-dependent mean square error estimators for vegetation height at different spatial scales. It should be noted that this analysis focused on height estimation rather than height change estimation because 3D data from DAP to be compared with ALS data were available just for one point in time. An operational methodology for change estimation could indeed exploit combined bi-temporal ALS and DAP data, but we wanted to quantify the effect of each of them on the precision of estimates without confounding the effects of the two 3D acquisition techniques.

Study Area
The study area is located in the municipality of Rollag in southern Norway (60 • 0 N 9 • 01 E, 910-950 m above sea level) ( Figure 1). The entire study was conducted within a 200 m × 600 m rectangle (12 ha). The work took place in the boreal-alpine tree line, which at this location was around 900-940 m above sea level. The main tree species in the trial area are Norway spruce (Picea abies (L.) Karst.), Scots pine (Pinus sylvestris L.), and mountain birch (Betula pubescens ssp. czerepanovii). The total stem density when the study area was established in 2006 was estimated to be 97 trees ha −1 , of which only 15 trees ha −1 were taller than 2 m [9].

Overview
This study comprised field data from two complementary datasets which both have been subject to analysis in previously published work but updated with new and original measurements for the purpose of the current study.
First, we selected and georeferenced individual trees that could be used as groundreference (1) for analysis of various types of remotely sensed data and the performance of such remotely sensed data to identify individual pioneer trees; (2) for detection of subtle changes at tree level over time using remotely sensed data; and (3) for development of methods for operational monitoring of vegetation changes over time by assistance of remotely sensed data, which subsequently could be adapted to monitoring over vast tracks of land. (4) An overarching purpose was to establish a long time series of individual tree data that could be used as reference for biological studies of changes in the tree cover in the boreal-alpine ecotone caused by anthropological drivers of change. The time series was established in 2006 and is still maintained, see details in Section 2.2.2. The individual tree dataset was subject to analysis under objectives #1 and #2.
Second, ALS data have to date been the primary source of remotely sensed data under study of pioneer trees. Because the sensitivity of ALS data for early detection of emerging pioneer trees depends on the accuracy of the digital terrain model (DTM) used for normalization of the ALS vegetation echoes, we collected ground reference points with known elevation and ground properties (e.g., terrain form and type of ground vegetation; see details in Section 2.2.3). The primary purpose was to assess systematic errors and accuracy of DTMs constructed from ALS data under different acquisition strategies, such as flying altitudes and pulse repetition frequencies. This dataset also contained valuable information on other vegetation than trees which complemented the individual tree dataset. It was subject to analysis under objective #2. , which were centers of the 25 m radius plots (gray circles) arranged along four sample lines; green is defined as forest according to the official N50 topographic map series; light yellow is above the tree line; the black triangle is the reference point (base station).

Individual Tree Data
The individual tree dataset was established in 2006 [9] and re-measured in 2012 [12] and 2017. In 2006, the point-centered quarter sampling method (PCQ; [27]) was used to select individual trees for the study. The tree sampling was conducted according to PCQ at 40 systematically distributed sample points within the 200 m × 600 m rectangle (Figure 1).
It should be noted that due to time constraints, only four points were measured on line #4, see Figure 1.
At each point, the closest tree in each of four height classes (<1 m, 1-2 m, 2-3 m, >3 m) within each of the four quadrants around the points defined according to the cardinal directions, i.e., the NE, SE, SW, and NW quadrants, was selected. Thus, a maximum of 16 trees were selected at each point. It should be noted that the selection of trees was restricted to a maximum distance from each sample point of 25 m [28].
The stem positions of the trees were recorded with a real-time differential global positioning system (GPS) and global navigation satellite system (GLONASS) receiver, with a local reference receiver for differential correction located within the study area ( Figure 1) at a national reference point of the Norwegian Mapping Authority. The expected accuracy was 3-4 cm [9]. For each tree, tree species, tree height, stem diameter at root collar, and crown diameter in two perpendicular directions (N-S and E-W) were recorded. In total, 342 trees were selected, ranging from 0.11 to 5.20 m in height. Details regarding the field work in 2006 can be found in Ref. [9]. The trees took many different forms, including distinct and solitary trees, groups of trees-for spruce often as krumholtz, and birch appearing as solitary trees as well as in the form of tall scrubby vegetation ( Figure 2). When the field work was repeated in 2012 and subsequently in 2017, the already defined 40 sample points were once again identified with real-time GPS + GLONASS. For each of the 40 points, the PCQ sampling was then conducted independently of the sampling in 2006, but according to the same protocol. Many of the trees selected in 2012 and 2017 were the same as those measured in 2006. During the 2017 campaign, which took place in the period 8 August to 21 September, we identified all trees that had been measured in 2006 and 2012 [12] and which were still alive, including those that were not selected for the 2017 PCQ sample. Thus, for the current study, we measured and analyzed all trees selected into the 2017 PCQ sample in addition to trees measured in 2006 and 2012 and which were alive in 2017. The individual tree recordings in 2017 followed the same protocol as in 2006. In total, 532 trees were recorded, ranging between 0.05 and 6.60 m in height ( Table 1). The field-recorded heights of the trees were designated h. The 2017 data have not been published before.

Ground Reference Points
The ground reference point (GRP) dataset was acquired in August 2010 [29]. A total of 440 GRPs were distributed at 5 m intervals in the N-S direction through the center points of the 40 PCQ plots using a measuring tape and a hand-held compass. Real-time GPS + GLONASS was used to record the coordinates of each GRP. For some points the positioning was unreliable due to poor radio link between the base and rover receivers. Thus, 426 of the initial 440 points were available for analysis ( Figure 1). For each GRP, three variables were recorded. They were "terrain form", "terrain surface", and "vegetation height". Vegetation height was recorded according to three mutually exclusive height classes, namely <0.10 m, 0.10-0.20 m, and >0.20 m. Terrain surface was recorded according to three mutually exclusive classes, namely "rock/bare", "lichen/heather", and "green vegetation". Heather comprised common heather (Calluna vulgaris), crowberry (Empetrum nigrum L.), cowberry (Vaccinium vitis-idaea), mountain heath (Phyllodoce caerulea), and alpine azalea (Loiseleuria procumbens), but not bilberry (Vaccinium myrtillus L.). The latter was classified as green vegetation.
In the analysis addressing objective #2, a model-dependent approach to inference was adopted, by which the height of trees and other vegetation was estimated for various domains of the study area. Under model-dependent inference, approximate unbiasedness of the estimators can only be assured if the model used for prediction is correctly specified for the domain of application. Since the vegetation in the study area is a mix of scattered trees of different species and other vegetation, the tree dataset alone would probably not warrant appropriate models for prediction of height of all vegetation, see detailed discussion in Ref. [12]. The GRPs with recorded vegetation height were considered complementary to the tree data for combined modelling of vegetation height for all vegetation, trees included.
Since the vegetation height was recorded in ordered classes only, we assigned a height value of 0.05 m to GRPs in the class <0.10 m. An exception was made for GRPs that were classified as "rock/bare" (see example in Figure 2E). They were assigned the value 0 m. A height value of 0.15 m was assigned to GRPs in the class 0.10-0.20 m. Because we did not know the upper height limit in the class >0.20 m, GRPs with vegetation height >0.20 m were discarded. In some cases, trees actually constituted the vegetation cover for points with height >0.20 m (see example in Figure 2C). Among the 426 recorded GRPs, 365 were subject to further analysis ( Table 2). In a similar way as for the trees, the vegetation height for the GRPs was designated h.

Laser Scanner Data
ALS data were acquired under leaf-on conditions using a fixed-wing aircraft. The acquisition took place on 18 June 2017 using a LMS-Q1560 laser scanner system (Riegel, Horn, Austria) and was part of the governmental effort to construct a new detailed terrain model for Norway. The study area was located within a 1169 km 2 ALS block for which ground control points were established across the entire block for calibration of the height of the laser measurements. The contracted minimum point density for the block expressed as number of first echoes per 10 m × 10 m cell tessellating the block was 5 points m −2 . The data satisfied this criterion within our study area. In fact, in certain parts of the 12 ha study area the point density was >25 points m −2 due to side overlap between adjacent, parallel strips and a single flight line perpendicular to the main direction of the scanned block. This dataset was used by the data vendor (TerraTec, Oslo, Norway) to produce the official national terrain model by classifying the points as ground and non-ground echoes using the progressive triangular irregular network (TIN) densification algorithm [30] in the TerraScan software [31].
The official national terrain model was used as terrain reference surface for the study. However, a harmonization of the point density was considered important because order statistics were used in the analysis of the tree and vegetation height. Order statistics, such as maximum height, are monotone increasing functions of number of points for a given target area [32]. In order to keep the point density stable across the study area, we discarded all data from the perpendicular flight line and from the overlap zone between adjacent, parallel strips. The resulting mean point density across the 12 ha area was reduced to 6.5 points m −2 for "first" and "single" echoes.
Normalized height values were computed for all "first" and "single" echoes relative to the official TIN by linear interpolation. Only "first" and "single" echoes with normalized height values were used in the subsequent analysis. All classified ground and non-ground points with negative normalized height values were assigned the value zero. All classified ground points were assumed to lie on the official terrain surface and where therefore assigned the value zero.

Unmanned Aerial Systems Image Data
UAS image data were acquired under leaf-on conditions using a eBee fixed-wing drone (senseFly Ltd, Cheseaux-Lausanne, Switzerland) weighing approximately 0.41 kg without payload [33]. The acquisition took place on 21 June 2017, three days after the ALS acquisition, using a Canon IXUS127 HS (Canon Inc., Tokyo, Japan) red, green, and blue camera producing three separate 16.1 megapixel images in the red (660 nm), green (520 nm), and blue (450 nm) wavelengths. The drone was equipped with an inertial measurement unit and an on-board Global Navigation Satellite System (GNSS) to control the flight parameters and provide rough positioning during flight operations [33]. The eBee flight plan was managed through senseFly's eMotion 2 software, ver. 2 [33], installed on a laptop computer. The longitudinal and lateral image overlaps were set to 90% and 80% respectively, although only a longitudinal overlap of 70% was achieved during the survey. The ground pixel resolution was set to 3.9 cm.
Prior to the image acquisition, the position of ten ground control points (GCPs) were determined and measured using the same RTK-based procedure as the one used to record positions of the GRPs. The GCP targets consisted of a set of 1 × 1 cross-shaped 4 cm × 46 cm timber planks painted orange to insure good contrast with the background vegetation.
The UAS images were processed in Agisoft PhotoScan Professional software, ver. 1.4.3 (Agisoft LLC, St.Petersburg, Russia), to produce a 3D point cloud [34]. The processing steps followed in the PhotoScan software together with the parameters used are described in Table 3. Table 3. Processing steps with corresponding parameters in Agisoft PhotoScan Professional software for the generation of 3D point cloud from UAS imagery.

Task Parameter
Alignment Accuracy: high a Generic preselection: yes a Reference preselection: yes a Key point limit: 40,000 a Tie point limit: 4000 a Adaptative camera model fitting: yes a Guided marker positioning Number of GCPs: 10 Depth maps and dense point cloud Quality: medium b Depth filtering: mild b a Parameters suggested in PhotoScan online tutorial [35]. b Parameters chosen using a trial and error approach.
After initial testing, an adaptive camera model fitting was used to perform the alignment. This function automatically selects the camera parameters to be included in the adjustment based on their reliability. The position of the GCPs were imported in the software to improve the estimates of the camera position and orientation. The GCP positions were manually refined and the camera alignment was optimized based on the GCPs to allow a more accurate model reconstruction. The average RMSEs associated with the estimated camera and GCP locations compared to the PhotoScan-estimated values were 0.92 m and 0.06 m, respectively. A dense point cloud was constructed using a medium quality parameter to reduce excessive processing time and a mild depth filtering parameter to remove outliers and reduce noise while allowing height variation between the 3D points. The point density of the resulting dense point cloud was around 50 points m −2 .
At this point we would like to clarify that the DAP methodology applied in this study is what is commonly known in the literature as structure-from-motion (SfM). Using an iterative least-squares solution, camera position and orientation, and scene geometry are simultaneously reconstructed by identification of matching features, or tie points, in multiple images. The output from SfM is fixed into a relative, not absolute, coordinate system [36]. The GCPs were used to transform the data to the absolute coordinate system adopted in this study. For the sake of simplicity, we refer to this methodology as DAP.
Normalized height values were computed for the DAP points relative to the official TIN by linear interpolation. Because the absolute height values of the DAP data were determined according to the elevation of the ten GCPs, we compared the elevation of the GCPs to the elevation of the official TIN. There was a mean difference in elevation of 0.055 m with a standard deviation of the differences of 0.028 m. The mean difference was subtracted from the normalized heights of all the DAP points. All DAP points with negative normalized height values after this subtraction were assigned the value zero.

Trees
Crown polygons for the recorded individual trees were constructed as ellipses around the recorded stem positions with the perpendicular crown diameter measurements (N-S, E-W) as minor and major axes. The tree crown polygons were laid atop the ALS dataset and the DAP dataset. Three different polygon height metrics were calculated for each of the two remotely sensed datasets. They were the maximum, mean, and the 90th percentile. The analysis conducted in this study revealed, however, that maximum height produced consistently greater accuracies in for example the tree height modeling and prediction. This is consistent with previous findings showing maximum height to be a strong predictor for tree height of small trees [10]. The fact that only a single point would be present for a large fraction of the small trees, at least in the ALS dataset (see examples in Ref. [9]), would exclude the use of, for example, deciles or moments of the height distributions, which are commonly used for modeling biophysical properties of larger forest trees. Only the results for the maximum values within each crown polygon were documented, and they were designated h ALSmax and h DAPmax .

Ground Reference Points
Circular polygons were constructed for each of the GRPs for which the vegetation height recorded in field was < 0.20 m ( Table 2). These polygons were laid atop the ALS and DAP datasets. When the GRP dataset was established in 2010, the variable "terrain form" was recorded within a circle with radius 1 m centered on the GRP [29]. However, the variable "vegetation height" was recorded for the point only without any further assessment of the vegetation height surrounding the point. In the current study, we chose a radius of 0.5 m for the circular polygon to which the recorded vegetating height was assigned. Even though we restricted the size of the polygon to a radius of 0.5 m, there was a risk of overhanging and non-recorded trees and bushes aside the GRP but with presence inside the polygon, which potentially could be represented by large positive height values in the remotely sensed point data and for which we had no field observations. We therefore inspected the ALS point data and the DAP data numerically and visually to identify such polygons where there would be a likely mismatch between the field-recorded vegetation height and the height in the remotely sensed data. There were 10 among the 365 polygons for which the ALS data had h ALSmax > 1.00 m when h ALSmax was defined in a similar way as for the tree polygons (Section 2.5.1). Because the point density of the ALS data was much smaller than for the DAP dataset, the DAP dataset contained a greater number of polygons with maximum heights > 0.20 m than the ALS dataset. Among the 365 polygons, 24 and 26 polygons in the ALS and DAP data, respectively, had maximum heights > 0.20 m. We decided to discard these polygons from both remotely sensed datasets with maximum height > 0.20 cm in either of the two datasets. Thus, 327 polygons were retained for the analysis addressed in objective #2. Because we did not have ground observations of the vegetation height within the polygons, there is uncertainty associated with the final data resulting from this data screening. For the retained polygons, h ALSmax ranged between 0 and 0.19 m with a mean value of 0.07 m. For DAP, the maximum value (h DAPmax ) was in the range 0-0.19 m with a mean value of 0.03 m.

Population Elements
Under objective #2, we compared the precision of vegetation height estimates across the chosen study area using the different remotely sensed 3D data. The 200 m × 600 m study area was tessellated into regular population elements of 1.5 m 2 in size. This size was a compromise between the size of the GRP polygons (0.79 m 2 ) and the tree polygons (2.10 m 2 ) subject to modeling under objective #2 (Section 2.7.1). The resulting 79,242 population elements constituted the overall population in a statistical sense (Section 2.7.2). The number of elements was slightly smaller than the theoretical size of 80,000 elements due to a small water body in the study area that was excluded from the population. The maximum point height was extracted for each individual element for each of the three 3D remotely sensed datasets.

Analysis-Objective #1
Under objective #1, we assessed and compared the performance of 3D remotely sensed data from ALS and from DAP for prediction of tree height of small pioneer trees and evaluated how tree size and tree species affected the predictive ability of the two types of 3D data. The main steps of the analysis are shown in Figure 3 for the sake of clarity and overview. A first step of the assessment was to analyse to what extent the different 3D data were sensitive to the small trees, i.e., if positive height values of the point clouds could be expected for a tree. Previous research (e.g., [9]) suggests that this will depend on factors such as tree height, size of the tree crown, tree species, the point density of the remotely sensed data which in the current study clearly differed between ALS and the DAP 3D data (Sections 2.3 and 2.4), and degree of laser pulse penetration into the tree crowns for ALS as opposed to a likely depiction of the outer surface of a tree crown with DAP data.
A logistic regression analysis with binary response supported this assessment. Among the 532 field-measured trees (Table 1), two trees had a substantially higher maximum height in the 3D remotely sensed datasets (1.51-2.86 m) than field-measured tree height. These two trees (trees #48 and #2100) had likely overhanging branches from taller, neighbouring trees and they were discarded from all subsequent analysis. They were both spruce trees. Fifteen trees with maximum heights in the remotely sensed datasets 0.20-0.97 cm greater than the corresponding field-measured tree heights were retained because we were unable to identify a specific reason for this pattern. We were thus rather conservative in the treatment of potential outliers. The logistic regression analysis was based on the remaining 530 trees (Table 4).
For each of the two remotely sensed datasets (ALS, DAP) every tree was classified as POSITIVE if the maximum height for the tree polygon (h ALSmax or h DAPmax ) had a positive value. If the maximum value was zero or the tree polygon did not contain any points for a given 3D remotely sensed dataset, the tree was classified as ZERO. The analysis was carried out in two steps. First, a general logistic regression model reflecting all effects mentioned above was fitted. This model of the probability of POSITIVE was formulated as follows: where π POSITIVE is the probability of maximum height of a tree polygon with a value greater than zero using observations from datasets (DATA) ALS and DAP. DATA DAP is a dummy variable for DAP (DATA DAP = 1 if DAP). Further, SP pine is a dummy variable for pine (SP pine = 1 if pine), SP birch is a dummy variable for birch (SP birch = 1 if birch), h (m) is the tree height measured in field, and A (m 2 ) is the elliptic tree crown area according to the field recordings of crown diameters. The betas (β 0 , β 1 , β 2 , β 3 , β 4 , β 5 ) are parameters to be estimated. Maximum-likelihood computation for fitting of the logistic model in Equation (1) was performed with the LOGISTIC procedure of the SAS package [37]. It should be noted that the reference in the model is the ALS dataset and the tree species spruce. Thus, the estimated parameters for the DATA and SP variables express differences relative to this reference (differences in intercept of the model). The effects of, for example, DAP relative to ALS will be expressed directly by the parameter estimate of the former variable. Finally, a Wald chi-square test was performed to test the null hypothesis that the parameter estimates for the two dummy variables for tree species were equal.
One of the results of the first step of the logistic regression analysis was that the effects of tree species on probability of detected trees differed significantly in the statistical sense between some of the species (p < 0.001, p = 0.037, and p = 0.080, respectively), see Table 7. On the other hand, the effect of dataset was not significant in the statistical sense (p = 0.733; Table 7). Further, both tree height and crown area were statistically significant (p < 0.001, Table 7).
Although some effects in the basic model in Equation (1) were significant and others not, some of the effects are likely confounded which may lead to incorrect interpretations. For example, the point density of the DAP point cloud was around 50 points m −2 whereas the corresponding density in the ALS data was 6.5 points m −2 . It is therefore reasonable that the area of a tree crown polygon is more critical for a crown polygon having a positive height value in the ALS data than in the DAP data. Likewise, tree species may affect the probability of positive height values differently in the two 3D remotely sensed datasets since laser pulses tend to penetrate the tree crowns before an echo is triggered while DAP may better capture the surface of a crown. Crowns of different species have different densities of biological matter (foliage and branches) and different shapes which may influence the point clouds for the two 3D remote sensing techniques differently. A more complex model was therefore formulated. In the model in Equation (1), it was assumed that the effect of dataset was similar for each individual tree species, i.e., that the different datasets only affected the intercept of the model. In addition to the basic effects accommodated by Equation (1), we allowed the effects of tree species to vary between the two 3D remotely sensed datasets. This was accommodated by introducing separate regression coefficients for tree species for the different datasets. Further, in the former model, it was assumed that the effect of dataset was constant across the entire range of tree heights and tree crown areas. In the second step of the analysis, we allowed the effects of dataset to vary according to the magnitude of the tree height and the tree crown area as well. This was accommodated by introducing separate regression coefficients for tree height and crown area for each individual dataset in the model: Similar to the model in Equation (1), the ALS dataset and the species spruce represent the reference in the model in Equation (2). Thus, the estimated parameter for the DATA variable (β 1 ) will express the overall difference in intercept relative to ALS, while parameters for the two SP variables (β 2 and β 3 ) will express the overall difference in intercept relative to spruce. The height and crown area parameter estimates (β 6 and β 8 ) will express the general effects of these two variables. The estimated parameters for the respective products of the DATA variable and the two species variables (β 4 and β 5 ), and the DATA variable and height and crown area (β 7 and β 9 ), will express differences in parameters for pine, birch, h, and A for DAP relative to ALS. Finally, a Wald chi-square test was performed to test the null hypothesis that the parameter estimates for pine and birch (β 2 and β 3 ) were equal. Likewise, Wald chi-square tests were performed to test the null hypotheses that the parameter estimates of the products of the DATA variable and the two SP variables (β 4 and β 5 ) were equal.
The second part of the analysis under objective #1 entailed modeling of tree height of the pioneer trees and evaluation of how tree size and tree species affect the predictive ability of the two types of the 3D remotely sensed data. Following a similar strategy as in the logistic regression analysis, we formulated a model with tree height observed in field as dependent variable and maximum height for each tree crown polygon in each of the 3D remotely sensed datasets and the factors to be evaluated as independent variables: where h max = h ALSmax when the dataset was ALS and h max = h DAPmax when the dataset was DAP. The other variables in the model were defined as above. The analysis was based on 389 of the trees for which h max ≥ 0 in both datasets ( Table 5). The least squares method for fitting the model was applied by using the REG procedure of the SAS statistical software package [37]. F-tests were performed to test the null hypotheses that (1) the parameter estimates for the two SP variables (β 2 and β 3 ) were equal and that (2) the parameter estimates of the products of the two SP variables and h max (β 6 and β 7 ) were equal. Finally, leave-one-out cross validation was adopted to assess the predictive ability of the two 3D datasets. However, the model in Equation (3) assumed the error variance to be the same for both 3D datasets. Separate models would be required if the error variances could be assumed to be different ( [38], p. 173). The cross validation was therefore performed for separate models constructed according to: where h max = h ALSmax when the model was constructed by using the ALS data. h max = h DAPmax when the model was constructed by using the DAP data.
In the cross validation of the two respective models, the prediction accuracy was assessed separately for different classes according to tree height and different tree species. The assessment was based on the differences between predicted and observed tree height for individual trees according to the statistics (1) mean difference, (2) standard deviation of the differences (Stdev), and root mean square error (RMSE). These statistics were also calculated across all trees for each individual 3D dataset and the null hypothesis of homogeneity of prediction variances among the two 3D datasets was tested by Levene's F-test [39] in the GLM procedure of the SAS package [37].

Analysis-Objective #2
To provide initial overview, the main steps of the analysis under objective #2 are shown in Figure 4.

Models for Vegetation Height
The first step of this analysis entailed construction of regression models used for prediction of vegetation height. These predictions were subsequently used to estimate mean vegetation height following model-dependent inferential principles (see details in Section 2.7.2). A critical assumption for model-dependent estimators to be approximately unbiased is that the model is correctly specified for the area of application. Misspecification of the model can lead to serious bias in the estimators [40].
Considerable research has been devoted to development of combinations of sampling designs and estimators that protect the inference from the adverse biasing effects of model misspecification [41]. A primary finding has been that model-dependent estimators tend to be biased unless the sample is balanced, i.e., the sample moments of the distribution of the independent variables equal the corresponding population moments [42]. In the current study, the dataset available for model construction was composed of the 389 tree polygons (Table 5) and the 327 polygons constructed for the 327 GRPs (see Section 2.5.2), i.e., 716 observations in total. The 327 GRP polygons were added to the tree dataset to better represent the vegetation structure in the study area. In the remote sensing community, criteria for characterizing the appropriateness of data and models for modeldependent inference have received little attention. A couple of exceptions are the studies by Refs. [12,43] who made explicit reference to the effects of sample imbalance on bias.
h max was the independent variable to be used in the models that were to be constructed. Prior to choosing a specific model form and model fitting technique, we constructed the distributions of h max for the sample of 716 observations and the 79,242 population elements that constituted the population and calculated the four first moments of the distributions where h max = h ALSmax when we used the ALS data and h max = h DAPmax when the DAP data were used. As is evident from the graphical presentation of the distributions (Figure 5), the population distributions were extremely skewed towards small values of h max whereas the samples had few observations in the lower end of the distributions. This is also evident in the calculated moments (Table 6). With a very large majority of population elements in the lower end of the distributions it is obvious that we should strive for models with appropriate prediction properties in that part of the population.  Non-linear models of the form y = β 0 x β 1 , where y is the dependent variable and x the independent variable, are commonly used to construct prediction models for biophysical parameters such as tree height and forest biomass with 3D data from ALS and DAP. One reason for choosing this model form is the fact that the predicted value of y will be zero when x is zero, which is a logical property in many applications. In our case, most of the population elements have a value of the independent variable equal to or close to zero. However, it is well documented that especially lasers tend to penetrate into tree and vegetation canopies before an echo is triggered. This is well illustrated for height of small trees in, for example, the study by Ref. [9], Figure 4. Therefore, forcing the model through origo will most likely result in a general tendency of under-prediction at the lower end of the range of the dependent variable and thus a biased estimator of vegetation height. The same would be the case with a simple linear zero-intercept model like y = β 1 x. Based on the 716 observations, we therefore chose to construct simple linear models of the form: where h max = h ALSmax when the model was constructed by using the ALS data. h max = h DAPmax when the model was constructed by using the DAP data. Now, given the large differences in population and sample distributions of h max , measures should be taken to reduce the risk of inappropriate models for the population in question and thus biased estimators. We chose to adopt a weighting scheme in the model fitting by which weights were assigned to each of the 716 sample observations in such a way that the weighted sample distributions of h max approximated the population distributions of h max . First, for each of the two remotely sensed datasets, we calculated the h max values corresponding to the nine percentiles p10, p20, . . . , p90 of the population distributions and formed ten equally large, ordered classes <p10, p10-p20, . . . , p80-p90, >p90. Then we assigned the 716 sample observations of each 3D dataset to these ten classes according to the h max value of each individual sample observation. Since each class in the population constitutes exactly 10% of the population, each class was given a total class weight of 0.1 (1/10). Further, each sample observation was given a class-specific weight depending on how many sample observations that were assigned to a particular class. For example, if a class contained 100 sample observations, each sample observation would have a within-class weight of 0.01 (1/100). Finally, we calculated individual weights for each sample observation by multiplying the class weight (0.1) by the within-class weight. This weighing scheme ensured that the sum of the weights was always equal to 1. It should be noted that for DAP, 25% of the population distribution had an h max value of zero. Thus, all sample observations falling in the three lowest classes (up to p30) where assigned the same weight. The weighted sample distributions of h max are illustrated graphically in Figure 5 and their moments are presented in Table 6. As is evident from the graphical presentation as well as the calculated moments, the weighted sample distributions were generally much more similar to the population distributions than the sample distributions in their original form.
The models were constructed according to Equation (5) and the weighting scheme outlined above using the least squares method as implemented in the lm function of the stats package [44]. White's and the Studentized Breusch-Pagan test statistics were calculated using the white_lm and breusch_pagan functions, respectively, of the skedastic package [45]. Both tests rejected the hypothesis of homoscedastic residuals for ALS as well as for DAP (p < 0.001). In the presence of heteroscedasticity, heteroscedasticity-consistent covariance matrix estimators were used, as recommended by Ref. [46]. The heteroscedasticityconsistent covariance matrix estimators of type HC 3 , presented by Ref. [47], were computed using the sandwich package [48,49] in R.
The constructed models were subsequently applied for prediction for every single 1.5-m 2 population element across the entire population. The predictions were performed for each of the two individual 3D remotely sensed datasets with h max of each element as predictor variable. The result was four prediction maps for vegetation height, two for each of the 3D remotely sensed datasets using the weighted and unweighted models, respectively. Two of these maps are presented in Figure 6.

Estimator for Vegetation Height
Based on the predictions, estimates of mean vegetation height were then produced for each of the 3D remotely sensed datasets and by using the weighted and unweighted models, respectively, and for different domains of the population (see details in Section 2.7.4). The general approach to estimation adopted in this study is known in the literature as the area-based approach. In the current study, we did not have a probabilistic sample that would have allowed design-based inference. Model-dependent estimation and inference was therefore adopted. Model-dependent inference has been applied frequently in recent years when estimating biomass, volume, and other biophysical parameters using remotely sensed data (e.g., [50][51][52][53]). An overview of the concept and a brief review of recent studies can be found in Ref. [54]. In the current study, we adopted the estimators in the way they were formulated by Ref. [12]. The context in the study by Ref. [12] was slightly more complex than in the current study. First, they estimated changes in height in bi-temporal data, not the height using single-date data. Second, they estimated height changes of trees in addition to height changes of all vegetation, which required two separate models-one for classification of trees and other (non-tree) vegetation and one for prediction of changes in height. This complicated in particular the estimation of the uncertainty.
This estimator can be used for smaller domains within the population as well.

Estimator of Mean Square Error
The term mean square error (MSE) rather than variance was used to characterize the uncertainty of the model-dependent estimator in Equation (6) because the modeldependent estimator of the population mean cannot be assured to be unbiased. For large areas, model-dependent MSE estimators will, in general, depend mainly on the uncertainty of the estimates of the model parameters (e.g., [55][56][57]). For small domains, an additional source of uncertainty must often be accounted for, namely the residual variance. Ref. [58] derived a model-dependent variance estimator that accounted for residual variance and incorporated a spatial autocorrelation structure of the residuals. Ref. [52] demonstrated with empirical data of timber volume from small forest stands and auxiliary data from DAP that ignoring the residual variance component may induce bias in the mean square error estimator. In the study on vegetation height change conducted in the current study area [12], the analysis suggested that " . . . most of the mean square error estimates (>95%) of the estimators will be accounted for by quantifying the variance attributable to the model parameter uncertainty". Nevertheless, it cannot be assumed that the magnitude and spatial structure of the residual variance of vegetation height predictions necessarily follow the same patterns as the height change predictions. In the current study, we therefore addressed all the mentioned components, i.e., (1) the variance due to uncertainty of the model parameters and the residual (2) variance and (3) spatial covariance components. These three components were treated as additive to obtain the total mean square error and will be described in detail below.
Under the model-dependent inferential framework, the variance due to uncertainty of the model parameters can be approximated either by using a closed-form formula based on a first-order Taylor series approximation (see e.g., [57]), or by using Monte Carlo simulations in the form of, for example, parametric bootstrap. Both approaches to estimation have been adopted in forest-and vegetation-related studies in recent years (e.g., [50,59]). Ref. [60] noted a few properties of parametric bootstrap that makes it attractive and demonstrated the technique with ALS data. In the current study, we chose to adopt this technique, as in Ref. [60]. In the following, it is assumed that the model parameter estimates of the predictive model (β in the model in Equation (5)) follow asymptotically multivariate normal distributions, i.e.:β where the expected value of the vector ofβ estimates is E[β] and Σβ is the heteroscedasticityconsistent estimates of the variance-covariance matrix ofβ. By sampling from the multivariate distribution in Equation (7), a large parametric bootstrap sample of random vectors β PB ∼ N β , Σβ was generated. The sample was denoted S PB , where S PB = {1, . . . , l, . . . , M}. This sample can be used to produce new predictions of h, according to Equation (5). Predictions of h were produced for all M random vectors β PB and for all N population elements. Thus, we obtained unique predictionsĥ PB,k,l for k ∈ U and l ∈ S PB . A parametric bootstrap variance estimator for the point estimator in Equation (6) is: where:ĥ and:ĥ Analytical estimators for residual variance under heteroscedasticity have been adopted in previous analysis of important parameters encountered in forest surveys, for example timber volume [52]. In the current study, every geographical domain subject to estimation had sample units (ground observations of vegetation height) that could be used to provide estimates of residual variance. Thus, for a particular domain with a sample S with n sample units, S = {1, . . . , k, . . . , n}, the residual variance for the point estimator for mean vegetation height (Equation (6)) was formulated as [12]: Residual covariance of substantial magnitude, as compared to the other sources of uncertainty when estimating forest resource parameters for small areas, has been encountered for shorter distances in several studies (e.g., [52]), while at greater distances-and consequently for larger areas-the residual covariance is often assumed or found to be negligible in magnitude (e.g., [51]). However, as noted above, Ref. [12] found the residual covariance to be negligible even for areas as small as 1.5 ha. Analytical ways of addressing the residual covariance have been demonstrated by e.g., Ref. [52].
Quantifying the spatial autocorrelation of the residuals is essential in the analysis of the residual covariance. Spatial correlation (ρ) is often estimated from the model prediction residuals by constructing a correlogram. Assuming that a correlogram has been fitted and that ρ can then be predicted to obtain predicted values of the correlation for all combinations of the N population elements in U, the residual covariance for the point estimator of mean vegetation height for a particular domain for which we had actual observations of the residuals, can be estimated by [12]: whereρ kl is the predicted value of the residual correlation between elements k and l in U and the correlogram was fitted using the residuals h k −ĥ k and h l −ĥ l .

Calculations
The first step of the analysis was to construct the regression models for the ALS and DAP data according to Equation (5) based on the 716 weighted sample observations. As noted above, for the sake of comparison of estimates we also constructed simple models according to Equation (5) without adopting the weighting scheme. These models are nevertheless not documented in the current article.
Once the regression models were constructed, we proceeded with the estimation. First, we estimated mean vegetation height across the entire study area, according to Equation (6), by using model predictions from the height models (Equation (5)).
Two cases of special interest were identified. First (Case A), we tessellated the study area into 1 ha cells that may serve as the primary mapping and monitoring units in, for example, a tree line monitoring program (see Figure 6B). This size was chosen for demonstration purposes, and this size can be changed to any size found meaningful for a particular application. 1 ha resolution may be found useful for some applications, while, for example, 100 m 2 resolution may be found relevant for other applications. We could even consider the primary resolution of 1.5 m 2 , but some level of aggregation is in many cases perhaps easier to comprehend and interpret. Another reason for choosing an area as large as 1 ha was that a restricted number of cells (12) would ease the interpretation of any potential differences in the properties of the estimates obtained from the two different 3D remotely sensed datasets.
Second (Case B), the current study area is subject to variations in wind exposure, temperature, soil depth, and moisture due to its location on and along a small ridge with variations in aspect. In particular, bedrock outcrops and mires are abundant in certain parts of the area. These factors tend to influence establishment of different types of vegetation, its height and density. We wanted to assess any potential differences in estimated height and precision for the two 3D remotely sensed datasets which could be attributed to different vegetation structures. The area was therefore delineated into four distinct sub-regions (see Figure 6A,B) based on visual interpretation of the orthophoto (natural colors). Thus, this delineation was independent of the height information in the ALS as well as the DAP point clouds. In the visual interpretation, we sought homogeneity within each polygon with respect to appearance of trees, greenness of the ground vegetation suggesting existence of bushes, and lighter areas suggesting lichens and bedrock outcrops.
Finally, we estimated the various MSE components for the entire study area and for each geographical domain (twelve 1 ha cells and four sub-regions) subject to analysis. First, we constructed separate correlograms for the residuals for each domain of interest. Correlations based on the observed residuals were graphed and then visually inspected. For natural phenomena, it is common to observe greater correlations at shorter distances and then declining correlations with increasing distances. Such spatial structures can be modeled, e.g., by some exponential model form. In the current dataset, no such spatial structure could be observed (see examples in Figure 7). Similar to the study on vegetation height change [12] linear regression models of the form: appeared to be suitable for vegetation height as well, where D is the distance between pairs of residuals. In total, 34 correlogram models were constructed according to Equation (13), i.e., separate models for each domain subject to estimation for each of the two 3D remotely sensed datasets. The results of the model construction revealed non-significant estimates (p > 0.05) of β 0 as well as for β 1 for 16 of the 17 models constructed for ALS and for 15 of the 17 models for DAP for β 0 and 14 of the 17 models for β 1 . Thus, the results suggested constant residual correlation equal to zero in 30 of the 34 cases. We therefore concluded that the residual covariance components of the estimates of MSE ( MSE) would be zero and could be ignored in these 30 cases. Even in the four remaining cases, the autocorrelation was negligible, even though the parameter estimates were statistically significantly different from zero. One of the three cases for DAP was the model constructed across the entire study area, which also had the smallest p-values (Figure 7, bottom). As is evident from visual inspection of the graphical presentation in the figure, there is hardly any structure in the spatial correlations of the residuals. For all the 34 cases, the residual covariance was therefore ignored.
In the characterization of the uncertainty of the various estimates of vegetation height, we proceeded with estimating MSE by ignoring the residual covariance and adding the two remaining variance components, accounting for model parameter uncertainty (Equation (8)) and residual variance (Equation (11)). Subsequently, we calculated the standard error of the mean estimate (SE; square root of MSE). The proportions of MSE that were attributed to the residual variance were also calculated.
In the bootstrap variance estimation used to characterize the model parameter uncertainty (Equation (8) Table 4, among the 530 trees subject to analysis, most of the trees greater than 1 m in height had positive maximum height values in the ALS as well as the DAP data. It is also evident that a greater portion of the trees lacked height observations in the ALS data than in the DAP data, which can be attributed to the greater point density of the DAP data than the ALS data. A surprising feature was noticed for some of the trees in the DAP data. For example, for a birch tree (tree #375) with h = 2.52 m and a crown area of A = 0.66 m 2 , the maximum height in the DAP data was h DAPmax = 0 whereas the corresponding height in the ALS data was h ALSmax = 2.23 m. Tree #375 appeared close by a pine tree (tree #1102) (gap between crowns of 0.2 m) with h = 3.00 m and crown width of 1.83 m 2 . Tree #1102 had h ALSmax = 2.33 m and h DAPmax = 0.38 m. Obviously, the DAP algorithm was unable to follow the outline of the crowns of these trees. Both trees are clearly visible to the human eye in the orthophoto and appears as a solitary cluster of two trees. This illustrates that the DAP algorithm sometimes produced spurious results.
The results of the regression analysis of probability of trees having a maximum height in the ALS and DAP datasets greater than zero (model in Equation (1)) suggested statistically significant differences among two of the tree species (p < 0.037; Table 7) while the effect of dataset (ALS versus DAP) was not significant (p = 0.733; Table 7), see Section 2.6. Both tree height and crown area were statistically significant (p < 0.001, Table 7). However, a more detailed regression analysis following the model in Equation (2), suggested that the statistical significance of many of the effects in the initial regression analysis following Equation (1) were misleading because effects were confounded. The detailed analysis showed that there is generally a higher probability of positive height values for the trees in the DAP data than in the ALS data (DATA DAP : p = 0.009; Table 7). There was hardly any effect of tree species at all, with p-values for the various parameter estimates for the tree species dummy variables and the Wald chi-square test ranging from p = 0.194 to p = 0.637 (Table 7). There was, however, a slightly negative and significant effect for the interaction of pine and DAP data (SP pine ·DATA DAP : p = 0.045; Table 7), indicating a smaller probability of pine trees having positive height values in the DAP data than in the ALS data. The pine trees in the study area have to a large extent more compact tree crowns than the two other species, which may lead to an early triggering of ALS echoes as the laser pulses penetrate into the tree crowns. The greater probability of positive height values for ALS for pine trees might therefore be attributed to better conditions for laser echo triggering in compact tree crowns.
The most striking result of the detailed analysis in Equation (2) was found for tree height and crown area. There was a very strong effect of crown area for ALS (A: p < 0.001; Table 7) with substantially greater probability of positive height values for the trees with increasing crown size, whereas this effect was strongly dampened for DAP expressed by the interaction of crown area and DAP data (A · DATA DAP : p < 0.001; Table 7). This is a very logical result and confirmed our expectations (see Section 2.6). Crown area is critical for even hitting the smaller trees with ALS pulses, given the point density in this study of 6.5 points m −2 , whereas this certainly is less critical for DAP, for which the point density was 50 points m −2 . Further, the height as such was not statistically significant for ALS (p = 0.540; Table 7), at least as long as crown area is included in the model. It should be noted that h and A were inter-correlated with Pearson r = 0.74. For DAP there was a tendency of greater probability of positive height values for the trees with increasing height than for ALS, but this effect was not significant in the statistical sense (h · DATA DAP : p = 0.052; Table 7).

Influence of Tree Size and Tree Species on Tree Height Modeling
The second part of the analysis under objective #1 entailed modeling of tree height of the pioneer trees and evaluation of how tree size and tree species affected the predictive ability of the two types of the 3D remotely sensed data using the 389 trees documented in Table 5. The regression analysis following the model in Equation (3) confirmed a statistically significant overall greater height of trees measured in field than by the two 3D remote sensing techniques as expressed by the positive intercept of the model (p < 0.001; Table 8). This tendency was smaller for DAP than for ALS as expressed by the negative sign of the dummy variable DATA DAP (p < 0.001; Table 8) but this difference between ALS and DAP was compensated for as the height of the trees increased. There was a strong effect of increasing differences between h and h max with increasing tree heights in general (p < 0.001; Table 8) with an additional effect with increasing tree height for DAP relative to ALS (h max · DATA DAP : p < 0.001; Table 8). The stronger underestimation of tree height for the DAP data than for the ALS data with increasing tree height can easily be seen even by visual inspection of the data ( Table 5). This is a somewhat surprising finding and contrary to the expectation that DAP may better capture the properties on the outer surface of a tree canopy of small trees than lasers, and that the DAP height points thus could resemble the surface of the crowns. This result is nevertheless consistent with recent findings reported by Ref. [26]. Whether this is a general phenomenon or other DAP algorithms/software packages and parameter settings would produce substantially different results is unknown and out of the scope of the current study. However, this might be a subject for further studies. The analysis further revealed that the effect of pine as opposed to spruce and birch was negative on modelled tree height, meaning that differences between h and h max for the two 3D remotely sensed datasets were greater for pine than for the two other species (p < 0.001; Table 8). Especially the difference between pine and birch was somewhat surprising (F-test: p < 0.001; Table 8) since the pine trees in most cases had compact and regularly shaped crowns which intuitively should be easier to handle using DAP and even result in early triggering of echoes in laser scanning, whereas the birch trees often constituted shrubby and more open vegetation forms (see illustrations in Figure 2). The separate analysis of ALS and DAP following the model in Equation (4) even showed that the specific effect associated with pine trees was more pronounced for DAP (p < 0.001; Table 8) than for ALS (p < 0.05; Table 8). This is counterintuitive, given the expectation that a more regular surface as represented by a pine tree crown would be better suited as an object for DAP.
Thus, to get a better understanding of the performance of the DAP, we inspected the data for individual pine trees in more detail by looking at their spatial context in the 3D DAP data as well as in the orthophoto. Among the 64 pine trees with h max = 0 m in the DAP data, 16 trees had h max > 0 m in the ALS data despite the much smaller point density in the ALS point cloud. These 16 trees had a mean tree height of 0.68 m (range: 0.16-1.09 m) and mean crown area of 0.36 m 2 (range: 0.02-0.93 m 2 ). Most of these trees appeared as solitary objects or as objects surround by only low vegetation. Thirteen of them were found in the open areas dominated by mosses and lichens in sub-region II, see orthophoto in Figure 6A. A reason for the DAP height values of zero for these trees might be an inability of the DAP algorithm to capture height variations when the objects have a small horizontal extension, i.e., a smoothing effect. As suggested below, the DAP algorithm may behave differently for larger groups of trees even if the height is not necessarily greater.
When the ability of ALS and DAP data to predict tree height by regression models was assessed by leave-one-out cross validation, separate models were fitted for ALS and DAP data respectively, following the model in Equation (2). The choice of fitting separate models was well justified by rejection of Levene's F-test of the null hypothesis of homogeneity of prediction variances between the two 3D datasets (p = 0.0032). The performance of the two models in terms of R 2 and RMSE were 0.90 and 0.83, and 0.36 m and 0.46 m for ALS and DAP, respectively ( Table 8). The RMSE of the model fitting corresponded closely with that of the cross validation of the same models (0.37 m and 0.47 m, Table 9). This is well in line with previous findings for small pioneer trees using ALS. For example, Ref. [18] reported an R 2 value of 0.86 and an RMSE value of 0.49 m in the same study area as the current study, but with data from another point in time and a laser scanning with greater point density. However, they used a simpler model with just a single independent variable (h max ).
For DAP, the current study is to our knowledge the first effort to model and predict height of individual pioneer trees in the boreal-alpine ecotone. Ref. [25] nevertheless modelled height of young forest plots with trees in a similar or somewhat greater height range (range: 0.5-13.0 m; mean: 2.5 m) than in the current study. They addressed mean height of plots and stands rather than individual trees. They reported similar or better prediction performance in terms of RMSE for DAP data than for ALS. However, their plots had a mean stem number of 5572 ha −1 as opposed to 97.2 ha −1 for the current study area [9]. A reason for a better performance of DAP for dense plots as opposed to more solitary individuals or clusters of trees might be an ability of the DAP algorithm to follow the surface of larger objects in the form of a forest plot with denser and more continuous canopy surface as opposed to solitary objects, which in our data tended to be smoothed out. Table 9. Results of leave-one-out cross validation of the models in Equation (4) ( Table 7) based on differences between predicted and observed tree height for the different 3D remotely sensed datasets, tree species, and tree height classes (n = 389).  The cross validation also revealed a tendency of greater standard deviation with increasing tree height, at least for DAP (Table 9). This heteroscedasticity should be mitigated by more sophisticated approaches to modeling when the models are to be used for estimation (see Section 2.7.1). For both ALS and DAP there was also a tendency of statistically significant over-prediction for small values of tree height and under-prediction of large values of tree height for spruce in particular (p < 0.001; Table 9).

Model Construction
Results for the models constructed for vegetation height based on the 716 sample observations are displayed in Table 10, which also includes the estimated heteroscedasticityconsistent variance-covariance matrices for the parameter estimates that were needed for the parametric bootstrap variance estimation. All parameter estimates for ALS as well as for DAP models using weighted as well as unweighted sample observations were highly significant (p < 0.001) and the models displayed an adequate fit to the data with R 2 values ranging from 0.80 to 0.91 and RMSE values ranging from 0.36 to 0.49 m. We are unaware of any previous studies that have constructed similar prediction models for vegetation height in the boreal-alpine ecotone. Model fit assessed by R 2 values and RMSE nevertheless corresponded well with models constructed for trees only, see e.g., Table 8 and Ref. [18].

Overall Vegetation Height Estimation and Inference
Based on the predicted vegetation height using the ALS and DAP data (see Figure 6C,D), the estimated mean vegetation height across the entire study area was 0.64 m (SE = 0.014 m) for ALS and 0.76 m (SE = 0.018 m) for DAP, respectively, when the prediction models were constructed with weighted sample observations (Table 11). An F-test revealed that the SE for DAP was significantly greater than for ALS (p < 0.001). Further, 95% confidence intervals of the mean vegetation height estimates clearly showed that the height estimate for DAP was significantly greater than for ALS, the difference in estimates being as great as 0.12 cm. This difference was surprisingly large, given that the same ground observations were used to construct the respective prediction models. This issue is further discussed in light of the different properties of the various sub-regions under Case B, see Section 3.2.3 below. Table 11. Estimates of mean vegetation height (ĥ; Equation (6)) for various domains of the study area and corresponding standard error estimates a (SE), and contribution of the residual variance component (Equation (11)) to overall mean square error estimates of mean vegetation height (var(ĥ) res ). a Standard error based on the sum of the variance due to parameter uncertainty (Equation (8)) and residual variance (Equation (11)). b Point estimates based on models constructed from unweighted sample observations.

ALS
A systematic difference on the order of 0.12 cm can in fact have a rather detrimental effect on change estimates. Ref. [12] estimated the overall change in vegetation height for our study area to be 0.16 m (SE = 0.02 m) for the time period 2006-2012 using ALS data for both points in time. Clearly, small but statistically significant change estimates can easily be confounded with systematic effects associated with choice of remote sensing technology used in the estimation.
It was also revealed that the mean vegetation height estimates were 0.03 to 0.04 m greater for ALS and DAP, respectively, when unweighted sample observations were used. Despite the fact that the differences in estimates using weighted and unweighted sample observations were not significant in the statistical sense, these differences were fairly constant across the entire study area-even for the smaller 1 ha cells of Case A representing a great variation in vegetation properties (Table 11).
The fact that there were significant differences between ALS and DAP estimates of mean vegetation height and between estimates based on weighted and unweighted sample observations, demonstrates that the model-dependent estimators are likely biased. The magnitude of the bias is however unknown, and based on the data and the analysis we can hardly tell which of the estimators that are the least biased. Simulations based on a known population would be needed to assess estimator bias. It is nevertheless reasonable to assume that weighting the sample observations during model construction to obtain more similar distributions of the predictor variables in the population and in the sample will provide models that are more apt for prediction purposes for the study area as a whole. Thus, the results illustrate that the composition of the sample is an important aspect of model-dependent inference, which often seems to be neglected by the remote sensing community.
Further, these results demonstrate that change estimation based on bi-temporal data with a combination of ALS and DAP may be challenging even when field data are available for model calibration at each point in time.
At the spatial scale of the entire study area (200 m × 600 m) the residual variance had a negligible effect on overall MSE estimates of mean vegetation height. The residual variance contributed only 0.85% and 0.96% to the overall MSE for ALS and DAP, respectively (Table 11). This is on the same order of magnitude as for mean change in vegetation height (0.4%) using ALS data in the very same study area [12]. Thus, residual variance can clearly be ignored at this spatial scale without any adverse effect on the statistical inference.

Domain-Specific Vegetation Height Estimation and Inference
The mean vegetation height estimates for the 1 ha cells (Case A) ranged from 0.30 m to 1.36 m for ALS (Table 11, Figure 6E) and from 0.36 m to 1.61 m for DAP (Table 11, Figure 6F) when using the weighted sample observations. The mean height estimates were generally higher for DAP for all cells as noted for the overall study area (Section 3.2.2). Despite varying differences in height estimates between ALS and DAP among cells, the two remotely sensed datasets captured the same trends in differences in vegetation height among the cells, as is evident in Figure 6E,F. This illustrates that for monitoring units of, for example, 1 ha in size, it should be possible to identify gradients in vegetation height with both remote sensing techniques, although the absolute value of the height estimates may suffer from systematic errors due to potential estimator bias, as discussed above. By assessing 95% confidence intervals of the respective cell-wise mean vegetation height estimates, it appeared that most cell estimates were significantly different from the others in pair-wise comparisons. As pointed out by Ref. [12], it should be noted, that when multiple statistical comparisons (tests) are performed simultaneously, one may wish to alter the level of significance to control the total Type I error (Bonferroni approach [62]). Especially for studies of trends in vegetation height over larger areas covering many monitoring units, this may be an important consideration.
The standard error estimates were generally smaller for ALS than for DAP, and the cell-wise estimates of SE also varied less among the cells for ALS than for DAP. The cellwise SE estimates ranged from 0.014 m to 0.018 m for ALS and from 0.019 m to 0.031 m for DAP (Table 11). At the 1 ha cell level the residual variance component accounted for 4.33% to 13.21% of the overall MSE estimate for ALS and 5.05% to 18.64% for DAP (Table 11). Thus, the contribution of the residual variance component was clearly much greater at 1 ha level than for the study area as whole (12 ha), as one would expect. Further, the contribution was also greater than at the intermediate spatial scale represented by the four sub-regions in Case B, see Figure 6A, which is also reasonable. Further, it should be noted that there did not seem to be a clear relationship between the height of the vegetation and the magnitude of the residual variance component. For ALS, for example, the greatest residual variance components were found for cells #2, #3, and #11 for which the vegetation height estimates were at the low as well as the high ends while the very smallest (cell #5) and very largest (cell # 1) height estimates showed smaller residual variance components. This may suggest that the relative magnitude of the residual variance may not be sensitive to the particular properties of the vegetation in the ecotone.
In an operational context, most spatial units subject to estimation will lack field observations and residuals cannot be estimated. The current results suggest that ignoring the residual variance component will lead to an underestimation of the overall MSE estimate by, say, 10-15% at a 1 ha level. The order of magnitude of this underestimation seems to be similar for ALS and DAP, and is somewhat greater than the~5% underestimation for change in vegetation height using ALS reported by Ref. [12].
The intention of Case B was to offer an opportunity to explore and understand how vegetation properties affected the performance of estimation based on the ALS and DAP data, respectively. As is evident in Table 11 and Figure 6G,H, the vegetation height estimates followed the same crude trends for ALS and DAP and the estimates were logical in the sense that they reflected the subjective criteria used when delineating the four sub-regions ( Figure 6A). Thus, it was rather surprising that the greatest differences in vegetation height estimates between ALS and DAP were found for one of the regions with low vegetation (sub-region I; difference of 0.20 m) and one of the regions with tall vegetation (sub-region IV; difference of 0.19 m). Likewise, the smallest differences in vegetation height estimates were found for one of the regions with low vegetation (sub-region II; difference of 0.06 m) and one of the regions with tall vegetation (sub-region IV; difference of 0.06 m). In fact, for these two latter regions the ALS-and DAP-based estimates were not significantly different in the statistical sense. We could not find any other properties of the sub-regions that could help explaining the differences in the estimates.
The observed patterns of the differences in estimates based on ALS and DAP therefore largely remain unexplained and further studies, preferably with a broader range of image acquisition settings, such as for example image overlap, and different parameter settings for the DAP would probably be useful.
During the analysis under objective #1 it was observed that solitary trees tended to be smoothed out in the DAP data (Sections 3.1.1 and 3.1.2). Similar findings have recently been reported even for large, scattered forest trees over large areas subject to operational forest inventory based on DAP and the area-based method [63], despite the generally good performance of DAP in denser and tall forests for estimation of a number of biophysical forest variables, including height (e.g., [20]). Along with the unexplained patterns of the height estimates when comparing ALS and DAP, this suggests that there is a critical need for more in-depth and fundamental understanding of the behavior of DAP in studies of trees and forests-from biomass-rich and productive forests in the lowlands to small trees in the boreal-alpine ecotone.

Conclusions
As one would expect, the results confirmed that the probability of obtaining positive height values for small individual pioneer trees will increase with increasing point density of the 3D remotely sensed data. Thus, DAP data with 10 times as great point density as the ALS data showed a significantly greater probability of obtaining positive height values for trees. Still, it was surprising and contrary to our expectation, that there was a more pronounced underestimation of tree height with DAP data than with ALS data. This tendency was stronger for pine than the other species, which was counterintuitive, given our expectation that DAP points clouds would capture the properties on the outer surface of a more regular object as represented by a pine tree crown and thus could resemble the surface of the crowns. We noticed that particularly solitary trees tended to be smoothed out. We were unable to explain these results, and general conclusions on the species effect could therefore hardly be drawn. We were unable to identify whether this was a general phenomenon, or a result of the particular DAP algorithm and parameter values adopted in this study. Although not documented in this article, it should be noted that numerous different settings for the image matching were tried in a preliminary phase of the analysis without any substantial impact on the results. It was revealed that tree height was predicted with significantly smaller residual variance when using ALS data, which gives greater promise for estimation of vegetation height with ALS than with DAP data.
The proposed method for estimation of vegetation height of small trees and other vegetation has the potential to produce very precise estimates. The analysis suggested that most of the mean square error estimates (>85%) of the estimators will be accounted for by quantifying the variance attributable to the model parameter uncertainty when the size of the target areas subject to estimation is as small as 1 ha. For larger areas, the model parameter uncertainty will account for an even greater portion of the total uncertainty, as demonstrated for the entire 12 ha study area (>99%). The most precise estimates were found for ALS. Statistically significant differences were found between ALS-based and DAP-based height estimates for certain sub-regions of the study area. Because we were unable to explain the causes of these differences, caution should be exercised regarding interpretation of differences in precision as well as systematic differences between use of ALS and DAP data for vegetation height estimation. Further investigations are needed to understand how DAP data in particular behave for small trees and low vegetation in the boreal-alpine ecotone under different image acquisition and DAP settings. This is of great importance for future operational monitoring of vegetation change in the ecotone since multi-temporal applications often will have to rely on combinations of 3D data from ALS and DAP from different points in time. Even small systematic effects of a particular technology on height estimates may compromise the validity of a monitoring system since change processes encountered in the boreal-alpine ecotone often are subtle and slow.  Data Availability Statement: The field data are not publicly available due to privacy of the private landowners. The DAP data are available on request from the corresponding author. The ALS data (3rd party data) are available at https://hoydedata.no/ (accessed on 22 June 2021).