Modeling structural traits of Aleppo pine (Pinus halepensis Mill.) forests with low-density LiDAR

ABSTRACT

Mediterranean forests of Aleppo pine (Pinus halepensis Mill.) have a crucial role in climate change, as they are extremely adaptive and provide valuable timber or carbon stocks. However, greater detail quantifying those attributes is needed: although National Forest Inventories are acceptable, continuous cover maps are normally lacking. Here, we use the public Spanish low-density LiDAR flights to model above-ground biomass, volume, tree density, basal area and dominant height of naturally regenerated Mediterranean Aleppo pine forests, comparing individual-tree detection and area-based approach. We found R² and RRMSE among 0.51–0.66 and 40–34% for above-ground biomass, 0.54–0.70 and 34–28% for volume, 0.23–0.45 and 33–28% for tree density, 0.48–0.62 and 32–27% for basal area, and 0.70–0.69 and 11–11% for dominant height. In all cases but dominant height, the area-based approach outperformed the individual-tree detection. Neither time difference between LiDAR flight and ground measurement or past land use affected the area-based approach models, yet the latter had a strong effect on observed productivity. The different definitions of dominant height were equivalent and did not influence the dominant height models. We believe these models, and their corresponding maps, will be a great asset for policymakers and different stakeholders for Aleppo pine forests throughout the Mediterranean basin.

KEYWORDS:

• LiDAR
• area-based approach
• individual-tree detection
• land historic use
• mediterranean pinewoods

Introduction

Aleppo pine (Pinus halepensis Mill.) is a Mediterranean forest species that is distributed over 36 million hectares across the Mediterranean basin (Caudullo et al., 2017), covering 6.8 million hectares (Farjon, 2010). This species is one the most representative forest species of the Mediterranean ecosystems (Mauri et al., 2016), being the most important forest species of Northern Africa and of great ecological importance in eastern Spain, covering 2.1 million hectares (MAGRAMA 2012). It is also widespread in southern France, Italy, and the coastal zones of Croatia and Bosnia-Herzegovina. On the eastern side of the Mediterranean, it is present in Turkey, Algeria, Tunisia, and Libya, where it shares distribution with Turkish or Calabrian pine (Pinus brutia Ten.), being able to hybridize with it (Dounavi et al., 2001). Despite having a reasonably good wood quality (Elaieb et al., 2017; Mòdol & Casals, 2012), it is normally deprecated in commercial forestry due to its irregular shape and small size (Praciak, 2013). Nonetheless, it has traditionally been used for mine props, railway sleepers, and telephone poles (Farjon, 2010), and its use is nowadays being boosted in the furniture market, packaging industry in pallet production (Elaieb et al., 2017) and even in newly engineered products such as oriented standard board (OSB) or cross laminated timber (CLT) (Brunet-Navarro et al., 2022). However, the bioenergy industry is still the most expanded use for this species (Lerma-Arce et al., 2021).

Pinus halepensis is a very plastic species with an outstanding adaptation capacity (Rojas-Briales et al., 2023). Drought resistant and extremely thermophilic, this pine species adapts well to poor soils, but also low rainfall and high temperatures, creating very resilient forest stands (Mauri et al., 2016) which is especially crucial under the current climate change scenarios. Additionally, due to the colonization of abandoned agricultural land or scrub-dominated natural sites, both naturally- (Daniels et al., 2018) and human-driven (Gil, 2002), the productive potential of this species has significantly increased. In this context, quantifying both the present and potential productivity and these forests’ stand attributes is a key aspect, regardless of the final aim: either for climate change mitigation (i.e. carbon stock quantification or wildfires risk mapping and management) or for productive management (i.e. volume stock or genetic improvement requirements).

Point-based retrieval of forest structural traits, such as the National Forest Inventories (NFIs), is a perfectly valid approach and a reliable source of information, however, multiple drawbacks need to be taken into consideration: (1) are very costly, physically demanding and time-consuming (Hall et al., 2005); (2) are carried out to provide national- and international-level statistics (Kangas et al., 2018) and therefore provide lower detailed analysis; (3) are normally performed at a multi-year interval (e.g. Spanish NFI every 10 years), and therefore the temporal analysis of the data is restricted to long-term analyses; and, directly related to the former (4) potential lack of consistency between different regional or temporal NFI programs. Another point-based source of information are Forest Management Inventories (FMIs), but these are normally generated based on field visual assessments and aerial orthophoto interpretation (Kangas et al., 2018), resulting in lower data quality. With the arrival of Light Detection and Ranging (LiDAR) technology, forest inventory has been in a state of constant redefinition (Hyyppä et al., 2008). LiDAR is considered the most accurate and expanded remote sensing technology to generate spatially explicit, wall-to-wall three-dimensional stand information applicable at strategic, tactic, and operational scales (Bergseng et al., 2015; Kangas et al., 2018; Tompalski et al., 2016). However, in Spain, the modeling and mapping of forest stand attributes is mostly restricted to more productive species of other regions of the Iberian Peninsula (e.g. Eucalyptus spp. or Pinus pinaster L. var. atlantica), with some exceptions (e.g. CREAF and ICGC, 2021; Domingo, Alonso Ponce, et al., 2019). Nevertheless, to improve climate change mitigation, wildfires prevention and management, or timber industry assessment, new efforts have to be made to extend those analyses to the remaining Iberian Peninsula.

There are two main LiDAR-based approaches to model forest stand attributes: individual-tree detection (ITD) and area-based approach (ABA) (Yu et al., 2010). In the latter method, several statistics calculated from the laser point clouds are obtained and used to generate models of forest stand attributes; whereas in the ITD method, individual trees are segmented from the laser point cloud, and treetops and their individual height values are located and used (either directly or as calibration proxies) to generate forest metrics models (Ke & Quackenbush, 2011). Nonetheless, there are several factors that affect the applicability and accuracy of these methodologies, being pulse density (pulses m⁻²), defined as the number of pulses sent by the sensor per unit area (Silva et al., 2017) and normally equated to the total number of first returns per unit area (Oh et al., 2022), a crucial aspect of the forest metric modeling (Coops et al., 2021). For instance, delineation of individual trees is possible when using small footprint (i.e. <1 m) LiDAR data (Sumnall et al., 2022), yet it has been reported that precise individual tree height estimations require >10 pulses m⁻² (Jakubowski et al., 2013). On the other hand, ABA has been shown to work with lower pulse densities (Guerra-Hernández et al., 2021) and has much lower computational power requirements to generate maps across the landscape (e.g. Domingo, Montealegre, et al., 2019; Fernández-Landa et al., 2018).

Although higher-density LiDAR point clouds rely less upon ground data to generate representations of stand structural traits (Kellner et al., 2019), the traditional workflow of validating remotely sensed models with field data from ground forest inventories is still a state-of-the-art procedure. Therefore, a common aspect for both methodologies is the time difference between the LiDAR data acquisition and the ground forest inventory. When LiDAR data is not harmonized in time with validation ground data, some noise can appear in the models, as the field-measured variable can be in a different stage than the observed one by the remote sensor. While this is something to be considered with high-density point clouds (either LiDAR or structure-from-motion) and fast-growing species, for Mediterranean species such as Pinus halepensis, with height growth rates of 0.25 to 0.4 m year⁻¹ for the most productive sites (Montero et al., 2001), this effect can be negligible (at least for some time), especially with low pulse densities and ABA approaches (Oh et al., 2022). Hence, multiple examples of ABA-based structural traits retrieval exist for different forest conditions (e.g. González-Ferreiro et al., 2012; Green et al., 2019).

Furthermore, the past land use and cover history normally has a non-negligible effect on the species’ growth and potential (Thom et al., 2018). In the Mediterranean basin, and specifically the Spanish Mediterranean forests, the rural exodus and land abandonment happened mostly in the 1950s and 1960s, generating a migration movement from the inland to the coastal areas, producing both a substantial cover increase in abandoned rural areas and the capitalization of the remaining forests (Delgado-Artés et al., 2022; Ribas-Costa & Dopazo, 2020). Nowadays, these forests behave differently, and the impact of the historical land cover in present Mediterranean forest attributes, yet commonly acknowledged, is complex to study due to the difficulty of obtaining good and representative data (Peñuelas et al., 2017). In this context, measuring, modeling, and mapping forest stand attributes (e.g. above-ground biomass, total volume, tree density, basal area, and dominant height – productivity), and analyzing the impact of past land use and cover on those can improve the understanding of global change on today’s forests and their ecosystem services, as well as help different stakeholders in their forest management decision processes.

Specifically, one important aspect of forest management is forest productivity. Site productivity is defined as the production that can be realized at a certain site with a given genotype, a specified management regime (Skovsgaard & Vanclay, 2008), and for a given environment, i.e. topography, soil characteristics, and climate (Subedi & Fox, 2016), and is normally understood as site index (SI), i.e. dominant height at a given age (Zhao et al., 2016). In Europe, dominant height is normally equated to the average height of the 100 thickest trees per hectare, although other definitions also exist (e.g. Weise, 1880 or Marschall, 1992). To assess forest productivity at a landscape scale, either single LiDAR acquisitions and forest age (e.g. Gopalakrishnan et al., 2019) or multi-temporal LiDAR acquisitions (e.g. Guerra-Hernández et al., 2021) are required. Discarded the multi-temporal LiDAR acquisition method due to the low-density point cloud restriction and the low growth rate of Pinus halepensis, having a LiDAR-derived map of dominant height appears as the most plausible option to achieve the landscape productivity assessment goal.

Accordingly, the specific objectives of this study are:

1. To assess the accuracy of freely available low-density LiDAR data to model Pinus halepensis forest attributes such as above-ground biomass or dominant height.

2. To compare two different methodologies, the area-based approach, and the individual-tree detection process.

3. To evaluate the effect of the following variables when modeling the previous forest attributes using the area-based approach methodology at a large scale:

   1. Time difference between the LiDAR data acquisition and the ground measurements.

   2. Pulse density of the LiDAR flights.

   3. Past land use and cover history.

   4. Dominant height definition.

4. To evaluate the effect of past land use and cover history in site productivity.

Methods

Spatial context

The study was located on the 571.6 km² island of Ibiza (Spain), situated in the western part of the Balearic archipelago, in the Mediterranean Sea (Figure 1). This island has been estimated to have 29,913.71 ha of forestland, out of which 23,246.04 ha are forests of Aleppo pine (Pinus halepensis Mill.) (MAGRAMA, 2012).

Figure 1. Locations of the ground data used in this study. Colors represent the number of plots per location, and the box summarizes the plot categories regarding historic land use. The Pinus halepensis Mill. coverage has been extracted from the IV National Forest Inventory Map National Forest Map 1:25,000 (2007–2017) (MAGRAMA 2012).

Data

Forest inventory data

A total of 113 ground plots were established in Aleppo pine forestland between October 2021 and July 2023 (Figure 1). The plots were 10-m radius and all live trees with a diameter at breast height (DBH) equal or higher than 7.5 cm inside the plot were measured. The DBH of each tree was recorded using a diameter tape (Weiss, Germany) and individual total height (HT) was measured using a laser rangefinder (Nikon, Japan). Age was also approximated in the field by coring the three thickest trees per plot using a tree core reader (Hagolf, Sweden). All the plots were categorized regarding historic land use, in (1) agriculture land, (2) scrubland, (3) sparse forestland, and (4) dense forestland by combining information from the field and interpretation of the 1956–1957 national orthophoto (National Plan of Aerial Orthophotography, PNOA) (Table 1). We recorded the plot center with a commercial handheld GPS device with a spatial accuracy of <6 m while co-registering it using the latest PNOA orthophotography available.

Table 1. Forest inventory summary of the measured and modeled variables, classified by historic land use and past cover.

Historic use	Dense forestland (N = 44)				Sparse forestland (N = 26)				Scrubland (N = 10)				Agricultural land (N = 33)
Statistic/variable	Min.	Max.	Mean	Std. Dev.	Min.	Max.	Mean	Std. Dev.	Min.	Max.	Mean	Std. Dev.	Min.	Max.	Mean	Std. Dev.
DBH (cm)	18.2	31.0	24.1	3.5	14.6	37.2	24.1	5.7	12.8	38.7	23.2	7.6	16.0	33.9	25.0	4.3
BA (m^2 ha^−1)	4.18	45.27	22.18	7.82	3.63	38.92	17.99	7.80	0.41	38.30	18.52	10.08	5.82	50.53	18.32	9.46
Hmean (m)	7.8	15.3	10.6	1.5	7.2	15.8	10.7	2.4	6.0	14.3	10.7	2.4	9.4	16.7	12.3	2.0
Hdom (m)	8.6	18.4	12.1	1.9	7.8	19.7	12.8	2.9	6.0	18.6	13.0	3.5	9.5	21.2	14.2	2.6
Age (years)	58	113	82	12	23	88	73	14	37	78	56	13	30	89	46	14
Density (trees ha^−1)	95	955	453	171	127	987	386	216	32	891	363	230	95	955	335	174
SI (m)	8.4	16.7	11.9	2.0	8.8	20.3	13.0	2.7	8.9	19.8	14.9	3.5	12.2	32.3	19.0	4.6
AGB (Mg ha^−1)	13.47	191.57	82.66	35.40	10.40	212.18	71.57	44.67	0.87	236.42	83.96	63.46	18.47	291.69	83.80	54.01
VOL (m^3 ha^−1)	18.32	229.09	104.89	41.42	17.04	219.93	89.31	45.21	2.31	240.63	98.58	62.92	22.63	303.36	99.11	56.32

Once the dataset was finished (N_trees = 1,460; N_plots = 113), we used the individual tree dataset to adjust a general DBH–HT curve using all the tree data. We based our final model on our own data exploration and on the work of Cabanillas-Saldaña (2010). Additionally, we also calculated the basal area per tree, the individual tree biomass, and volume by using Ruiz-Peinado et al. (2011) allometric equations (Equations 1(7) $\mathit{SI}=H_{\mathit{dom}}\cdot {\left(\frac{1-e^{-b\cdot A_0}}{1-e^{-\mathit{b}\cdot \mathit{A}}}\right)}^{\frac{1}{c}}$(7) –5(11) $\mathit{d}=4.7217+1.2407\cdot \mathit{h}+0.0392\cdot h^2$(11) , Table 2) and the IV Spanish National Forest Inventory (MAGRAMA, 2012) volumetric equations (Equation 6, Table 2).

Table 2. Pinus halepensis Mill. allometric and volumetric equations.

Item	Units	Expression	Equation	Reference
Stem	kg	$W_s=0.0139\cdot d^2\text{cdoth}$	(1(7) $\mathit{SI}=H_{\mathit{dom}}\cdot {\left(\frac{1-e^{-b\cdot A_0}}{1-e^{-\mathit{b}\cdot \mathit{A}}}\right)}^{\frac{1}{c}}$(7) )	Ruiz-Peinado et al. (2011)
Thick branches	kg	If d ≤27.5 cm then Z = 0; If d >27.5 cm then Z = 1 $W_{\mathit{b}7}=\left[3.926\cdot \left(\mathit{d}-27.5\right)\right]\text{cdotZ}$	(2(7) $\mathit{SI}=H_{\mathit{dom}}\cdot {\left(\frac{1-e^{-b\cdot A_0}}{1-e^{-\mathit{b}\cdot \mathit{A}}}\right)}^{\frac{1}{c}}$(7) )	Ruiz-Peinado et al. (2011)
Medium branches	kg	$W_{\mathit{b}2-7}=4.257+0.00506\cdot d^2\text{cdoth}-0.0722\text{cdotdcdoth}$	(3(9) $\mathit{RMSE}=\sqrt{\frac{\sum _{\mathit{i}=1}^n{\left(y_i-{\hat{y}}_{\mathit{i}}\right)}^2}{\mathit{n}}}$(9) )	Ruiz-Peinado et al. (2011)
Thin branches + needles	kg	$W_{\mathit{b}2+\mathit{n}}=6.197+0.00932\cdot d^2\text{cdoth}-0.0686\text{cdotdcdoth}$	(4(10) $\mathit{RRMSE}=\frac{\mathit{RMSE}}{\bar{y}}$(10) )	Ruiz-Peinado et al. (2011)
Above-ground Biomass	kg	$\mathit{AGB}=W_s+W_{\mathit{b}2-7}+W_{\mathit{b}2+\mathit{n}}$	(5(11) $\mathit{d}=4.7217+1.2407\cdot \mathit{h}+0.0392\cdot h^2$(11) )	Ruiz-Peinado et al. (2011)
Total Volume	dm^3	$\mathit{VOL}=47.45+0.0000257\cdot {\left(\mathit{d}\cdot 10\right)}^2\text{cdoth}$	(6(11) $\mathit{d}=4.7217+1.2407\cdot \mathit{h}+0.0392\cdot h^2$(11) )	MAGRAMA (2012)

Biomass yields in kilograms, volume in cubic decimeters. Diameter at breast height (d) in cm and total height (h) in meters. Equations 1(7) $\mathit{SI}=H_{\mathit{dom}}\cdot {\left(\frac{1-e^{-b\cdot A_0}}{1-e^{-\mathit{b}\cdot \mathit{A}}}\right)}^{\frac{1}{c}}$(7) to 5(11) $\mathit{d}=4.7217+1.2407\cdot \mathit{h}+0.0392\cdot h^2$(11) extracted from Ruiz-Peinado et al. (2011), Equation 6(11) $\mathit{d}=4.7217+1.2407\cdot \mathit{h}+0.0392\cdot h^2$(11) extracted from MAGRAMA (2012). For allometric equations, thick branches are of more than 7 cm in diameter, medium branches are between 7 and 2 cm in diameter, and thin branches are below 2 cm in diameter.

After that, we collapsed all the single-tree variables per plot, obtaining the target variables:

• Dominant height: Hdom. Five different definitions of dominant height were tested:

  • domHtws: The average height of the 20% thickest trees per plot (Weise, 1880).

  • domHtassb: The average height of the 100 thickest trees per hectare (adaptation from Assman, 1970), in this case, 3 trees per plot. However, to account for bias in the dominant height definition when plots are smaller than 5,000 m² (Bouchon, 1994; Pardé & Rennolls, 1978), we calculated H_dom as the average height of the 2 thickest trees per each 314.15 m² plot.

  • domHtmsc: The average height of the 20% tallest trees per plot (Marschal, 1992).

  • domHt3tt: The average height of the 3 tallest trees per plot.

  • domHtpct85: The 85th percentile of the individual tree heights in the plots.

• Basal area: BA (m² ha⁻¹). As the sum of the individual tree basal areas divided by the plot area.

• Tree density: TPH (trees ha⁻¹). As the total number of stems divided by the plot area.

• Above-ground Biomass: AGB (Mg ha⁻¹). As the sum of the individual tree above-ground biomasses divided by the plot area.

• Total volume with bark: VOL (m³ ha⁻¹). As the sum of the individual tree volumes divided by the plot area.

• Age: A (years). As the average age of the three thickest trees in the plot.

We used the SI models for Pinus halepensis by Montero et al. (2001) to generate the SI equation as a function of observed dominant height and age, transforming the anamorphic curves by Montero et al. (2001) to an Algebraic Difference Approach (ADA) model as in Trim et al. (2020). The final model is (Equation 7(7) $\mathit{SI}=H_{\mathit{dom}}\cdot {\left(\frac{1-e^{-b\cdot A_0}}{1-e^{-\mathit{b}\cdot \mathit{A}}}\right)}^{\frac{1}{c}}$(7) ).

(7) $\mathit{SI}=H_{\mathit{dom}}\cdot {\left(\frac{1-e^{-b\cdot A_0}}{1-e^{-\mathit{b}\cdot \mathit{A}}}\right)}^{\frac{1}{c}}$(7)

In which $\mathit{SI}$ is the site index value (in meters) at the reference or base age $A_0$, in this case of 80 years (Montero et al., 2001), $H_{\mathit{dom}}$ is the observed dominant height per plot (in meters), $\mathit{A}$ is the measured age per plot (in years), and $\mathit{b}$ and $\mathit{c}$ are the model coefficients (0.0203954 and 1.046295, respectively).

LiDAR data

LiDAR data was obtained from the second nationwide coverage acquisition (PNOA, 2016–2019). Specifically, the data that was used for this study was acquired between November and December of 2019. All PNOA data is publicly available through the web of PNOA at CNIG (Centro Nacional de Información Geográfica) online servers, from which LiDAR swaths can be downloaded in 2-km side square tiles in a compressed LAS file format (LAZ) (https://centrodedescargas.cnig.es/CentroDescargas/index.jsp#). The scanning sensor was a LEICA ALS80 (Leica, Germany), with an XY accuracy of 0.3 m, a Z accuracy of 0.15 m, and a nominal laser pulse density of 1 pulse m⁻². Data was delivered in the ETRS89 UTM Zone 31N coordinate system (EPSG:25831).

Despite following two different approaches in the LiDAR processing, the pre-processing steps were common, and consisted in the following workflow: (1) tile download and merge to obtain the complete coverage for a given location using LASMergePro function from LAStools (Isenburg, 2020); (2) clipping of the LAS file to the location extent using LASClip function from LAStools (Isenburg, 2020); (2) statistical filtering implementation to remove noise based on the distance to neighbor points (Rusu et al., 2008), algorithm implemented using the Filters.outlier function of the Point Data Abstraction Library – PDAL (PDAL contributors, 2020) performed with Python 3.8.13 (https://www.python.org/), and (3) a ground reclassification algorithm implementation using the LASGround function from LAStools (Isenburg, 2020) to separate ground returns from vegetation. Finally, in step (4) we conducted a normalization of the point cloud to obtain above-ground heights using the PDAL (PDAL contributors, 2020) filters.hag_nn function. This function calculates an average ground height weighted by the distance of each ground point from the non-ground points. The height above-ground (HAG) is the difference between the z value of the non-ground point and the interpolated ground height.

LiDAR methodologies

Area-based approach

The Area-Based Approach (ABA) method consisted in clipping the LiDAR data around each of the ground plots. Ground plots were matched with the LiDAR data by geographical intersection, and a 20-m side square (pixel) was buffered from the plot center and then used to clip the LiDAR point cloud, obtaining the LiDAR data that corresponded to the plot. After this, statistical descriptors were obtained. We exclusively used the first return per pulse, as they are more related to ground-based height variables, and the models generated with only first returns have been proved to be more transferable (Bater et al., 2011; Cao et al., 2014; Nord-Larsen & Riis-Nielsen, 2010). We also excluded the returns with a lower HAG value than 2 m to exclude both ground and understory returns (Kankare et al., 2013; Torre-Tojal et al., 2022). Once excluded, the following statistics were computed for each voxel: mean, median, standard deviation, variance, coefficient of variation, maximum, minimum, skewness, kurtosis, percentiles (5th, 10th, 25th, 85th, 90th, 85th, 90th, 95th, 99th), percentage of returns above a certain height threshold (1 m, 2 m, 5 m, 10 m, 15 m, 20 m), and vertical canopy cover at both 2 m and at the average crown height, calculated as FCI – First Echo Cover Index (Equation 8(8) $\mathit{FCI}=\frac{\sum \mathit{Singl}{\mathit{e}}_{\mathit{canopy}}-\sum \mathit{Firs}{t}_{\mathit{canopy}}}{\sum \mathit{Singl}{e}_{\mathit{all}}-\sum \mathit{Firs}{t}_{\mathit{all}}}$(8) ) from Korhonen et al. (2011), which compares the proportion of first and single canopy echoes. Table 3 shows all the computed variables.

(8) $\mathit{FCI}=\frac{\sum \mathit{Singl}{\mathit{e}}_{\mathit{canopy}}-\sum \mathit{Firs}{t}_{\mathit{canopy}}}{\sum \mathit{Singl}{e}_{\mathit{all}}-\sum \mathit{Firs}{t}_{\mathit{all}}}$(8)

Table 3. Names of the variables computed for the ABA methodology.

avgHTagrn1	Mean	pct85HTagrn1	85th percentile
medianHTagrn1	Median	pct90HTagrn1	90th percentile
stdHTagrn1	Standard deviation	pct95HTagrn1	95th percentile
varHTagrn1	Variance	pct99HTagrn1	99th percentile
vcHTagrn1	Variation coefficient	pte1rn1	Percentage of first returns above 1 m
maxHTagrn1	Maximum	pte2rn1	Percentage of first returns above 2 m
minHTagrn1	Minimum	pte5rn1	Percentage of first returns above 5 m
skHTagrn1	Skewness	pte10rn1	Percentage of first returns above 10 m
ktHTagrn1	Kurtosis	pte15rn1	Percentage of first returns above 15 m
pct5HTagrn1	5th percentile	pte20rn1	Percentage of first returns above 20 m
pct10HTagrn1	10th percentile	fci2	First echo cover index at 2 m
pct25HTagrn1	25th percentile	fciavgHt	First echo cover index at average crown height
pct75HTagrn1	75th percentile	•	•

All the variables are computed for the first returns above ground (i.e. height above ground ≥ 2).

Individual-tree detection

For the individual-tree detection (ITD) method, we followed a low-density canopy height model (CHM) approach based on the adaptation of previous methodologies (Hui et al., 2022; Yu et al., 2010) that consisted in the following steps:

1. Rasterization of the normalized point cloud with a 1-m pixel size, using the maximum HAG value within each cell.

2. Morphological filtering to remove “salt-and-pepper noise” by applying a median filter to the CHM.

3. Gaussian filtering to remove high-frequency noise and to emphasize local maxima and crown contours.

4. Local maxima filter to identify the treetops and their individual height.

We applied the previous workflow to all the ground plot locations, geographically intersecting the ground plots with the LiDAR point clouds and buffering their centroids into a 20-m side square pixel. All the trees that were within or overlaid the pixel were clipped and then several variables were calculated at a per-plot basis (Table 4).

1. Dominant height: we obtained statistical descriptors of the individual heights, i.e. mean, median, percentiles (75th, 80th, 85th, 90th, 95th, 99th), and maximum.

2. Tree density: using the total number of treetops per plot.

3. Basal area: we used our DBH-HT curve to estimate the individual diameter of each tree, from which we estimated a per-tree value of basal area. Then we obtained the total basal area per plot.

4. Above-ground biomass: we also used our DBH-HT curve to estimate the individual diameter of each tree and used Equations 1(7) $\mathit{SI}=H_{\mathit{dom}}\cdot {\left(\frac{1-e^{-b\cdot A_0}}{1-e^{-\mathit{b}\cdot \mathit{A}}}\right)}^{\frac{1}{c}}$(7) –5(11) $\mathit{d}=4.7217+1.2407\cdot \mathit{h}+0.0392\cdot h^2$(11) to compute the aerial biomass per tree. Then all the individual tree AGBs in the plot were summed.

5. Volume: again, we used our DBH-HT to estimate the individual diameter of each tree and used Equation 6(11) $\mathit{d}=4.7217+1.2407\cdot \mathit{h}+0.0392\cdot h^2$(11) to compute the volume per tree. Then all the individual tree VOLs in the plot were summed.

Table 4. Names of the variables computed for the ITD methodology.

proxyTPH	Number of detected treetops	Htpct90	90th percentile of individual heights
avgHt	Mean height of the trees	Htpct95	95th percentile of individual heights
medianHt	Median height of the trees	Htpct99	99th percentile of individual heights
Htpct75	75th percentile of individual heights	Htmax	Maximum individual height
Htpct80	80th percentile of individual heights	proxyBA	proxy of BA
Htpct85	85th percentile of individual heights	proxyAGB	proxy of AGB
–	–	proxyVOL	proxy of VOL

Statistical analyses

As a first approximation to the data, we performed a set of non-parametric models built with random forest (RF) using the “ranger” R package (Wright and Ziegler, 2017). We then performed a visual inspection of the data and used ordinary least squares (OLS) regressions to model the target variables using the most explanatory covariates. We performed an internal validation approach as in Steyerberg et al. (2001) via a 500-repetition bootstrapping method, using the “rms” R package (Harrell, 2023) all in R 4.2.3. statistical software (https://www.r-project.org/). In ABA, to compare among different models built with different variables, the predictive capability was assessed via the coefficient of determination (R²), the root mean square error, the relative root mean square error (RMSE and RRMSE, Equations 9(9) $\mathit{RMSE}=\sqrt{\frac{\sum _{\mathit{i}=1}^n{\left(y_i-{\hat{y}}_{\mathit{i}}\right)}^2}{\mathit{n}}}$(9) and 10(10) $\mathit{RRMSE}=\frac{\mathit{RMSE}}{\bar{y}}$(10) ), and the Akaike Information Criterion (AIC).

(9) $\mathit{RMSE}=\sqrt{\frac{\sum _{\mathit{i}=1}^n{\left(y_i-{\hat{y}}_{\mathit{i}}\right)}^2}{\mathit{n}}}$(9)

(10) $\mathit{RRMSE}=\frac{\mathit{RMSE}}{\bar{y}}$(10)

Where $y_i$ is the observed value, ${\hat{y}}_{\mathit{i}}$ is the predicted value for each of the plots $\mathit{i}$, and $\bar{y}$ is the mean of all observed values. In ITD, the dominant height was compared to different percentiles and the best fit was selected. Then, the LiDAR-derived variables were used as proxies and input covariates to calibrate the OLS regression models. Again, the predictive capability was assessed via R², RMSE, and RRMSE. Given that no substantial increase in predictive performance of the models was found when comparing OLS to RF, we chose to continue the analysis with OLS for interpretability and simplification.

Once the models were built and validated, we assessed the impact of several variables on their performance through a workflow (Figure 2). First, we studied the models’ residuals against the time difference between the LiDAR flights and the ground measurements and the average pulse density of each plot. Both regression analysis and Tukey’s Honestly Significant Difference (HSD) test were performed between the residuals and each of the studied variables. Additionally, we analyzed the interaction effect between the most explanatory variable for each model and the past land use and cover. To do so, we ran OLS between the most explanatory variables and analyzed the coefficients for each combination. Finally, we compared the effect of past land use and history on the observed dominant height and SI via a Tukey’s HSD test.

Figure 2. General overview of the methodologies applied in this study. In blue, are the method steps that involve ground data, in green are the steps that involve LiDAR data and in orange are the steps that involve both.

Results

The DBH–HT relation

The best function to approximate diameter at breast height from total height was Equation 11.(11) $\mathit{d}=4.7217+1.2407\cdot \mathit{h}+0.0392\cdot h^2$(11)

(11) $\mathit{d}=4.7217+1.2407\cdot \mathit{h}+0.0392\cdot h^2$(11)

Where $\mathit{d}$ is diameter at breast height in centimeters and $\mathit{h}$ is total height in meters. This equation had an R² of 0.371, an RMSE of 7.3 cm, and a RRMSE of 31.61% (outputs from a 10-fold cross-validation procedure).

Model building and methodology comparison

Area-based approach

After analyzing the effect of the covariates in each of the models, we found that the best combinations were the ones presented in Table 5, ABA section. Dominant height was strongly correlated to the 95th percentile of the height distribution of the returns, whereas basal area was related to the 25th percentile. Adding more variables to these models resulted in increased complexity, risk of multicollinearity, and reduction of the model goodness of fit (higher AIC), and therefore were not considered. The most explanatory covariates for both above-ground biomass and volume were the average height return and the percentage of first returns above 5 m. Despite the latter having a non-negligible effect, it had a lower predictive power than the other covariate in the model. Tree density was found to be correlated to the first echo cover index at average return height and to the variation coefficient of the heights, and both had a similar importance in the model.

Table 5. Expressions to predict forest attributes for both LiDAR-based methodologies.

Method.	Forest Attribute	Expression	Núm.
Area-Based Approach	AGB (Mg ha^−1)	$\mathit{AGB}=65.2339-21.5444\cdot \mathit{avgHTagrn}1+1.9932\cdot \mathit{avgHTagrn}{1}^2+0.6983\cdot \mathit{pte}5\mathit{rn}1$	11
	VOL (m^3 ha^−1)	$\mathit{VOL}=63.4626-21.2720\cdot \mathit{avgHTagrn}1+1.8986\cdot \mathit{avgHTagrn}{1}^2+1.1192\cdot \mathit{pte}5\mathit{rn}1$	12
	TPH (trees ha^−1)	$\mathit{TPH}=299.9573-5.6318\cdot \mathit{vcHTagrn}1+603.3975\cdot \mathit{fciavgHT}$	13
	BA (m^2 ha^−1)	$\mathit{BA}=-5.2010+3.3765\cdot \mathit{pct}25\mathit{HTagrn}1$	14
	Hdom (m)	$H_{\mathit{dom}}=0.5694+0.9952\cdot \mathit{pct}95\mathit{HTagrn}1$	15
Individual-Tree Detection	AGB (Mg ha^−1)	$\mathit{AGB}=34.3894+1.6520\cdot \mathit{proxyAGB}$	16
	VOL (m^3 ha^−1)	$\mathit{VOL}=39.1917+1.6896\cdot \mathit{proxyVOL}$	17
	TPH (trees ha^−1)	$\mathit{TPH}=137.5623+1.6524\cdot \mathit{proxyTPH}$	18
	BA (m^2 ha^−1)	$\mathit{BA}=8.7911+1.5835\cdot \mathit{proxyBA}$	19
	Hdom (m)	$H_{\mathit{dom}}=1.1977+0.9077\cdot \mathit{HTpct}95$	20

Tree density model has been calculated for N = 106 plots, when measured TPH ≥ 750; RRMSE of the dominant height model has been obtained using the most extended definition of dominant height (domHtassb).

The best fit was for the volume model, with a R² of 0.691 and a RRMSE of 27.66%. After volume, the second highest R² was for the dominant height model, which had a RRMSE of 11.41%. The model for basal area had the lowest correlation, with an R² of 0.623, yet it had lower RRMSE, 27.05%. Above-ground biomass had the highest RRMSE, 33.86%, despite having an R² of 0.655. In the case of the tree density model, data with more than 750 trees ha⁻¹ (N = 7, 6.2% of the total data) had to be excluded from further analysis as the behavior of the explanatory covariates did not explain those outliers, which had a mean tree density of 887 trees ha⁻¹ (i.e. almost 2.5 times more than the average of the remaining 94% of the data). Regardless of the former adjustment, the model achieved the lowest fit, with an R² of 0.449, but had a similar error to the other models, with an RRMSE of 27.89%. Table 6 summarizes all the results and Figure 3 shows two examples of the maps generated.

Figure 3. Maps of basal area (1), dominant height (2), volume (3) and above-ground biomass (4) for a forested area of the southwestern part of Ibiza Island. The maps have been generated using the models shown in Table 5 (ABA) and the resolution is 20 m.

Table 6. Goodness of fit in terms of R², RMSE, and RRMSE for all the models shown in Table 5 for both LiDAR-based methodologies.

Methodology	Forest Attribute	R^2	RMSE (m)	RRMSE (%)
Area-Based Approach	AGB (Mg ha^−1)	0.655	27.28	33.86
	VOL (m^3 ha^−1)	0.691	27.40	27.66
	TPH (trees ha^−1)	0.449	110	27.87
	BA (m^2 ha^−1)	0.623	5.35	27.05
	Hdom (m)	0.685	1.48	11.41
Individual-Tree Detection	AGB (Mg ha^−1)	0.508	32.55	40.40
	VOL (m^3 ha^−1)	0.539	33.40	33.71
	TPH (trees ha^−1)	0.226	130	33.00
	BA (m^2 ha^−1)	0.479	6.28	31.76
	Hdom (m)	0.704	1.44	11.10

Values are the result of a 500-repetition bootstrapping validation procedure. Tree density model has been calculated for N = 106 plots, when measured TPH ≥ 750; RRMSE of the dominant height model has been obtained using the most extended definition of dominant height (domHtassb).

Individual-tree detection

The main finding of the ITD method was a strong underprediction of the individually located treetops against the real trees in the plot, finding an average underprediction of 255 trees ha⁻¹, with a standard deviation of 174 trees ha⁻¹. Accordingly, the proxies of all the target variables were in all cases of lower value than the observed ones, and therefore rather than using them as direct estimates, we used them as covariates to calibrate the final models. We found almost all the ITD-based models worse in predictive capability than the ABA ones. Except the dominant height model, which had the best fit of all models with a R² of 0.704 and a RRMSE of 11.10%, the rest of the models had a R² under 0.6, and a RRMSE over 30%. In the case of above-ground biomass, the RRMSE of the model was 40.40%.

However, when analyzing the mismatch between the observed tree density and the treetops identified in the ITD method, we found that when limiting the number of plots to a higher average diameter at breast height (plot average DHB >25 cm, N = 43), the correlation between the two variables increased significantly (R² = 0.422 and p-value <0.0001).

Assessing the variability of the model’s predictions

Given the accuracy decrease in the ITD methodology, the substantial increase in computational power and our aim of mapping landscape at a large scale, we chose the ABA models to assess their applicability and to map the forest structure attributes.

The first analysis was performed to check the effect of time difference between ground and LiDAR data acquisitions (Figure 4, left side). This time difference ranged from 22 months to 47 months, but we did not find evidence of effect on the residuals. None of the models showed any effect of time lag on the model residual distribution, both in the Tukey’s HSD mean comparison nor in the regression analysis (p-value >0.311 for all models). Similarly, there was no effect of pulse density on any of the models (Figure 4, right side), except for the dominant height model groups (1.5, 2], flagged with moderate evidence of being different than (1.0, 1.5] (p-value = 0.025) and (2.0, 2.5] (p-value = 0.044). However, the regression analysis showed that overall, there was null effect of pulse density on the model’s performance (p-value >0.369 for all models).

Figure 4. Performance of the models of above-ground biomass (1), total volume (2), trees per hectare (3), basal area (4), and dominant height (5) regarding time lag (left, in blue), and pulse density (right, in green). Residuals are calculated as the observed value minus the predicted value. The whiskers represent the 5th and 95th percentiles of the distribution. The white numbers correspond to the median values, and the black numbers above each box correspond to the number of observations in each group. The tree density model has been calculated for N = 106 plots, when measured TPH ≥ 750; the dominant height represents the average height of the 100 thickest trees per hectare (domHtassb). The fitted resolution is 20 m.

Land history and dominant height definition effect

We also assessed if there was an effect of past land use and cover history on the performance of the models. We performed the analysis with the variables shown in Table 5 but as a linear input, to avoid overprediction and misbehavior of the individual models, due to low sample size for each level.

The relation between the average return height and AGB (Figure 5) was not affected by the past land use (p-value >0.474), except for dense forestland. Sites with past land use of dense forestland had on average 66 Mg ha⁻¹ more than other past land uses (p-value <0.052), yet the model slope was lower: an increase in 1 m of the average return height translated to 7 Mg ha⁻¹ less of AGB than the same increase for past agricultural land (p-value = 0.034), what overall implied a negligible interaction effect. Furthermore, we did not find any substantial improvement in the models’ goodness of fit when including the historic land use. The effect was similar in the VOL models (Figure 5), but the interaction evidence was weaker (p-value >0.093). Regarding TPH, very weak evidence of history effect on the model was found for historic scrubland cover (p-value >0.108), but probably due to low sample size (N = 10). For BA, again no effect was found (p-value >0.368).

Figure 5. Behavior of the most explanatory covariate of above-ground biomass (AGB, 1) and volume (VOL, 2) (avgHtagrn1) by land cover historic use. Regression is fitted using a linear model to prevent overfitting due to limited sample size of each level. The fitted resolution is 20 m.

In terms of dominant height (using the average height of the 100 thickest trees per hectare or domHTassb), the analysis showed moderate evidence that scrubland (p-value = 0.039) had a lower intercept value than all other land histories, i.e. the slope of the dominant height models for historic scrublands was steeper than for historic dense forestlands (Figure 6). This effect shades when building the regression with all the levels of historic land use at the same time (see slope coefficient in Equation 15).

Figure 6. Behavior of the most explanatory percentile of height (pct95HTagrn1) for the five different dominant height definitions tested, and by land cover historic use. R², RMSE, and RRMSE correspond to all data. The fitted resolution is 20 m.

Furthermore, the different definitions of dominant height perform equally (Figure 6): regardless of which one was used, the performance of the LiDAR-derived 95th percentile was found to be the best fit, and its behavior was equivalent among the five definitions of dominant height tested. However, despite finding no significant differences among them, it is true that the dominant height definitions based on tree height rules alone, such as the average height of the three tallest trees per plot (domHT3tt) or the average height of the 20% tallest trees per plot (domHTmsc) normally outperformed the ones based on stem thickness rules, such as the average height of the 100 thickest trees per hectare (domHTassb) or the average height of the 20% thickest trees per plot (domHTws). Looking at the combined effect of dominant height definition and historic land use, the lowest RRMSE is always achieved with domHT3tt: 7.8% for dense forestland; 11.1% for sparse forestland; 8.7% for scrubland and 9.1% for agricultural land. Nevertheless, when comparing all the definitions via a Tukey’s HSD test, all were considered equal (p-value >0.305).

Finally, when looking at the effect of historic land use and cover on forest productivity (in terms of SI), we observed strong evidence (p-value <0.01) that past agricultural land was the most productive, with 4 m more site index than historic scrublands, 6 m more than sparse forest, and over 7 m more than dense forestland. Accordingly, we found the ranking of historic land use effect on productivity to be agricultural land, scrubland, sparse forestland and dense forestland (Table 7).

Table 7. Tukey’s HSD test results (alpha 0.05) for the analysis of past use and land history effect on site index in meters for a base age of 80 years.

Historic land use and cover	Letter report of significance			Site index base age 80 years mean (m)
Agricultural land	A			19.0
Scrubland		B		14.9
Sparse forestland		B	C	13.0
Dense forestland			C	11.9

Levels not connected by the same letter are significantly different.

Discussion

Methodology comparison

In this study, we compared two LiDAR-based methodologies to model forest attributes and found the area-based approach (ABA) to outperform the individual-tree detection (ITD) in terms of predictive capability. However, in all the forest attributes, models’ accuracies were between 11% (for dominant height both in ABA and ITD) and (40% for above-ground biomass in ITD). These results are in line with previous studies where for both methodologies RRMSEs have been reported between 6% and 67% (e.g. Adhikari et al., 2023; Fekety et al., 2018; Kankare et al., 2013) for different stand attributes. Yet more complex machine learning algorithms can lower the relative errors (e.g. Torre-Tojal et al., 2022), the risk of losing extrapolation power requires a trade-off (Hennigar et al., 2016; Jeong et al., 2016). Moreover, notwithstanding that ground-based forest inventories may offer enhanced precision, as less than 15% relative error should be expected (Silva et al., 2017), remote sensing techniques facilitate the assessment of significantly larger geographical extents while reducing the surveying effort. It is in such cases where traditional inventory methods require an excessively high surveying intensity that remote sensing-based approaches are specially needed, and it is in that scenario that ABA outperforms ITD.

ABA models performed better than ITD mainly because the canopy height model-based individual location process was not able to detect all the trees inside the plot. This omission error was due to: (1) the low pulse density of the point cloud and (2) the physical structure of the Pinus halepensis forests. Despite the former is a commonly described limitation of publicly available LiDAR data (Sánchez, 2018), efforts are being made to overcome it (e.g. 3rd PNOA flight, 2019–2025, with 5 pulses m⁻²). However, a slight increase in pulse density will still not easily allow to perform ITD to Mediterranean Pinus halepensis forests. Low-density (<10 pulses m⁻²) ITD methods normally rely on raster-based approaches, which have the disadvantage of losing the 3D information of the point cloud, causing errors from extraction, interpolation, and smoothing procedures (Zhen et al., 2016). More complex ITD approaches normally require higher pulse density for multi-scale estimation methods for crown width (Hui et al., 2022). In addition, Pinus halepensis trees have very globularly-shaped, heterogeneous crowns and tend to cluster, often occluding each other. This hinders the location of the real tree-tops when applying a local maxima filter. Additionally, these forests have a multilayered structure, with a thick and complex understory layer (Alvarez, Gracia, and Retana 2012), and normally the heterogeneous terrain morphology (Pausas et al., 2004) can add more complication to the process. All these factors cause ITD methods to miss trees (Adhikari et al., 2023; Chen et al., 2018), and therefore produce biased – and normally underestimated – outputs. However, the increased predictive performance when modeling plots with average DBH >25 cm could suggest that the ITD method works better in more mature (i.e. more dominant and codominant trees) forests.

In addition, our ITD methodology relied on a pre-established DBH – HT relation, as in Kankare et al. (2013), which introduced variability. Relationships between DBH and HT can be affected by surrounding tree density or local conditions (Sumida et al., 2013; Zhang et al., 2004) and more complex relations have shown better performance (Özçelik et al., 2018). However, we opted for a general and simpler relation, as we were aiming for a landscape approach. While this general approach introduced another source of error to the ITD methodology, the goodness of fit of the chosen model was similar to other equations tested for the same species in the Mediterranean region of Aragón, Spain (Cabanillas-Saldaña, 2010). To minimize the impact of the biased estimates in the predictions, we used them as covariates for calibrated linear regressions, and as suggested by Kankare et al. (2013).

Area-based approach modeling

Regarding the nature of the ABA models, we found that the average height return had a very strong effect both on above-ground biomass and volume, like other previous studies (Laurin et al., 2020; Li et al., 2015). We included the percentage of first returns above 5 m as a descriptor of vertical structure of the forest to improve the models, as it is reported to be closely related to above-ground biomass (Toraño Caicoya et al., 2015). The final models for above-ground biomass presented errors that were similar to previous studies (e.g. Laurin et al., 2020; L. Li et al., 2015), yet lower than for the same species in the Iberian Peninsula (Domingo et al., 2018; Domingo, Montealegre, et al., 2019). Return height coefficient of variation, which is related to forest vertical structure, improved the prediction of number of stems (Sumida et al., 2013). After height variability, also found important by A. Li et al. (2017), the other covariate chosen was the first echo cover index at the average height return to build the model of stem density. Similarly, Pearse et al. (2019) also found that for both stem density and basal area, a descriptor of canopy closure was important. The errors for stem density under different forest conditions have been reported between 24% and 67%, even when modeled with more complex approaches (Cao et al., 2019; Fekety et al., 2018). Regarding basal area, the 25th percentile of the height returns was selected as the unique explanatory variable, variable also found key in Domingo et al. (2018). After visually inspecting the data, we found that in most of the plots, the 25th percentile corresponded to the height threshold prior to where most of the returns were in higher parts of the canopy. In other words, this variable could be understood as a proxy of height to live crown, which could explain the relation with basal area (Yang & Huang, 2018). Regarding the modeling errors, they were also in line with previous studies (e.g. Cao et al., 2019; Pearse et al., 2019). Finally, we used a percentile-based approach for dominant height modeling. Given that dominant height is better correlated to higher percentile, we used the 95th percentile and got a similar error to other point-cloud -based analyses (e.g. Gopalakrishnan et al., 2019; Ritz et al., 2022). While different definitions of dominant height can affect the posterior modeling (García & Batho, 2005; Sharma et al., 2002), we did not find any substantial differences.

Assessing the variability of the model’s predictions

We did not find evidence that the time difference between the ground measurements and the LiDAR acquisition affected the model residuals. This result is supported by previous studies in different regions (e.g. González-Ferreiro et al., 2012; Marino et al., 2022; Mauro et al., 2021), where a 3-year time difference is considered negligible, supporting the relative independence of the forest structural traits predictions from time difference in data acquisition. While some work has shown some effect due to time lag between ground and LiDAR measurements (Oh et al., 2022), these differences are normally in faster-growing species than Pinus halepensis. Likewise, pulse densities ranging from 0.7 to 6 pulses m⁻² did not affect the model residuals, in line with other work (Pearse et al., 2019; Strunk et al., 2012).

The models performed equally for each historical use, making it generally applicable at the landscape scale. However, we found differences in observed site productivity due to past land covers where historic agricultural land (normally terraced land with deeper soil), had a much higher SI than other historic land covers. Mediterranean abandoned terraced soils are generally more productive as (1) they were allocated in the most suitable areas for agriculture in the first place (Abadie et al., 2018), (2) previous agricultural/livestock practices improved their chemical, physical, and fertility properties, as long as no soil erosion took place (Stanchi et al., 2012; Zornoza et al., 2009), and (3) they were normally much deeper (Djuma et al., 2020). It is logical to expect that agricultural-related benefits are transferred to forest productivity, as suggested by Vance (2000), sometimes after the soil experiments an edaphic maturation process (Arnaez et al., 2011). This carry-over effect of historic use on current forest attributes has also been found in the Amazon rainforests (Vieira Wandelli & Martin Fearnside, 2015) and the South American Mediterranean climates context for above-ground biomass (Hernández et al., 2016), but also in the European Mediterranean environment for post-fire regeneration (Pausas et al., 2004) or plant species composition and diversity (Kouba et al., 2015).

Limitations and future directions

The first limitation of the methodology presented in this study is the scale of applicability. These ABA models have been developed based on a 20-m size grid matching 10-m radius plots, and this area mismatch can shift the predictions to underestimation, but also help to reduce the prediction error (Packalen et al., 2019). The GPS error, even though normally reported to be lower than 7 m for commercial hand-held devices (Gopalakrishnan et al., 2019), can potentially add some noise. However, if long enough occupation times are waited, these errors can be substantially reduced (Tomaštík & Everett, 2023). Though spatial autocorrelation with neighboring pixels improves estimates in more homogeneous environments, such as forest plantations (Prior et al., 2022), this effect can be lower in a more heterogeneous environment such as the Pinus halepensis forests. Therefore, we encourage the use of this models for a landscape scale, to obtain average values for larger areas, rather than to apply them in high resolution applications. Nevertheless, we believe that this methodology is simple enough to extrapolate to other similar environments of the natural habitat of this species.

Secondly, these models have been generated using publicly available LiDAR data, with the aim of developing a widely applicable tool for stakeholders, both in time and space. We have shown that there is no effect on the models’ performance regarding time lag between the ground and LiDAR data, making them very flexible, as a single LiDAR flight can be used for over 3 years to model forest attributes. However, it is still unknown whether two consecutive LiDAR acquisitions (every 5 years for PNOA LiDAR), will be able to capture the actual variation in the forest properties, at least when following the ABA methodology presented here.

Regarding the analysis of land history and past land cover, the orthophoto interpretation process is known to be challenging (Delgado-Artés et al., 2022). Also, while the anthropic effect on the landscape is well known, the segmentation of the past land uses, and especially their spatial representation is somewhat complex. For instance, Ibiza’s Mediterranean forests were subject to very intense forest uses (Ribas-Costa & Dopazo, 2020) such as timber, pine tar or lime production, and the visual interpretation of natural vs. human driven historic sparse forest land can sometimes be misleading. To partially overcome this issue, we included on-the-ground interpretation of the plots prior to their categorization, yet discrepancies could still occur.

Finally, in this study, we modeled dominant height as the preliminary step to generate a forest productivity map of Pinus halepensis forests in the Mediterranean basin. We performed the productivity analysis regarding historical land use, so we could better understand and interpret Pinus halepensis productivity potential.

Conclusion

We have shown that PNOA LiDAR can be used for large-scale forests’ assessment, presenting a methodology widely applicable in time and space. We provide a tool to estimate Pinus halepensis stands’ structural traits that are accurate enough to generate reliable maps of above-ground biomass, volume stock, tree density, basal area, and dominant height. These maps will benefit not only private landowners by being able to quantify the potential of their forest (i.e. operational scale) but also for public administration and land-use policy makers, as it can be used for landscape (i.e. tactic and strategic scale) assessment of Pinus halepensis forests’ potential in terms of both productive and protective forestry. For instance, mapping above-ground biomass will be useful to improve current estimations of carbon, a key aspect in climate change mitigation; mapping volume stock will imply a boost on the local-based Mediterranean timber industry, opening new management opportunities and allowing precision forestry; and mapping dominant height will allow, eventually, to map forest productivity, crucial not only to provide a current interpretation of the forest, but also a measure of its potential, a concept closely related to future landscape understanding, and therefore enhancing better management.

Author contribution

VARC conceptualized the study, carried out the fieldwork, performed the analyses and wrote the manuscript. RC contributed to the manuscript revision and reviewed the language. AGG supervised VARC, contributed to the manuscript revision and helped with statistical analyses.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The PNOA LiDAR and orthophotography data is openly available at https://centrodedescargas.cnig.es/CentroDescargas/index.jsp. Forest inventory data is property of the landowners and managed by TERRAPI WORLD S.L. and will be available from the corresponding author (VARC) upon reasonable request.

Additional information

Funding

This work was supported by the Doctoral INPhINIT Retaining fellowship LCF/BQ/DR22/11950023 from ‘la Caixa’ Foundation [ID 100010434] awarded to VARC in June 2022. Additional support was provided by AG through grant PID2021-127241OB-I00, funded by MCIN/AEI/10.13039/501100011033, and by “ERDF A way of making Europe”; NSF I/UCRC Center for Advanced Forestry Systems [Award #1916552], the Forest Productivity Cooperative, the Department of Forestry and Environmental Resources at NC State University, and the company TERRAPI WORLD S.L. plus the landowners who shared the inventory data.