ESTIMATION OF OPTIMAL CROWN COVERAGE AND CANOPY SHAPE FOR SHADOW ESTIMATION ON TROPICAL MOIST BROADLEAF FOREST

Shadow fraction is essential for improving the estimation of gross primary production, but it is difficult to be observed by satellite due to the diurnal variations. Therefore, it is necessary to estimate the 3D model with physical parameters by simulating virtual forest reflectance. In this study, we aim to estimate the optimal combination of canopy shape and Crown Coverage (CC) through simulating virtual forests reflectance. First, satellite-derived Tree Height (TH) and CC for virtual forests were compared with the ones obtained by Canopy Hight Model (CHM). Second, virtual forests with different CC and canopy shapes were created, and the reflectance and shadow fraction were simulated. The canopy shape used were cylinder, ellipsoid, half-ellipsoid, and inverted half-ellipsoid. Finally, the simulated reflectance and shadow fraction were validated with Sentinel-2 reflectance and shadow fraction from voxel model. Our results show that the mean TH is 15±2m, and the CC was increased from 10% to 60% in 10% intervals. TH and CC obtained from the satellite had the Root Mean Square Error (RMSE) of 5m and 40%. Ellipsoid with 20% CC shows the lowest RMSE and the smallest discrepancy for shadow fractions at the same sun position. However, other combinations were more accurate in estimating mean daily shadow fraction. This would be caused by only one image adopted in validation, which could be improved by using multi-season images in the future.


INTRODUCTION
Gross Primary Production (GPP) is one of the important components of carbon flux (Beer et al., 2010). Satellite remote sensing is useful for estimating GPP on a global scale, and the MODerate resolution Imaging Spectroradiometer (MODIS) GPP product (Running et al., 2004) has been widely used. For GPP estimation, the theoretical model devised by Monteith (1972) using absorbed photosynthetically active radiation and Light Use Efficiency (LUE) to estimate photosynthetic production has been used, but large uncertainties still remain in the global estimation. One of the uncertainties is the complexity of the tree forest 3D structure. Shadows are created by the geometric relationship between forest structure and sun position. Since shadow is a major determinant of leaf surface conductance (Brooks et al., 1997), it has been observed that shaded leaves are not light saturated (Knohl and Baldoccohi, 2008) and that diffuse light shows higher LUE (Kanniah et al., 2013). Therefore, both shade and sun leaves contribute to GPP, but their ratio is different (Gu et al., 2002). Therefore, incorporating their proportions in the GPP estimation will improve the accuracy. Yan et al. (2017) estimated GPP with higher accuracy than before by using the clearance index (He et al., 2013) to estimate the ratio of diffuse to direct radiance, and the clumping index to take into account the sunlit and shaded leaves. Clumping index, however, is given a certain value for each tree species (Chen et al., 1999;He et al., 2013), but the shadow ratio varies with different structural parameters such as tree height and density even if the tree species is the same. In recent years, in addition to the clumping index, vegetation indices related to shadow fraction (Shadow index; Ono et al., 2015, Shadow Eliminated Vegetation Index;Jiang et al., 2019, Normalized Difference Canopy * Corresponding author Shadow Index; Xu et al., 2019) have been developed. It is expected that these indices used to estimate the shadow fraction independent of the tree species. However, since the ratio of direct to diffuse radiation is diurnal variation due to the shading effect of the forest structure (Li et al., 1992), it is still questionable to estimate GPP using only the sun-shade component of satellite observation time. Hilker et al. (2008) reported that adding the mean daily shadow fraction calculated from a 3D model of the forest to the LUE improved the explanatory power. However, the method using satellite images can only observe at a certain time, so the mean daily shadow fraction cannot be calculated. One possible solution to this problem is to estimate the 3D shape of the forest and use it to estimate the mean daily shadow fraction.
In recent years, structural parameters such as Tree Height (TH) (Simard et al., 2011) and tree density (Crowther et al., 2015) have become globally available. Yang et al. (2017) developed a statistical approach to generate forest structure based on TH and tree density to create a highly accurate essential climate variables retrieval algorithm using the Monte Carlo-based Radiative Transfer (MCRT) models. The MODIS reflectance simulated by this approach was high accurate than the theoretical error of the atmospheric correction product. However, they used field measurement values, so it is expected that the error will increase when used only satellite-derived product is used because the product include errors. Therefore, when creating a virtual forest using only products, optimization of structural parameters is necessary. In addition, it is unclear whether the MCRT model can reproduce actual shadows, since shade and sun leaves are determined by random. Therefore, it is desirable to use a geometric optics model (Fan et al., 2015) to separate sunlit leaves from shaded leaves, but this model requires many parameters such as leaf shape and leaf density, some of which are difficult to estimate globally. On the other hand, Fujiwrara and Takeuchi (2020) developed a simple method for simulating reflectance using the shadow and leaf reflectance as parameters. The accuracy of the reflectance simulated by this method is highly dependent on that of the input shadow fraction. In other words, it was expected that the shadow ratio could also be estimated by selecting the optimal parameter that has the smallest error to the reflectance observed by the satellite from the combination of structural parameters obtained from the products.
An important parameter for creating a virtual forest is the canopy shape. This parameter needs to be assigned explicitly. In MCRT, the canopy shape is usually triangular pyramidal for coniferous forests and ellipsoid for broadleaf forests (Pisek and Chen, 2009;Ligot et al., 2016). Nelson et al. (1997) used several canopy geometries to simulate the airborne laser profiling response. The results showed that the different roughness of each canopy shape had a significant effect on the estimates. Since roughness also affects the shadow fraction, errors in shadow and reflectance due to differences in canopy shape should be evaluated and the optimal canopy shape should be selected. In particular, the difference is expected to be larger in forests where various canopy shapes are assumed, such as tropical forests.
The objective was to estimate the optimal structural parameters for shadow fraction estimation by comparing simulated reflectance and satellite-observed one using a virtual forests. In this study, we tackled the following three tasks: 1) to create several virtual forests using satellite-derived structural parameters, 2) to simulate the reflectance and compare with observed reflectance by satellite sensors, and 3) to simulate the shadow fraction and compare with that calculated from the 3D model. The target satellite sensor was Sentinel-2, because using a sensor with lower spatial resolution would increase the spatial heterogeneity of the forest. The structural parameters of interest were TH, Crown Coverage (CC) and crown shape. TH and CC were obtained from Global Forest Canopy Height (GFCH) 2019 (Potapov et al., 2021) and Global 2010 Tree Cover (GTC) (Hansen et al., 2013), respectively. The spatial resolution of both products is 30m and GTC stores the CC of trees over 5 meters and GFCH stores the tree height in 1 meter units. However, creating a virtual forest by changing all combinations of TH, CC, and canopy shape would increase the computational cost. Therefore, we assumed that the value of TH did not have significant variation in the target forest. In other words, we created a virtual forest with different combinations of CC and canopy shape. The target forest was a tropical forest in Myanmar.

Flow of this study
The flow of this study is shown in Figure 1. First, the TH and CC inputs to the virtual forest were determined by sampling, respectively. Those values were compared to the Canopy Hight Model (CHM) generated by Unmanned Aerial Vehicle (UAV) -Structure from Motion (SfM). Note that the measured values are not used to correct TH and CC. The four tree canopy shapes used are cylinder, ellipsoid, half-ellipsoid, and inverted halfellipsoid, as shown in Figure 2. Second, a virtual forest is created based on those structural parameters, and reflectance and shadow are simulated. The canopy shape is not mixed, and the tree structure is determined based on the allometric equations. The canopy structure is TH, Canopy Length (CL), Canopy Radius (CR) and Diameter at Breast Height (DBH). Third, the simulated reflectance was evaluated using the Sentinel-2 Bottom of Atmosphere (BOA) reflectance. It was evaluated by Root Mean Square Error (RMSE) to the average value of each band. A test of normality was performed for all bands of sampled Sentinel-2 reflectance, and bands with p-values greater than 0.05 were used for evaluation. Finally, the simulated shadow fraction was evaluated by that calculated by the voxel model. The shadow fraction was evaluated by the difference of the mean values.

Study forest
The target forest is located on the Yangon Technological University (16.8 • N, 96.1 • E), Myanmar. There were no mountains or other high elevations around, and the terrain was flat. The climate zone of the region is tropical monsoon climate. (Myanmar Information Management Unit, 2021). Several species were mixed together, and it was difficult to obtain allometric equations for each species. Therefore, we use the allometric equations developed by Antin et al. (2013) in their study of tropical rainforests in western India. The model was recalibrated to estimate DBH, CR and CL using TH.

Creation of validation data for tree height, crown coverage and shadow ratio
The target forest was observed by UAV on 21 September 2018. The camera zenith angle is 0 degrees and the altitude from the ground is about 100 meters. The point cloud was generated by SfM based on this data. The UAV used was Phantom 4 pro, and the software used for SfM was Pix4D. Geometric corrections were performed using GPS log data. Digital Surface Model (DSM) and Digital Terrain Model (DTM) were created respectively, and the CHM was generated by subtracting DTM from DSM. CHM is used for TH and CC validation. The spatial resolution of CHM was 30 m. Since GTC is the percentage of tree cover above 5 m, data with CHM above 5 m were used to validate CC. The shadow ratio was calculated by converting the point cloud into a voxel model. The voxel size was determined so that each voxel has at least one point.

Sentinel-2
), acquired on September 30, 2019. Although about one year has passed since the UAV observation date, the forest condition is considered to be close because it is the same season. In this study, the normality of the reflectance of bands 2, 3, 4, 5, 6, 7, 8, 8A, 11, and 12 obtained by the sample was investigated, and the bands with p-values greater than 0.05 were used.

Setting of tree height and canopy coverage
40 points were visually sampled for CC and TH, respectively (0% coverage was excluded). The obtained CC values were converted to values in 10% increments, assuming that a few percent difference in CC does not affect reflectance or shading. The Weibull distribution (Siipilehto 2006) is usually used for TH, but field measurements are required. Since the purpose of this study is to use only products, we did not assume that distribution. A normal distribution was assumed in this study. If the p-value of the sampled TH is greater than 0.05, the TH of the forest is determined using the normal distribution. If the p-value is less than 0.05, a normal distribution with the mean as the most frequent value and 2m as the standard deviation is used.

Create virtual forest
The trees used were represented by voxel models. The size should match the voxel size of the UAV-SfM. The size of the forest was assumed to be 30 m square with a flat ground surface. The virtual forest is randomly arranged with trees to meet the coverage determined by sampling. It was assumed that if trees were placed completely randomly, the distance between adjacent trees may become very near in some cases.Therefore, tree locations were randomly placed to meet a canopy overlap of less than 10%, same as Yang et al. (2017). 40 forest patterns were created for the same CC and canopy shape. Even the canopy shape was different, the tree positions were the same in the forests with the same pattern.

Reflectance simulation model
The method developed by Fujiwara and Takeuchi (2020) was used to simulate the reflectance. This method calculates the shielding factor for direct and diffuse radiation, defined as Cast Shadow (CS) and Self Cast Shadow (SCS), respectively, using a voxel model and the following equation. In this method, the leaf inclination is assumed 0 degrees, and multiple reflections and atmospheric effects are not considered. The Simple Model of the Atmospheric Radiative Transfer of Sunshine (SMARTS) code, version 2.9.5, was used to calculate the direct and diffuse irradiance.

Leaf and understory reflectance
The spectral reflectance data for the simulation was selected from the ECOSTRESS Spectral Library (version 1.0). Prosopis articulata, a species also found in tropical South Asia, was used for leaf reflectance and Avena fatua for understory reflectance. The wavelength interval is 1 nm.

Point cloud generation
The average density of the generated point cloud was 10.38 points/m 3 . The geometric correction was done by referring to the GPS log data without using ground control points. Therefore, geometric accuracy was not calculated. As mentioned in Section 2.3, there must be at least one point in a voxel, so a voxel size of 50 cm was used.

Comparison of TH and CC obtained by product and measurement
The TH and CC of the 40 sampled points were compared with those generated from the CHM. Figure 3 shows the comparison results of TH. The TH calculated using CHM ranged from 2.5 m to 13 m, whereas the TH obtained from GFCH ranged from 7.5 m to 17.5 m. The GFCH is higher than measured TH. The RMSE was 5.4 m and the maximum error was about 11 m, which was an overestimation. Figure 4 shows the results of the CC comparison. Most of the measured CC values ranged from 60% to 100%, but the CC values obtained from GTC ranged from 20% to 60%. The RMSE was 40.3 % and the maximum error was about 70%, which was an underestimated.

Reflectance simulation
Since the p-value of TH obtained by sampling was 0.02<0.05, the mean value of TH used for the virtual forest was the most frequent value of 15 m, as described in Section 2.5. The minimum and maximum CC values were 10% and 60%, respectively. 24 patterns of virtual forests were therefore created, with 6 different CC (10%, 20%, 30%, 40%, 50%, 60%) and 4 different canopy shapes (cylinder, ellipsoid, half-ellipsoid, and inverted half-ellipsoid). For each pattern, 40 virtual forests were created to simulate the reflectance.   Figure 5 shows the virtual forests created using a cylinder and a half-ellipsoid for a crown shape when CC is 60%. Green indicates the canopy, brown indicates the trunk, and black indicates the forest floor. The height and position of the trees are the same. Shadow images under the same sun position are also shown. It can be seen that the half-ellipsoid shape cause more shadows than the cylinder. Figure 6 shows the difference in the simulated reflectance when the CC is fixed at 60% and the tree canopy shape is changed. In the simulation method used (eq. 4), the reflectance is inversely related to the percentage of shadow. This means that for the same CC, the inverted halfellipsoid has the lowest shadow fraction, followed by the cylinder, ellipsoid and half-ellipsoid. Comparing the reflectance between inverted half-ellipsoid and the half-ellipsoid, the difference was about 1.3 times.
On the other hand, the Figure 7 shows the difference in the sim-ulated reflectance when the crown shape is fixed to an ellipsoid and the CC is changed from 0% to 60%. The reflectance decreased as CC increased from 10% to 40%, but the change from 40% to 60% was small. The reflectance at 10% CC differs from that at 50% by about 1.3 times. It can be seen that the variation value of reflectance becomes smaller as CC increases. Figure 5. An example of a virtual forest using a cylinder and a half-ellipsoid as crown shapes. The height and position of the trees are the same in both forests, but the shadow ratio differs due to the difference in crown shape. Figure 6. Reflectance difference due to canopy shape when CC is 60%. The difference beteween an inverted half-ellipsoid and a half-ellipsoid is about 1.3 times.

Evaluation of simulated reflectance
As described in Section 2.1, the reflectance of the normality band is used to evaluate the result of reflectance simulation. The p-values of the sampled Sentinel-2 pixels and the simulated reflectance at each band were examined, and both were greater than the 0.05 significance level in band 6, 7, 8, 8A, and 11. The sample size is 40.  Table 2 shows the difference between the shadow fraction calculated using the voxel model and that calculated from the virtual forest at the same sun position. The solar azimuth angle is 136 degrees and the zenith angle is 26 degrees. The average shadow fraction, calculated using the voxel model, is 0.27. The results showed that 0.017 was the smallest difference, the crown shape was half-ellipsoid and the CC was 20%. This combination was consistent with the second lowest RMSE combination in Table 1. The difference estimated from a virtual forest with a half-ellipsoid canopy shape and 60% CC, however, is 0.021, a difference of only 0.004. Furthermore, when cylinders or inverted half-ellipsoids are used for the crown shape, the difference is almost always smaller than 0.1 regardless of CC. For all crown shapes, the difference was the largest when CC is 50%. Table 3 shows the difference of the mean daily shadow fraction calculated using the voxel model and virtual forest. The mean daily shadow fraction is the average of the shadows calculated from the sun position at 8:00, 10:00, 12:00, 14:00, and 16:00, and the result calculated using the voxel model is 0.34. The smallest difference was 0.003, the canopy shape was halfellipsoid, and the CC was 50%. There was no tendency for the difference to change with respect to CC. Regardless of the canopy shape, the difference was about 0.2 in the virtual forest with 10% CC. In virtual forests where the crown shape is ellipsoid or half-ellipsoid and the CC is larger than 30%, the difference is smaller than 0.1.  Table 3. The difference between the mean daily shadow fraction estimated from the voxel model and that estimated from the virtual forest. The average mean shadow fraction, calculated using the voxel model, is 0.34. Potapov et al. (2021) compared the accuracy of GFCH with the tree heights estimated by The Global Ecosystem Dynamics Investigation (GEDI) lidar and airborne laser scanner. The results showed that the overall RMSRs were 6.6 m and 9.07 m, respectively. They also investigated the accuracy at each 1 m strata, and reported an underestimation of about 2 m for the present tree canopy coverage. Although we did not conduct a survey using the same method as theirs, the RMSE for the TH calculated by CHM was 5.4 m (Figure 3), which was better than their result. However, it is an overestimation not the underestimation. This is probably due to the training data used to generate the GFCH was limited in the tropical forest due to cloud cover (Potapov et al., 2021). We used the most frequent value of 15 m for the average TH of the virtual forest, but it turned out to be higher than the actual forest. In the future, it will be necessary to determine the tree height by referring to the location of the training data used for the GFCH and the error in each strata.

DISCUSSION
The CC estimated using GTC is about 40% (Figure 4) lower than that calculated from CHM. GTC provides CC for trees taller than 5 m, and the height was calculated using the LiDAR (footprint size is 70 m) of NASA's GLAS (Geoscience Laser Altimetry System) instrument aboard the IceSat-1 satellite (Hansen et al., 2013). This product has been used to estimate CC using Landsat or high spatial resolution images such as Quickbird as training data, but it has been reported that there is high uncertainty in the estimates. The CC used for the virtual forest used a value between the minimum and maximum values, but in the future, a higher CC value should be used.
It can be seen from Figure 5 and Figure 6 that a shape with a rounded top of the tree canopy produces fewer shadows and higher reflectance than a flat shape. This indicates that a round top has larger area to block sunlight. The fact that the reflectance is different even when the upper part of the canopy is the same, as in the case of an ellipsoid and a half-ellipse, or an inverted half-ellipse and a cylinder, indicates that the lower part of the canopy also affects the shadow ratio. From Figure 7, it can be seen that the reflectance decreases as the CC increases from 10% to 40%. This is due to the increase in the area where sunlight is shielded by the increase in the number of trees. On the other hand, from 40% to 60%, the reflectance is almost constant, indicating that there is little change in the percentage of area that blocks sunlight. In this study, we assumed that the TH of the virtual forest follows a normal distribution with a mean value of 15 m and a standard deviation of 2 m. As the standard deviation becomes larger, the roughness is expected to become even larger and the percentage of shadows larger. On the contrary, as the standard deviation becomes smaller, the roughness is expected to become smaller and the percentage of shadow becomes smaller. Table 1 shows that the virtual forest with an ellipsoid canopy shape and CC is 20% has the smallest RMSE, but the difference with other canopy shapes when CC is 20% was small. Therefore, changing the reflectance of the leaves and understory used may change the optimal canopy shape. On the other hand, as CC increased, the difference in RMSE due to differences in canopy shape became larger. Nelson et al. (1997) suggested using ellipsoids when the tree canopy shape could not be identified, but when CC was high, it was expected that the error would be more significant if an appropriate tree canopy shape was not used. As shown in Table 2, the smallest difference between the shadow fraction estimated from the virtual forest and voxel model at the same sun position was 0.017. In that case, the suitable canopy shape was half-ellipsoid and the CC was 20%, but there were several virtual forests where the difference was less than 0.05. Even when the virtual forest with the lowest RMSE in Table 1 ( canopy shape is ellipsoid and CC is 20%) was used, the difference was 0.033 in Table 2. This indicates that the optimal CC and canopy shape derived from the reflectance simulation can be used to estimate the shadows fraction of the forest. However, in Figure 4, CC derived from CHM ranges from 60% to 100%, indicating a large difference from the simulation results. As shown in Figure 7, in our simulations, the reflectance decreased with higher CC because of the increase in shadows. However, it was assumed that even at high CC, if TH variation is small, roughness would be smaller and shadows would be less. In this study, the standard deviation of TH in the virtual forest was fixed at 2 m, so the forest in such a case was not considered. As shown in Table 3, the virtual forest with the smallest difference in daily shadow fraction had a half-ellipsoid crown shape and a CC was 50%. The virtual forest with the lowest RMSE in Table 1 had a difference of 0.149, which was larger than the difference in Table 2. The reason seems to be that the sun position was limited because only one scene Sentinel-2 image was used for the validation. Therefore, we expect to improve the accuracy by adding images from different seasons and other satellite images such as Landsat and Aster to the validation.

CONCLUSION
The objective of this study was to estimate the optimal CC and canopy shapes for shadow fraction estimation through reflectance simulation. CHM generated from UAV-SfM, Sentinel-2 BOA reflectance and voxel models were used to validate TH and CC, reflectance simulation, and shadow fraction, respectively. RMSE of TH and CC obtained from each product was 5m and 40% to those value calculated from CHM, respectively. TH tended to overestimate and CC tended to underestimate. The average TH of the virtual forest was 15 m with a standard deviation of 2 m. The CC was increased from 10% to 60% in 10% intervals. A total of 24 patterns of virtual forests were created, with 40 forests per pattern. When the CC was fixed and only the crown shape was changed, the highest reflectance was obtained for the inverted half-ellipsoid and the lowest for the half-ellipsoid. When the crown shape was fixed and the CC was changed, the reflectance decreased from 10% to 40%. The combination of canopy shape and CC that resulted in the lowest RMSE for Sentinel-2 BOA reflectance was inverted halfellipsoid and 20%. This combination was also 0.03 difference to the shadow fraction calculated from the voxel model. However, CC calculated from the CHM ranged from 60% to 100%, so the difference with CC estimated was large. The reason for this is that even when CC is high, a small variation in TH may result in a less shadow, which has not been taken into account in this study. In estimating mean daily shadow fraction, other combinations were more suitable than the best combination to Sentinel-2 reflectance. This is due to the reflectance simulation, where only one scene Sentinel-2 image was used for validation. In this target forest, by optimizing CC and canopy shape through reflectance simulation, it was possible to estimate the shadow fraction at the time of observation by the satellite sensor, but CC was different from real one, and it was found that the variation of TH also needed to be optimized. In addition, the use of multiple time period images and other satellite images is expected to improve the estimation accuracy of the mean daily shadow fraction.