*Article* **Characterizing the Error and Bias of Remotely Sensed LAI Products: An Example for Tropical and Subtropical Evergreen Forests in South China**

**Yuan Zhao 1,2, Xiaoqiu Chen 1,\*, Thomas Luke Smallman 2,3, Sophie Flack-Prain 2,3, David T. Milodowski 2,3 and Mathew Williams 2,3**


Received: 31 August 2020; Accepted: 21 September 2020; Published: 23 September 2020

**Abstract:** Leaf area is a key parameter underpinning ecosystem carbon, water and energy exchanges via photosynthesis, transpiration and absorption of radiation, from local to global scales. Satellite-based Earth Observation (EO) can provide estimates of leaf area index (LAI) with global coverage and high temporal frequency. However, the error and bias contained within these EO products and their variation in time and across spatial resolutions remain poorly understood. Here, we used nearly 8000 in situ measurements of LAI from six forest environments in southern China to evaluate the magnitude, uncertainty, and dynamics of three widely used EO LAI products. The finer spatial resolution GEOV3 PROBA-V 300 m LAI product best estimates the observed LAI from a multi-site dataset (*R*<sup>2</sup> = 0.45, bias = <sup>−</sup>0.54 m2 m<sup>−</sup>2, RMSE = 1.21 m2 m<sup>−</sup>2) and importantly captures canopy dynamics well, including the amplitude and phase. The GEOV2 PROBA-V 1 km LAI product performed the next best (*R*<sup>2</sup> <sup>=</sup> 0.36, bias <sup>=</sup> <sup>−</sup>2.04 m2 <sup>m</sup><sup>−</sup>2, RMSE <sup>=</sup> 2.32 m2 <sup>m</sup><sup>−</sup>2) followed by MODIS 500 m LAI (*R*<sup>2</sup> = 0.20, bias = <sup>−</sup>1.47 m2 m−2, RMSE = 2.29 m2 m−2). The MODIS 500 m product did not capture the temporal dynamics observed in situ across southern China. The uncertainties estimated by each of the EO products are substantially smaller (3–5 times) than the observed bias for EO products against in situ measurements. Thus, reported product uncertainties are substantially underestimated and do not fully account for their total uncertainty. Overall, our analysis indicates that both the retrieval algorithm and spatial resolution play an important role in accurately estimating LAI for the dense canopy forests in Southern China. When constraining models of the carbon cycle and other ecosystem processes are run, studies should assume that current EO product LAI uncertainty estimates underestimate their true uncertainty value.

**Keywords:** remotely sensed LAI; field measured LAI; validation; magnitude; uncertainty; temporal dynamics

#### **1. Introduction**

The Leaf Area Index (LAI), the total leaf area per unit ground area, is a key biophysical variable playing an important role in global carbon, water, and energy cycles [1,2]. As such, it acts as an important parameter for several applications, such as land surface models [3], ecological models [4], and yield prediction models [5]. The amount of leaf area has a first-order control on photosynthesis, transpiration, and absorption of radiation, varying in both space and time. LAI seasonal dynamics provide information about phenological processes of canopy development, senescence, and plant

traits [6]. Therefore, retrieving LAI over large areas and having a good knowledge of their yearly variations, errors, and bias is extremely important. Such information is central to accurately estimating primary productivity, understanding land surface-atmosphere exchanges, and detecting the response of terrestrial vegetation to climate change [7]. It is also beneficial for a large remote sensing community because it provides insights for the interpretation and correct usage of the LAI maps. Satellite-based Earth Observation (EO) offers the opportunity to retrieve information on LAI which is global in coverage and at increasing temporal and spatial resolution [8–12]. Robust estimates of uncertainties associated wth EO LAI estimates are needed to ensure their appropriate integration within further analyses such as data assimilation [6,13]. However, such robust evaluations are rarely possible due to the scarcity of in situ data at appropriate temporal and spatial resolutions [14,15].

LAI product uncertainty information represents the performance of the products' algorithm and reflects the uncertainties in the input data, model imperfections, and the inversion process [16–18], which could be called the theoretical uncertainties [19]. Evaluation studies typically assess the theoretical uncertainty of LAI based on the standard from Global Climate Observing System (GCOS) [20], which sets a target for absolute LAI uncertainty range to be within 0.5 m2 m−<sup>2</sup> and a relative uncertainty of less than 20%. However, few studies validate the uncertainty estimates of EO LAI using in situ data [14] at appropriate scales. A key challenge for the validation process is the mismatch in the scale of satellite products (typically > 300 m) compared to that of the in situ data (<10 m). Heterogeneity of LAI at scales finer than the satellite resolution [21] means that simple comparison between measurements at different scales can be problematic [22].

The validation of LAI at temporal scales includes a more important component, the LAI temporal dynamics, besides the traditional absolute value checking [23]. The temporal dynamics of LAI time series contain useful ecological information (e.g., phenological infomation), which help study the relationship between plant phenology and climate [24] or assist in land cover classification [25]. However, the temporal dimension has been neglected for most of the LAI validation work [23]. To validate the seasonal variations of LAI products, higher field sampling rates are need, especially for the Evergreen Broadleaf Forests (EBF) across several regions (Africa, Eurasia, South America, Australia and Asia) [14].

Despite the important role played by tropical forests in regulating the global carbon [26,27] and energy balance [28], tropical phenology is highly uncertain. Large-scale monitoring of tropical forest LAI with satellites is hindered by several challenges, including signal saturation [23], poor observation conditions (cloud-aerosol contamination, etc.) [29], coarse spatial and temporal resolutions [30], and scarcity of both ground data and detailed validation tests [31–33]. Uncertainties are therefore particularly high in these regions [34]. Some progress has been made in the Amazon [35–39] and African forests [40–42]; however, the Southern China region has typically been neglected in remote sensing phenology studies [43–46]. Consequently, the phenological character of forests in this region remains uncertain, with additional complexity driven by fragmentation [47], high species diversity, and complex topography [48]. The lack of detailed validation studies in Southern China means that there is little information regarding the suitability of global EO LAI products for the forests in this region, which play a major role in the carbon sink across China, and span the climate gradient from the subtropics to tropics [49,50].

In this study, we take advantage of a network of 6 sites with ~8000 ground-based measurements of LAI for forests located across the tropical and subtropical Southern China region to fully evaluate three satellite-based remote sensed LAI products with global coverage (Moderate Resolution Imaging Spectroradiometer (MODIS) 500 m and PROBA-V GEOV2 1 km and V3 300 m). Specifically, we address the following questions:


To answer these questions, we evaluate errors across multiple LAI retrieval algorithms and spatial resolutions and their associated temporal dynamics. We use additional fine resolution satellite data to evaluate the heterogeneity of canopy cover at the scales intermediate between in situ data and the satellite products at each location. Finally, we discuss the importance of robust error characterization for LAI products in the context of constraining models of the terrestrial carbon cycle and other ecosystem processes. We aim to characterize and understand errors and bias in the EO LAI products but do not make retrievals over the entire region. Whilst our scientific questions are focused on the Southern China region, the approach used can be applied elsewhere if in situ data is available. Assessing the uncertainties in LAI products through comparison with in situ measurements, i.e., direct validation is critical for their proper use in land surface models [51,52]. A better understanding of the uncertainties embedded in current LAI products will improve the assimilation of the LAI into land surface modeling studies [53,54] by providing a more robust error weighting [55].

#### **2. Materials and Methods**

#### *2.1. Study Area*

Site-level LAI observational data are collected from CERN (Chinese Ecosystem Research Network) [56]. The six sites included are Ban Na Forest (BNF), Ai Lao Forest (ALF), Gong Ga Forest (GGF), Ding Hu Forest (DHF), He Shan Forest (HSF), and Shen Nong Forest (SNF). These sites are spread across southern China between 21.9◦ and 31.3◦N, covering major plant functional types from northern subtropical zones to the northern tropical zones in this region (Figure 1). The regional climate is characterized by a wet, warm summer and a dry, mild winter [57]. The site-specific Koeppen climate classifications [58] are Cwa (Temperate zone warm summer; dry winter) for ALF, BNF, HSF, and SNF. Cfa (Temperate zone hot summer; no dry season) for DHF and Dwb (Cold zone dry winter; cold summer) for GGF. The mean annual temperature ranges between 4.2 and 21.8 ◦C. Annual precipitation ranges between 1506 mm and 2175 mm. The sites vary in elevation (above sea level) between 70 and 3160 m (Table 1). Forest types include the tropical seasonal rainforests, mixed coniferous forest, evergreen and deciduous mixed broad-leaved forest. The sites used support publically forests which are natural, except for forests at HSF and DHF, which support planted forests over 40 years old. All forests have at least three vertical layers: woody, shrub, and herbs. The tree density (trees per hectare) varies widely across sites. The highest tree density reaches over 7000 and the lowest is around 600 (Table S1). The range of the mean diameter breath heights is 5–20 cm (Table S1). The dominant species composition differs between sites but mainly consists of the evergreen species (Table S1). An image of forests measured at each site is included in Figure S1.

**Figure 1.** Map of the study sites in Southern China.


**Table 1.** Background information for each site. MAT: mean annual temperature; AP: annual precipitation.

#### *2.2. Field LAI Data Measurements*

In each of the six sites, there are one to four subplots at which LAI is estimated, and each subplot area ranged from 400–10,000 m2 (Figure S1). Field measurements of LAI were made each month using an LAI-2000 [59]. Field LAI measurement processing is consistent for each site and is based on the unified data collection and quality control protocols specified for CERN [60]. The sampling positions are distributed evenly and diagonally within the plot. The horizontal distance between each sampling points or plot boundary are to exceed 15 m and 10 m separately to avoid overlap sampling and reduce marginal effects. LAI-2000 were used to scan the canopy twice a day (8:00 am and 16:00 pm) to get the three-layer LAI (woody, shrub, and herbs) at each sampling point at the measurement day. The position for each sampling measurements and scan height is fixed. The median number of monthly measurements made at each site ranged from 6 to 80 (Table S2). BNF is the most intensively surveyed site, with between 203 and 360 measurements per year; therefore, BNF is particularly suited to evaluating the temporal dynamics of EO LAI products. For the other sites, only a subset of years and growth season months are measured (Table S2). A total of 795, 9 ground measurements from 2005 to 2017 were acquired across all the sites (Table S2).

#### *2.3. Earth Observation LAI Estimates*

In this study we evaluate three global EO LAI products. MODIS 500m product (MOD15A2H) provides estimates of LAI at 500 m spatial resolution with an 8-day interval from 2000-present [10]. GEOV2 1 km [11] and GEOV3 300 m LAI [12] were derived from the SPOT/VEGETATION sensor data at a 10-day interval, at spatial resolutions of 1/112◦ and 1/336◦, respectively (approximately 1-km and 300 m at the equator), in the Plate Carrée projection. For the MODIS 500 m LAI product, we consider only the products derived from the main algorithm, which is based on the use of Look Up Tables (LUTs) built for six different plant functional types, with simulations from a three-dimensional radiative transfer model [61]. For the GEOV2 1 km and GEOV3 300 m LAI products, LAI is estimated using Neural Networks(NNTs) applied on Top-of-Canopy (TOC) input reflectance in the red, near-infrared, and shortwave infrared bands, at 1km resolution and 300m, respectively [11,12]. The period of MODIS 500 m LAI products and GEOV2 1 km LAI products is from 2005 to 2017 and 2014 to 2017 for the GEOV3 300 m LAI product (released from 2014).

The definition of uncertainty information varies between LAI products. For MODIS LAI, the standard deviation is used to measure the uncertainty of pixel LAI values, which is calculated over all acceptable solutions of a look-up table (LUT) retrieval method [17]. For GEOV2 1 km and GEOV3 300 m LAI products, the uncertainties are computed as the RMSE between the final decadal value and the daily NNTs estimates in the compositing period [12,62]. In this study, they are both treated as the LAI products' uncertainty. The quality assessment information is used to clean the LAI and inform uncertainty information at the pixel level for all of these LAI products. For MODIS 500 m products, data effected by cloud shadow, internal cloud masks, or aerosol are removed and only main retrievals

were used in the analysis. For Copernicus LAI data (GEOV2 1 km and GEOV3 300 m), pixels that were filled or interpolated, or out of LAI range, were removed.

#### *2.4. Mapping Heterogeneity of Canopy Properties*

To evaluate the upscaling of field LAI measurements to satellite products, Landsat Level-2 Surface Reflectance products including Enhanced Vegetation Index (EVI) and Normalized Difference Vegetation Index (NDVI) derived from Landsat 5 (TM), Landsat 7 (ETM+) [63] and Landsat 8 (OLI) [64] scenes are also analyzed. EVI is calculated from red (R), near-infrared (NIR) and blue band (B) based on Equation (1) and incorporates parameters to adjust for canopy background, atmospheric resistance.

$$EVI = 2.5 \ast \left( \left( \text{NIR} - \text{R} \right) / \left( \text{NIR} + 6 \ast \text{R} - 7 \ast \text{B} + 1 \right) \right) \tag{1}$$

NDVI is calculated as a ratio between the R and NIR values in a traditional fashion (Equation (2)).

$$NDVI = \begin{pmatrix} NIR - R \end{pmatrix} / (NIR + R) \tag{2}$$

Landsat 5 (TM), Landsat 7 (ETM+), and Landsat 8 (OLI) Level-2 surface reflectance products were generated, using the Land Surface Reflectance Code (LaSRC) and the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) algorithms [65,66]. The LEDAPS and LaSRC surface reflectance algorithms correct for the temporally, spatially, and spectrally varying scattering and absorbing effects of atmospheric gases and aerosols, which is essential to derive the Earth's reflectance surface values. Landsat Spectral Indices are generated at 30 m spatial resolution on a Universal Transverse Mercator (UTM) mapping grid. The Spectral Indices could be obtained from the USGS Earth Resources Observation and Science (EROS) on-demand processing system.

The GTOPO 30 1 km elevation map [67] in the China region with Landsat scene footprints at each site is presented in Figure S2. In each of the six sites, there are one to four sampling subplots. The location of each site with the sampling subplots is shown on the Landsat scenes with path and row information (Figure S3) and on the 90 m SARTM elevation map [68], respectively (Figure S4). There are in total 778 Landsat 5–7 and Landsat 8 scenes with good quality (collection category is Tier 1) and cloud cover smaller than 20% of the total scene area from the year of 2005 to 2017 for all of these forests. Because there is no quality assessment information for the Landsat spectral indices, we use the pixel quality assurance band (pixel\_qa) instead to obtain quality information. Pixels with cloud shadow, aerosols, etc., are filtered out, and we keep pixels with the best quality.

Landsat EVI and NDVI were extracted at the field sampling plots spatial level (20~100 m). Based on the Landsat scene quality information, we cleaned and kept just good quality Landsat pixel data. Then, daily EVI and NDVI were retrieved from the area-weighted averaged Landsat pixels inside the field sampling plot. We averaged the daily EVI and NDVI values by month. Finally, the monthly EVI and NDVI at the field sampling plots level were regressed with corresponding monthly in situ LAI values to calculate the statistics (R2, etc.).

We investigated whether forest cover fraction in different LAI products' spatial resolution levels influenced product bias and error. We use the 30 m resolution GlobeLand30 land cover product for 2010 [69] to extract the forest cover fraction inside the relevant pixels for each of the LAI products at each site. We classify the land cover for each 30 m pixel on different spatial scales (LAI products' scale: 300 m, 500 m, and 1 km) at each site and plot. We only kept the 30 m pixels located entirely inside the LAI pixels. The total number of 30 m pixels and pixels classified as forests were both counted and summed to determine the total area and the forest areas at different spatial scales. Then we can obtain the forests cover proportion at different spatial scales for each plot and site. Finally, the forests' cover proportion was compared with the LAI products' bias at each site.

#### *2.5. Statistical Analysis*

Daily satellite LAI and uncertainty information were retrieved from the three LAI products at the location of each subplot of the six sites (Figure S1). If the field subplots overlapped with multiple pixels for a given LAI product, then the satellite-based value was estimated using the weighted sum of the contributing pixels, weighted by the area of overlap (Equation (3)). Data were cleaned based on the products' quality assessment information, before aggregation into monthly mean values, and comparison with the mean monthly ground-based LAI measurements at each site and plot. We calculate R2, bias, and RMSE between the field and satellite LAI time series to evaluate the accuracy and error of the three LAI products.

$$LAI\_R = \sum\_{i=1}^{n} \frac{Ai}{A} \ast Pi\_i \tag{3}$$

where *LAIR* is the extracted remote sensing LAI value; *n* represents the total number of the contributing pixels which overlapped with the field plots; *A* is the total area for the plots, *Ai* is the overlapped area between pixel and plots; *Pi* is the remote sensing pixel value.

LAI products' uncertainty information was evaluated using the ratio of LAI products' uncertainty (LAIU) with the bias against the field measurements, *r*, as follows (Equation (4)),

$$r = \frac{LAI\_{U}}{|Bias|} \,\prime \tag{4}$$

where the bias (Equation (5)) is defined as:

$$Bias = LAI\_{product} - LAI\_{field} \tag{5}$$

The ratio, *r*, can be used to evaluate whether or not the product uncertainty is robust. When r ≥ 1, the bias in the LAI product is within the reported uncertainty range, indicating that the uncertainty estimate for LAI is robust; conversely, if r < 1, the reported uncertainty fails to capture the observed bias and therefore the uncertainty estimate is not robust. We calculate this ratio *r* for each month at different plots and sites.

#### *2.6. Calculation of Amplitude, Phases and Periods for LAI Time Series*

We use metaCycle [70] to calculate the periodic characteristics (periods, phases, and amplitude) for the monthly field, MODIS, GEOV2 and V3 LAI time series. This method can calculate the periodic characteristics for time series with missing values [70], which frequently occurred for the satellite LAI times series. MetaCycle incorporates three methods: ARSER (Autoregressive spectral analysis) [71], JTK\_CYCLE(Jonckheere-Terpstra-Kendall algorithm) [72], and Lomb-Scargle [73], for periodic signal detection, and it could output integrated analysis. The integrated period from MetaCycle is an arithmetic mean value of multiple periods, while phase integration based on the mean of circular quantities calculated from the above three mentioned methods.

MetaCycle recalculates the amplitude with the following model (Equation (6)):

$$Yi \;= \; B + TRE \star \left( ti - \frac{\sum\_{i=1}^{n} ti}{n} \right) + A \star \cos \left( 2 \star \pi \star \frac{ti - PHA}{PER} \right) \tag{6}$$

where *Yi* is the observed value at time *ti*; *B* represents the baseline value (mean value) of the LAI time series; *TRE* is the trend level of the time-series profile; *A* is the amplitude of the waveform. *PER* and *PHA* are the integrated period and phase, respectively. *PER* represented the periods of the LAI time series; here, it is 12 months. *PHA* represented the phases of the LAI time series at the period of 12 months. "*i*" is the time points of the time series. "*n*" represents the total number of the time points of the time series. In this model, only *B*, *TRE*, and *A* are unknown parameters and can be calculated using ordinary least squares. Fisher's method is implemented in MetaCycle for integrating multiple *p*-values. Relative amplitude (rAMP), which is the *A*/*B*, can be used to compare the amplitudes between time series with different baselines levels.

In this study, each LAI time series is decomposed into the baseline value (mean value), the trend, and the wavelength with specific amplitude, period, and phase. The integrated phases can be used to compare whether the LAI time series is lead, lag, or synchronous. The amplitude can be used to quantify the magnitude of seasonal fluctuations for LAI time series and the relative amplitude can be used to compare the LAI seasonal fluctuations from different sites and data source. We apply this method to the monthly averaged field LAI time series at BNF, and the monthly averaged retrieved satellite LAI at each site because only BNF has frequent enough field measurements over multiple months from the year 2005 to 2017.

#### *2.7. Software*

Remote-sensing data processing was carried out by Rv3.4.0 [74] using the package "raster" [75]. Amplitude, phases, and periods of LAI time series were calculated using the functions "meta2d" in the package 'MetaCycle' [70] in R. Figures and maps were produced in R and ArcGIS 10.4.

#### **3. Results**

#### *3.1. Accuracy of the LAI Products Magnitude*

Satellite derived mean LAI estimates were typically lower than the field measurements, which varied between 3.2 and 5.6 m<sup>2</sup> m−<sup>2</sup> (Table 2). GEOV3 300 m estimates were 10% lower than field measurements (mean = 4.50 m<sup>2</sup> m−2), MODIS 500m were 30% lower (mean = 3.6 m<sup>2</sup> m−<sup>2</sup> and GEOV2 1 km were 40% lower (mean = 3.07 m2 m−2) (Figure 2; Table 2). LAI mean and relative uncertainty increased with product resolution: GEOV3 300m (0.22, 7%) < MODIS 500 m (0.37, 10%) < GEOV2 1 km (0.56, 19%).

**Figure 2.** Time series of the LAI for the field, GEOV3 300 m, GEOV2 1 km, and the MODIS 500 m for different sites. The timespan of the LAI time series is from 2005 to 2017. Each field data point represented the mean value of the measurements from the sample subplots in each month. "N" denoted the number of total sample measurements for each site, including the repeat measurements in each month at each subplot.



Field mean LAI was best captured by GEOV3 300 m (*R*<sup>2</sup> = 0.45; Figure 3), followed by GEOV2 1 km (*R*<sup>2</sup> = 0.36), and MODIS 500 m (*R*<sup>2</sup> = 0.20) for all of these six forests in this region. This performance varied between forests and all of these products showed the highest *R*<sup>2</sup> (Table 2) for the forests at SNF, which showed more LAI seasonal variation (Figure 2). GEOV3 300 m LAI also had the lowest RMSE and bias (RMSE = 1.21, bias = −0.41) compared with GEOV2 1 km LAI (RMSE = 2.32, bias = −2.04) and MODIS 500m LAI (RMSE = 2.29, bias = −1.47) (Figure 3, Table 2).

**Figure 3.** Regression between the field LAI and the satellite retrieved LAI. Each point represents the monthly LAI field measurement against the satellite LAI retrieval. The blue solid lines represent the regression between the field LAI against the satellite LAI for all the sites combined on a monthly scale. Point colors and shapes represent different sites.

#### *3.2. Robustness of the LAI Products Uncertainty*

The magnitude of the reported EO LAI uncertainty is typically smaller than the bias to in situ estimates (Figure 4), here shown as the ratio of uncertainty and reported bias (*r*), and this result is consistent across forests (Figure S5, Table 2). GEOV3 300 m has the smallest proportion of EO uncertainty estimates below the observed bias (86%) while MODIS 500m LAI and GEOV2 1 km have 87% and 91% of uncertainty estimates below the reported bias. The median ratio indicated similar average skill between both GEOV3 300 m (0.26) and GEOV2 1km LAI products (0.28) and MODIS 500 m LAI products (0.21) (Figure 4, Table 2).

**Figure 4.** Distribution of the ratio of EO uncertainty estimates to bias, based on a comparison against in situ LAI measurements. Dashed lines represent the median values and the color represents different EO products. The vertical black line is the reference line where the ratio is equal to 1. A ratio larger than 1 indicates the uncertainty estimate is larger than the bias and is therefore considered robust.

#### *3.3. Validation of the LAI Products Temporal Dynamics*

The metaCycle analysis showed that a statistically significant period of 12 months (*p* < 0.05) is present in the monthly field LAI time series at BNF and in all monthly satellite LAI time series at all sites, except for MODIS 500 m LAI, which showed no significant 12-month periodic characteristics at BNF. Overall, LAI retrieved from GEOV3 300 m showed very close values of the base, phase, amplitude and relative amplitude with the field LAI at BNF (Figure 5). LAI retrieved from GEOV2 1km showed a lower base and amplitude, similar phase, but larger relative amplitude compared with GEOV3 300 m LAI for most sites and the field at BNF (Figure 5). LAI retrieved from MODIS 500 m product showed the lower base values, but had a high variation of phases with one underestimated phase for the forest at ALF compared with the other two LAI products (Figure 5b). The relative amplitude of LAI retrieved from MODIS 500 m is lower for the forest at GGF and SNF, but higher at HSF compared with the other two LAI products (Figure 5c).

**Figure 5.** Periodic characteristics of the LAI time series from 2005 to 2017 for different sites. (**a**): Baseline values of the LAI time series; (**b**): Phase of the LAI time series; (**c**): Amplitude of the LAI time series; (**d**): Relative amplitude of the LAI time series. Point colors and shapes represent different sites. A detailed description of the calculation can be found in the Methods section of this paper.

#### *3.4. Evaluation of Landcover Heterogeneity*

The Landsat data provide a more spatially consistent comparison of surface reflectance against field data of LAI. The results (Figure S6) show that there are significant correlations between NDVI/EVI and LAI across sites, but the effect sizes are small (*R*<sup>2</sup> = 0.1). Thus, the relationships between remotely-sensed NDVI/EVI at the fine resolution and field-measured LAI are less significant than those between LAI products at the coarse resolution and field-measured LAI. The forest cover proportion at the pixel level varied between 38%–100% (Table S3). The landcover data show that LAI retrievals were degraded by the heterogeneity of landcover within product pixels. There was a clear pattern of LAI bias increasing with reductions in the proportions of forest cover at the pixel scale (Figure S7).

#### **4. Discussion**

The performance of three EO LAI products was evaluated at six forests across southern China. GEOV3 300 m LAI showed the best fit (*R*<sup>2</sup> = 0.45) and the smallest bias (bias = <sup>−</sup>0.54) (Figure 3) and the most similar LAI time series seasonal dynamic variation (phase and amplitude) compared with in situ LAI measurements (Figure 5). For the GEOV2 1 km LAI product, its performance had larger bias (Figure 3) and seasonality but similar phases of the LAI time series with the in situ observations compared with GEOV3 300 m LAI products (Figures 3 and 5).

The different performance of the two Copernicus products is likely due to the mismatch of the scaling between the EO pixel size and the field site. LAI is strongly non-linearly related to reflectance, making its estimation from remote sensing observations scale-dependent [76,77]. In contrast, the core operational algorithm (neural network techniques), data filtering and smoothing processes are similar for these two products. There are differences in the method used for temporal compositing, where temporal smoothing and gap filling using a climatology are used for GEOV2 1 km and

interpolation applied in GEOV3 300 m) [78]. Differences in the applied gap-filling approach between these two products do not impact our conclusion that resolution is the primary driver of performance improvement at 300 m relative to 1 km, as all gap-filled and interpolated retrievals were removed in our study. The in situ measurements were usually conducted in protected mature forest areas with high plant densities, complex canopy stratums, and rich species diversities (Figure S1, Table S1). However, at coarser spatial resolutions, pixels integrate heterogeneity over a greater diversity of landscapes, habitats, and species, distributed across a variety of stages of growth and succession [79]. Thus, fine resolution GEOV3 300 m LAI product showed lower bias and similar seasonal dynamic variation to the in-situ measurements compared to the higher bias for the GEOV2 1 km LAI product.

For Collection 6 500 m MODIS LAI, our results indicate that this product does not capture the dynamics observed in situ, irrespective of the accuracy of estimated in situ LAI (Figure 3). Furthermore, it does not adequately capture the LAI seasonal dynamics (Figure 5). This result is consistent with two validation studies which both showed MODIS LAI had the poorest performance for the evergreen forests in the south of China [80,81]. The MODIS LAI values for tropical evergreen forests are severely impacted by atmospheric conditions, especially clouds during the growing season (around 42% data are influenced by the cloud in this study, Figure S8), which lead to strong noise in the input reflectance data and affect the retrieval [78]. Additionally, the reflectance saturation usually happened in dense canopies and the main algorithm is sensitive to uncertainties in atmospheric correction, particularly when red and NIR BRFs are saturated [82,83]. This means the reflectance does not contain sufficient information to estimate the LAI value [84] and leads to the instability of the LAI retrievals [14]. All of these may result in the poor performance of representation of evergreen forests in southern China for the MODIS 500 m LAI products. Furthermore, the product does not adequately capture LAI periodic characteristics in comparison to the field data and showed temporal inconsistency with GEOV2 and V3 products (Figure 5). Jiang (2017) also found the large temporal inconsistency between existing global LAI products at a longer time scale [85]. Cammalleri (2019) found GEOV2 fAPAR showed a systematic overestimation of the fAPAR anomalies compared with the MODIS fAPAR and proposed a two-step harmonization procedure to remove this discrepancy [86]. However, the homogenization may alter the magnitude of the original fAPAR time series in an undesirable way. These results highlight the need to validate the temporal consistency between different satellite products and explore more solutions to deal with such inconsistencies. In addition, geolocation uncertainty due to the spatial mismatch could also have influenced validation reliability, As our study sites are not homogeneous and sampling area is consistently smaller than the pixel of remote sensing data, mean or median remotely-sensed LAI values of surrounding 3 × 3 array of pixels [15,87,88] cannot be used for validation. To solve this problem, future field measurements should be conducted at a larger spatial scale and across more homogenous habitats.

Overall, the absolute and relative uncertainty of LAI products tends to be smaller for the fine resolution LAI products (Table 2). Collection 300 m GEOV3 has the smallest absolute uncertainty at around 0.22. Collection1 km GEOV2 has the largest absolute uncertainty at around 0.56. Uncertainty magnitude for Collection 6 500 m MODIS LAI is 0.37. This performance is consistent with the LAI relative uncertainty (Table 2). If compared with the absolute uncertainty requirements (±0.5) and relative uncertainty requirements (20%) set by the GCOS [12], all three products (Collection 300 m GEOV3, 0.22, 5%; Collection 6 500 m MODIS LAI, 0.37, 10%; Collection 1 km GEOV2 0.56, 18%) appear satisfactory. This is consistent with two global studies comparing different EO LAI products (MODIS, CYCLOPES, and GLOBCARBON) at two coarse spatial resolutions (5 km and 1 km). More pixels can meet the absolute uncertainty and relative uncertainty requirements at the spatial resolution of 1 km compared with the 5 km [19,34]. All of these results indicated that the LAI product uncertainty is scale-dependent. The algorithms will produce smaller uncertainty estimated for the LAI values at finer spatial scales. This could make the uncertainty estimates of the fine-resolution products become very conservative.

Changes in LAI error with spatial resolution could be due to vegetation heterogeneity, especially for the forests and scale-dependent reflectance values. However, it is still unclear how product uncertainty changes over different spatial scales. For MODIS, the product uncertainty is the standard deviation over all acceptable solutions of a look-up table (LUT) retrieval method [17]; for Collection 300 m GEOV3 and 1 km GEOV2 LAI products, the quantitative uncertainties (LAI\_ERR) are computed as the RMSE between the final dekadal value and the daily NNTs estimates in the compositing period [12,62]. Standardizing different LAI products to the same spatial scale and averaging the LAI and uncertainty for different pixels during the upscaling and resampling process [19,89] could all generate errors [90]. When dealing with different spatial scales, the land surface properties (e.g., land cover proportion, soil properties) could change significantly. Thus, LAI retrieval algorithms are based on inherent empirical assumptions on the distribution of their parameters, which can depart significantly from the actual canopy and soil characteristics [14].

The ratio (*r*, in Equation (4)) of LAI product uncertainty to the bias of in situ measurements is dramatically lower than the reference value for all of these LAI products (Figure 4, Table 2). This result indicated that none of the three LAI uncertainty products tested here were able to provide a robust estimate for the in-situ measurements. The most promising product is the Collection 300 m GEOV3 LAI; however, while its accuracy is the best, its uncertainty is still too conservative to capture the offset from co-located field LAI measurements. This weakness reflects the complexity of this biome type and highlights that the algorithm used to generate the uncertainty information for these LAI products needs to be improved. For areas of high uncertainties, data producers may need to refine the algorithms and verify the information using additional data sources, particularly in situ data [34]. However, with respect to field data collection, more accurate measurements are also needed to test the results, such as destructive methods or using the litterfall traps. Direct measurement methods obtaining 'true' LAI values could also act as a reference to correct the bias and errors in indirect measurement techniques [31].

The accuracy of LAI data retrieved from remote sensing products is very important for a range of earth system models, for which LAI is a key internal variable. LAI uncertainty is particularly critical for model-data fusion in terrestrial carbon cycle studies [91]. The uncertainty produced for the existing remote sensing LAI products in the analysis would further lead to biases in vegetation productivity estimates [92–94]. For instance, the saturation of satellite-derived VIs generally results in an underestimation of GPP or NPP in areas with dense vegetation. Using finer resolution remote sensing data in carbon modeling studies could be an effective way to reduce the uncertainty in the estimates of C fluxes and/or stocks, but more robust errors are critical [95].

Uncertainty estimates of the observations are as important as the observations themselves because together they co-determine the outcome of the assimilation systems [96]. Bayesian approaches are commonly used data assimilation systems and the influence of observations are weighted by observational uncertainties. Uncertainties in data are critical as they often determine the outcome of analyses and forecasts [97]. The assimilation process requires clear error quantification for LAI to resolve model results and limit biases [98]. Besides, the presence of bias in a data stream will limit the utility of using multiple observation types in an assimilation framework. Therefore, it is imperative to characterize the error in the observations and understand better the error associated with the direct measurements of LAI to landscape pixel [99]. Our results indicated that none of these LAI uncertainty products are robust and therefore users should consider inflating the existing LAI product uncertainty when using these data sources in data assimilation frameworks [98,100,101] or in other studies for which the uncertainty of the product plays an important role.

One way of solving this spatial mismatch problem is to upscale the field measurements based on its relationship with the plant index retrieved from high-resolution images (e.g., 10–30 m resolution) [102], or using airborne Light Detection and Ranging (LiDAR) [103,104]. However, we failed to find a good relationship between the Landsat based EVI or NDVI and the field measurements (Figure S5). LAI-VI relationships are limited by a number of factors, including vegetation type, sun-surface-sensor

geometry, leaf chlorophyll content, background reflectance, and atmospheric quality [31]. Geolocation uncertainties due to spatial mismatch and insufficient overlap data between in situ measurements and good quality Landsat scenes partly lead to this result. Saturation effects occur for both field measurement and spectral plant indices for very dense canopies, which also increases the challenge of upscaling field LAI to generate high-resolution LAI reference maps. In addition, the temporal mismatch between field and Landsat data also contributes to the poor relationship between field and Landsat data. This result may be also associated with the complex forest environments (such as terrain, tree densities, canopy structures, and species diversities) in southern China (Figure S4, Table S1). Thus, their LAI cannot be accurately modeled based on simple spectral indices, even with relatively high resolution (30 m) satellite data. More complex upscaling approaches or those utilizing very high-resolution sensors (e.g., LIDAR) to bridge the scaling gap may yield better results [102,105,106]. In particular, there is an exciting potential for unmanned aerial vehicles to provide critical upscaling information from field to satellite [21].

When using EO LAI products in the absence of reliable validation data, we recommend that the provided uncertainties are treated cautiously. Remote-sensing LAI products require improved algorithms with particular attention given to generating robust uncertainty estimates. New cloud and aerosol detection techniques based on time series and spatial analysis may help to improve the uncertainty products and LAI estimated for evergreen forests [107]. In addition, the new remote-sensing systems such as Global Ecosystem Dynamics Investigation (GEDI) [108], National Aeronautics and Space Administration's (NASA's) Ice, Cloud, and land Elevation Satellite-2 (ICESat-2) and Synthetic aperture radar (SAR) provide promising ways to solve cloud coverage issues [109]. Given the importance of data uncertainties (e.g., integration within terrestrial carbon modeling frameworks [92–94]), there is an urgent need for robust uncertainty characterization, based on links to in situ observations. The findings in this study are constrained to southern China. Future studies could explore the LAI and error products for other regions, or focus on the seasons where products show the largest uncertainties [19].

#### **5. Conclusions**

In this study, we validated three global LAI products and their uncertainty products using ~8000 in situ measurements of LAI from subtropical and tropical forests in the southern China region. The finest resolution product, Collection 300 m GEOV3, performed best, with the lowest RMSE and the lowest bias. It also best captured the temporal dynamics observed in the in-situ dataset, including the magnitude, amplitude, and phases of the LAI time series. Collection 6 500 m MODIS LAI showed the poorest performance and importantly it did not capture the temporal dynamics observed in situ for the forests in this specific region. Critically, for all three LAI products, the accompanying uncertainties were all far smaller than the bias compared to the in situ measurements, indicating that each product uncertainty estimate is not robust. Given the importance of LAI uncertainty in the context of constraining models of carbon cycling and other ecosystem processes, users should use the product's uncertainty with caution and consider inflating the existing LAI product uncertainty when using it within data assimilation frameworks or other studies in which the uncertainty plays an important role.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2072-4292/12/19/3122/s1, Table S1: Detailed forest community information in the southern China region, Table S2: The detailed information for field LAI measurements at each site, Table S3: Distribution of the ratio of EO uncertainty estimates to bias, based on a comparison against in situ LAI measurements, Figure S1: Forest appearance for each forest in this study in the southern China region, Figure S2: TOPO30 DEM (1 km) in the China region. Figure S3: Location of the sampling subplots and the GEOV2 1km LAI pixels showed on the Landsat EVI scenes for different sites, Figure S4: Location of the sampling subplots and the GEOV2 1 km LAI pixels showed on the SARTM 90 m DEM images for different sites. Figure S5: Frequency distribution of the ratio of LAI products uncertainty with LAI products' bias against the in situ measurements for different sites. Figure S6: Regression analyses between the Landsat EVI/NDVI values and the field observed LAI values at sampling subplots, Figure S7: Regression analyses between the proportion of forests cover (LAI pixel level) and the remotely sensed LAI bias against the field measured LAI for different LAI products, Figure S8: Summary of the dates for the MODIS LAI pixels which were influenced by the cloud in this study.

**Author Contributions:** Conceptualization, Y.Z. and M.W.; methodology, Y.Z., M.W., T.L.S., and D.T.M.; formal analysis, T.L.S. and S.F.-P. and D.T.M.; writing—original draft preparation, Y.Z. and M.W.; writing—review and editing, M.W., X.C., T.L.S., D.T.M., S.F.-P., and Y.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China under grant number 41771049 and NCEO, CSSP Newton Fund, and China Scholarship Council.

**Acknowledgments:** The authors thank all the workers in the Chinese Ecosystem Research Network for field data collection and database establishment.

**Conflicts of Interest:** The authors declare no conflict of interest.

## **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Analysis of the Spatial Di**ff**erences in Canopy Height Models from UAV LiDAR and Photogrammetry**

## **Qingwang Liu 1,2,**†**, Liyong Fu 1,3,**†**, Qiao Chen 1,**†**, Guangxing Wang 4,\*, Peng Luo 1, Ram P. Sharma 5, Peng He 6, Mei Li 1, Mengxi Wang <sup>1</sup> and Guangshuang Duan 1,7**


Received: 19 July 2020; Accepted: 1 September 2020; Published: 6 September 2020

**Abstract:** Forest canopy height is one of the most important spatial characteristics for forest resource inventories and forest ecosystem modeling. Light detection and ranging (LiDAR) can be used to accurately detect canopy surface and terrain information from the backscattering signals of laser pulses, while photogrammetry tends to accurately depict the canopy surface envelope. The spatial differences between the canopy surfaces estimated by LiDAR and photogrammetry have not been investigated in depth. Thus, this study aims to assess LiDAR and photogrammetry point clouds and analyze the spatial differences in canopy heights. The study site is located in the Jigongshan National Nature Reserve of Henan Province, Central China. Six data sets, including one LiDAR data set and five photogrammetry data sets captured from an unmanned aerial vehicle (UAV), were used to estimate the forest canopy heights. Three spatial distribution descriptors, namely, the effective cell ratio (ECR), point cloud homogeneity (PCH) and point cloud redundancy (PCR), were developed to assess the LiDAR and photogrammetry point clouds in the grid. The ordinary neighbor (ON) and constrained neighbor (CN) interpolation algorithms were used to fill void cells in digital surface models (DSMs) and canopy height models (CHMs). The CN algorithm could be used to distinguish small and large holes in the CHMs. The optimal spatial resolution was analyzed according to the ECR changes of DSMs or CHMs resulting from the CN algorithms. Large negative and positive variations were observed between the LiDAR and photogrammetry canopy heights. The stratified mean difference in canopy heights increased gradually from negative to positive when the canopy heights were greater than 3 m, which means that photogrammetry tends to overestimate low canopy heights and underestimate high canopy heights. The CN interpolation algorithm achieved smaller relative root mean square errors than the ON interpolation algorithm. This article provides an operational method for the spatial assessment of point clouds and suggests that the variations between LiDAR and photogrammetry CHMs should be considered when modeling forest parameters.

**Keywords:** digital surface model; digital terrain model; canopy height model; constrained neighbor interpolation; ordinary neighbor interpolation; point cloud density; stereo imagery

#### **1. Introduction**

Forest structure information is a prerequisite for forest resource inventories and forest ecosystem modeling [1–3]. Various techniques and methods can be used to obtain such information. Light detection and ranging (LiDAR) and photogrammetry have the ability to depict the three-dimensional canopy structure of forests and therefore can be used to monitor structural changes over time [4–8]. Unmanned aerial vehicle (UAV) LiDAR, hereafter called LiDAR, can be used to accurately measure the spatial distributions of forest canopies with a high density of point clouds [9–13]. Laser pulses can penetrate the gaps between branches and leaves of tree crowns and detect middle parts or areas under forest [14–17]. LiDAR point clouds are usually classified as vegetation, ground and other objects, which are used to generate digital surface models (DSMs) and digital terrain models (DTMs) using interpolation algorithms [10,18,19]. Generally, a DSM is created using the maximum algorithm in a regular grid at the given spatial resolution, which is determined according to the point cloud density [20–22]. A canopy height model (CHM) depicts the variations in the forest canopy height above the terrain, and these height variations can be determined by subtracting the DTM from a DSM [3,10,14,23–27]. Many algorithms that estimate the parameters at the individual tree or sample plot level are developed based on a CHM [14–16,23,26–32].

UAV photogrammetry, hereafter called photogrammetry, usually acquires images with high areas of overlap (usually 60–90% along-track and 30–60% across-track), which are processed to build the spatial structure of forest canopies using stereo imagery algorithms, such as structure from motion (SfM) [25,32–37] or semi-global matching [24,38–41]. Dense point clouds reconstructed from image pairs are used to create a DSM with a similar algorithm used for LiDAR point clouds [23,42]. The distribution of dense point clouds is substantially affected by the image along- and across-track overlap, flying height, terrain, and other factors [43–46]. It is quite difficult to reconstruct the terrain under dense forest due to mutual obscuration among tree crowns. A photogrammetry CHM (P-CHM) can be generated as the difference between a photogrammetry DSM (P-DSM) and LiDAR DTM (L-DTM) [23,33,35,47]. Forest attributes can be precisely estimated using the CHMs of either LiDAR or photogrammetry [23,25,48–50].

The observation geometry of LiDAR is obviously different from that of photogrammetry [23,46]. LiDAR can directly measure the spatial positions of branches and leaves at different heights, while photogrammetry obtains the surface envelope of the forest canopy using a stereo imagery algorithm. Thus, canopy heights estimated by LiDAR and photogrammetry show inherent spatial differences, and these differences have not been discussed in depth in the literature. The spatial resolutions of LiDAR and photogrammetry CHMs vary from centimetres to metres based on the different point cloud densities [9,14,15,23,25–27,37]. Detailed canopy structures tend to be suppressed with a coarse spatial resolution. The densities of LiDAR and photogrammetry point clouds are unevenly distributed due to various factors, such as flying height and crown shadows, and these variations produce varied points within different grid cells at a specified spatial resolution. A cell with one or more points is referred to as an effective cell (valid cell), while a cell without any points is referred to as a void cell (no data, missing or blank cell). The void cells can form variable holes in the grid. These holes are different from canopy gaps (canopy openings) created by the snapping and falling of trees, the impacts of insects or pathogens, or the mortality of single trees or small groups of trees; moreover, these gaps will eventually be closed in the ecological processes of forest ecosystems [14,23].

The size of the holes in a grid can be expressed as the number of continuous void cells, which varies at differing spatial resolutions. The same-sized hole will have more void cells at a fine resolution than at a coarse resolution. Some void cells in a hole will have effective neighboring cells, which are referred to as outer void cells, while the other void cells within a hole are referred to as inner void cells. The outer and inner void cells can be estimated or interpolated from the effective cells by using the ordinary neighbor (ON) interpolation algorithm. The interpolated inner void cells are subject to more uncertainties than the interpolated outer void cells considering the spatial relevance [9,14,51]. The number of effective neighboring cells can also affect the uncertainties of interpolated cells. The interpolated cells with many effective neighboring cells will have higher confidence than those with few effective neighboring cells. The ON algorithm should be constrained to obtain more reliable values with more effective neighboring cells.

The distribution of effective or void cells in a grid can reflect the homogeneity of the point cloud to some degree. Certain questions remain regarding the distribution of point clouds and estimated canopy heights, such as (1) how can the distribution of point clouds for a given spatial resolution be quantitatively described, (2) what spatial resolution is optimal for the specified point cloud, and (3) are spatial differences of the canopy heights depicted by LiDAR and photogrammetry affected by different spatial resolutions?

The purpose of this study was to assess the spatial distribution of point clouds and compare the differences between CHMs derived by LiDAR and photogrammetry. Three descriptors of the spatial distribution were introduced, namely, the effective cell ratio (ECR), point cloud homogeneity (PCH) and point cloud redundancy (PCR), to characterize the point clouds and the derived CHMs in the grid. The ON and constrained neighbor (CN) interpolation algorithms were used to fill the void cells of CHMs, and the CN algorithm was used to distinguish small and large holes. LiDAR CHMs were used as references to assess the photogrammetry CHMs that were created at different flying heights and with different image overlaps. This article will provide an operational method for point cloud assessment and spatial distribution analysis.

#### **2. Materials and Methods**

#### *2.1. Study Site*

The study site (114◦05 E, 31◦52 N) is located in the Jigongshan National Nature Reserve of Henan Province, Central China (Figure 1). This area belongs to a transition zone from a subtropical to temperate climate zone. The area covered is approximately 340 by 360 m, and the terrain elevation ranges from approximately 115 to 198 m. The forest is dominated by deciduous broad-leaved trees, including sawtooth oak (*Quercus acutissima* Carruth.), Chinese cork oak (*Quercus variabilis* Blume), Chinese sweet gum or Formosan gum (*Liquidambar formosana* Hance), and bald cypress (*Taxodium distichum* (L.) Rich.). Abundant shrubs are observed in the understory layer. This site has been used for studies on atmospheric nitrogen deposition in forest ecosystems [52,53].

**Figure 1.** (**a**) Location of the study site in the Jigongshan Natural Reserve, Xinyang, Henan Province, China; and (**b**) true color map of the study site (red dots indicate the centres of ground plots).

The ground plots used were square with a size of 25 <sup>×</sup> 25 m (625 m2), which were set up on the basis of the existing studies [52,53], and there were a total of 28 ground plots. The parameters of individual trees with a diameter at breast height (DBH) ≥5 cm, including the DBH (cm), tree height (m), crown radius in four directions (m) and species, were determined in August 2017 (Table A1). The number of stems within the ground plots varied from 30 to 96. Lorey's heights of the ground plots ranged from 15.7 to 31.2 m.

#### *2.2. UAV Data Sets*

The LiDAR data and photogrammetry images were acquired in August 2017. The LiDAR system used in this study is a Velodyne laser scanning system (VLP-16) with a high-precision global navigation satellite system and inertial measurement unit (GNSS and IMU) mounted on an eight-rotor aircraft [54,55]. The laser sensor has 16 channels with maximum scan angles of 30◦ and 360◦ in the along- and across-track directions, respectively. More detailed specifications of the LiDAR and photogrammetry system are shown in Table 1. The flying heights of the aircraft varied from 50 m to 55 m above the terrain or take-off point of the UAV at a speed of approximately 4.8 m s<sup>−</sup>1.


**Table 1.** Light detection and ranging (LiDAR) and photogrammetry system specifications.

The LiDAR data acquired on 10 August 2017 did not cover the entire study area because the travel distances of some laser pulses exceeded the maximum range (120 m) of the LiDAR system. The areas with elevations below 130 m had almost no data, and they were mainly located in the west and southwest parts of the study site. The second LiDAR flight was designed to obtain supplementary point cloud data for the areas without data on 11 August 2017. The point clouds of the two flights were matched and combined as a single LiDAR data set (denoted as L-D55). Detailed information on the LiDAR data set is shown in Table A2. The point cloud density of LiDAR data varied from 0 to 1757 points m<sup>−</sup>2, with a mean of 168 points m−<sup>2</sup> in the grid at a 1.0 m resolution (Figure A1).

The images were captured using a Canon EOS 5DS camera with a high-precision GNSS mounted on a six-rotor aircraft (see Table 1 for detailed specifications) [56,57], and the lens had a focal length of 35 mm. The size of the image consisted of 8688 × 5792 pixels, which corresponds to the size of the sensor of 36.0 × 24.0 mm. The photogrammetry system was operated at different heights varying from 80 to 300 m (the ground sample distance (GSD) varied from 1 cm to 4 cm) with image overlaps varying from 64% to 84% to obtain five data sets (denoted as P1-D300, P2-D150, P3-D150, P4-D150 and P5-D80) covering the same area as shown in Table A2. P1-D300 was acquired at a flying height of 300 m with a corresponding GSD of 4 cm. P2-D150, P3-D150 and P4-D150 were obtained at 150 m with a GSD of 2 cm using different across-track overlaps and flying speeds. P5-D80 was obtained at a flying height of 80 m with a corresponding GSD of 1 cm. The images were captured with an exposure time of 1/800 s or 1/1000 s on cloudy or sunny days using an automatic aperture and an ISO of 160. The images were processed to build a dense point cloud of the forest structure using the SfM method [58]. Several steps were involved to generate high-precision dense point clouds. The images were aligned with high accuracy using the camera positions obtained by static differential processing of high-precision GNSS data. The camera and lens parameters, including the focal length, principal point coordinates, affinity and skew (non-orthogonality) transformation coefficients, radial distortion coefficients, and tangential distortion coefficients, were optimized to minimize the camera position errors using the optimize camera alignment tool [58]. The dense points were built with high quality and mild depth filtering to maintain detailed characteristics of the forest canopies. The point cloud density of the photogrammetry data sets varied from 1 to 12,761 points m<sup>−</sup>2, with means ranging from 257 to 2562 points m−<sup>2</sup> in the grid at a 1.0 m resolution (Figure A1).

#### *2.3. Data Processing*

The LiDAR point cloud was classified as ground, vegetation, building, and noise using TerraSolid [59]. The noise points were each carefully identified by visual evaluation. The ground algorithm was used to separate ground points from other objects with a maximum building size of 20 m, a terrain angle of 88◦, an iteration angle of 6◦ to the plane and an iteration distance of 1.4 m to the plane. The ground points were modified by visual evaluation to eliminate false ground points. The points of artificial objects, such as water towers and buildings, were manually recognized. The remaining points were classified as vegetation after other objects were excluded. The L-DTM was generated from the ground points by a triangulated irregular network interpolation algorithm with spatial resolutions of 0.1, 0.2, 0.5, and 1.0 m [9,23,35,60,61]. The original LiDAR DSM (L-DSM) was created from the vegetation and ground points by using maximum algorithms with spatial resolutions of 0.1, 0.2, 0.5, and 1.0 m [9,14,28,61]. The original LiDAR CHM (L-CHM) was produced based on the difference between the original L-DSM and L-DTM. The void cells of the original L-DSM were partially interpolated by the CN interpolation algorithm to generate the interpolated L-DSM. The interpolated L-CHM was produced based on the difference between the interpolated L-DSM and L-DTM.

The displacements between the LiDAR and photogrammetry point clouds were adjusted using the reference points extracted from the buildings and open ground area. The horizontal displacement was calculated based on roof ridges and building edges. The vertical displacement was computed according to the ground control points after horizontal displacement was applied. The aligned data sets were generated after the vertical displacement was applied. The photogrammetry dense point cloud was used to generate the original P-DSM by using the maximum algorithm. The interpolated P-DSM was created by the CN interpolation algorithm. The original P-CHM and the interpolated P-CHM were created from the corresponding P-DSMs normalized by the L-DTM [23–25,32,33,37,48]. The DSMs and CHMs from LiDAR and photogrammetry were masked using the boundaries of the buildings and water towers to maintain consistency. The boundaries of the buildings were digitized according to the building points of the point cloud. Seventeen water towers with heights of approximately 38 m were distributed across the study site. The centres of the water towers were determined and used to create circular boundaries by buffering with a radius of 1.2 m.

#### *2.4. Descriptors of Spatial Distribution*

This study used three descriptors of the spatial distribution, namely, the ECR, PCH and PCR, to quantitatively depict the characteristics of the point cloud in the grid at the given spatial resolution. The original DSMs or CHMs obtained from LiDAR and photogrammetry may have void cells due to complex conditions, such as the laser scanning patterns, flight attitudes, and sunshine conditions [14,31,51]. The ECR is defined as the proportion of effective cells to total cells (Equation (1)) and used to depict the distribution patterns of effective and void cells. The ECR reflects the uneven characteristics of point clouds to some degree. A small ECR means that many points are clustered in a small area.

$$ECR = \frac{N\_E}{N} \tag{1}$$

where *ECR* is the effective cell ratio, *NE* is the number of effective cells, and *N* is the number of total cells, including effective and void cells. If the *ECR* equals 1, then the point cloud is evenly distributed over each cell of the grid.

Different densities of points may occur in grid cells for different data sets with the same ECR. To explain the homogeneity of the point cloud distributed in the grid, the PCH is defined in Equation (2) by introducing the number of points per cell, which can be calculated as the product of the point cloud density per square metre and the cell area.

$$PCH = \begin{cases} \left(\frac{ECR}{D\_C}\right)^{1/D\_C} & \text{, if } D\_C < 1\\ \left(\frac{ECR^{D\_C}}{ECR^{D\_C}}\right) & \text{, if } D\_C \ge 1 \end{cases} \tag{2}$$

where *PCH* is the point cloud homogeneity, *DC* is the mean number of points per cell by *DC* = *n*/*N*, *n* is the number of total points, and *N* is the number of total cells. If *DC* equals 1, then one point is located in each cell and the *PCH* is only determined by effective cells. If the *ECR* equals 1, then the *PCH* will be 1, regardless of how many points are located in each cell.

The PCR depicts the richness of points in each cell as defined in Equation (3). A lower PCH corresponds to a higher PCR, which indicates that some cells hold more unnecessary points. An evenly distributed data set should have a higher PCH and lower PCR.

$$PCR = 1 - \frac{ECR}{D\_{\odot}} \tag{3}$$

where *PCR* is the point cloud redundancy, *ECR* is the effective cell ratio, and *DC* is the mean number of points per cell. If the *ECR* and *DC* both equal 1, then the points are ideally distributed evenly in each cell and the *PCR* equals 0. If the *DC* is much greater than the *ECR*, then the *PCR* approximately approaches 1.

The ECR, PCH, and PCR of the point cloud could also be calculated based on the interpolated DSM or CHM. The interpolation introduced some degree of variation in these descriptors of the spatial distribution.

#### *2.5. Constrained Neighbor Interpolation*

The original DSM or CHM had void cells in the grid at a specified spatial resolution due to the uneven distribution of point clouds. The void cells can be filled using the effective neighboring cells. The number of effective neighboring cells (i.e., not including void neighboring cells) for one void cell varies from 1 to 8, which affects the spatial variation in the filled cells [9,14]. The CN interpolation algorithm is used to calculate the interpolated values for the void cells (Equation (4)).

$$\mathbf{C}\_{i} = \frac{\sum\_{j=1}^{k} \text{NC}\_{j}}{k}, \ k \le 8, \text{ if } k \ge Q \tag{4}$$

where *Ci* (m) is the interpolated value of the *i*th void cell, *NCj* (m) is the value of the *j*th effective neighboring cell of the *i*th void cell, k is the number of effective neighboring cells of the *i*th void cell, and *Q* is the threshold of effective neighboring cells of the void cell. The *i*th void cell will be interpolated if the number of its effective neighboring cells is not less than the threshold (*Q*).

The value of the threshold *Q* varies from 1 to 8. The void cell will be interpolated if there is at least one effective neighboring cells; that is, *Q* = 1. The ON interpolation algorithm will similarly interpolate the void cells with any effective neighboring cells, which can be regarded as the special case (*Q* = 1) of the CN interpolation algorithm. The CN and ON interpolation algorithms use an iterative process to interpolate the void cells. Some void cells with effective neighboring cells can be interpolated in the first loop, while other void cells may be interpolated in the second or more loops.

The void cells of different holes have a varied number of effective neighboring cells (Figure 2). For example, the maximum number of effective neighboring cells of the void cells in Figure 2a is 4, which will not be interpolated if the threshold *Q* equals 5. Eight void cells will be interpolated in the first loop and 4 void cells interpolated in the second loop if the threshold *Q* equals 4. In practice, any holes could be interpolated by iterative loops if the threshold *Q* is less than 5. The runs of iterative loops will vary from several times to dozens of times depending on the number of void cells. The CN algorithm will continue to loop until every hole is either eliminated or reduced to a specific size as depicted in Figure 2. The double width cross pattern of 12 cells would be the minimum hole if *Q* equals 5 under normal circumstances in Figure 2a. Another rarely occurring case would be that the interleaved chess pattern in Figure 2e could not be interpolated. The small holes in Figure 2b–d will not be interpolated if the threshold *Q* equals 6, 7, or 8, while those holes can be filled if the threshold *Q* equals 5. Therefore, the threshold of *Q* = 5 is chosen for filling void cells in this study.

**Figure 2.** Minimum holes left by the constrained neighbor (CN) interpolation algorithm with different threshold values (*Q*) (void cells are white, and effective cells are grey): (**a**) *Q* = 5; (**b**) *Q* = 6; (**c**) *Q* = 7; (**d**) *Q* = 8; and (**e**) *Q* ≥ 5.

The hole in Figure 2a is larger than the other holes in Figure 2, and it is referred to as the diagnostic hole, which is used as the criteria to distinguish small and large holes. Small holes have one or more void cells, and the numbers of columns and rows of void cells were not less than 4. Large holes have twelve or more void cells, and both the numbers of columns and rows were equal to or greater than 4. The actual area of the diagnostic hole was determined based on the spatial resolution, for example, the area was 0.12 or 0.48 m2 with a spatial resolution of 0.1 or 0.2 m, respectively.

#### *2.6. Analysis of Spatial Resolution*

The spatial resolution will affect the number of points in each cell of the original DSM or CHM. A fine spatial resolution tends to result in many void cells in the grid, while a coarse spatial resolution will discard many points in the effective cells of the grid for a given point cloud, which represents a compromise for determining the spatial resolution [9,14].

The spatial distribution characteristics of point clouds depend on the spatial resolution. The ECR will increase as the number of void cells decreases at a coarse spatial resolution. The interpolated void cells by the CN algorithm will also, to some degree, cause the ECR to increase. The variation in the ECR can be used to analyze the optimal spatial resolution for the DSM or CHM, in which as much information about the point cloud is retained as possible so that the void cells can be confidently interpolated in the grid. Series of ECR pairs of original DSMs and interpolated DSMs by the CN algorithm can be calculated at different spatial resolutions. The difference between the ECR of the interpolated DSM and that of the original DSM is referred to as the ECR pair difference, which is used as an indicator to select the spatial resolution. The ECR pair differences might decrease from fine to coarse spatial resolution when the void cells of the grid are interpolated by the CN algorithm. The optimal spatial resolution could be determined as the corresponding spatial resolution if the ECR pair differences do not become greater than the specified threshold (for example, 0.10). The series of

spatial resolutions used as candidates is 0.1 m, 0.2 m, 0.5 m, and 1.0 m in this study. The optimal spatial resolution is considered relative to the other spatial resolution candidates.

#### *2.7. Data Set Assessment*

The estimated canopy heights obtained from the LiDAR data set were used as references to compare the estimated canopy heights obtained from five different photogrammetry data sets with spatial resolutions of 0.1, 0.2, 0.5, and 1.0 m. The overall mean difference (*d*) (Equation (5)) and mean absolute difference (|*d*|) (Equation (6)) were calculated for this purpose [14,23]. The canopy heights were divided into bins at 1.0 m intervals and used as a stratified variable to calculate a series of mean differences for analyzing the potential trends of canopy heights.

$$\overline{d} = \frac{\sum d\_i}{N} = \frac{\sum LH\_i - PH\_i}{N},\tag{5}$$

$$\left|\overline{\overline{d}d}\right| = \frac{\sum \left|d\_i - \overline{d}\right|}{N},\tag{6}$$

where *d* (m) is the overall mean difference, *di* (m) is the canopy height difference between LiDAR and photogrammetry for the *i th* effective cell, *LHi* (m) and *PHi* (m) are the canopy heights of LiDAR and photogrammetry for the *i th* effective cell (i = 1, ... , N), N is the total number of effective cells, and <sup>|</sup>*d*<sup>|</sup> (m) is the mean absolute difference.

The original and interpolated CHMs obtained from the LiDAR data set (response) were linearly regressed and analyzed with the corresponding CHMs of the five photogrammetry data sets (predictor). Three widely used statistical criteria, namely, the coefficient of determination (*R*2) (7), root mean square error (*RMSE*) (8), and relative *RMSE* (*rRMSE*) (9), were used to assess the model accuracy [41]. The *RMSE* combines the mean error and error variance to provide a robust measure of overall accuracy [62].

$$R^2 = 1 - \frac{\sum \left( LH\_i - L\hat{H}\_i \right)^2}{\sum \left( LH\_i - \overline{LH} \right)^2},\tag{7}$$

$$RMSE = \sqrt{\overline{\varepsilon}^2 + \sigma\_{\varepsilon}^2} \tag{8}$$

$$rRMSE = \frac{RMSE}{\overline{LH}},\tag{9}$$

where *LHi* and *LH*ˆ *<sup>i</sup>* are the reference and estimated forest canopy heights (m) for the *i*th effective cell (i = 1, ... , N), N is the total number of the reference canopy heights, *LH* is the mean of the reference canopy heights, *e* is the mean error calculated by *e* = - *ei*/N = - *LHi* <sup>−</sup> *LH*<sup>ˆ</sup> *<sup>i</sup>* /N, *ei* is the *i*th error and σ2 *<sup>e</sup>* is the error variance calculated by σ<sup>2</sup> *<sup>e</sup>* = -(*ei* − *e*) 2 /(N − 1).

#### **3. Results**

#### *3.1. CHMs of LiDAR and Photogrammetry*

The original L-DSM of the L-D55 data set with a spatial resolution of 0.1 m had many holes within and between crowns. The interpolation of the L-DSM by the CN algorithm filled all the small holes with areas less than 1.2 m2. The hole-free L-DSM was generated by the ON algorithm. The original and interpolated L-CHMs were calculated based on the corresponding L-DSMs and L-DTM. The small holes of the original L-CHM with a spatial resolution of 0.1 m (Figure 3a) were interpolated by the CN algorithm (Figure 3b). All of the holes of the original L-CHM were filled by the ON algorithm except for the holes representing the water tower and building areas (Figure 3c).

**Figure 3.** Original canopy height models (CHMs) and interpolated CHMs by the constrained neighbor (CN) and ordinary neighbor (ON) algorithms from the LiDAR and photogrammetry data sets with a spatial resolution of 0.1 m (no data areas are white): (**a**) original CHM of L-D55; (**b**) CHM of L-D55 interpolated by the CN algorithm; (**c**) CHM of L-D55 interpolated by the ON algorithm; (**d**) original CHM of P1-D300; (**e**) CHM of P1-D300 interpolated by the CN algorithm; (**f**) CHM of P1-D300 interpolated by the ON algorithm; (**g**) original CHM of P3-D150; (**h**) CHM of P3-D150 interpolated by the CN algorithm; (**i**) CHM of P3-D150 interpolated by the ON algorithm; (**j**) original CHM of P1-D80; (**k**) CHM of P1-D80 interpolated by the CN algorithm; and (**l**) CHM of P1-D80 interpolated by the ON algorithm. The color describes the lowest to highest height from blue to green. The CHMs of P2-D150 and P4-D150 are similar to that of P3-D150 at the same flying height of 150 m and thus are not shown here.

The P-DSMs of the five photogrammetry data sets were processed similarly to those of the LiDAR data set to generate original P-DSMs and those interpolated by the CN and ON algorithms, and the corresponding P-CHMs were generated by subtracting the L-DTM (Figure 3d–l). The original P-CHM of the P1-D300 data set also had some small holes within crowns and large holes between the crowns, and the distributions of the holes were similar to the distribution of the holes in the L-CHM. The original P-CHMs of the P2-D150, P3-D150 and P4-D150 data sets had few small holes within crowns but large holes existed between the crowns. The original P-CHM of the P5-D80 data set had some large holes within the crowns and much larger holes between the crowns than any other data set. The large holes were mainly distributed in the shadow area of tree crowns.

#### *3.2. Spatial Distribution of Point Clouds*

The spatial distribution of the point clouds was described using the ECR, PCH and PCR. The calculated ECRs, PCHs and PCRs of the original and interpolated L-CHMs with spatial resolutions of 0.1, 0.2, 0.5, and 1.0 m were shown in Figure 4. The ECRs of the original and interpolated L-CHMs increased as the spatial resolution changed from 0.1 to 1.0 m. The ECRs of the original L-CHM with a spatial resolution of 1.0 m were still less than 1.0, which indicated that some areas had sparsely distributed points due to weak ability of backscattering. The CHMs interpolated by the ON algorithm achieved the ideal ECR of 1.0 for any spatial resolution, although the method introduced great spatial variations and decreased the accuracy of the estimated canopy heights.

**Figure 4.** Effective cell ratios (ECRs), point cloud homogeneities (PCHs), and point cloud redundancies (PCRs) of the original CHMs and interpolated CHMs by the constrained neighbor (CN) algorithm from LiDAR and photogrammetry data sets: (**a**) ECRs of the original CHMs; (**b**) ECRs of the CHMs interpolated by the CN algorithm; (**c**) PCHs of the original CHMs; (**d**) PCHs of the CHMs interpolated by the CN algorithm; (**e**) PCRs of the original CHMs; and (**f**) PCRs of the CHMs interpolated by the CN algorithm.

The ECRs, PCHs and PCRs of P-CHMs are shown in Figure 4. All the ECRs of the original P-CHMs were greater than those of the L-CHM at spatial resolutions less than 1.0 m, which indicated that the photogrammetry data sets had fewer holes than the LiDAR data set. The ECRs of the original and interpolated P-CHMs also increased as the spatial resolution decreased.

The PCH of the original L-CHM was less than that of the original P-CHMs of most photogrammetry data sets except for P-D80 with a spatial resolution of 0.1 m, which showed that the dense point cloud of photogrammetry was more evenly distributed than the LiDAR point cloud. The low PCH of the original P-CHM of the P-D80 data set was caused by many large holes due to the failed reconstruction of the dense point cloud. The PCHs of the original L- and P-CHMs decreased continuously as the spatial resolution decreased. The PCRs of the original P-CHMs were all higher than those of the L-CHM, indicating that there were more redundant points in the photogrammetry data sets.

The difference between the ECRs of the original L-CHM and interpolated L-CHM by CN at the spatial resolution of 0.1 m was the largest (0.19) compared with those at other spatial resolutions. The ECR pair differences decreased as the spatial resolution decreased. The optimal spatial resolution of L-CHM was 0.2 m if the threshold of the ECR pair difference was set as 0.10. The optimal spatial resolution of P-CHMs was 0.1 m when the same threshold of the ECR pair difference as that of L-CHM was used.

#### *3.3. Di*ff*erences between L-CHMs and P-CHMs*

The differences between the original and interpolated L- and P-CHMs were calculated on a cell-by-cell basis. The void cells of the holes might not have corresponding cells among different CHMs. These cells were ignored if there was a void cell in either of the CHM pairs used to calculate the differences. The mean differences between the original L-CHM and original P-CHMs with a spatial resolution of 0.1 m varied from −0.1 to −0.5 m, while the corresponding mean absolute differences varied from 0.9 to 1.1 m (Table A3). This result indicated that the original L-CHM was lower overall than the original P-CHMs and large positive and negative differences occurred within the original CHMs of different data sets. The mean differences between the original L-CHM and P-CHMs increased from negative to positive as the spatial resolution became coarser because low values within cells were ignored when CHMs with a coarse spatial resolution were generated using the maximum algorithm. The mean absolute differences between the original L-CHM and original P-CHMs decreased as the spatial resolution increased, which indicated that the spatial variation in canopy heights was suppressed at a coarse spatial resolution.

The mean differences and mean absolute differences between the LiDAR and photogrammetry data sets slightly increased when the CHMs were interpolated by the CN algorithm compared with the original CHMs. The mean differences and mean absolute differences further increased using the CHMs that were interpolated by the ON algorithm compared with the CHMs interpolated by the CN algorithm.

The stratified mean differences between the original L-CHM and P-CHMs within the 1.0 m bin of the canopy heights were calculated to analyze the spatial variations among canopy heights (Figure 5). Three mark heights of mean differences were observed along the canopy heights, including the minimum mark height, negative mark height and positive mark height. The minimum mark height of the mean difference between CHMs at a spatial resolution of 0.1 m was the 3.0 m bin of the canopy height corresponding to a minimum mean difference of 3.3 m. All the mean differences between CHMs with a spatial resolution of 0.1 m were negative when the canopy height was less than the 21.0 m bin and were positive when the canopy height was greater than the 26.0 m bin, which were referred to as negative and positive mark heights, respectively. The mean differences increased from negative to positive as the canopy height increased above the 3.0 m bin. This result indicated that photogrammetry tends to overestimate the low canopy heights and underestimate the high canopy heights compared to those obtained from LiDAR data.

**Figure 5.** Stratified mean differences between the original LiDAR CHM and photogrammetry CHMs within 1.0 m bins of the canopy heights: (**a**) mean differences between CHMs at a spatial resolution of 0.1 m; (**b**) mean differences between CHMs at a spatial resolution of 0.2 m; (**c**) mean differences between CHMs at a spatial resolution of 0.5 m; and (**d**) mean differences between CHMs at a spatial resolution of 1.0 m.

The minimum mean differences increased from −3.3 to −0.6 m as the spatial resolution changed from 0.1 to 1.0 m. The negative and positive mark heights both shifted from high bins to low bins as the spatial resolution became coarser.

#### *3.4. Correlation between L-CHMs and P-CHMs*

The original and interpolated L-CHMs (responses) were regressed with the corresponding original and interpolated P-CHMs (predictors) (see Figure 6). High correlations were observed between the LiDAR and photogrammetry data sets at different spatial resolutions as shown by the *R*2, which varied from 0.87 to 0.98. The original L-CHMs were better correlated with P-CHMs at a coarse spatial resolution than a fine spatial resolution. The *RMSE* between the original L-CHMs and P-CHMs decreased as the spatial resolution changed from 0.1 to 1.0 m.

**Figure 6.** Coefficient of determination (*R*2), *RMSE*, and relative *RMSE* between the CHMs from LiDAR and photogrammetry data sets: (**a**) *R*<sup>2</sup> between the original CHMs; (**b**) *R*<sup>2</sup> between the CHMs interpolated by the CN algorithm; (**c**) *R*<sup>2</sup> between the CHMs interpolated by the ON algorithm; (**d**) *RMSE* between the original CHMs; (**e**) *RMSE* between the CHMs interpolated by the CN algorithm; (**f**) *RMSE* between the CHMs interpolated by the ON algorithm; (**g**) relative *RMSE* between the original CHMs; (**h**) relative *RMSE* between the CHMs interpolated by the CN algorithm; and (**i**) relative *RMSE* between the CHMs interpolated by the ON algorithm.

The L-CHMs and P-CHMs interpolated by the CN algorithm had slight effects on the *R*<sup>2</sup> and root mean square errors (RMSEs) compared with those between the original L-CHM and P-CHMs. The RMSEs between the interpolated L-CHM and P-CHMs by the ON algorithm were higher than those between the interpolated CHMs by the CN algorithm, which was expected. The relative RMSEs between the interpolated L-CHM and P-CHMs by the ON algorithm increased by approximately 3.1% compared with those between the original L-CHM and P-CHMs with a spatial resolution of 0.1 m. This result indicated that the ON interpolation algorithm would introduce obvious spatial variation into the CHM at fine spatial resolution. The relative RMSE decreased as the spatial resolution became coarser.

#### **4. Discussion**

#### *4.1. Spatial Distribution of Point Clouds*

The spatial distribution of point clouds is usually uneven across a survey area due to various conditions, such as the LiDAR or photogrammetry system configuration, flying height, crown shadows, or terrain conditions [43–45,63]. At different locations, some points are closely clustered while other points are sparsely distributed. The mean point density can describe only the overall characteristics of point clouds and does not reflect the features of uneven distributions. The ECR, PCH and PCR are used to describe the spatial distribution of point clouds based on regular grids.

The spatial resolution determines the number of points located within each cell of the grid and affects the values of these spatial descriptors. The coverage, homogeneity and redundancy of point clouds are dependent on the spatial resolution of the grid. The coverage of the point clouds will increase at a coarse spatial resolution if the ECR is less than 1.0, the homogeneity of the point clouds will decrease at a coarse spatial resolution, while the redundancy of the point clouds will increase at a coarse spatial resolution.

The relative correlation of the spatial distribution of different data sets might remain constant. For example, P3-D150 had the best coverage, L-D55 had the worst coverage, P1-D300 was the most evenly distributed data set, P5-D80 had the worst homogeneity, L-D55 had the lowest redundancy, and P5-D80 had the most redundant points, and these observations were all based on CHMs with the same spatial resolution (Figure 4). The ideal data set should have the maximum coverage, best homogeneity and lowest redundancy. No data set fulfilled these criteria in this study case.

#### *4.2. E*ff*ects of the CN and ON Algorithms on the CHM*

The void cells of the CHM could be filled by a neighbor interpolation algorithm. The ON algorithm continuously interpolates void cells from the outside cells to inside cells. On the other hand, the CN algorithm interpolates only void cells that meet the criteria of neighboring cells. For example, if the threshold of effective neighboring cells (*Q*) is 5, the void cells along the edges of a large hole will not be filled while the void cells at the corners will be conditionally filled by the CN algorithm.

The interpolated void cells of the CHM by the CN algorithm varied as the spatial resolution changed and produced different ECR values (Figure 4). The ON algorithm interpolated all void cells and obtained the same ECR at all spatial resolutions. If the dimension of the spatial resolution was much smaller than the distance of the point cloud, then the number of the effective neighboring cells of void cells was smaller than the threshold of effective neighboring cells (*Q*) and no void cells were interpolated by the CN algorithm. If the dimension of the spatial resolution was larger than the size of the maximum hole, then the ECR was equal to 1.0. The CN algorithm had no effect on the ECR of the CHM according to these two mentioned cases.

The CN algorithm had a weak effect on the mean absolute differences between the L-CHM and P-CHMs at the spatial resolution of 0.1 m. The mean absolute differences between the original L-CHM and original P-CHMs at the spatial resolution of 0.1 increased by 0.2 m compared with that between the L-CHM and P-CHMs interpolated by the ON algorithm (Table A3). The effects of the CN and ON algorithms on the mean absolute differences both became weaker at a coarser spatial resolution.

The relative RMSEs between the L-CHMs and P-CHMs interpolated by the ON algorithm were larger than those between the CHMs interpolated by the CN algorithm at all spatial resolutions. The CN algorithm had the ability to control the spatial variation introduced by the ON algorithm. Thus, the CN algorithm is recommended for interpolating the void cells of a CHM rather than the ON algorithm.

#### *4.3. Distinguishing Small and Large Holes within a CHM*

The holes within a CHM with a specified spatial resolution consist of one or several continuous void cells. Small holes in an L-CHM are usually caused by the fluctuation of the UAV platform, mutual obscuration of crowns, gaps within and between crowns, lost echoes of laser pulses, etc. Large holes in an L-CHM are mainly affected by the large incident angles on crowns and the loss of many backscattering signals. The influencing factors of small and large holes of the P-CHM include the reconstruction algorithm of dense point clouds, the quality of images, illumination conditions, texture of crowns, pits of crowns, gap between crowns, and swinging of crowns due to wind. Large holes in

the P-CHM are mostly due to weak texture and poor lighting (extremely bright or dim). The light conditions of different photogrammetry data sets are very different and will cause different shadows to affect the spatial distribution of the reconstructed point cloud.

The coverage of holes is the ratio of void cells to all cells expressed as 1-ECR. The original CHM with a fine spatial resolution might include small and large holes, while the CHM interpolated by the CN algorithm would include only large holes. For example, the area covered by holes of the original L-CHM with a spatial resolution of 0.1 m was 38% of the whole area. The large holes (≥0.12 m2) of the interpolated L-CHM occupied 19% of the whole area. This result indicated that 50% of the holes were large holes. The percent of large holes (≥0.48 m2) decreased to 40% of the holes when the spatial resolution of the original L-CHM was 0.2 m. The minimum size of large holes further increased, and the percent of large holes further decreased when the spatial resolution of the L-CHM became coarser.

#### *4.4. Optimal Spatial Resolution of the CHM*

The spatial resolution of a CHM is usually determined by considering the mean point cloud density. Many void cells will exist within the CHM at this spatial resolution if the point cloud is unevenly distributed across the survey area. Small holes of the CHM can be interpolated by the CN algorithm. The interpolated cells reflect the potential space that would become effective cells at a coarse spatial resolution. The optimal spatial resolution may be determined by the potential space, which is referred to as the potential space criteria. The potential space is calculated as the difference between the ECRs of the original CHM and the CHM interpolated by the CN algorithm.

The optimal spatial resolution can also be determined by the ECR of the CHM interpolated by the CN algorithm, which is referred to as the ECR criteria. If the ECR approaches 1 and is greater than the specified threshold (for example, 0.90), then the corresponding spatial resolution will be regarded as optimal. The ECR criteria are easily affected by very large holes, whose sizes are far greater than the mean space of the point cloud. Additional criteria can be used to select the optimal spatial resolution. We recommend potential space criteria that will not be affected by large holes. Moreover, the optimal spatial resolution can be simply determined by the ECR pair differences among predefined spatial resolutions, which should be calculated by automatic methods in future studies.

#### *4.5. E*ff*ects of Flying Configuration on the CHM*

Five photogrammetry data sets at three different flying heights were used in this study. P1-D300 had the highest flying height and a highly efficient capability of data acquisition when compared to the other data sets collected at lower flying heights. The optimal spatial resolution of P1-D300 was 0.1 m if the threshold of the difference between the ECRs of the original CHM and the CHM interpolated by the CN algorithm was 0.1. Therefore, an optimal flying height of approximately 300 m was recommended to generate a CHM at a spatial resolution of 0.1 m using the similar photogrammetry system in this study.

The P5-D80 data set had the lowest flying height among the photogrammetry data sets, although the ECRs of the original P-CHM of P5-D80 were smaller than those of the other photogrammetry data sets. Very large holes exist within the bright crown area in the original P-CHM of the P5-D80 data set. This finding indicated that some tree crowns failed to be reconstructed using images with very high spatial resolution.

Different shadows occurred in the photogrammetry data sets with different light conditions, and these differences would affect the distribution patterns of void cells in the P-CHMs. The overall effects of light conditions on the reconstructed point cloud and the P-CHMs were weaker than the effects of the flying height in this study. The ECRs of the original P-CHMs with a spatial resolution of 0.1 m had slight differences under different light conditions at the same flying height of 150 m. The differences among the ECRs of the P-CHMs with a spatial resolution of 0.1 m at flying heights of 80, 150 and 300 m were greater than those among the ECRs of the P-CHMs at the same flying height of 150 m.

The mean differences between the combined pairs of photogrammetry data sets were smaller and fluctuated from −0.5 m to 0.5 m when the canopy heights were less than 30 m. Obviously large variations occurred when the canopy heights were greater than 30 m. This variation might be due to the swinging of leaves and twigs of tall trees affected by wind. The mean differences between the combined pairs of the P2-D150, P3-D150 and P4-D150 data sets at the same flying heights were smaller than those of other pairs at different flying heights for all canopy height bins.

The mean differences between the P2-D150 and P4-D150 data sets with different image overlaps varied from 0 m to 0.1 m, which meant that image overlap (>68%) had a minor effect on the estimated canopy height. However, the point density of the photogrammetry data set was substantially affected by image overlap (Figure A1). The mean differences between the P2-D150 and P3-D150 data sets with different flying speeds varied from −0.2 to 0 m, indicating that canopy height was only slightly affected by flying speed (<8ms<sup>−</sup>1).

#### *4.6. E*ff*ects of Gaps within and between Crowns on the CHM*

The CHMs obtained from LiDAR and photogrammetry data sets for the given spatial resolution were affected by the gaps within and between crowns [14,23,27]. Small gaps within crowns could be depicted in the L-CHM with fine spatial resolution, while the P-CHM tended to ignore such small gaps [23]. A number of deep gaps were observed between the tall tree crowns, which were much darker than the surrounding crowns under the variable illumination conditions. It was difficult to reconstruct the spatial structure in these deep gaps [14,38].

The spatial resolution will affect the height values of gaps within and between crowns. The cells of a CHM are often calculated as the highest value if there are numerous points within the corresponding cell. This condition will suppress the low values of gaps and promote high values of crowns, thus reducing the spatial variation in the canopy heights. LiDAR can provide more detailed structural information than photogrammetry on a fine scale and more details will be lost on a coarse scale. The CHMs between LiDAR and photogrammetry tend to have more consistency as the spatial resolution decreases.

#### **5. Conclusions**

This study compared the canopy heights obtained from UAV LiDAR and UAV photogrammetry and interpolated by two spatial interpolation algorithms CN and ON. The comparison was conducted based on three proposed spatial distribution descriptors of point clouds: the ECR, PCH and PCR quantifying the unevenness, homogeneity and redundancy characteristics, respectively, of point clouds in the grid at a given spatial resolution. The stratified mean differences revealed that there existed an inherent trend between the estimated canopy heights from LiDAR and photogrammetry, which changed from negative to positive as the canopy heights increased. The LiDAR CHM strongly correlated with the photogrammetry CHM. More importantly, the CN algorithm had the ability to distinguish small and large holes and determine the optimal spatial resolution according to the ECR pair differences, while the ON algorithm did not have this ability. Compared with the ON algorithm, the CN algorithm apparently reduced the spatial variation in the CHM, led to smaller RMSE values, and could be recommended to obtain reliable CHM values.

Some large holes in the CHMs interpolated by the CN interpolation algorithm still occurred at very high point densities, which were typically distributed around the deep gaps between tall tree crowns. Precisely measuring such deep gaps is quite challenging. Photogrammetry tends to ignore gaps within the crowns, overestimate low canopy heights and underestimate high canopy heights; these issues need to be considered while canopy cover, canopy closure, and forest stand height are modelled. Overall, this article provides an operational method for the spatial assessment of point clouds and suggests that the differences between LiDAR and photogrammetry derived CHMs should be considered when forest parameters are estimated. In particular, further study is necessary to enhance understanding the quality measures of the point cloud spatial distribution, optimal spatial

resolution, small and large holes within CHMs, gaps within and between crowns, and their effects on estimation of forest parameters using additional data sets.

**Author Contributions:** Conceptualization, Q.L. and L.F.; data curation, Q.C., P.L., P.H., M.L. and M.W.; formal analysis, M.L., M.W. and G.D.; funding acquisition, Q.L. and L.F.; methodology, Q.L. and L.F.; supervision, G.W.; visualization, Q.C., P.L. and P.H.; writing—original draft, Q.L.; writing—review and editing, G.W. and R.P.S. All authors contributed to interpreting results and the improvement of the article. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was partially funded by the National Key R&D Program of China under grant number 2017YFD0600904, the Central Public-interest Scientific Institution Basal Research Fund (Grant No. CAFYBB2019QD003 and CAFYBB2016SZ003).

**Acknowledgments:** We thank Xincheng Shi, Denglong Ha, Yixi Ju, Qiuyan Wang, et al. for assisting field work. We thank the support of the Jigongshan National Nature Reserve, Xinyang, Henan Province, China and the Xinyang Normal University, Xinyang, Henan Province, China for experiments. We thank the anonymous reviewers for their constructive comments and recommendations, which we used to significantly improve our article.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A**

The statistics of the plot parameters are shown in Table A1.

**Table A1.** The statistics of the plot parameters.


The acquisition information of the LiDAR and photogrammetry data sets is shown in Table A2.


**Table A2.** LiDAR and photogrammetry data sets.

<sup>1</sup> The flying height is the relative height above the take-off point (194 m above sea level) of the UAV, which is near the highest point of the study site (198 m above sea level); <sup>2</sup> LiDAR flight with a height of 50 m above the terrain other than the take-off point of the UAV.

The calculated mean differences and mean absolute differences between the L-CHMs and P-CHMs are listed in Table A3.

The point density distributions of the LiDAR and photogrammetry data sets are shown in Figure A1.


**Table A3.** The mean differences and mean absolute differences between the CHMs of LiDAR and photogrammetry data sets.

<sup>1</sup> CN means the constrained neighbor interpolation algorithm; <sup>2</sup> ON means the ordinary neighbor interpolation algorithm.

**Figure A1.** Point density of the 1 m resolution grid (no data areas are white; the grid is stretched with a percent clip of 0.5): (**a**) L-D55 point cloud with a mean of 168 points m<sup>−</sup>2; (**b**) P1-D300 point cloud with a mean of 303 points m<sup>−</sup>2; (**c**) P2-D150 point cloud with a mean of 1319 points m<sup>−</sup>2; (**d**) P3-D150 point cloud with a mean of 1191 points m<sup>−</sup>2; (**e**) P4-D150 point cloud with a mean of 1050 points m<sup>−</sup>2; (**f**) P5-D80 point cloud with a mean of 2529 points m<sup>−</sup>2.

## **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
