Extracting Canopy Closure by the CHM-Based and SHP-Based Methods with a Hemispherical FOV from UAV-LiDAR Data in a Poplar Plantation

: Canopy closure (CC), a useful biophysical parameter for forest structure, is an important indicator of forest resource and biodiversity. Light Detection and Ranging (LiDAR) data has been widely studied recently for forest ecosystems to obtain the three-dimensional (3D) structure of the forests. The components of the Unmanned Aerial Vehicle LiDAR (UAV-LiDAR) are similar to those of the airborne LiDAR, but with higher pulse density, which reveals more detailed vertical structures. Hemispherical photography (HP) had proven to be an effective method for estimating CC, but it was still time-consuming and limited in large forests. Thus, we used UAV-LiDAR data with a canopy-height-model-based (CHM-based) method and a synthetic-hemispherical-photography-based (SHP-based) method to extract CC from a pure poplar plantation in this study. The performance of the CC extraction methods based on an angular viewpoint was validated by the results of HP. The results showed that the CHM-based method had a high accuracy in a 45 ◦ zenith angle range with a 0.5 m pixel size and a larger radius (i.e., k = 2; R 2 = 0.751, RMSE = 0.053), and the accuracy declined rapidly in zenith angles of 60 ◦ and 75 ◦ ( R 2 = 0.707, 0.490; RMSE = 0.053, 0.066). In addition, the CHM-based method showed an underestimate for leaf-off deciduous trees with low CC. The SHP-based method also had a high accuracy in a 45 ◦ zenith angle range, and its accuracy was stable in three zenith angle ranges ( R 2 : 0.688, 0.674, 0.601 and RMSE = 0.059, 0.056, 0.058 for a 45 ◦ , 60 ◦ and 75 ◦ zenith angle range, respectively). There was a similar trend of CC change in HP and SHP results with the zenith angle range increase, but there was no signiﬁcant change with the zenith angle range increase in the CHM-based method, which revealed that it was insensitive to the changes of angular CC compared to the SHP-based method. However, the accuracy of both methods showed differences in plantations with different ages, which had a slight underestimate for 8-year-old plantations and an overestimate for plantations with 17 and 20 years. Our research provided a reference for CC estimation from a point-based angular viewpoint and for monitoring the understory light conditions of plantations.


Introduction
Plantations, as a major component of forest ecosystems, play an important role in biodiversity conservation, climate change and economic development [1,2]. Poplar plantations provide a large amount of timber and forest products in China [3]. So, it is essential to investigate the resources of poplar plantation timely and accurately [4].
The canopy structure affects different species in forest ecosystems by affecting light levels [5]. Canopy closure (CC), as a useful biophysical parameter for forest canopy structure, is widely used as an index for carbon fluxes, forest production, biodiversity and ecosystem functions by affecting energy transmission and microclimate in forest the DBH estimation (R 2 = 0.98) [56]. Liu et al. (2018) evaluated the capability of the UAV-LiDAR system for estimating forest structural attributes and analyzed the effects of point cloud densities. They found that it was robust for estimating forest attributes when the point cloud density is higher than 16 pts·m −2 [57]. Wu et al. (2019) used UAV-LiDAR data to estimate the canopy cover of a ginkgo-planted forest with three methods and found that the CHM-based method with a 0.5 m resolution had a high accuracy (R 2 = 0.92) [54]. Thus, the UAV-LiDAR data with high density point cloud are highly advantageous in describing the forest canopy structure, which has a high potential for an accurate CC extraction.
Although there have been many studies using LiDAR data for CC extraction, few of them evaluate the method based on CHM and SHP comprehensively, such as exploring the effects of different FOVs on different methods and explaining the scope of applicability of different methods. Moreover, previous studies rarely focused on the CC estimation in forest plantations based on UAV-LiDAR data. Therefore, this study aims to evaluate an efficient method to extract CC directly from UAV-LiDAR data. To achieve this goal, the specific objectives are: (1) to estimate CC in three zenith angle ranges, of 45 • , 60 • and 75 • , based on CHM and SHP data from UAV-LiDAR data; (2) to validate the results of the two CC extraction methods by the classification results from HP data; (3) to explore the effect of different resolutions and stand ages on CHM-based and SHP-based methods.

Study Area
The study area is located in Dongtai City, Jiangsu Province ( Figure 1). The area is flat, with a range of elevation from 11 to 14 m. It falls in the climate region of the transition zone between the subtropical and warm temperate area [58]. The annual mean temperature is 15.4 • C, the annual rainfall is 1494.0 mm and the average relative humidity is 76.0% [58]. The main soil type is desalted meadow soil, and the soil texture is sandy loam, alkaline (pH = 8.2). The main tree species include poplar (Populus deltoids), fir (Metasequoia glyptostroboides), ginkgo trees (Ginkgo biloba L.). Among all the three species, poplar has the largest plantation area, as well as with a variety of stand ages.  [52][53][54][55]. Brede et al. (2017) compared terrestrial LiDAR and UAV-LiDAR for estimating forest canopy heights and diameter at breast heights (DBH), and their results showed a strong correlation of two LiDAR data in the DBH estimation ( 2 = 0.98) [56]. Liu et al. (2018) evaluated the capability of the UAV-LiDAR system for estimating forest structural attributes and analyzed the effects of point cloud densities. They found that it was robust for estimating forest attributes when the point cloud density is higher than 16 pts·m −2 [57]. Wu et al. (2019) used UAV-LiDAR data to estimate the canopy cover of a ginkgo-planted forest with three methods and found that the CHM-based method with a 0.5 m resolution had a high accuracy ( 2 = 0.92) [54]. Thus, the UAV-LiDAR data with high density point cloud are highly advantageous in describing the forest canopy structure, which has a high potential for an accurate CC extraction. Although there have been many studies using LiDAR data for CC extraction, few of them evaluate the method based on CHM and SHP comprehensively, such as exploring the effects of different FOVs on different methods and explaining the scope of applicability of different methods. Moreover, previous studies rarely focused on the CC estimation in forest plantations based on UAV-LiDAR data. Therefore, this study aims to evaluate an efficient method to extract CC directly from UAV-LiDAR data. To achieve this goal, the specific objectives are: 1) to estimate CC in three zenith angle ranges, of 45°, 60° and 75°, based on CHM and SHP data from UAV-LiDAR data; 2) to validate the results of the two CC extraction methods by the classification results from HP data; 3) to explore the effect of different resolutions and stand ages on CHM-based and SHP-based methods.

Study Area
The study area is located in Dongtai City, Jiangsu Province ( Figure 1). The area is flat, with a range of elevation from 11 to 14 m. It falls in the climate region of the transition zone between the subtropical and warm temperate area [58]. The annual mean temperature is 15.4 °C , the annual rainfall is 1494.0 mm and the average relative humidity is 76.0% [58]. The main soil type is desalted meadow soil, and the soil texture is sandy loam, alkaline (pH = 8.2). The main tree species include poplar (Populus deltoids), fir (Metasequoia glyptostroboides), ginkgo trees (Ginkgo biloba L.). Among all the three species, poplar has the largest plantation area, as well as with a variety of stand ages.

Field Design and Field Data Collection
The field work was conducted in May 2021. According to the difference in the size of poplar sublots, two different scale plots were set (60 × 60 m, 30 × 30 m; Figure 1). The plots contain varying stand ages, planting spacing and canopy structural characteristics (Table 1). A total of 29 square sample plots were measured, including 20 small plots (i.e., the total area is 900 m 2 ) and 9 large plots (i.e., the total area is 3600 m 2 ). Each of the plots was divided into several small blocks by a 10 × 10 m grid, and the center of each small block was the location for taking HP (Figure 1). In total, 570 HP were collected from 29 plots. All the HP were taken horizontally using a Canon M50 with a LAOWA CF 4 mm F2.8 circular fisheye lens. The camera was fixed on a tripod and each photo was taken 1.4 m above the ground. The HP was taken before sunrise, after sunrise or during daytime at an overcast sky to ensure a low light condition. All photos were taken in auto-exposure mode, and the shooting location was avoided too close to the trees. The shooting location was recorded below the tripod using a HUACE T10 real-time-kinematic (RTK) equipment (centimeter accuracy), and corrected with high precision real-time signals received from continuously operating reference stations (CORS). The recorded location of the global positioning system (GPS) points was regarded as the shooting location of HP. The UAV-LiDAR data were collected via a Velodyne VLP-16 LiDAR sensor carried by a six-rotor DJI M600 PRO UAV. The parameters of the LiDAR sensor are as follows: wavelength 903 nm, vertical scanning angle ±15 • , horizontal scanning angle 360 • (a range of ±70 • was reserved), pulse frequency 30 kHz. The flight parameters are the following: flight altitude 70 m, flight speed 3.6 m·s −1 , flight interval 60 m. The point density of different sample plots was not the same because of the difference of flight strip numbers (Table 1). In addition, the forest planting spacing, the average tree height and the average branch height of each plot were collected for the subsequent data processing and analysis (Table 1).

Methodology
After the pre-processing of UAV-LiDAR data, including georeferencing, strip alignment, merging, de-noising and the interpolation of point clouds with different types, CHM and a normalized point cloud were generated. Then the shooting locations of HP data were used to delineate the FOV in CHM and normalized point cloud data, and the CC of different plots was extracted by height based on the FOV-CHM. SHP data were generated using a 3D polar coordinate conversion method to extract CC from all the normalized point cloud data directly. Finally, the CC extracted by the UAV-LiDAR data of different resolutions was validated by the results from HP data in different FOVs. The age effects of poplar plantations on CC extraction from UAV-LiDAR point cloud data were also explored in this study. In addition, a new classification model was developed to extract CC from HP data in this study as the validation data to compare the CC extraction result from UAV-LiDAR point cloud data using the CHM-based and SHP-based methods. An overview of the flowchart for CC extraction is shown in Figure 2.

Preprocessing of UAV-LiDAR Point Cloud Data
The raw UAV-LiDAR point cloud coordinate was transformed based on UAV station GPS data, inertial measurement unit (IMU) data and base station data. Each strip was adjusted by a surface matching method based on the iterative closest point (ICP) algorithm [59]. After merging all the strips, the noised points were removed depending on the maximum distance from the point to its neighboring points in the LiDAR360 software. The ground points were classified by an improved progressive triangulated irregular network (TIN) densification filtering algorithm [60]. Then, the ground points were interpolated to generate the digital elevation model (DEM) by the TIN algorithm, and the digital surface model (DSM) was generated by the same interpolation method. Then, we subtracted the DEM from the DSM to generate CHM of different resolutions. We also normalized the height of the denoising point cloud based on the DEM to generate normalized point cloud data in preparation for the subsequent generation of SHP data.

Preprocessing of HP and FOV Delineation
Each photo was checked for quality, and the ones with overexposure and underexposure were discarded. We cropped all the photos and reduced the effect of image local highlights in the Adobe Photoshop CC 2019 software. The resolution of the processed photo was 3124 × 3124 pixels. The directed upward HP had a nearly polar projection, so the azimuth angle and zenith angle ranges of the photo could be divided according to the center of the photo [20]. Considering the difference of CC extraction results in different zenith angle ranges [42,61], we calculated the CC from three zenith angle ranges, of 45°, 60° and 75°. The different zenith angle ranges represented the different FOVs and radius in the image (Equation 1).
where is the radius of the entire FOV in the photograph, and 90° represents the theoretical maximum zenith angle.
is the zenith angle, which in this study is 45°, 60° and 75°. is the new radius with different FOVs in HP.

CC Extraction by the CHM-based Method with a Hemispherical FOV
Compared to the estimation of canopy cover using CHM, it is a challenge to estimate the CC using CHM by a hemispherical FOV. In this study, the delineated FOV was used as the analysis window to calculate CC based on the CHM data generated by the prepro-

Preprocessing of UAV-LiDAR Point Cloud Data
The raw UAV-LiDAR point cloud coordinate was transformed based on UAV station GPS data, inertial measurement unit (IMU) data and base station data. Each strip was adjusted by a surface matching method based on the iterative closest point (ICP) algorithm [59]. After merging all the strips, the noised points were removed depending on the maximum distance from the point to its neighboring points in the LiDAR360 software. The ground points were classified by an improved progressive triangulated irregular network (TIN) densification filtering algorithm [60]. Then, the ground points were interpolated to generate the digital elevation model (DEM) by the TIN algorithm, and the digital surface model (DSM) was generated by the same interpolation method. Then, we subtracted the DEM from the DSM to generate CHM of different resolutions. We also normalized the height of the denoising point cloud based on the DEM to generate normalized point cloud data in preparation for the subsequent generation of SHP data.

Preprocessing of HP and FOV Delineation
Each photo was checked for quality, and the ones with overexposure and underexposure were discarded. We cropped all the photos and reduced the effect of image local highlights in the Adobe Photoshop CC 2019 software. The resolution of the processed photo was 3124 × 3124 pixels. The directed upward HP had a nearly polar projection, so the azimuth angle and zenith angle ranges of the photo could be divided according to the center of the photo [20]. Considering the difference of CC extraction results in different zenith angle ranges [42,61], we calculated the CC from three zenith angle ranges, of 45 • , 60 • and 75 • . The different zenith angle ranges represented the different FOVs and radius in the image (Equation (1)).
where R is the radius of the entire FOV in the photograph, and 90 • represents the theoretical maximum zenith angle. θ is the zenith angle, which in this study is 45 • , 60 • and 75 • . r is the new radius with different FOVs in HP.

CC Extraction by the CHM-Based Method with a Hemispherical FOV
Compared to the estimation of canopy cover using CHM, it is a challenge to estimate the CC using CHM by a hemispherical FOV. In this study, the delineated FOV was used as the analysis window to calculate CC based on the CHM data generated by the preprocessed point cloud data. Because the FOV for field measured CC has the shape of an inverted cone, the horizontal radius of the analysis window should be generated based on the tree height [38]. However, the analysis window radius in the CHM-based method was not clear in previous studies [35,62]. Parent et al. (2014) put it forward using the height multiplied by the tangent of the zenith as the CHM radius range [37], but we found the radius of this model too large to fit the actual results of HP data in our study. Thus, we improved this FOV-delineated model based on the average height of the sample plot and the zenith angle of interest (Equation (2)).
where r is the radius of the CHM analysis window; k is the distance coefficient to adjust the radius (i.e., the radius that with k is 2 and 3 in the three zenith angle ranges); h is the average height of a sample plot, and θ is the zenith angle of interest. CC extraction from CHM images is a method based on the ratio of extracted canopy pixels to total pixels. The separation of canopy pixels and non-canopy pixels usually used a fixed height threshold (i.e., 2 or 3 m [28]). By using this threshold, the point cloud of the low plants had little effect on the generation of CHM data, so a CHM image with a height greater than this threshold is assumed as canopy pixels ( Figure 3). However, the spatial resolution of CHM images might have influences on CC extraction. Thus, it is worth to explore whether there is an optimum spatial resolution to extract CC with the CHM-based model. Thus, we generated CHM with three spatial resolution, of 0.5, 2 and 5 m, for a sensitivity analysis. The CHM images were generated by a preprocessed point cloud data, and the different spatial resolutions were generated through an interpolation algorithm (see the description in Section 2.3.1). All the steps of the CHM-based method were written as a python script (see Supplementary Materials) to enhance the capacity of automatic batch-processing of this method. cessed point cloud data. Because the FOV for field measured CC has the shape of an inverted cone, the horizontal radius of the analysis window should be generated based on the tree height [38]. However, the analysis window radius in the CHM-based method was not clear in previous studies [35,62]. Parent et al. (2014) put it forward using the height multiplied by the tangent of the zenith as the CHM radius range [37], but we found the radius of this model too large to fit the actual results of HP data in our study. Thus, we improved this FOV-delineated model based on the average height of the sample plot and the zenith angle of interest (Equation 2).
where is the radius of the CHM analysis window; is the distance coefficient to adjust the radius (i.e., the radius that with is 2 and 3 in the three zenith angle ranges); ℎ is the average height of a sample plot, and is the zenith angle of interest. CC extraction from CHM images is a method based on the ratio of extracted canopy pixels to total pixels. The separation of canopy pixels and non-canopy pixels usually used a fixed height threshold (i.e., 2 or 3 meters [28]). By using this threshold, the point cloud of the low plants had little effect on the generation of CHM data, so a CHM image with a height greater than this threshold is assumed as canopy pixels ( Figure 3). However, the spatial resolution of CHM images might have influences on CC extraction. Thus, it is worth to explore whether there is an optimum spatial resolution to extract CC with the CHM-based model. Thus, we generated CHM with three spatial resolution, of 0.5, 2 and 5 m, for a sensitivity analysis. The CHM images were generated by a preprocessed point cloud data, and the different spatial resolutions were generated through an interpolation algorithm (see the description in section 2.3.1). All the steps of the CHM-based method were written as a python script (see Supplementary Materials) to enhance the capacity of automatic batch-processing of this method. with a spatial resolution of 0.5 m generated from the UAV-LiDAR point cloud data collected in the same location of (a): the green circle and the purple circle represent the analysis window range when is 2 and 3, respectively; (c) a height threshold of 2 m was applied to (b): the black and grey pixels in (c) are canopy pixels and non-canopy pixels, respectively.

Estimation of CC with a SHP-based Method
In order to extract CC from a point-based angular viewpoint, we transformed the coordinate system of the normalized point cloud data. The traditional Cartesian coordinates (x, y, z) were converted into polar coordinates ( , , ), and all the point cloud data were printed in the polar coordinate system to simulate the photos taken by a circular fisheye lens [44] where is the zenith angle in the polar coordination; is the azimuth in the polar coordination; is the distance between point cloud and the shooting location; , and are (c) a height threshold of 2 m was applied to (b): the black and grey pixels in (c) are canopy pixels and non-canopy pixels, respectively.

Estimation of CC with a SHP-Based Method
In order to extract CC from a point-based angular viewpoint, we transformed the coordinate system of the normalized point cloud data. The traditional Cartesian coordinates (x, y, z) were converted into polar coordinates (r, θ, ϕ), and all the point cloud data were printed in the polar coordinate system to simulate the photos taken by a circular fisheye lens [44] (Equation (3)). where θ is the zenith angle in the polar coordination; ϕ is the azimuth in the polar coordination; r is the distance between point cloud and the shooting location; x, y and z are the 3D properties of the point cloud; x 0 and y 0 are the GPS-coordinates of the shooting location; z 0 is the height of the camera set (a constant of 1.4 m in this study).
The maximum horizontal distance between the transformed point cloud and origin (i.e., shooting location) would have an effect on the results [44]. A large distance would ensure the correctness of the results [44], but the huge data of the point cloud would reduce the efficiency of the computing process. Considering that the maximum zenith angle in our study was 75 • , we used the maximum tree height multiplied by a tangent of 75 • as a theoretical distance range of the point cloud data in this study (i.e., 80 m). To eliminate the influence of the different point cloud density, we used the point cloud thinning algorithm to transform all the data into the same point density (50 pts·m −2 ). Then, the average branch height of each plot was used as the threshold to mask out the point clouds of the main tree trunks. The zenith angle information of the converted point cloud was used to correspond to the results from photos in three FOVs. Finally, we generated the SHP data, which correspond to the HP data ( Figure 4).
Remote Sens. 2021, 13, x. https://doi.org/10.3390/xxxxx www.mdpi.com/journal/remotesensing duce the efficiency of the computing process. Considering that the maximum zenith angle in our study was 75°, we used the maximum tree height multiplied by a tangent of 75° as a theoretical distance range of the point cloud data in this study (i.e., 80 m). To eliminate the influence of the different point cloud density, we used the point cloud thinning algorithm to transform all the data into the same point density (50 pts·m -2 ). Then, the average branch height of each plot was used as the threshold to mask out the point clouds of the main tree trunks. The zenith angle information of the converted point cloud was used to correspond to the results from photos in three FOVs. Finally, we generated the SHP data, which correspond to the HP data ( Figure 4). In order to extract the optimum value of CC with this point cloud density, we need to establish different grid resolutions based on the zenith and azimuth angles of all the point clouds data. Too many grids would make the projected area of the point clouds too small in the overall grid; on the contrary, too few grids would make it too large. So, we selected three different resolutions, of 1, 1.5 and 2, to extract the CC based on the range of zenith angles and azimuth angles. For example, resolution 1 represents the use of 1 degree to split the gird composed of the zenith angle (the maximum zenith angle is 90°) and the azimuth angle (the maximum azimuth angle is 360°), and it has 360 × 90 grids when the zenith angle is 90°. Similarly, resolution 1.5 represents 240 × 60 grids in a zenith angle of 90°, and resolution 2 refers to 180 × 45 grids in a zenith angle of 90°. Then, the CC was calculated by the ratio of the number of grids with existing point clouds to the total number of grids. We used this model instead of the method based on a binarized image to avoid the subsequent image processing. All the steps of the SHP-based method, including the coordinate system transformation and the CC automatic extraction, was written as a python script based on the version of Python 3.8 (see Supplementary Materials).  In order to extract the optimum value of CC with this point cloud density, we need to establish different grid resolutions based on the zenith and azimuth angles of all the point clouds data. Too many grids would make the projected area of the point clouds too small in the overall grid; on the contrary, too few grids would make it too large. So, we selected three different resolutions, of 1, 1.5 and 2, to extract the CC based on the range of zenith angles and azimuth angles. For example, resolution 1 represents the use of 1 degree to split the gird composed of the zenith angle (the maximum zenith angle is 90 • ) and the azimuth angle (the maximum azimuth angle is 360 • ), and it has 360 × 90 grids when the zenith angle is 90 • . Similarly, resolution 1.5 represents 240 × 60 grids in a zenith angle of 90 • , and resolution 2 refers to 180 × 45 grids in a zenith angle of 90 • . Then, the CC was calculated by the ratio of the number of grids with existing point clouds to the total number of grids. We used this model instead of the method based on a binarized image to avoid the subsequent image processing. All the steps of the SHP-based method, including the coordinate system transformation and the CC automatic extraction, was written as a python script based on the version of Python 3.8 (see Supplementary Materials).

A New Semi-Automated Classification Method for CC Extraction from HP
We developed a new semi-automatic classification method based on the morphology image processing and threshold classification. According to the spectral characteristics of the tree trunks in the three bands of RGB in the images, we used an index (R > G and R > B) to extract most of the trunks. Then, the extracted portion was morphologically processed to preserve the main trunk, this process included one mode filtering, two minimum filtering and two maximum filtering. The sky and canopy were then distinguished after the main trunk components were masked out from the images. The blue band was often used as the optimal channel for classifying sky and canopy pixels [35,63]. However, a fixed threshold (G < 200) based on the green band was used to distinguish canopy pixels from the sky pixels due to the influence of bright blue leaves under the influence of light in the blue band. Then, CC was calculated based on the ratio of classified canopy pixels to total pixels (Equation (4)). All the steps were written as a python script (see Supplementary Materials) to enhance the capacity of automatic batch-processing of this method.
where CC means canopy closure; canopy pixel were canopy pixels classified from the new semi-automatic method; total pixel refers to total pixels of the preprocessed HP, including the pixels of sky, tree trunk and tree canopy. Two random photos were selected from each sample plot to perform the accuracy assessment of the semi-automatic classification method of HP. Twenty random points were generated in each photo for accuracy assessment. The confusion matrix F-score was selected as the indicator of accuracy assessment (Equations (5)- (7)).
where r is the recall of the results, p is the precision of the results and F is an index of the overall accuracy. TP is the number of canopy pixels classified correctly by the model. FN is the number of canopy pixels classified incorrectly by the model. FP is the number of other pixels which were extracted as canopy by the model.

Validation and Accuracy Assessment of CC Extraction from UAV-LiDAR Point Cloud Data
The CC extraction results from HP were used to validate the CC extraction from the UAV-LiDAR point cloud data based on the CHM-based method by a hemispherical FOV and the SHP-based method. The accuracy of the CC estimation was evaluated with the coefficient of determination (R 2 ) and root mean squared error (RMSE), and the significance test was also conducted (Equations (8) and (9)). where n is the number of all hemispherical photos, y is the field measured CC, y i is the average of the field measured CC,ŷ i is the estimation value by the model.

The Extracted CC from UAV-LiDAR Data in Different Poplar Plots
Along with the increase of the FOV, the result range of extracted CC from UAV-LiDAR data by both CHM-based and SHP-based methods becomes narrower (Figure 5a), which indicates that the CHM-based and SHP-based methods with smaller FOV are more sensitive to extract CC. CC extraction from HP, the validation data in this study, had a narrower data range than the CC extracted from UAV-LiDAR (Figure 5a). The extracted CC results by HP classification and the SHP-based method increased according to the increase of the zenith angle range from 45 • to 75 • , while the extracted results by the CHM-based method are not sensitive to the FOV change (Figure 5a). The performance of CC extraction based on the three methods was different in the poplar plantations with different stand ages (Figure 5b-d). The range of the extracted value of CC by the three methods was relatively consistent in the 45 • zenith angle range comparing to that in the 60 • and 75 • zenith angle ranges (Figure 5b). Comparing to the CC extraction from the HP classification, the CC extraction by the CHM-based method had an underestimation trend along with the increasing FOV (Figure 5c,d).

of 18
Along with the increase of the FOV, the result range of extracted CC from UAV-Li-DAR data by both CHM-based and SHP-based methods becomes narrower (Figure 5a), which indicates that the CHM-based and SHP-based methods with smaller FOV are more sensitive to extract CC. CC extraction from HP, the validation data in this study, had a narrower data range than the CC extracted from UAV-LiDAR (Figure 5a). The extracted CC results by HP classification and the SHP-based method increased according to the increase of the zenith angle range from 45° to 75°, while the extracted results by the CHMbased method are not sensitive to the FOV change (Figure 5a). The performance of CC extraction based on the three methods was different in the poplar plantations with different stand ages (Figure 5b-d). The range of the extracted value of CC by the three methods was relatively consistent in the 45° zenith angle range comparing to that in the 60° and 75° zenith angle ranges (Figure 5b). Comparing to the CC extraction from the HP classification, the CC extraction by the CHM-based method had an underestimation trend along with the increasing FOV (Figure 5c,d).

Validation of CC Extraction by the CHM-based Method with the Extraction Results from HP
The relationship between the CC extracted from HP and the CC extracted by the CHM-based model fit the broke-line relationship instead of a simple linear relationship ( Figure 6). It indicated that the CC extracted from HP might have a saturation issue in the range of a high CC value, comparing to the CC extract by the CHM-based method with a hemispherical FOV. The validation model had different responses for the CHM-based

Validation of CC Extraction by the CHM-Based Method with the Extraction Results from HP
The relationship between the CC extracted from HP and the CC extracted by the CHM-based model fit the broke-line relationship instead of a simple linear relationship ( Figure 6). It indicated that the CC extracted from HP might have a saturation issue in the range of a high CC value, comparing to the CC extract by the CHM-based method with a hemispherical FOV. The validation model had different responses for the CHM-based method with different FOVs (Figure 6). Among all the three validation models (i.e., 45 • , 60 • and 75 • zenith angle ranges; Figure 6), the model with a 45 • zenith angle range had the best performance (R 2 = 0.751, RMSE = 0.053; Figure 6a)   Different spatial resolutions of CHM images also had an impact on the accuracy of the model (Table 2). Along with the decreasing CHM pixel size, the accuracy of the model gradually declined for all the FOVs. The model accuracy declined faster from the 0.5 m to 2.0 m and from the 2.0 to 5.0 m spatial resolutions in a 45° zenith angle range than in the 60° and 75° zenith angle ranges (Table 2). In other words, the model accuracy was relatively stable in the 60° and 75° zenith angle ranges (Table 2). Table 2. Accuracy assessments of CHM-based models of varying spatial resolutions (the range coefficient for the radius is 2). The different radius of the CHM analysis window (i.e., different range coefficient (k) for radius) also had an influence on the model's accuracy (Table 3). There was an obvious decline of model accuracy in the 45° zenith angle range, along with the reduction of the radius range (i.e., the range coefficient for radius (k) from 2 to 3). A slight decline of accuracy was observed in the 60° zenith angle range when k changes from 2 to 3. However, the model accuracy had a great improvement for the validation model with the 75° zenith angle range (i.e., R 2 increased from 0.490 to 0.589; RMSE decreased from 0.066 to 0.059). Table 3. Accuracy assessments of the CHM-based models of varying areas (the resolution of CHM is 0.5 m).

Coefficient
Zenith  Different spatial resolutions of CHM images also had an impact on the accuracy of the model (Table 2). Along with the decreasing CHM pixel size, the accuracy of the model gradually declined for all the FOVs. The model accuracy declined faster from the 0.5 m to 2.0 m and from the 2.0 to 5.0 m spatial resolutions in a 45 • zenith angle range than in the 60 • and 75 • zenith angle ranges (Table 2). In other words, the model accuracy was relatively stable in the 60 • and 75 • zenith angle ranges (Table 2). Table 2. Accuracy assessments of CHM-based models of varying spatial resolutions (the range coefficient k for the radius is 2). The different radius of the CHM analysis window (i.e., different range coefficient (k) for radius) also had an influence on the model's accuracy (Table 3). There was an obvious decline of model accuracy in the 45 • zenith angle range, along with the reduction of the radius range (i.e., the range coefficient for radius (k) from 2 to 3). A slight decline of accuracy was observed in the 60 • zenith angle range when k changes from 2 to 3. However, the model accuracy had a great improvement for the validation model with the 75 • zenith angle range (i.e., R 2 increased from 0.490 to 0.589; RMSE decreased from 0.066 to 0.059).

Validation of CC Extraction by the SHP-Based Method with the Extraction Results from HP
Comparing to the CC extraction from HP, the SHP-based models also had a good performance in the CC extraction from UAV-LiDAR data (Figure 7). The validation models for the relationship between the CC extraction with the SHP-based method and the CC extraction from HP fit the broke-line relationship better than a simple linear relationship (Figure 7), which also indicated that the CC extraction from HP might have a saturation issue in the range of the high CC value, comparing to the extraction results with the SHPbased method in all three FOVs. The grid resolution has slight effects on the model accuracy. The SHP-based method with a 1.5 grid resolution showed the best fit along with the 1:1 line before the segmented point (Figure 7b Figure 7e). However, the model in a 75 • zenith angle range had a slightly lower performance (R 2 = 0.601; RMSE = 0.064; the segmented point is 0.546; Figure 7h).

The Advantages and Disadvantages of CC Extraction from UAV-LiDAR Data by CHM-Based and SHP-Based Methods with a Hemispherical FOV
The extracted CC based on the CHM and SHP data had a good relationship and similar regression trends with the CC extraction from HP classification. Another similarity between the two methods using UAV-LiDAR data was that they both have the best performance at the zenith angle of 45 • compared to that of 60 • and 75 • (Figures 5-7). The accuracy of the two models decreases with the increase of FOV, which was more obvious in the CHM-based model (Figures 6 and 7). One possible reason was that the CHM-based method does not correspond well to a ground-based hemispherical FOV for a large zenith angle [11,37]. A previous study also showed that the CHM-based method might have a better description for the spatial averages of CC [44]. This could be confirmed with the CC results of the CHM-based model with the zenith angle range increase (Figure 5a), which might be the main reason why the accuracy of the CHM-based model continues to decrease with the increase of the FOV. Another interesting phenomenon was that the SHP-based model's result was usually consistent with the 1:1 line before the segmented point, while the CHM result showed an underestimation (Figures 6 and 7). In addition to the reasons above, another reason was that the CHM-based method would underestimate CC for leaf-off deciduous trees (i.e., "17-L" in Figure 5) because the pulses of UAV-LiDAR were more likely to pass between the branches [37]. Overall, the SHP-based model describes the change of the ground hemisphere FOV better than the CHM-based method based on UAV-LiDAR data.
Remote Sens. 2021, 13, x. https://doi.org/10.3390/xxxxx www.mdpi.com/journal/remotesensing estimation in the range of lower CC values before the segmented point with a grid resolution of 1.0 and 2.0, respectively (Figure 7). With the optimum grid resolution of 1.5, the validation model in a 45° zenith angle range had the best performance (R 2 = 0.688; RMSE = 0.065; the segmented point is 0.474; Figure 7b), and the model in a 60° zenith angle range had a similar result (R 2 = 0.674; RMSE = 0.062; the segmented point is 0.463; Figure 7e). However, the model in a 75° zenith angle range had a slightly lower performance (R 2 = 0.601; RMSE = 0.064; the segmented point is 0.546; Figure 7h).

The Advantages and Disadvantages of CC extraction from UAV-LiDAR Data by CHMbased and SHP-based Methods with a Hemispherical FOV
The extracted CC based on the CHM and SHP data had a good relationship and similar regression trends with the CC extraction from HP classification. Another similarity Although the SHP-based model was relatively stable, its accuracy slightly decreased with the increase of the FOV (Figure 7), which might be related to the coverage range of the point cloud and the geographic conditions of poplar plantations [41]. We used the algorithm of SHP to meet the range of the zenith angle as much as possible, but the range of the laser pulse is limited and it is hard to detecting the objects far from the shooting location. Moreover, the point cloud data had been normalized to eliminate the impact of the terrain, which cannot be achieved in the HP classification. Other factors, including understory vegetation, the unevenness of foliage lighting, the noise effects of UAV-LiDAR and other differences between HP and LiDAR data might be magnified in the range of a high zenith angle [64].
The spatial resolution had some impact on the accuracy of the CHM-based model, but the CC extraction results were still robust even when we increase the spatial resolution to 5 m ( Table 2). The accuracy of the CHM-based model had a pronounced decline with the decrease of spatial resolution at a low zenith angle range (i.e., 45 • ). The gap fraction of the canopy was more sensitive to the CC extraction with the range of the low zenith angle, and the small canopy gap might not be the primary factor in CC extraction with a lager zenith angle range [37]. However, the SHP-based method showed a better robustness with the grid resolution change (Figure 7). On the one hand, the grid resolution interval of the SHP-based method was not very large. On the other hand, the SHP-based method is more accurate for the description of the canopy gap in a FOV range [48].
Our radius model of the CHM-based method had different CC extraction results in different zenith angle ranges, and the different distance coefficient (k) provided a probable reference for it. The fixed radius range was not suitable for our study because of the distinguishing difference of tree height between the plots. A diverse FOV was also a big challenge to the accurate radius model. To improve the accuracy of the CHM-based model, it is necessary to find a better radius model or improve the CHM-based algorithm in subsequent research. However, the range of area in the SHP-based method was simple, and the only need was to determine the theoretical maximum horizontal radius at the maximum zenith angle to obtain the same result as the ground FOV [13,44,45], which was verified by the visualization results of the SHP algorithm ( Figure 4). The point cloud density's unevenness in different plots might have an effect on the results [57], but this impact was eliminated by thinning the point clouds of all sample plots.

The Accuracy of the Validation Data (CC Extraction from HP)
The validation data were extracted from HP in three zenith angle ranges ( Figure 8). The distribution of CC in different FOVs was different from the histogram of the results (Figure 8c), and this distribution pattern appeared both in high and low CC conditions (Figure 8a,b). Overall, the extracted CC from HP is larger when the zenith angle is larger, which is consistent with Fiala et al.'s results [65]. The accuracy of the HP method was robust in three zenith angle ranges ( Table 4 [65]. The accuracy of the HP method was robust in three zenith angle ranges ( Table 4). The overall accuracy of all plots was highest in the zenith angle of 45° (F = 0.967), followed by the zenith angle of 75° (F = 0.942) and the zenith angle of 60° (F = 0.938).
In the study, we improved the accuracy of CC extraction in HP by removing the main trunk of poplar trees because the tree trunks account for a certain percentage (Figure 8a,b). In addition, the UAV-LiDAR system could describe the tree crown well, but had difficulty in describing the trunk [66], so it is necessary to eliminate the trunk's influence in HP. HP is essentially a two-dimensional passive sensor in a hemisphere FOV [14], which is different from the UAV-LiDAR point cloud data. Thus, there is a certain underestimate after removing the trunks in HP. Moreover, HP data and UAV-LiDAR data are sensitive to different factors (i.e., sensor differences and sensitivity to the environment) in the CC extraction. In addition to the difference in sensor parameters, the accuracy of the CC results from HP is more sensitive to environmental factors, such as radiometric condition, terrain, wind, unevenness of foliage lighting and the choice of subjective threshold, which would increase this uncertainty [14,41].

Stand Age Effects on the Model Accuracy for CC Extraction
The residual boxplot distribution results were similar in the young poplar plantations (i.e., 8, 11 and 14 years old) in the 45° zenith angle range. A significant overestimation of the CC extraction from UAV-LiDAR data was observed in the old plantations (i.e., 17 and 20 years old; Figure 9a). However, by the CHM-based models, the results appeared an underestimation in the young plantations and an overestimation in the old plantations  In the study, we improved the accuracy of CC extraction in HP by removing the main trunk of poplar trees because the tree trunks account for a certain percentage (Figure 8a,b). In addition, the UAV-LiDAR system could describe the tree crown well, but had difficulty in describing the trunk [66], so it is necessary to eliminate the trunk's influence in HP. HP is essentially a two-dimensional passive sensor in a hemisphere FOV [14], which is different from the UAV-LiDAR point cloud data. Thus, there is a certain underestimate after removing the trunks in HP. Moreover, HP data and UAV-LiDAR data are sensitive to different factors (i.e., sensor differences and sensitivity to the environment) in the CC extraction. In addition to the difference in sensor parameters, the accuracy of the CC results from HP is more sensitive to environmental factors, such as radiometric condition, terrain, wind, unevenness of foliage lighting and the choice of subjective threshold, which would increase this uncertainty [14,41].

Stand Age Effects on the Model Accuracy for CC Extraction
The residual boxplot distribution results were similar in the young poplar plantations (i.e., 8, 11 and 14 years old) in the 45 • zenith angle range. A significant overestimation of the CC extraction from UAV-LiDAR data was observed in the old plantations (i.e., 17 and 20 years old; Figure 9a). However, by the CHM-based models, the results appeared an underestimation in the young plantations and an overestimation in the old plantations with the 60 • and 75 • zenith angle ranges (Figure 9b,c). This phenomenon also appeared in the results of the SHP-based method, but the underestimation was slighter in the 8-year-old plantations and more obvious in the 20-year-old plantations (Figure 9b,c). However, along with the increase of the zenith angle range, there was an obvious increase of the CC value extracted by the HP and SHP-based methods, while there was no significant change by the CHM-based method ( Table 5). The characteristics of the CHM-based method, which had a better description for the spatial averages of CC, cause this phenomenon (see the description in Section 4.1 for a detailed discussion).  Table 5). The characteristics of the CHM-based method, which had a better description for the spatial averages of CC, cause this phenomenon (see the description in Section 4.1 for a detailed discussion).  The distinguished differences in the extraction of CC in young poplar plantations and old poplar plantations had a great impact on the accuracy of the model. The variety of tree heights in different ages plots is the main reason for the difference in CC extractions, and different structures in the poplar age plots would affect the extraction accuracy as well [13,67]. The longer distance between the camera lens and the tree canopy lead to a  The distinguished differences in the extraction of CC in young poplar plantations and old poplar plantations had a great impact on the accuracy of the model. The variety of tree heights in different ages plots is the main reason for the difference in CC extractions, and different structures in the poplar age plots would affect the extraction accuracy as well [13,67]. The longer distance between the camera lens and the tree canopy lead to a low resolution of the canopy, which brings errors in the CC extraction of HP [61,68]. Another reason is the difference of light conditions between plots, because we found there was a slight exposure in HP of these old poplar plantations, especially at the larger zenith angle. Thus, this might be due to the underestimation of CC from HP instead of the overestimation of CC from the UAV-LiDAR point cloud data in old poplar plantations, because the photographs are more sensitive to changes of the environment [41,64]. The age effects caused by the limitation of the camera sensor led to a saturation tendency in the model. However, the regression relationship of the model under the same age plot is theoretically linear. One way of eliminating age effects is to consider regression models for each age plot, but it might not be effective in plantations with uniform planting density, because the narrow range of CC often leads to an insignificant regression relationship in the model [69]. Another way is to find a balance coefficient to normalize all plots, but it is also a challenge due to the difficulty of quantifying complex environmental conditions. Further research is needed to assess the age effects and how to eliminate them.
Overall, the UAV-LiDAR system can obtain reliable forest CC information over extensive areas. This is also a good way to store the 3D information of the forest canopy, in which there is an agreement of the extracted CC results with the hemispherical photos taken in the same location. The automatic CC extraction from UAV-LiDAR data would benefit the description of understory light conditions as well, which is consistent with the finding, in the previous study, that CC extraction is promising for assessing regional solar radiation metrics [70]. Moreover, CC often have a strong relationship with various forest parameters. Thus, our results and methods of CC extraction can become good references for further studies on the extraction of other forest parameters in plantation forests.

Conclusions
Both CHM-based and SHP-based methods show a good performance when extracting CC from UAV-LiDAR data in poplar plantations with a uniform planting density. The CHM-based method had the highest accuracy in a 45 • zenith angle range (R 2 = 0.751, RMSE = 0.053), and the accuracy showed a decline when choosing a zenith angle of 60 • and 75 • (R 2 = 0.707, 0.490; RMSE = 0.053, 0.066). The higher resolution of pixel size enhances the performance of CC extraction by the CHM-based model, and a suitable radius of the CHM analysis window is another factor which can mitigate the bias in the model. The accuracy of the SHP-based method was stable in three zenith angle ranges (R 2 : 0.688, 0.674, 0.601 and RMSE = 0.059, 0.056, 0.058 for 45 • , 60 • and 75 • zenith angle ranges), and it showed a better consistency with the 1:1 line before the segmented point. Comparing to the HP results, the CC extracted from UAV-LiDAR data showed an overestimation in old poplar plantation plots due to the difference between camera and LiDAR sensors. Compared to the HP results, the SHP-based method was more sensitive to the FOV changes than the CHM-based method. Based on our research, the SHP-based model had a better performance in CC extraction than the CHM-based method in different FOVs. Our research proved that the UAV-LiDAR data are useful for CC estimation in forest plantation, which provided a reasonable reference for monitoring regional understory lighting conditions in the future.
Author Contributions: Y.P. contributed by developing the methodology, writing the python script, analyzing the data and writing the manuscript. D.X. contributed by coming up with the initial ideas, guiding the research directions and revising the manuscript. H.W. and D.A. contributed to the field work design, field data collection and data analysis. X.X. contributed by providing information of the study area and interpreting the research results. All authors have read and agreed to the published version of the manuscript.