Using GEDI lidar data and airborne laser scanning to assess height growth dynamics in fast-growing species:a showcase in Spain

Juan Guerra-Hernández and Adrián Pascual

Abstract

Keywords: Forest inventory,Productivity, Multi-temporal LiDAR, Full-waveform (FW), Satellite mapping

Introduction

Forest inventories employing airborne laser scanning(ALS) data have become common in many countries(Wehr and Lohr 1999; Lefsky et al. 2002; Næsset et al.2004;White et al.2016).The ALS-based forest inventory methods are typically categorized into two groups: the area-based approach (ABA) (Næsset 2002; Guerra-Hernández et al. 2016) and individual tree detection(ITD) (Hyyppä et al. 2001). So far, most operational forest inventories framed under ALS data utilization have been designed under the ABA method (Maltamo et al.2014). Under the ABA method, stand attribute models are fitted with sample plots using metrics calculated from ALS data and then these models are used to predict stand attributes wall-to-wall for whole inventory area, typically using square grid cells whose area exactly matches the size of the sampling plots as inventory units(Næsset 2002). As a result, the elevation of tree vegetation or the predicted forest inventory variables can be expressed as a spatially continuous raster map. The science and the implementation of ALS-based forest inventory is solid, evolving fast and mainstream worldwide(Maltamo et al. 2014).

During the last decades, nationwide laser scanning has helped forest managers to integrate LiDAR in forest inventory (Nilsson et al. 2017). The impact of ALS data is strong on forest decision-making when data from National Forest Inventories is estimated with remote sensing sources under good co-registration goodness(Pascual et al. 2019). The number of countries that have been scanned at least once relying on ALS technology has increased and in many areas forest dynamics and natural disturbances can be estimated using multitemporal LiDAR (Hopkinson et al. 2008). The time-gap between ALS surveys, the growth rates of the species and forest management actions challenge the management of fast-growing plantations (McRoberts et al.2018). As a showcase, the rotation age for Eucalyptus stands in the Northwest of Spain ranges from 13 to 20 years while the time gap between the last two nationwide ALS surveys for the area was half the rotation (7 years). Therefore, capturing the fast growth rates under commercial silviculture can benefit from independent and flexible sources of 3D data operated on a shorter interval basis compared to ALS nationwide surveys(Silva et al.2017).The launch of the NASA’s Global Ecosystem Dynamics Investigation (GEDI) satellite mission in 2018 has paved the way to the integration of global, multitemporal and full-waveform laser data into operational forestry (Duncanson et al. 2020). The estimation of forest structure using GEDI laser statistics expands the frontiers of forest monitoring and support a change of paradigm on the spatio-temporal evaluation of forest productivity (Bontemps and Bouriaud 2013; Coops et al.2018).

The GEDI missions is continuously scanning over the earth’s land surfaces since 2019 using optical Light Detection and Ranging (LiDAR) to monitor tropical and temperate forests ecosystems (Dubayah et al. 2020a).The billions of GEDI laser shots are expanding the frontiers of forest monitoring worldwide especially in remote forest areas (Silva et al. 2017; Duncanson et al. 2020).The GEDI laser instrument generates total of 8 beam ground tracks that are spaced approximately 600 m apart on cross-track direction. Each beam track consists of ～25 m footprint samples spaced every 60 m along the ground track. The GEDI data science products include arrays of laser waveform metrics for each footprint sample, including canopy height and profile metrics on above-ground height or ground elevation (Hancock et al.2019;Dubayah et al.2020a).The understating and development of solutions based on GEDI data is a promising research topic (Duncanson et al. 2020).One of the applications of the GEDI technology is the enhancement of large-scale forests maps generated with ground-based National Forest Inventories (NFI).

The aim of the investigation was to test the capability of GEDI to support ALS scanning when it comes to describe dynamics in fast-growing forest ecosystems. The 2.5-year time gap between ALS surveys and the first GEDI orbits in the Northwest of Spain reflects the scenario in which GEDI laser shots can be used to update forest structure and detect spatio-temporal changes.This investigation is one of the first attempts to validate the capability of GEDI to detect spatial and temporal changes using nationwide ALS surveys as baseline. We tested the capability of GEDI data science products under challenging fast-growing commercial forest ecosystems and the small-scale forestry carried out in Galicia, a mosaic of private and public forest patches (2 ha on average).

The objectives of this study were:i)to use GEDI to detect disturbances at regional level, ii) to investigate,evaluate and explain the source of uncertainty when interpreting GEDI footprint ground shots and iii) to validate the height growth dynamics detected from ALS to GEDI compared to ground NFI data. The paper emphasizes on the filtering mechanisms and geospatial operations required to combined ALS data and GEDI data.The GEDI laser data will become a reference to monitor forest dynamics worldwide once researchers and practitioners fully understand the technology and its integration under operational forest monitoring.

Material and methods

Training area

Fig.1 Location of the training area in Spain in the region of Galicia.The extent of the three forest ecosystems assessed in the province of Lugo is displayed using the Forest Map of Spain(FMS)

The training area is the province of Lugo, located within the region of Galicia in the Northwest of Spain(Fig. 1). The landscapes are characterized by rugged orography and mild temperatures with low thermal oscillation between winter and summer and frequent rainfalls. The area is a rural region and forestry is one of the principal economic activities. The conditions of the region boost the high productivity of the main tree commercial species in Galicia and even Spain: Eucalyptus globulus Labill., Pinus pinaster Ait.and Pinus radiata D. Don. Galicia is one of the most important regions in Spain in terms of forestry production. To overview the economics, around 6.5 million m3are harvested every year and the turnover is steadily growing due to the intensive use of forest biomass as supply of the pulp industry, a major driver of the economy in Galicia. Around 30% of the 1.4 million hectares of forest land comprises pure and even-aged pine stands, mainly maritime pine (P. pinaster) (271,281 ha) and radiata pine (P. radiata) (96,(916,71777 ha) (MAGRAMA 2018). The high productivity of P. pinaster and Pinus radiata forests, for instance,allows profitable timber production with short rotation lengths. The integration of laser scanning technology in the region is becoming mainstream in operational forest inventory (Gonçalves-Seco et al.2011; González-Ferreiro et al. 2012). However, the short rotation and the fast dynamics under commercial forestry have constrained the impact of nationwide ALS surveys on the decision-making arena. The use of 3D data collected on shorter schedules can support the industry while giving forest owners of small-scale forest properties the chance to monitor remotely located assets. The fragmentation of the landscape, the fast growth of forest species and the silviculture turn Galicia into a laboratory to test the measuring capability of GEDI laser technology.

The Forest map of Spain

The last version of the Forest Map of Spain (FMS) at 25-m resolution was used to select the forest areas corresponding to the most important commercial species in Galicia (MAPA 2018). The FMS contains wall-to-wall information at polygon level extracted from the National Forest Inventory (NFI) including species occupation, forest stage of development or forest cover, which was assessed by photointerpretation of aerial RGB imagery from the PNOA project (Plan Nacional de Ortofotografía Aérea) and ground NFI data. The integration of data science products from GEDI and ALS surveys can improve the description of FMS units (polygons). The scope of the study was narrowed to a bounding box of 73.2 km width and 140.6 km long (Easting 588,902.2-661,214.47, Northing 4,695,137.4-4,835,977.0 UTM)covering around a million hectares, which is 2% of the country and almost the full extent of the Lugo Province.The research evaluated the three main fast-growing species in Spain: P. pinaster, P. radiata and Eucalyptus spp.The final scale of the assessment after filtering data reached a total of 211,346 ha in three groups: 47,820 ha(P. pinaster), 72,620 ha (P. radiata) and 90,906 ha (Eucalyptus spp.).

Airborne laser scanning data

The ALS surveys are part of the PNOA project (MATA 2019), which allows the use of nationwide laser data for free. The nominal pulse density of the 2015-2017 ALS surveys, was 1.0 pts.·m−2on average with a measuring accuracy of 20 cm in the vertical profile. The scheduled calendar for the flightpaths was available. The scale of the LiDAR processing was regional (Galicia) and the FMS polygons were used to account for the spatial distribution of the three species across the region. The raw LiDAR data files retrieved from the PNOA server comprised 282 GB of data. The licensed software Lastools(Isenburg 2020) was used to process the point clouds and compute 3D laser statistics using Intel®Core i9-7900X workstation with 64 GB of RAM and 10 cores. As a guide for the reader, we recommend the papers from Wehr and Lohr (1999) and Lefsky et al. (2002) to understand the principles of ALS and our previous study on ALS processing steps using large-scale 3D point clouds(Pascual et al. 2020). Briefly, lasheight were used to normalize the classified point cloud of ALS echoes. The lascanopy command was used to extract the 99th height percentile and CC (canopy cover) metrics from the normalized point cloud using a circle of 25 m of diameter for the center of each GEDI footprint ground shot, while lasclip was used to extract the normalized point cloud for the extent of FMS polygons. Then, lascanopy was also used to generate a 25-m raster grid of the 99th height percentile for the training area. The height_cutoff and cover_cutoff parameters were set to disregard ALS echoes below 2 m when computing ALS statistics. The 25-m grid raster for the 99th height percentile (ALSH99)was generated for the three species within the extent of the bounding box created to retrieve GEDI data. The grid scale of 25 m matched the ground footprint of GEDI laser shots (Dubayah et al. 2020b) to avoid the bias resolution when comparing laser statistics.

Filtering operations were very important along the research project as they were applied to process the ALS surveys, to process the GEDI laser shots and to combine both while acknowledging the edge effect of the irregular boundaries of FMS polygons and the resolution dependence of the computed statistics. The notation for filters along the manuscript follows their chronological implementation in the research workflow.Shots GEDI not completely contained within FMS polygons in the training area were removed to avoid edge effects(Filter#1).Statistics from the 3D point clouds were computed for each grid cell. Growth dynamics under multi-temporal laser scanning are better detected when using the highest height percentiles on ALS echoes distribution(Hopkinson et al. 2008; Socha et al. 2020). The proper benchmarking between GEDI and ALS statistics should be carried out using ALS point cloud versus GEDI FW to avoid resolution dependance and the averaging effect on ALS statistics when using raster pixels (Packalen et al. 2019). From the center points of the valid GEDI laser shots considering the 25-m footprint, we clip the 3D point cloud using only the ALS echoes within the ground GEDI footprints.Hereby,the 99th height percentile(ALSH99)from the ALS data was computed for each GEDI shot.

The GEDI satellite LiDAR

The bounding box was used to select the GEDI beams within the 1-million-ha training area. The rGEDI package(Silva et al. 2020) was used to retrieve and process the GEDI data under the 3.6 version of the R statistical software (R Core Team 2020). There were 31 GEDI orbit tracks in h5 format available to date March 15th, 2020(Fig. 1). The GEDI laser shots within the bounding box considered in the research were collected in 2019 between April 21st and August 3rd. The GEDI laser equipment uses 3 lasers to produce 8 laser beams creating laser shots every 60 m along the longitudinal track in the ground.The distance in the cross-track direction for two consecutive longitudinal tracks is 600 m for a given orbit. When combining GEDI shots from different acquisition dates and orbits, the cross-track separation between tracks decreased. For instance, the study area contained distances in the cross-track direction down to 60 m when we merged the multi-temporal GEDI shots collected in ten different days along the specified time interval. To date,the stage of development of the GEDI mission restricts the possibility of wall-to-wall mapping across our training area (Hofton and Blair 2020) for which the level of landscape fragmentation is high(Fig.2).

The use of GEDI Level 3 to produce wall-to-wall estimates should be carefully assessed when evaluating small-scale properties as occurs in Galicia. Hereby, the authors framed the research using the elevation statistics in Level 2A data science product (Dubayah et al. 2020b).The full-waveform (FW) for each GEDI laser shot was processed to compute the following statistics for each of the 25-m footprint ground shot. The focus of using Level 2A was to compute the same height percentiles as with the ALS data. The relative height was computed along 100 unitary intervals from the minimum to the maximum elevation (RH100), and RH99was the selected statistic on which to benchmark ALS data.

Geospatial operations and statistics with GEDI, ALS and the FMS

Fig.2 Partial overview of the P. pinaster forests in the training area showing the aerial satellite image and canopy height model derived with airborne laser data in the background.The ground beams of the GEDI shots(GEDI_level2A)intersecting forest areas are presented

Fig.3 Spatial layout of Filters#2-8 use to disregard GEDI shots displayed sequentially. a Filter #2,b Filter#3, c Filter#4-6 and d Filter #7-8.The maps correspond to forest areas classified as P.pinaster. The raster map for ALSH99 (99th height percentile,LP99 in the legends)derived from ALS surveying is shown over the 2017 orthophoto

The benchmark between ALS and GEDI was carried out within the boundaries of FMS polygons intersecting GEDI ground tracks (Fig. 3a). For a given GEDI laser shot, we derived the RH99from the 25-m footprint ground shots (GEDI99), the 99th percentile of the ALS echoes height distribution within the extent of the 25-m diameter samples (ALS99). From the initial set of 288,529 laser shots within the bounding, we applied the following filters:

1) Selection of GEDI shots completely contained within the species-specific FMS polygons(Filter #2).

2) Selection of GEDI shots located further than 30 from the edge of the ALS-based raster map boundaries (Filter #3).

3) Selection of GEDI shots whose GEDI99was above 70 m (Filter #4).

4) The column quality flag was used to disregard all GEDI shots classified as 0 meaning (Filter#5) the technical and quality attributes for the given shot number were not according to standards. See detail of interpretation of L2A Quality Flag in Dubayah et al. (2020b). A quality flag value of 1 indicates the laser shot meets criteria based on energy,sensitivity,amplitude, and real-time surface tracking quality.

5) All pairwise values of GEDI99and ALSH99ranged between 2 and 70 m to remove outliers (Filter #6).The upper bound was set based on the ground NFI tree height measurements across Galicia for the three species.

6) We disregard GEDI shots for which the canopy cover (CC)associated ALSH99was below 70% to reduce the impact multi-layered and sparse forest conditions on the GEDI laser reflectance homogeneity (Filter#7).

7) The average height increment (ALSH99-GEDI99)and its standard deviation was computed for FMS polygons comprising 3 or more GEDI shots within their boundaries. Preliminary tests showed the importance of this filter to remove outliers for very large positive increments (Filter #8).

The final selection of GEDI shots were used to compute height differences from GEDI to ALS statistics corresponding to the ground extent of GEDI shots. As a showcase, we display the geospatial operations from Filter #2 to Filter #6 for a sample of P. pinaster stands in the training area (Fig. 3).

The estimation of height growth dynamics required from the exact date for which ALS data was collected for all GEDI points. For the set of valid GEDI shots, we computed the difference in days and in years between ALS survey acquisition retrieved from PNOA server(2015-2017) and time of the corresponding GEDI FW was registered (Dubayah et al. 2020b). The difference between GEDI99and ALS99_3Dwas calculated for each GEDI shot available after Filter #8. The assessed mean and standard deviation values for height increments at FMS polygon level were converted to annual height increments for better interpretation. The information for the time intervals is presented in Table 1.

The sign of the difference was used to classified GEDI points as candidate metrics for harvests (negative) and natural dynamics (positive). The mean value and the standard deviation of height increments computed from ALS to GEDI was computed at FMS polygon level. The average values were compared to ground validation data and the proposed thresholds to cut off very large positive increments. The tendency and dispersion of the height increments was assessed in preliminary tests toward the effect of slope, the area of the FMS polygons and goodness of the DEM (Digital Elevation Model) derived from both laser sources considering ALS as the reference value. The results were visually displayed to confirm the capability of GEDI lasers to detect harvests (i.e., negative values higher than 5 m per year), fast growing dynamics(positive increments) while the transition region denoted as “uncertainty” was carefully inspected using satellite imagery to confirm the occurrence of cuttings and addressing the effect of intra-polygon variability in the FMS. The dataset was enlarged by adding polygon area,elevation, slope in degrees and the near distance from each GEDI point to the boundary of the FMS polygon(greater than 25 m after applying Filter #3). The geospatial operations were iteratively implemented in the R statistical software (R Core Team 2020). The captured laser-based increments were validated using an independent set of high-quality ground NFI observations.

Validation data for ground elevation

For the final DEM surface generation, the ALS ground points were converted to a TIN surface raster using 2-m raster DEM. A memory-efficient-streaming technology under three parallel processes was computed using the blast2dem command available in Lastools (Isenburg 2020). For ground ALS, the average of the 2-m raster grid cells was computed at 25-m diameter footprint level, while the elev_lowestmode attribute in GEDI Level 2A provided the ground elevation of GEDI laser shots.The vertical accuracy of ground elevation was evaluated using the differences between ground ALS and ground elevation GEDI (elev_lowestmode) within the 25-mdiameter footprint. The following statistics was calculated: correlation coefficient of Pearson (r), the Root Mean Squared Error (RMSEz, Eq. 1), mean of the differences referred to as biasz (Eq. 2), and relative bias(rbiasz) (Eq. 3).

Table 1 Summary information for the three species on time intervals between ALS adquisition date and GEDI laser shooting computed for the set GEDI shots analysed

Table 2 Summary of the 256 samples from National Forest Inventory of Spain (NFI) located in Lugo Province and used for validation

where n is the number of GEDI shoots, yiis the mean elevation of ALS raster surface model for the 25-m diameter footprint GEDI shot i and ^yiis the ground elevation estimation from the 25-m footprint GEDI defined by the elev_lowestmode attribute in the GEDI Level 2A.

The accuracy of ground elevation at footprint GEDI level was also tested via linear regression model using DEMALSas variable dependent as ground truth value(Eq. 1). To determine the accuracy of ground elevation at footprint GEDI the following statistics were computed from the model: model efficiency (adjR2, Eq. 2), the overall root mean square error (RMSE, Eq. 3), the relative root mean square error (rRMSE, Eq. 4).

where n is the number of GEDI shoots as the total number of observations used to fit the linear regression,yiis the mean elevation of ALS raster surface model from the 25-m footprint GEDI shot i and ^yiis the estimated ground elevation from the fitted linear regression model using elev_lowestmode the GEDI Level 2A as dependent variable.

Validation data for height growth dynamics

The rotation of the NFI in the Northern and most Atlantic provinces of Spain is set to 5 years, half or a third of the normal rotation for a province (10-15 years). The last two rounds of the NFI in Galicia are denoted as 4th and 5th NFI, respectively. We used the permanent NFI samples located in the province of Lugo for the three species to compute the annual height increment using ground tree measurements. The Spanish NFI protocol for the tree measurements can be found in Álvarez-González et al. (2014) and in previous research by the authors (Pascual et al. 2021; Guerra-Hernández et al.2021). The ground validation of heights increments for Eucalyptus globulus can be substantially biased with NFI data so we restricted the validation to P. pinaster and P.radiata using 198 and 60 plots, respectively. Data for the 4th NFI data for P. pinaster and P. radiata were measured between March and September 2009. The measurement of the 5th NFI was between May and October 2018. The difference in dominant height between the inventories ranged between 0.01 and 1.471 m per year in NFI samples for P. radiata, while the upper growth bound for P. pinaster samples ranged from 0.01 to 1.1 m. The maximum height increments were used to filter GEDI ground shots whose height dynamics (i.e.,ALSH99subtracted from GEDI99) was strictly positive.The maximum values for the ground data were validated using studies in Galicia for the species (Diéguez-Aranda et al. 2009, 2012). The validation data is presented in Table 2 and the flowchart of the proposed methodology is shown in Fig.4.

Fig.4 Flowchart of the research showing the collection and integration of laser data from both airborne and GEDI satellite laser into the operational forest inventory to estimate height growth dynamics

Results

Comparison of ALS and FW lidar-derived ground elevation

Fig.5 a Scatterplots of LF and ALS-derived DEM elevation at level of footprint GEDI from Filter #5(quality flag value=1,n=9875,Filter #5),b Boxplot of elevation differences between GEDI (elev_lowestmode)and DEMALS_2m considering terrain slope classes at GEDI footprint level.The upper and lower whiskers extend from to the highest and lowest value respectively within the 1.5 times the interquartile range

The FW ground elevation and ALS lidar-derived DEM(2 m resolution) data was strongly correlated (r=0.99;n=9875, Fig. 5a). In terms of RMSEz, the vertical accuracy across study area was 4.48 m (rRMSEz=0.79%)using ALS data as ground truth value. Mean differences between FW-derived elevation (elev_lowestmode) and ALS-derived elevation (denoted as DEMALS_2m) was 0.18 m (biasz=0.18 m) at footprint GEDI level. FW ground elevation estimations were positively biased meaning FW predominantly underestimated ground elevation when compared with ALS-based DEM. Linear regression improved the accuracy of FW ground height estimations from ALS data in terms of RMSE (RMSE=0.23 m and rRMSE=0.041%). Although slope effects also increased the dispersion values of differences between FW and observed ALS ground elevation (Fig. 5b), we did not observe a clear patter and substantial improvements in our results at FMS polygon level adding a filter based on steep conditions (Fig. 5b).

Filtering and detection of fast-growing areas and harvests

The use of filters significantly constrained the pool of valid GEDI shots: after Filter #2 the total set comprised 56,235 shots:18,269 for Eucalyptus spp.,13,527 for P.radiata and 24,439 for P.pinaster.After Filter#8,the final pool of valid GEDI shots decreased by 92.3%, reaching a total of 3566 observations:1218 for Eucalyptus spp.(16,313 ha),1726 for P. radiata (16,710 ha) and 622 for P. pinaster (7445 ha).The array of FMS polygons reached 384 observations.

The negative height increments computed from the 2015-2017 ALS surveys to the 2019 GEDI scanning showed similar patterns for the three cases (Fig. 6). Negative increments classified as “harvest” corresponded to tree elevations above 15 m in ALS and those cases corresponded to low height values(i.e.,early stages of stand development)in 2019 when computing GEDI shot information. The negative annual height threshold of 5 m was sufficient to detect homogenous FMS polygons entirely harvested(Fig. 6). The mean annual height increment calculated for FMS polygons was displayed using the individual GEDI ground beams to confirm the results in those cases where the harvest coincided with the 2017 orthophotomap available for Galicia in the PNOA repository.The spatial layout of harvested polygons was displayed for the three species(Fig.7).

The positive height dynamics computed from ALS surveys to GEDI scanning at FMS polygon level were, on average, higher comparing the mean values but more similar when comparing median values, especially for P.radiata. The variability of annual growth estimates with the ALS-GEDI method was higher than NFI ground observations (Fig. 8). The median and maximum annual height increments in NFI plots were 0.56 (SD=0.37)and 1.47 m for P. radiata, while for P. pinaster the values were 0.44 (SD=0.24) and 1.11 m, respectively.Under the ALS-GEDI approach the values were 0.62(SD=0.42) and 2.93 m for P. radiata, and 0.61 (SD=0.63) and 2.55 m for P. pinaster, respectively. The values for Eucalyptus spp. were 0.69 and 2.55 m, respectively.The proportion of FMS polygons classified as “growth”decreased by a third when adding Filters #7 and #8.Same as for harvests, we present examples of FM polygons classified as “growth” (Fig. 9).

Uncertainty and scaling when interpreting height growth dynamics

Fig.6 Mean height difference between 99th height percentiles compute from GEDI shots and airborne laser scanning data. The mean values were computed using polygons in the Forest Map of Spain(FMS).The inside-polygon GEDI shots were classified considering the above-ground height derived from GEDI(left),from ALS data(right)and the sign of the height difference.a,b: Eucalyptus; c, d: P.pinaster; e,f: P.radiata

The “uncertain” region as presented in Fig. 6 comprised the greater number of non-positive height increments. We provide a graphical description of the most common sources of error when calculating mean annual height increment. The most recent aerial image from PNOA supported the verification process (Fig. 10) regarding the following identified sources:

1) Variability of forest structure withing FMS polygons, especially for Eucalyptus and P. radiata from intensive thinning,partial clear-cuts and irregular planting.

2) No-adjustment of FMS polygons to harvesting operations and real forest edges.

3) Erratic behavior of variable GEDI99due to the positioning error of shots in sparse forest areas showing presence of gaps and discontinuities between forest patches, which was often the case for P.pinaster stands and P.radiata.

Discussion

The combined use of ALS-based forest maps and the timely GEDI laser data to estimate growth dynamics certainly depended on the spatial co-registration between laser sources. The research focused on the ability to aggregate, scale and re-use existing ALS-based raster data for the assessment of fast height growth dynamics and intensive management at a scale never collected before.The multi-temporal assessment of height dynamics in the ALS-GEDI approach arose the strengths and limitations of data fusioning. It is well accepted that ALS-based forest inventory is solid (Maltamo et al.2014). Therefore, the challenge comes from the integration of GEDI laser data to support forest mapping missions. There are important caveats in the interpretation of this work regarding data filtering, spatial co-registration or the singularities of commercial forestry and forest ownership in the area. In the paper,we highlighted the geospatial operations and filters required to combine GEDI data and large-scale ALSbased mapping.

Fig.7 Showcases of FMS polygons classified as harvests from the height dynamics between ALS(2015-2017)and GEDI scanning(2019).a,b Clear-cuts in Eucalyptus stands;c Unevenly distributed planting for P.pinaster; d Clear-cut for P.pinaster; e Intra-polygon variation on management decisions; f P.radiata thinned stand.The 2017 orthophotomap is the background

Calibration of ground level

The nominal accuracy of horizontal co-registration for GEDI observations is between 8 and 10 m(Luthcke et al. 2020) while in the vertical the accuracy is close to 50 cm (section 4.31 in Dubayah et al.2020b). The vertical accuracy from GEDI sensor (filter criteria quality flag=1) was 4.45 m in terms of RMSEz with a positive bias of 0.18 m (biasz)considering as ground truth the DEM derived from the ALS surveys. Our values were similar to what Hancock et al. (2019) reported (bias of less 0.22 m and RMSE of 5.7 m using ALS as baseline) and higher than values Silva et al. (2018) showed (RMSE=1.64 m, bias=0.29, n=2987, spatial resolution=25 m). The consistent pattern in the GEDI FW sensor for ground elevation was under-estimation, which likely impacted our calculations by means of annual non-positive increments as documented by Hancock et al. (2019):underestimation of upper RH metrics at low canopy cover (see Fig. 7d in Hancock et al. 2019). The integration of GEDI data will benefit from improved horizontal positioning accuracy and locally calibrated algorithms (Luthcke et al. 2000, 2005).

Fig.8 Average annual height growth for Eucalyptus spp. (Ec),P.pinaster(Pp)and P.radiata(Pr)estimated for the ALS-GEDI approach. Mean(dot)and median(bar line)values were computed using GEDI shots positioning within the polygons of the Forest Map of Spain.The results for NFI data were computed using the dominant height values at plot-level in Lugo for the species in the 4th and 5th NFI

The processing of ALS data is always oriented to rescale the elevation of many laser echoes to ground level.The way ground bare land was calculated differed from the ALS surveys to GEDI data. In the first case, the multiple multi-layered ALS echoes per laser pulse allowed the estimation of the DEM using multiple ground points,while for the case of GEDI the estimation of all laser statistics relied on one FW laser pulse per beam. The estimation of ground level from satellite LiDAR is challenging (Duncanson et al. 2020) and the process becomes more complex in sparse forest ecosystems as we observed for the case of sparse P. pinaster stands and low-canopy cover conditions. We used the description of the GEDI geolocation procedure (see Section 6 Luthcke et al. 2000) to interpret our analyses. As previous studies reported (Hyde et al. 2005; Silva et al. 2018),negative increments were attributable to the spatial configuration of canopy elements in footprints located across the transition areas between non-forest and forest areas and multi-layered forests (Fig. 11). The accuracy of ground finder is a significant source of error in FW laser scanning to estimate bare-ground elevation before calculating RH metrics (i.e., above-ground heights). Despite disregarding the addition of a slope filter based on the DTM assessment in this study, slope can affect ground and above-ground height estimates in other ecosystems covering wide altitudinal gradientes (Hofton et al. 2002;Silva et al. 2018).

Height growth dynamics and the training area

The use of Galicia as training area favored the detection of fast-growing dynamics and the occurrence of scattered disturbances. The time gap of 2-3 years between the ALS surveys in Galicia and first GEDI shots was enough to ensure height increments beyond the precision of the laser sensors (Luthcke et al. 2020). Harvests were detected and their evaluation could have been simplified by accounting for the absolute values computed for each GEDI shot. The distinction between harvests and GEDI shots classified as uncertain was according the annual threshold of 5 m. The forest area affected by cuttings was greater but heights increments were computed at FMS polygon-level, the intra-polygon variation of forest attributes and the large size of FMS minimized the detected total harvested area. A representative FMS polygon was a mosaic of small forest patches for which management, ownership and species composition could be diverse. Consequently, many of the FMS polygons allocated in the uncertain region comprised harvested area. One way to detect those conflictive cases is to evaluate the homogeneity of height annual increments within the boundaries.

Fig.9 Areas classified as fast-growing when accounting for height dynamics using laser data.a,b Even aged and homogenous polygons are presented for Eucalyptus spp.;c, d P.pinaster stands;e mature P.radiata plantation;f early-stage P.radiata plantation

The comparison of FMS polygons attending to the variability of the mean annual increments was a good proxy to detect variability across the landscape and detect those management units above a certain variability threshold on which harvests must have occurred above a certain height. The results showed, for instance, how large Eucalyptus plantations can be as homogeneous as small pine stands, and it makes sense if we consider a large plantation as the sum of small homogeneous forest units. The recognition of landscape fragmentation in the FMS was not as accurate as desirable, and that is the main reason for the extent of the “uncertain” level in our results(Fig. 12). The reasoning is supported by previous studies that documented the issue under multitemporal observations based on ALS surveys (Hopkinson et al. 2008).

Fig.10 Showcases of non-positive height dynamics from ALS to GEDI statistics in Eucalyptus areas focusing on canopy structure.Low (blue)and dense canopy cover(yellow,CC>70%)cases are displayed with the corresponding both ALS point clouds and simulated GEDI full waveformusing gediWFSimulator command from rGEDI package(Hancock et al.2019)

The results confirm that the accuracy of FW estimates of forest growth depends on the complexity of the horizontal and vertical structural diversity of the vegetation.From the earned experience and extensive visual inspection attending to the spatial layout of GEDI shots within FMS polygons, the presented approach is robust for homogeneous ecosystems on which treatment prescriptions have not impacted horizontal tree spacing. In this sense, the annual rates of lidar estimated growth from more homogeneous P. radiata plantations were close to(but slightly higher than) the observed field growth estimate for the 2009 to 2018 period of NFI plots. However,due to the P. pinaster forest structural diversity of Galicia (even-aged and uneven-aged forest), annual forest growth estimates were more variable than in homogeneous P. radiata plantations. The observed noise on sparse, naturalized, and more multi-layered P. pinaster forests could be explained by the estimation of bareground elevation from the GEDI full-waveform laser and the associated variability when computing ALS statistics in the presence of vegetation gaps (Hyde et al. 2005;Pascual et al.2019).The growth dynamics captured from GEDI to ALS surveying were slightly higher (median values) to the NFI plots used as the validation data, although the ground-based evidence reported in Galicia showed faster growth rates (Diéguez-Aranda et al. 2012).The impact of laser precision and co-registration requires the determination of some upper bounds for the height increments and, in this research, we benefit from the support of extensive studies on growth and yield modeling in the region to identify that upper bound to filter “growth” from “outliers”.

Fig.11 Showcases of non-positive height dynamics in P.radiata polygons with GEDI footprints located in the proximity of forest edges,discontinuities and low canopy cover

Fig.12 Showcases of areas classified as uncertainty when accounting for height dynamics between ALS and GEDI data and the spatial configuration of ground footprints in the transition areas between non-forest and in heterogenous patches in terms of canopy cover

Co-registration of ALS and GEDI data

The technical specifications of the GEDI laser system and official recommendations suggest considering between 8 and 10 m as the precision accuracy to determine the easting and northing the GEDI ground beams(Luthcke et al. 2020) (GEDI manual). The high degree of fragmentation in the landscape and the small size of forest assets conflicted with the nominal co-registration goodness of GEDI shots. Validation studies are being encouraged from the development team to improve the co-registration performance of the GEDI laser mission.The 30-m buffering distance towards the edge FMS polygons used in this study could have been higher but the impact on the computations might have been not as relevant as expected due to the heterogeneity of FMS polygons. However, by restricting the analyses to GEDI shots of >70% of forest cover and using the 99th height percentile to calculate height growth differences we certainly controlled the impact of edge effects. On the other hand, a slope filter was not used in our study. We did not observe substantial improvements in our results at Lugo province since the number of GEDI shots with high slopes values was low in our dataset. However, this type of filters could be restrictive in other forest areas topographically more complex with steep slopes. The reduction of 92.3% on the valid GEDI shots was important but necessary to isolate valid GEDI shots from the factors of uncertainty regarding positioning and ALS referencing. The constraint of data use is not relevant as the mission is collecting data on a weekly basis so the almost 300,000 GEDI shots initially evaluated could turn into 2 million by the end of 2020 while reducing the between-track distance. The reason for not using GEDI L3 products (i.e., gridded statistics on heights, forest cover or ecological indexes) was the high between-track distance for the small-scale and dispersed forest properties in the region. We evaluated the capability of GEDI to detect changes, but the accuracy of the measurements could be tested if ALS and GEDI were collected simultaneously. In future studies, Extremadura region (Spain)will be tested including the GEDI Level 3 as: i) the scale of forest management units is greater, ii) growth rates are lower and iii) forest management is not as intense as in commercial forestry in Galicia. The assimilation of GEDI data would benefit from parametrization. We hypothesize with the utility of allowing the user to rescale the estimates for upper RH percentiles based on accurate ALS-based ground information. The presented results and ideas can help GEDI scientists to improve and adjust the data processing mechanism.

Conclusions

The use of GEDI laser data provides valuable insights for forest industry operations especially when accounting for rapid spatio-temporal changes in forest structure high-quality ALS laser scanning data as baseline. Harvest and fast-growing forest are detected as well as the variability of forest attributes. The ALS-supported filtering methods described in this paper bring relevant insights to implement GEDI mission data, which is still operative. The correction of edges and the behavior of the laser signal in sparse forest conditions are important factor to consider and mitigate. Biome-specific algorithms to calibrate GEDI data science products with dense ALS data can upgrade the quality of the ground DEM and easy the assimilation of GEDI statistics under the novel and multi-temporal conceptualization of forest productivity.

Acknowledgements

We gratefully acknowledge Vicente Sandoval and Elena Robla from National Forest Inventory Department for supplying the inventory databases NFI data and the latest version of FMS. We also thanks to i) Dr. Carlos Alberto Silva(University of Maryland) for the relevant remarks, suggestions and his advice using rGEDI package and ii) Juan Gabriel Álvarez González (Department of Agroforestry Engineering, University of Santiago de Compostela) for his support with NFI database.

Authors’contributions

The two authors were involved in all research stages from experiment conceptualization to scientific reporting. The author(s) read and approved the final manuscript.

Funding

This work was partially supported by‘National Programme for the Promotion of Talent and Its Employability’ of the Ministry of Economy, Industry, and Competitiveness (Torres-Quevedo program) via postdoctoral PTQ2018-010043 to Dr. Juan Guerra Hernández.The authors also thank to Forest Research Centre, a research unit funded by Fundação para a Ciência e a Tecnologia I.P. (FCT), Portugal (UIDB/00239/2020). Arizona State University,USDA Forest Service and the Asner Lab supported Dr. Adrián Pascual in the final stages of the research.

Availability of data and materials

The research data and materials are available from the corresponding author on reasonable request.

Ethics approval and consent to participate

Not applicable.

Consent for publication

We express our consent to publish the paper in the journal.

Competing interests

The authors declare that they have no competing interests.

Author details

13edata,Centro de iniciativas empresariais, Fundación CEL. O Palomar s/n,27004 Lugo, Spain.2Forest Research Centre, School of Agriculture, University of Lisbon, Tapada da Ajuda, 1349-017 Lisbon, Portugal.3Center for Global Discovery and Conservation Science, Arizona State University, Hilo, HA 96720,USA.

Received: 2 July 2020 Accepted: 1 February 2021