Forest stand height predicted from ICESat-2 ATLAS data for forest inventory and comparison to airborne laser scanning metrics

The study includes the entire Estonia except its western islands (Figure 1), covering an area of 40,000 km². Estonia is located in north-east Europe with its centre coordinates found approximately at 25.6° E, 58.8° N. The ground surface topography is flat with ground elevation rising from the north and west towards the south and east (Arold, 2005). Recognizable hills with local relative height greater than a few dozens of metres can be found in the south and south-east Estonia (Arold, 2005).

Figure 1.

Density of ICESat-2 ATLAS measurements and years of airborne laser scanning (ALS) data. Labels for the Estonian ALS data are “S” for summer and “K” for spring.

Joonis 1. ICESat-2 ATLAS mõõtmiste paiknemine ja Maa-ameti aerolaserskaneerimise aastad. Märgised näitavad suvist ja kevadist aerolidari andmestikku.

The forests in Estonia are dominantly composed of evergreen coniferous Norway spruce (Picea abies (L.) H. Karst.), Scots pine (Pinus sylvestris L.), and broadleaf deciduous European aspen (Populus tremula L.), silver birch (Betula pendula Roth), downy birch (B. pubescens Ehrh.), black alder (Alnus glutinosa (L.) Gaertn.), and grey alder (A. incana (L.) Moench). A diverse nomenclature of soil types is present in forest land ranging from mineral soils (Leptosols, Cambisols, Luvisols, Rectisols, Luvisols) to Histosols (Kõlli et al., 2004). Abundant grasses and bushes can be found below forest stands growing on fertile soils.

According to the Estonian National Forest Inventory the forest cover in Estonia is 53.5% (Valgepea et al., 2023). Mature forest height remains usually below 35 m (Tappo, 1982). Multi-layer forests are common with Norway spruce and other shade-tolerant species growing below the upper canopy. The average stand size is 1.24 ha, with larger stands found more frequently in state forests (Valgepea et al., 2023). The forests usually develop a dense canopy cover which is modified with thinnings and maintenance fellings in managed forests. Retention forest management is dominant (Gustafsson et al., 2012) with retention trees left growing on regeneration felling areas, including so called clear-cuts. The maximum allowed clear-cut area is 2–7 ha depending on the forest type (Forest Act, 2023) and this yields and maintains a diverse forest landscape. The share of protected forests with limited or entirely prohibited management interventions is 30.3% (Raudsaar et al., 2023).

The Estonian Forest Register is maintained by the Environment Agency and is a database of forest management inventory data. The database contains tabular data and stand border maps. The database was used to identify forest stands that are sufficiently large to contain at least two ICESat-2 ATLAS observations. The stand boundaries were shrunk by 20 m to exclude edge effects and ICESat-2 ATLAS pixels that extend over different stands. The forest register was used to assign a data flag indicating dominance of evergreen coniferous tree species in ICESat-2 ATLAS pixel locations. Only the ICESat-2 ATLAS pixels were used in the following analyses that had been inventoried in the year of ATLAS measurements.

ATLAS is the only instrument for Earth observation on board ICESat-2. ATLAS is a photon counting laser scanner that emits green light (532 nm) pulses at 10 kHz frequency (Neumann et al., 2019). The pulses are narrow (<1.5 ns) and have an about 17 m footprint size on the ground at 500 km average orbital altitude of ICESat-2. The combination of pulse frequency and spacecraft velocity of about 7 km s⁻¹ yields a ground sampling distance of 0.7 m on the ground.

The laser pulse energy is split into six exiting laser beams from which three strong beams have an energy ratio of about 1:4 compared to weak beams. The distance on ground of weak and strong beams in a pair is approximately 90 m and the three footprints of pairs have an about 3.3 km separation distance in the across track direction (Neuenschwander & Pitts, 2019). At the equator, the distance between the reference ground tracks (the imaginary lines through the six-beam patterns) is 28.8 km. The areas between the reference ground tracks are covered by ICESat-2 ATLAS by the changing off-nadir pointing angle of the satellite after every 91 days (Markus et al., 2017).

The orientation of ICESat-2 is changed 180° about every 9 months to direct radiators away from the sun and to maximize illumination of the solar panels (Brunt et al., 2019) and this operation swaps the strong and weak beams in tracks GT1L, GT1R, GT2L, GT2R, GT3L, and GT3R because the light path inside the ATLAS is fixed.

The canopy height product is quantified within discrete 100-meter segments. Concurrently, canopy height values are derived for smaller 20-meter segments (pixels) from the same measurement dataset. The geographic coordinates corresponding to each canopy height measurement are centred within the 100-meter segment from which the measurement originates.

In our analysis, ICESat-2 data (version 006) were acquired from the NASA Earth-data website for Estonia during the years 2019–2022, focusing on the summer and early autumn months (June, July, August, and September). Our examination concentrated on the juxtaposition of the canopy height data obtained at 20-meter segments with the information contained within the Forest Registry. Only ICESat-2 data points situated at a minimum of 20 metres inward from the perimeter of forest compartments were chosen for analysis. Furthermore, the ICESat measurements and the inventory dates recorded in the Forest Registry had to originate from the same calendar year. To ensure adequate spatial representation, a minimum of two ICESat pixels per compartment were required for inclusion in the analyses. Adhering to these criteria, the mean ICESat-derived canopy height per parcel was computed and juxtaposed with the corresponding height recorded in the Forest Registry.

The ALS measurements are done by the Estonian Land Board using Riegl VQ-1560i (RIEGL Laser Measurement Systems GmbH, Riedenburgstraße 48, A-3580 Horn, Austria) operated in discrete return registration mode. In this study we used data from the years 2018 until 2022 (Figure 1). The point clouds have a density of 0.8 m⁻² in summertime measurements and 2.0 m⁻² in spring measurements which result in a lower flight altitude. The ALS pulse footprint in nadir direction is 0.56 m in summertime data. The springtime ALS measurements start after snowmelt and are carried out for one quarter of Estonia. Usually, the measurements are completed before final leaf unfolding of the trees. The summertime ALS measurements are carried out after final leaf unfolding when the springtime project is completed.

The ALS data processing was done using FUSION/LDV (McGaughey, 2020) tools. The ground surface elevation model was obtained from the Estonian Land Board. All pulse returns with a height over ground greater than 1 m were used for height distribution metrics. First returns were used to calculate canopy cover (Lang, 2010; Jennings et al., 1999) (K_ALS, %) at a 2 m reference height over ground. All ALS metrics were calculated for 20 m that extended over each dataset shown in Figure 1. The constructed map pixels are about the same size as the ICESat-2 ATLAS ATL08 pixel on the ground. The ATL08 pixel boundaries were transformed into the Estonian base map coordinate system EPSG:3301. The ZonalStatistics plugin in QGIS 3.18 was used to extract ALS metrics for each ATL08 pixel.

Preliminary analysis showed that the best matching ALS metrics to ATL08 forest height from ALS was the 95^th percentile (H_ALS) of point cloud height distribution. The finding was in concordance with conclusions made by Neuenschwander & Pitts (2019). However, before further seeking the possible influential factors that may be significant for the ATL08 forest height (H_ICESat) determination, the combined ALS and ATL08 dataset had to be cleaned from outliers. Filters were imposed on the data-set to remove pixels: (1) where ATL08 forest height was greater than 35 m, (2) that were located within forest stands where the K_ALS and H_ALS decreased substantially between ALS measurements (ALS data from 2018 was used for this filter), (3) that were located in stands where the range of H_ALS was greater than 8 m, and finally (4) that were in stands with less than one remaining ATL08 pixel. About 100 pixels were removed manually according to scatterplots of data. The outliers were mostly found in stands with larger retention trees and low canopy cover.

After removal of outliers the dataset contained 12,711 observations from summers in 2019–2022 with H_ICESat mean value in the range of 15–17 m which in turn was very similar to the mean forest height (H_FI) estimated from the forest inventory data (Table 1). The relative standard deviation around the mean values was 30–40% of mean height which is explained by the diversity of stands in the sample. The determination coefficient was 0.85 for a linear model fitted for 3,065 stands between stand mean H_ICESat and H_FI with the model residual error remaining 2.7 m. The average H_ICESat was about 0.3 m greater than H_FI, and the difference may be a combined result from forest height definitions and the underlying methods in both datasets.

In the analysis, subsets of observations were drawn from the database with the condition of equality of the ALS measurement year and ICESat-2 ATLAS measurement years. Since only summertime data was used from the ICESat-2 ATLAS ATL08 dataset, the measurement month information was not used in data partitioning.

We assumed linear relationships between H_ICESat and H_ALS. Following the observations made by Moudrý et al. (2022) about the influence of beam energy on canopy height predicted from ICESat-2 ATLAS data, we first analysed the dependence of the difference H_ALS-H_ICESat on the measurement track. With the help of box-plots (Figure 2) it was confirmed that the beam energy in combination with the track may have a small systematic influence, and track label (six levels) and also the beam energy level factor (according to the sc orient parameter in ATL08) were included into models. While canopy cover may have an influence on the H_ICESat according to Neuenschwander & Pitts (2019) the ALS-based K_ALS was included into models. The factor of coniferous evergreen species dominance on pixels was included to account for deeper pulse penetration into the canopy in case of springtime ALS measurements. Because forest canopy structure may depend on growth conditions, we used a factor discriminating between deep peat soils from other soils. We used a threshold of 15% of organic carbon content in soil imposed to the dataset published by Kmoch et al. (2021). The threshold separated the deep Histosols from other soils well when results were checked against other soil maps in Estonia and forest site type codes in the forest inventory database.

Figure 2.

Examples of the influence of ICESat-2 ATLAS scanning track on height difference of canopy surface height calculated from airborne laser scanning (ALS) data and ICESat-2 data ATLAS. The legend shows the ALS measurement time and number of observations. Labels for ALS data are “S” for summer and “K” for spring. Only strong impulses are shown for 2019; in 2020 and 2022 individual tracks contained only strong or weak impulses.

Joonis 2. Näiteid satelliidi ICESat-2 ATLAS skaneerimisraja mõjust satelliidilt ja lennukilt laserskanneriga (ALS) mõõdetud taimkatte pinna kõrguse erinevusele. Legendil on toodud ALS mõõtmisaeg ja ICESat-2 ATLAS pikslite arv. 2019. aasta suve andmetes on näidatud ainult tugevad impulsid, 2020. ja 2022. aasta andmetes sisaldasid rajad ainult nõrku või tugevaid impulsse.

Model equations (1) – (7) were constructed where b_x are estimated parameters, S_species is a factor for evergreen coniferous dominance, P_soil is an indicator of deep peat soil, T_ICESat.trac is the factor for ATLAS track, T_ICESat.beam is a binary factor separating strong and weak energy beams in the ATL08 dataset, and ε is the model residual error. The model parameters were estimated with least squares fitting in R version 4.2.2 (R Core Team, 2022). The determination coefficient, model residual error and significance of parameters were used to assess the impact of factors that may explain the relationships between the observed H_ALS and H_ICESat. (1) HICESat=b0+b1HALS+ε, {H_{ICESat}} = {b_0} + {b_1}{H_{ALS}} + \varepsilon , (2) HICESat=b0+b1HALS+b2KALS+ε, {H_{ICESat}} = {b_0} + {b_1}{H_{ALS}} + {b_2}{K_{ALS}} + \varepsilon , (3) HICESat=b0+b1HALS+b2KALS+b3Sspecies+ε, {H_{ICESat}} = {b_0} + {b_1}{H_{ALS}} + {b_2}{K_{ALS}} + {b_3}{S_{species}} + \varepsilon , (4) HICESat=b0+b1HALS+b2KALS+b3Sspecies+b4Psoil+ε, {H_{ICESat}} = {b_0} + {b_1}{H_{ALS}} + {b_2}{K_{ALS}} + {b_3}{S_{species}} + {b_4}{P_{soil}} + \varepsilon , (5) HICESat=b0+b1HALS+b2VICESat.beam+ε, {H_{ICESat}} = {b_0} + {b_1}{H_{ALS}} + {b_2}{V_{ICESat.\rm{beam}}} + \varepsilon , (6) HICESat=b0+b1HALS+b2KALS+b3VICESat.beam+ε, {H_{ICESat}} = {b_0} + {b_1}{H_{ALS}} + {b_2}{K_{ALS}} + {b_3}{V_{ICESat.{\rm{beam}}}} + \varepsilon , (7) HICESat=b0+b1HALS+b2KALS+b3Sspecies+b4TICESat.track+ε {H_{ICESat}} = {b_0} + {b_1}{H_{ALS}} + {b_2}{K_{ALS}} + {b_3}{S_{species}} + {b_4}{T_{ICESat.track}} + \varepsilon

On average, the ground surface elevation in ICESat-2 ATLAS ATL08 (Z_ICESat) was 18.92 m greater (standard deviation 1.21) than the Estonian digital ground surface model (Z_ALS). The systematic difference is caused by geodetic datums, because there was almost a functional linear relationship (8) ZALS=−19.33+1.006ZICESat+ε, {Z_{ALS}} = - 19.33 + 1.006{Z_{ICESat}} + \varepsilon , with the residual standard error of 1.2 m and determination coefficient R²=0.9982 at 3476 degrees of freedom. However, scatterplots (Figure 3) revealed that there are systematically deviating clusters from the main relationship.

Figure 3.

Difference between the Estonian digital terrain elevation model (DEM) and the model used for ICESat-2 ATLAS ATL08 as a function of longitude and latitude (expressed with colour).

Joonis 3. Eesti ja ICESat-2 ATLAS ATL08 maapinna kõrgusmudelite erinevus sõltuvalt geograafilisest asukohast (värv vastab laiuskraadile).

We found that the geographic latitude λ and longitude φ of pixel location explained most of the clustering shown in Figure 3 when inserted into the regression model (9) ZALS=259.9+0.9984ZICESat−14.45φ−5.051λ+0.2579φλ+ε. {Z_{ALS}} = 259.9 + 0.9984{Z_{ICESat}} - 14.45\varphi - 5.051\lambda + 0.2579\varphi \lambda + \varepsilon .

The model (9) had R²=0.9997 and a much decreased residual error of 0.51 m. Other variables like soil type, the share of coniferous trees, and forest canopy cover were not significant if added additionally into the model (9). This residual error indicates that ground elevation errors in the ICESat-2 ATLAS ATL08 dataset do introduce some uncertainty into the canopy height predictions. The error is comparable to grass layer and shrub layer height in forest under-storey. The reasons for outliers in Figure 3 were not explained by the variables in our datasets. For comparison, the Z_ICESat was converted to the reference system of Z_ALS that is based on the EST-GEOID2017 model (Ellmann et al., 2019). We used a converter published on the Estonian Land Board web page. The Z_ICESat remained on average 15 cm higher than Z_ALS with a standard deviation of 41 cm. This is about the same uncertainty as shown by the model (9).

The airborne lidar data-based canopy height explained 89–95% of the variability in canopy height prediction in the ICESat-2 ATLAS ATL08 dataset according to the determination coefficient of linear models (Table A1.1 in Appendix 1). The residual error remained between 1.6 m and 2.0 m. Surprisingly, the additional variables that described forest canopy cover, site conditions and ICESat-2 ATLAS measurement geometry had only marginal power to explain the rest of variability. Even if parameter estimates for the additional model variables were frequently significant (Table A1.2 in Appendix 1), the increase in the model determination coefficient was marginal and model residual error decreased less than 10–20 cm. Neither the visual inspection of scatterplots (Figure 4) nor inclusion of pixel location data into models revealed any dependence of the variability in the geographic location of pixels as it was found in the case of ground elevation models. The unexplained residual scatter is probably related to the small number of registered backscattered photons (0–4) per emitted pulse (Neuenschwander et al., 2022).

Figure 4.

Examples of the 95^th percentile of airborne laser scanning (ALS) point cloud height distribution and ICESat-2 ATLAS canopy surface height by year and growth season when ALS was done. The symbol colour corresponds to the ICESat-2 pixel location longitude.

Joonis 4. Lennukilt (ALS) ja kosmoselidariga ICESat-2 ATLAS tehtud lasers-kaneerimise andme-test arvutatud puistu pinnamudeli võrdlus ALS mõõtmisaja järgi. Sümboli värv kirjeldab ICESat-2 piksli asukoha geograafilist idapikkust.

Canopy cover was almost always a significant variable if added alone or with others into the models. The share of evergreen coniferous trees in stand species composition was also significant in most of the models. The indicator of peat soil was significant, but the significance was much less compared to canopy cover (Table A1.2 in Appendix 1). The beam energy indicator variables were significant for measurements from the year 2020 and later. Most likely, the residual error is caused by the uncertainties in the ICESat-2 ATLAS pixel location coordinates, noise in the photon counting, and errors in determining ground surface elevation. Part of the unexplained residual variability is related probably to the optical properties and structure of trees combined with the spatial location pattern of the trees on the ICESat-2 ATLAS pixel.