Correction of Canopy Shadow Effects on Reflectance in an Evergreen Conifer Forest Using a 3D Point Cloud

: Correction of spectral reﬂectance for shadow and topography in optical remote sensing data is challenging. Here, we corrected for canopy self-shadowing in an evergreen conifer forest in mountainous terrain using a three-dimensional (3D) point cloud. In our approach, the surface was modeled from structure-from-motion processed images provided by an unmanned aerial vehicle; then, the relationship between the observed spectral reﬂectance of the Sentinel-2A / B multispectral instrument, and the simulated sunlit fraction (the percentage of the solar-illuminated area within a Sentinel-2 pixel) was determined based on a ray-tracing scheme using a 3D point cloud. Scene-based and pixel-based linear regressions were applied to remove canopy-shadow and topographic e ﬀ ects from satellite-observed reﬂectance. Scene-based regression resulted in large seasonal changes that caused overcorrection in winter. Pixel-based regression generated stable seasonal proﬁles in both the red and near-infrared reﬂectance values over the conifer canopy; however, excessively smoothed seasonal changes were implemented over deciduous vegetation. In both correction methods, the reﬂection of incident light by the canopy likely improved the correction accuracy by decreasing the contrast between illuminated and shaded pixels in summer and increasing it in winter. The results were also visually compared with those from the Sun-Canopy-Sensor with C (SCS + C) method and the Sentinel-2 Level-2A product. This research demonstrates the e ﬀ ectiveness of using a 3D point cloud to correct for self-shadowing and topographic e ﬀ ects on remotely sensed reﬂectance.


Introduction
Spectral reflectance is a crucial information source in optical remote sensing that is used to extract the status of, and changes in, biogeophysical parameters, land cover, and surface energy budgets. Numerous processing methods have been developed in an attempt to derive accurate, detailed information from remotely sensed spectral reflectance data; however, the ability to resolve specific features is limited, as the reflectance of a target differs depending on its position and the relative positions of the sun and the sensor (i.e., the sun-target-sensor geometry). Major causes of reflectance variation are shadow and terrain, which are closely related because terrain affects shadow coverage. Asner and Warner [1] showed that shadow accounted for up to 30-50% of the variance in red and near-infrared (NIR) spectral reflectance data in the Amazon tropical forest. Vanonckelen et al. [2] indicated that topographic correction of reflectance data improved the accuracy of land cover classification, especially in forest areas. main target of the study due to its seasonally stable reflectance and structure combined with complex terrain for canopy self-shadowing analysis. Satellite data from the Sentinel-2 MSI were adopted due to their 10-m spatial and 5-day temporal resolutions.
The remainder of this paper is organized as follows. Section 2 describes the study site and the data collected. Two correction methods, scene-based and pixel-based linear regression, are described in Section 3 along with the calculation of the sunlit fraction using a 3D point cloud. In Section 4, the corrected results are demonstrated and compared with those obtained using the SCS+C method and the Level-2A product of Sentinel-2. The advantages and improvements of this approach are discussed in Section 5. Section 6 presents a summary of our results and concluding remarks.

Research Site
The research site was a conifer forest in Mochida-danchi artificial forest, located in Tosa, Kochi Prefecture, Japan, where 60-year-old Japanese cedar trees (Cryptomeria japonica) were growing at the time of data acquisition. Japanese cedar is one of the most popular timber trees and has the largest planting area (44% of artificial forest and 18% of the domestic forest areas) in Japan [39]. Tree heights in the research plots were about 25 m, and the average density ranged from 390-420 trees/ha per plot [40]. Figure 1 shows the location and an orthoimage of the research site. A logging road passed through the center of the target plot, and there were small open spaces due to harvesting and windfallen trees. From a broader perspective, this site is located on a north-facing slope; however, there are small south-and north-facing slopes in the local terrain (e.g., the digital surface model shown in Figure 5). The white frame in the orthoimage corresponds to a 10-m-resolution Sentinel-2 pixel that was used for the analysis. The area is 82,500 m 2 ; it corresponds to 825 pixels of Scentinel-2. Notably, the area outside the frame is excluded due to the low density of the 3D point cloud. . The right image is the orthoimage of the site generated by unmanned aerial vehicle/structure-from-motion (UAV/SfM) data. The white frame indicates the area used for the analysis, which corresponds to 10-m-resolution Sentinel-2 multispectral instrument (MSI) data.

Three-Dimensional Point Cloud
The site was observed by a commercially available UAV, a Phantom 4 Pro (https://www.dji.com/), on 8-9 December 2018; at this time of the year, cedar trees have leaves, but deciduous vegetation and herbaceous vegetation are leafless. Video was recorded for about 30 minutes (15 min × 2 flights) at a flight altitude of approximately 120 m above the ground. To limit the matching error in the SfM process, the nose of the UAV was pointed in a fixed direction. Images were captured from the video files every 30 frames (1 second [s]) and were input to the SfM software, Figure 1. Left image shows the location of the research site (Kochi Prefecture, Japan). The right image is the orthoimage of the site generated by unmanned aerial vehicle/structure-from-motion (UAV/SfM) data. The white frame indicates the area used for the analysis, which corresponds to 10-m-resolution Sentinel-2 multispectral instrument (MSI) data.

Three-Dimensional Point Cloud
The site was observed by a commercially available UAV, a Phantom 4 Pro (https://www.dji.com/), on 8-9 December 2018; at this time of the year, cedar trees have leaves, but deciduous vegetation and herbaceous vegetation are leafless. Video was recorded for about 30 minutes (15 min × 2 flights) at a flight altitude of approximately 120 m above the ground. To limit the matching error in the SfM process, the nose of the UAV was pointed in a fixed direction. Images were captured from the video files every 30 frames (1 second [s]) and were input to the SfM software, Pix4Dmapper (https://www.pix4d.com/). Remote Sens. 2020, 12, 2178 4 of 21 In total, 1726 images were used to generate the 3D point cloud. Noise reduction filtering was not applied after point cloud generation.
Seven ground control points were installed along the logging road to match the geographic position of the point cloud to that in the Sentinel-2 MSI image. During the UAV flight, 1-m-diameter white circular cloths were placed on the control points. After the flight, the geographic coordinates (latitude, longitude, and ellipsoidal height) of the control points were measured using a differential Global Navigation Satellite System (GNSS) receiver. Based on these coordinates, the point cloud was converted into a Japan Plane Rectangular IV projection using a 3D affine transformation. The root mean square errors associated with the georeferensing results obtained using seven control points were 0.12 m and 0.13 m in the horizontal and vertical directions, respectively. An ortho-rectified image with 10-cm spatial resolution was also generated (Figure 1).

Sentinel-2 MSI
The Sentinel-2A/2B MSI was adopted for this study due to its spatial and temporal resolutions. Despite our small target area (410 m × 330 m), the 10 m resolution of MSI bands 2, 3, 4, and 8 produced images showing sufficient detail of the forest stand and other structures such as logging roads, within a reasonable number of sample pixels. We also found that a 5-day temporal resolution was favorable for investigating seasonal changes in vegetation and illumination. Multitemporal data were downloaded from the Copernicus Open Access Hub (https://scihub.copernicus.eu/). To balance the seasonal data, we intended to use two scenes for every month in 2019. However, due to the difficulty in obtaining sufficient data, especially owing to frequent cloud cover in summer, we used additional data from the previous two years, 2017 and 2018. Table 1 shows the observation dates and the associated satellite sources (satellites A/B of Sentinel-2). Twenty scenes were used in total (only one scene was available in July and one in September, and no data were available in June). Bands 3, 4, and 8 (corresponding to green, red, and NIR spectral regions, respectively) were geometrically rectified to the same projection as the 3D point cloud with 10-m spatial resolution. We used the Level-1C product [41], so no atmospheric correction was applied.

Methods
Shadow has a significant effect on remotely sensed radiance, so it is a first-order term with respect to controlling radiance [42]. Our basic assumption in this study was that satellite-observed reflectance values have a linear relationship with the percentage of an area illuminated in a pixel (sunlit fraction). Two kinds of regression were adopted: scene-based and pixel-based. Scene-based regression uses the pixels in a single scene as a statistical linear regression sample; in pixel-based correction, time series reflectance with respect to individual pixels is used for linear regression analysis. Below, we discuss the use of a 3D point cloud to calculate the sunlit fraction, as well as scene-based and pixel-based correction. In addition, the DEM (SCS+C correction) and Sentinel-2 level-2A product are discussed, as these were included in a comparison analysis of the reflectance results.

Calculation of the Sunlit Fraction
The sunlit fraction (the percentage of solar-illuminated area within a Sentinel-2 pixel) at a particular observation time was simulated based on a ray-tracing scheme using the 3D point cloud. During this process, the point cloud was used to calculate the shadow area directly, without generating a canopy surface model. Thus, the influence of the point cloud noise was very small, even without the noise-reduction filtering of the SfM-generated point cloud. A conceptual diagram is shown in Figure 2. We divided individual 10-m-resolution Sentinel-2 pixels into 400 (20 × 20) sub-pixels with 0.5-m resolution. Binary judgement of illumination/shade was carried out on a sub-pixel basis, as described below: 1.
"Mean_Sun_Angle" was extracted from the MSI metafile to obtain the solar position (zenith and azimuth angles) at individual observation times.

2.
The sensor position was assumed to be at nadir, as our research site is located almost at the center of the MSI observation swath.

3.
A unique radius was assigned to each point in the 3D point cloud. Here, a "point" is regarded as a sphere for ray tracing. A 0.1-m radius was adopted after several trials.

4.
In our approach, the line-of-sight vector started at the satellite position and was directed toward the center of individual sub-pixels, with 3D coordinates calculated for the first intersection of the vector and the sphere (i.e., a point in the 3D point cloud).

5.
The vector was reflected from the intersection in the solar direction. 6.
A sub-pixel was categorized as shaded or illuminated by determining whether the reflected vector intersected with another sphere.
In step 4, the intersection of the line-of-sight vector and point cloud was evaluated based on whether the distance between the line-of-sight vector and the center of the point cloud was smaller than the radius. Because the satellite was positioned at the nadir, the distance could be calculated from the horizontal coordinates (east, north). The same method was applied in step 6, after rotating the point cloud to orient the solar direction to the nadir. Finally, the sunlit fraction of the 10-m-resolution Sentinel-2 pixel was calculated by counting the number of illuminated sub-pixels among the 400 sub-pixels considered. Because the sub-pixels near the edge of the site did not contain any spheres, the pixels used in the analysis had to have at least one sphere within at least 360 sub-pixels (90%). It should be noted that the optical properties of the 3D point cloud, such as color or brightness, were not used for the sunlit fraction simulation.
Our research site is shaded by mountainous terrain in winter, as it is located on the northern slope of Mt Kuishi. We focused on canopy self-shadowing in this study; therefore, the terrain-induced shadow area was excluded from the regression by calculating the terrain shadow using the 10-m-resolution DEM provided by the Geospatial Information Authority (GSI) of Japan (http://www.gsi.go.jp/index.html).
Remote Sens. 2020, 12, x FOR PEER REVIEW 5 of 21 During this process, the point cloud was used to calculate the shadow area directly, without generating a canopy surface model. Thus, the influence of the point cloud noise was very small, even without the noise-reduction filtering of the SfM-generated point cloud. A conceptual diagram is shown in Figure 2. We divided individual 10-m-resolution Sentinel-2 pixels into 400 (20 × 20) subpixels with 0.5-m resolution. Binary judgement of illumination/shade was carried out on a sub-pixel basis, as described below: 1. "Mean_Sun_Angle" was extracted from the MSI metafile to obtain the solar position (zenith and azimuth angles) at individual observation times. 2. The sensor position was assumed to be at nadir, as our research site is located almost at the center of the MSI observation swath. 3. A unique radius was assigned to each point in the 3D point cloud. Here, a "point" is regarded as a sphere for ray tracing. A 0.1-m radius was adopted after several trials. 4. In our approach, the line-of-sight vector started at the satellite position and was directed toward the center of individual sub-pixels, with 3D coordinates calculated for the first intersection of the vector and the sphere (i.e., a point in the 3D point cloud). 5. The vector was reflected from the intersection in the solar direction. 6. A sub-pixel was categorized as shaded or illuminated by determining whether the reflected vector intersected with another sphere.
In step 4, the intersection of the line-of-sight vector and point cloud was evaluated based on whether the distance between the line-of-sight vector and the center of the point cloud was smaller than the radius. Because the satellite was positioned at the nadir, the distance could be calculated from the horizontal coordinates (east, north). The same method was applied in step 6, after rotating the point cloud to orient the solar direction to the nadir. Finally, the sunlit fraction of the 10-mresolution Sentinel-2 pixel was calculated by counting the number of illuminated sub-pixels among the 400 sub-pixels considered. Because the sub-pixels near the edge of the site did not contain any spheres, the pixels used in the analysis had to have at least one sphere within at least 360 sub-pixels (90%). It should be noted that the optical properties of the 3D point cloud, such as color or brightness, were not used for the sunlit fraction simulation.
Our research site is shaded by mountainous terrain in winter, as it is located on the northern slope of Mt Kuishi. We focused on canopy self-shadowing in this study; therefore, the terrain-induced shadow area was excluded from the regression by calculating the terrain shadow using the 10-mresolution DEM provided by the Geospatial Information Authority (GSI) of Japan (http://www.gsi.go.jp/index.html). Examples of this terrain-induced shadow are shown in Figure 6.

Scene-Based Correction
Scene-based correction uses the pixel data of a single scene as a statistical linear regression sample. The observed reflectance showed a linear relationship with the sunlit fraction ( Figure 3a). After linear regression, the reflectance of an individual pixel can be expressed as follows:

Scene-Based Correction
Scene-based correction uses the pixel data of a single scene as a statistical linear regression sample. The observed reflectance showed a linear relationship with the sunlit fraction ( Figure 3a). After linear regression, the reflectance of an individual pixel can be expressed as follows: where ρ (l,p) is the observed reflectance at the pixel position (l, p) in a Sentinel-2 image, R (l,p) is the sunlit fraction at the same pixel position, G and O are the gain and offset, respectively, of the linear regression, and ε (l,p) is the regression residual. The regression residuals represent local variation in land cover reflectance (Figure 3b). The regression residuals are preserved following the correction. Corrected reflectance is made by adding this regression residual to the extrapolated reflectance at a 1.0 sunlit fraction (Figure 3c).
Remote Sens. 2020, 12, x FOR PEER REVIEW 6 of 21 where ( , ) is the observed reflectance at the pixel position ( , ) in a Sentinel-2 image, ( , ) is the sunlit fraction at the same pixel position, and are the gain and offset, respectively, of the linear regression, and ( , ) is the regression residual. The regression residuals represent local variation in land cover reflectance (Figure 3b). The regression residuals are preserved following the correction. Corrected reflectance is made by adding this regression residual to the extrapolated reflectance at a 1.0 sunlit fraction (Figure 3c).

Pixel-Based Correction
In pixel-based correction, linear regression analysis of the sunlit fraction and reflectance was carried out pixel by pixel using multitemporal data. This method was used because Japanese cedar shows little seasonal variation in its reflectance, although the leaves change to brownish green for several weeks in early spring. Figure 4a shows the time series of the sunlit fraction and reflectance for one pixel over the cedar canopy; the reflectance changes seasonally along the linear regression line. The regression model for individual pixels is given below: where is the observed reflectance at time , is the sunlit fraction at , and are the respective gain and offset of the linear regression, and is the regression residual at . Figure 4b shows the time series of the observed reflectance, along with the estimated reflectance values derived by the regression. Here, we assume that these regression residuals are related to the seasonal changes in the cedar trees. These regression residuals were used in the correction for canopy shadow, similar to the scene-based correction approach. Figure 4c shows the corrected reflectance values obtained by adding individual regression residuals to the estimated reflectance values derived from regression coefficients and a sunlit fraction of 1.0.

Pixel-Based Correction
In pixel-based correction, linear regression analysis of the sunlit fraction and reflectance was carried out pixel by pixel using multitemporal data. This method was used because Japanese cedar shows little seasonal variation in its reflectance, although the leaves change to brownish green for several weeks in early spring. Figure 4a shows the time series of the sunlit fraction and reflectance for one pixel over the cedar canopy; the reflectance changes seasonally along the linear regression line. The regression model for individual pixels is given below: where ρ t is the observed reflectance at time t, R t is the sunlit fraction at t, G and O are the respective gain and offset of the linear regression, and ε t is the regression residual at t. Figure 4b shows the time series of the observed reflectance, along with the estimated reflectance values derived by the regression.
Here, we assume that these regression residuals are related to the seasonal changes in the cedar trees. These regression residuals were used in the correction for canopy shadow, similar to the scene-based correction approach. Figure 4c shows the corrected reflectance values obtained by adding individual regression residuals to the estimated reflectance values derived from regression coefficients and a sunlit fraction of 1.0. The maximum sample size of the regression was 20; however, this number decreased if a particular pixel was shaded by terrain in winter. The analysis results showed that 42% of the analyzed pixels had the maximum number of samples, with 88% having more than 15 samples. The minimum sample size was 12.

Comparison with SCS+C Correction
SCS+C correction [5] was adopted for comparison with existing topographic correction methods, although the topographic techniques used a DEM as opposed to a point cloud. SCS+C correction is based on scene-based linear regression of reflectance values with respect to the angle of incidence (i.e., the angle between the sun direction and the normal vector of the inclined surface). Following Soenen and Peddle [5], the SCS+C correction can be expressed as where is the normalized reflectance, is the uncorrected reflectance, is the terrain slope, is the solar zenith angle, is the incidence angle. Parameter is given by where and are derived from linear regression of the uncorrected reflectance and the incidence angle, as follows: In Equation (3), the denominator represents the variation in radiance from the inclined surface with respect to the angle between the surface normal and solar direction. Parameter C compensates indirect solar illumination, which contributes a small amount of brightness due to atmospheric scattering, even under a shadow. The cosine of the solar zenith angle (cos ) in the numerator functions to convert reflectance into radiance, eliminating the effects of the solar zenith angle in the calculation of top-of-atmosphere reflectance. The cosine of the slope angle (cos ) is adopted for the SCS correction to consider the increase in the area of instantaneous field of view (IFOV) of the sensor due to the slope surface. For the calculation of slope angle and incidence angle, we used the same DEM as that for the terrain-shadow calculation mentioned above.

Sentinel-2 MSI Level-2A Product
Topographic and atmospheric corrections were applied to the Level-2A product of Sentinel-2 MSI [15][16][17]. We also used this product for corrected reflectance comparisons; however, based on the availability of Level-2A data (available after 12/15/2018 for our research site at the time of the analysis), 15 scenes among 20 observation dates (listed in Table 1) were used for the comparison. These were The maximum sample size of the regression was 20; however, this number decreased if a particular pixel was shaded by terrain in winter. The analysis results showed that 42% of the analyzed pixels had the maximum number of samples, with 88% having more than 15 samples. The minimum sample size was 12.

Comparison with SCS+C Correction
SCS+C correction [5] was adopted for comparison with existing topographic correction methods, although the topographic techniques used a DEM as opposed to a point cloud. SCS+C correction is based on scene-based linear regression of reflectance values with respect to the angle of incidence (i.e., the angle between the sun direction and the normal vector of the inclined surface). Following Soenen and Peddle [5], the SCS+C correction can be expressed as where L n is the normalized reflectance, L is the uncorrected reflectance, α is the terrain slope, θ is the solar zenith angle, i is the incidence angle. Parameter C is given by where a and b are derived from linear regression of the uncorrected reflectance and the incidence angle, as follows: In Equation (3), the denominator represents the variation in radiance from the inclined surface with respect to the angle between the surface normal and solar direction. Parameter C compensates indirect solar illumination, which contributes a small amount of brightness due to atmospheric scattering, even under a shadow. The cosine of the solar zenith angle (cos θ) in the numerator functions to convert reflectance into radiance, eliminating the effects of the solar zenith angle in the calculation of top-of-atmosphere reflectance. The cosine of the slope angle (cos α) is adopted for the SCS correction to consider the increase in the area of instantaneous field of view (IFOV) of the sensor due to the slope surface. For the calculation of slope angle and incidence angle, we used the same DEM as that for the terrain-shadow calculation mentioned above.

Sentinel-2 MSI Level-2A Product
Topographic and atmospheric corrections were applied to the Level-2A product of Sentinel-2 MSI [15][16][17]. We also used this product for corrected reflectance comparisons; however, based on the  Table 1) were used for the comparison. These were downloaded from the Copernicus Open Access Hub and geometrically corrected to represent the same projection as the other data. Figure 5 shows the 3D point cloud and digital surface model (DSM) generated by SfM. The 3D point cloud provides an oblique view of a 30 × 30-m 2 area (3 × 3 Sentinel-2 pixels). The canopy assumes a conical shape composed of needle leaf clumps of sub-meter size. The canopy height in the target plot has lower variation because all trees were planted in 1956 and 1957 [40], and a dense canopy evenly covers the area, with the exception of several small gaps, as shown in the DSM. Figure 5 shows the 3D point cloud and digital surface model (DSM) generated by SfM. The 3D point cloud provides an oblique view of a 30 × 30-m 2 area (3 × 3 Sentinel-2 pixels). The canopy assumes a conical shape composed of needle leaf clumps of sub-meter size. The canopy height in the target plot has lower variation because all trees were planted in 1956 and 1957 [40], and a dense canopy evenly covers the area, with the exception of several small gaps, as shown in the DSM.

Point Cloud and Sunlit Fraction
The top images of Figure 6 show six examples of the Sentinel-2/MSI top-of-the-atmosphere reflectance. In January and November, terrain-induced shadows (described in Section 3.1) were cast onto some parts of the research area. These terrain-induced shadows in Figure 6 are overlapped with white frames that were calculated using the DEM. In March, May, and September, north-facing slopes were darker than south-facing slopes, although the difference was less than that in winter. In July, the difference in reflectance between slope directions was much smaller due to the higher solar elevation. In contrast to the cedar canopy, the reflectance values of the logging road were high in band 8 (NIR) because it was fully covered by grass-type vegetation. Small gaps also showed higher NIR reflectance of low deciduous broad-leaved trees and understory vegetation. The bottom of Figure 6 shows the sunlit fraction simulated by the 3D point cloud. The appearances of brightness were similar to the Sentinel-2 reflectance; that is, the contrast was higher in winter, and south-facing slopes were bright throughout the year. Because the terrain-induced shadow was not calculated, the brightness in the white frames differed from satellite reflectance; however, these areas were eliminated from the regression analysis.   The top images of Figure 6 show six examples of the Sentinel-2/MSI top-of-the-atmosphere reflectance. In January and November, terrain-induced shadows (described in Section 3.1) were cast onto some parts of the research area. These terrain-induced shadows in Figure 6 are overlapped with white frames that were calculated using the DEM. In March, May, and September, north-facing slopes were darker than south-facing slopes, although the difference was less than that in winter. In July, the difference in reflectance between slope directions was much smaller due to the higher solar elevation. In contrast to the cedar canopy, the reflectance values of the logging road were high in band 8 (NIR) because it was fully covered by grass-type vegetation. Small gaps also showed higher NIR reflectance of low deciduous broad-leaved trees and understory vegetation. The bottom of Figure 6 shows the sunlit fraction simulated by the 3D point cloud. The appearances of brightness were similar to the Sentinel-2 reflectance; that is, the contrast was higher in winter, and south-facing slopes were bright throughout the year. Because the terrain-induced shadow was not calculated, the brightness in the white frames differed from satellite reflectance; however, these areas were eliminated from the regression analysis.  Figure 7 shows examples of reflectance values corrected by the three correction methods (scenebased, pixel-based, and SCS+C), together with the uncorrected observed reflectance values. The corrected pixels are inside the white frame (reflectance in the area outside the frame corresponds to the observed reflectance). In scene-based correction, the reflectance values were significantly higher in winter (January and November) than in other seasons. In addition, differences due to variation in slope direction and canopy cover (evergreen cedar canopy versus other deciduous vegetation) were smaller, especially in November. In summer, the corrected scene was similar to the observed scene. In this correction, differences in reflectance values resulting from slope direction were reduced in summer, but remained in spring through autumn.    Figure 7 shows examples of reflectance values corrected by the three correction methods (scene-based, pixel-based, and SCS+C), together with the uncorrected observed reflectance values. The corrected pixels are inside the white frame (reflectance in the area outside the frame corresponds to the observed reflectance). In scene-based correction, the reflectance values were significantly higher in winter (January and November) than in other seasons. In addition, differences due to variation in slope direction and canopy cover (evergreen cedar canopy versus other deciduous vegetation) were smaller, especially in November. In summer, the corrected scene was similar to the observed scene. In this correction, differences in reflectance values resulting from slope direction were reduced in summer, but remained in spring through autumn.  Figure 7 shows examples of reflectance values corrected by the three correction methods (scenebased, pixel-based, and SCS+C), together with the uncorrected observed reflectance values. The corrected pixels are inside the white frame (reflectance in the area outside the frame corresponds to the observed reflectance). In scene-based correction, the reflectance values were significantly higher in winter (January and November) than in other seasons. In addition, differences due to variation in slope direction and canopy cover (evergreen cedar canopy versus other deciduous vegetation) were smaller, especially in November. In summer, the corrected scene was similar to the observed scene. In this correction, differences in reflectance values resulting from slope direction were reduced in summer, but remained in spring through autumn.   In pixel-based correction, scenes showed similar reflectance values throughout the year, and there was no apparent difference with respect to slope direction. Given that some pixels were shaded by terrain in winter, the reflectance values could be estimated by regression using another season. Additionally, the logging road and deciduous trees showed higher reflectance values than the cedar canopy in all three bands (green, red, and NIR).

Corrected Reflectance
In the SCS+C correction, gray pixels were not corrected because the solar incident angle on the sloped terrain was greater than 90 degrees. Scenes from May to September looked similar to those using scene-based correction; differences with respect to slope were corrected more than when using the scene-based method. However, winter scenes had a different appearance, with high local contrast due to the terrain slope calculated using the 10-m resolution DEM. Figure 8 shows time series of the 10th, 50th, and 90th percentile points of observed and corrected reflectance values in bands 4 (red) and 8 (NIR). The observed reflectance values were higher in summer in both bands, although there was much less seasonal change in the red reflectance values. The gradual increase in red reflectance values from October is most likely due to the seasonal change in the cedar trees, as it is not synchronized with solar elevation; cedar leaves are brown from winter to early spring in this area due to low temperatures. In contrast, NIR reflectance during the same season was stable. The ranges between the 10th and 90th percentiles for bands 4 and 8 were larger in winter due to the significant contrast between illuminated and shaded pixels.
From Figure 8, a summary of some of the key characteristics of the three correction schemes is given below: • Scene-based correction generated much higher reflectance values in winter than in summer. In pixel-based correction, scenes showed similar reflectance values throughout the year, and there was no apparent difference with respect to slope direction. Given that some pixels were shaded by terrain in winter, the reflectance values could be estimated by regression using another season. Additionally, the logging road and deciduous trees showed higher reflectance values than the cedar canopy in all three bands (green, red, and NIR).
In the SCS+C correction, gray pixels were not corrected because the solar incident angle on the sloped terrain was greater than 90 degrees. Scenes from May to September looked similar to those using scene-based correction; differences with respect to slope were corrected more than when using the scene-based method. However, winter scenes had a different appearance, with high local contrast due to the terrain slope calculated using the 10-m resolution DEM. Figure 8 shows time series of the 10th, 50th, and 90th percentile points of observed and corrected reflectance values in bands 4 (red) and 8 (NIR). The observed reflectance values were higher in summer in both bands, although there was much less seasonal change in the red reflectance values. The gradual increase in red reflectance values from October is most likely due to the seasonal change in the cedar trees, as it is not synchronized with solar elevation; cedar leaves are brown from winter to early spring in this area due to low temperatures. In contrast, NIR reflectance during the same season was stable. The ranges between the 10th and 90th percentiles for bands 4 and 8 were larger in winter due to the significant contrast between illuminated and shaded pixels.
From Figure 8, a summary of some of the key characteristics of the three correction schemes is given below: • Scene-based correction generated much higher reflectance values in winter than in summer.

Scene-Based Correction
As shown in Figure 8, scene-based correction generated much higher reflectance values in winter. Notably, the corrected NIR-band reflectance in winter (0.4) was almost double that in summer (0.2). Figure 9a shows all regression lines in the NIR band (band 8). Two scatter plots are shown in the figure, for December and July, corresponding to the highest (0.362) and lowest (−0.029) regression gains, respectively. The large change in reflectance values versus a unit sunlit fraction in winter (plotted by a circle) resulted in the higher regression gain. On the other hand, no clear change in reflectance values was evident in July; hence, a lower regression gain was recorded. This trend is the

Scene-Based Correction
As shown in Figure 8, scene-based correction generated much higher reflectance values in winter. Notably, the corrected NIR-band reflectance in winter (0.4) was almost double that in summer (0.2). Figure 9a shows all regression lines in the NIR band (band 8). Two scatter plots are shown in the figure, for December and July, corresponding to the highest (0.362) and lowest (−0.029) regression gains, respectively. The large change in reflectance values versus a unit sunlit fraction in winter (plotted by a circle) resulted in the higher regression gain. On the other hand, no clear change in reflectance values was evident in July; hence, a lower regression gain was recorded. This trend is the same as the results for July and October reported by Gu and Gillespie [4], although they reflect the relationship between the DEM-based solar incidence angle and shortwave infrared reflectance. In addition to this higher gain, the lower range of the sunlit fraction in winter (about 0.0-0.4) required large extrapolation to the 1.0 sunlit fraction in this correction. Thus, corrected reflectance values were higher in winter.
Remote Sens. 2020, 12, x FOR PEER REVIEW 11 of 21 same as the results for July and October reported by Gu and Gillespie [4], although they reflect the relationship between the DEM-based solar incidence angle and shortwave infrared reflectance. In addition to this higher gain, the lower range of the sunlit fraction in winter (about 0.0-0.4) required large extrapolation to the 1.0 sunlit fraction in this correction. Thus, corrected reflectance values were higher in winter.
(a) (b) (c)  Figure 9b shows the time series of regression gains for the three bands. The gains were synchronized with solar elevation, with the highest (lowest) gain at the winter (summer) solstice. Figure 9c shows the linear relationship between the regression gains and the solar zenith angle. The gains in bands 3 and 8 showed sharp changes at both ends of the solar zenith angle. In band 4, restrained change was evident at both ends, implying that the accuracy of this correction method could be improved by considering reflection by a sunlit canopy. Figure 10 shows the gain and offset of pixel-based linear regression for the three Sentinel-2/MSI bands. The gains in bands 4 and 8 of the cedar canopy on the north-facing slope are larger than those of the south-facing slope; thus, the north-facing canopy had higher reflectance sensitivity to the sunlit fraction. On the other hand, the intercept of the south-facing slope was high; i.e., it showed high reflectance even for a zero sunlit fraction.   Figure 9b shows the time series of regression gains for the three bands. The gains were synchronized with solar elevation, with the highest (lowest) gain at the winter (summer) solstice. Figure 9c shows the linear relationship between the regression gains and the solar zenith angle. The gains in bands 3 and 8 showed sharp changes at both ends of the solar zenith angle. In band 4, restrained change was evident at both ends, implying that the accuracy of this correction method could be improved by considering reflection by a sunlit canopy. Figure 10 shows the gain and offset of pixel-based linear regression for the three Sentinel-2/MSI bands. The gains in bands 4 and 8 of the cedar canopy on the north-facing slope are larger than those of the south-facing slope; thus, the north-facing canopy had higher reflectance sensitivity to the sunlit fraction. On the other hand, the intercept of the south-facing slope was high; i.e., it showed high reflectance even for a zero sunlit fraction.

Pixel-Based Correction
Remote Sens. 2020, 12, x FOR PEER REVIEW 11 of 21 same as the results for July and October reported by Gu and Gillespie [4], although they reflect the relationship between the DEM-based solar incidence angle and shortwave infrared reflectance. In addition to this higher gain, the lower range of the sunlit fraction in winter (about 0.0-0.4) required large extrapolation to the 1.0 sunlit fraction in this correction. Thus, corrected reflectance values were higher in winter.
(a) (b) (c)  Figure 9b shows the time series of regression gains for the three bands. The gains were synchronized with solar elevation, with the highest (lowest) gain at the winter (summer) solstice. Figure 9c shows the linear relationship between the regression gains and the solar zenith angle. The gains in bands 3 and 8 showed sharp changes at both ends of the solar zenith angle. In band 4, restrained change was evident at both ends, implying that the accuracy of this correction method could be improved by considering reflection by a sunlit canopy. Figure 10 shows the gain and offset of pixel-based linear regression for the three Sentinel-2/MSI bands. The gains in bands 4 and 8 of the cedar canopy on the north-facing slope are larger than those of the south-facing slope; thus, the north-facing canopy had higher reflectance sensitivity to the sunlit fraction. On the other hand, the intercept of the south-facing slope was high; i.e., it showed high reflectance even for a zero sunlit fraction.   Figure 11 shows scatterplots for four sample pixels: pixels 1 and 2 are on the north-facing slope and pixels 3 and 4 are on the south-facing slope (pixel locations are shown in Figure 5b). On the north-facing slope, the larger gain and smaller intercept were due to wide ranges of seasonal change in both the sunlit fraction and reflectance values. In contrast, the small gain and large intercepts for the south-facing slope are the result of the lower sensitivity of the reflectance to small changes in the sunlit fraction at higher values. Although the trends differed depending on the slope, almost the same reflectance values were observed for the highest sunlit fraction in summer, for example, about 0.20-0.25 reflectance in the case of band 8. This means that the decrease in reflectance during winter was restrained on the sunward slope due to canopy reflection. the sunlit fraction and reflectance values. In contrast, the small gain and large intercepts for the southfacing slope are the result of the lower sensitivity of the reflectance to small changes in the sunlit fraction at higher values. Although the trends differed depending on the slope, almost the same reflectance values were observed for the highest sunlit fraction in summer, for example, about 0.20-0.25 reflectance in the case of band 8. This means that the decrease in reflectance during winter was restrained on the sunward slope due to canopy reflection. Figure 12 shows the coefficient of determination and root mean square of the residual (RMSR) of the linear regression. The coefficient of determination (the proportion of the variance in reflectance that is predictable by the sunlit fraction) was higher for the north-facing slope; it was almost a value of 1 in band 8. This demonstrates the trend in the regression gain shown in top panels in Figures 11. Dynamic seasonal changes in both the reflectance and sunlit fraction on the north-facing slope resulted in a higher coefficient of determination. On the other hand, the regression gain was almost zero for the south-facing slope, i.e., the reflectance was less predictable due to flat scatterplots. Although the coefficient of determination was lower (almost a value of 0 in bands 4 and 8), the regression residual was smaller on the south-facing slope. Thus, the relationships were linear (but flat) with respect to the sunlit fraction, with a similar magnitude of regression error.   Figure 12 shows the coefficient of determination and root mean square of the residual (RMSR) of the linear regression. The coefficient of determination (the proportion of the variance in reflectance that is predictable by the sunlit fraction) was higher for the north-facing slope; it was almost a value of 1 in band 8. This demonstrates the trend in the regression gain shown in top panels in Figure 11. Dynamic seasonal changes in both the reflectance and sunlit fraction on the north-facing slope resulted in a higher coefficient of determination. On the other hand, the regression gain was almost zero for the south-facing slope, i.e., the reflectance was less predictable due to flat scatterplots. Although the coefficient of determination was lower (almost a value of 0 in bands 4 and 8), the regression residual was smaller on the south-facing slope. Thus, the relationships were linear (but flat) with respect to the sunlit fraction, with a similar magnitude of regression error. Figure 13 shows the time series of observed and corrected NIR reflectance values of four sample pixels over the cedar canopy. There were large differences in the observed values between north-facing and south-facing slopes as the solar elevation declined during the winter. Because there is only a small seasonal change in cedar trees and all pixels are selected from the same plot, this difference was due mostly to the illumination conditions associated with canopy shadow and terrain. On the other hand, corrected reflectance values showed a significantly smaller seasonal change, although the reflectance displayed large variations among pixels. The reflectance values of north-facing pixels (pixels 1 and 2) increased in summer due to a larger regression gain, as shown in Figure 11. In particular, the corrected reflectance values of point 1 were approximately 0.1 higher than those of the other three pixels, as the extrapolated reflectance values were higher for a 1.0 sunlit fraction. A significant change was not evident in south-facing pixels due to the small change in the reflectance versus the sunlit fraction. Seasonal change was considered in this correction using regression residuals; however, the change was very small (i.e., noise).

Coefficient of determination
Root mean square residual  Figure 13 shows the time series of observed and corrected NIR reflectance values of four sample pixels over the cedar canopy. There were large differences in the observed values between northfacing and south-facing slopes as the solar elevation declined during the winter. Because there is only a small seasonal change in cedar trees and all pixels are selected from the same plot, this difference was due mostly to the illumination conditions associated with canopy shadow and terrain. On the other hand, corrected reflectance values showed a significantly smaller seasonal change, although the reflectance displayed large variations among pixels. The reflectance values of north-facing pixels (pixels 1 and 2) increased in summer due to a larger regression gain, as shown in Figure 11. In particular, the corrected reflectance values of point 1 were approximately 0.1 higher than those of the other three pixels, as the extrapolated reflectance values were higher for a 1.0 sunlit fraction. A significant change was not evident in south-facing pixels due to the small change in the reflectance versus the sunlit fraction. Seasonal change was considered in this correction using regression residuals; however, the change was very small (i.e., noise).    Figure 13 shows the time series of observed and corrected NIR reflectance values of four sample pixels over the cedar canopy. There were large differences in the observed values between northfacing and south-facing slopes as the solar elevation declined during the winter. Because there is only a small seasonal change in cedar trees and all pixels are selected from the same plot, this difference was due mostly to the illumination conditions associated with canopy shadow and terrain. On the other hand, corrected reflectance values showed a significantly smaller seasonal change, although the reflectance displayed large variations among pixels. The reflectance values of north-facing pixels (pixels 1 and 2) increased in summer due to a larger regression gain, as shown in Figure 11. In particular, the corrected reflectance values of point 1 were approximately 0.1 higher than those of the other three pixels, as the extrapolated reflectance values were higher for a 1.0 sunlit fraction. A significant change was not evident in south-facing pixels due to the small change in the reflectance versus the sunlit fraction. Seasonal change was considered in this correction using regression residuals; however, the change was very small (i.e., noise).  Although we focused on corrections for the conifer canopy, we also provide results for some deciduous vegetation types. Pixel-based regression assumes a small seasonal change in the cedar tree. Thus, the results are unrealistic for deciduous vegetation, such as the grass on the logging road or deciduous trees. In the case of the logging road, the trends in the regression gain and offset differed in the north-south and east-west directions ( Figure 10). This was due to the relative position of the sun. The north-south direction was less affected by shadows because there was no canopy on the south side; thus, the seasonal change was mainly due to the vegetation. However, the east-west direction had canopies on the south side that were likely responsible for the large seasonal change in reflectance. It should be noted that the RMSR was larger in both directions of the logging road and for other deciduous vegetation ( Figure 12); that is, the seasonal changes in the reflectance were not clearly linear with respect to the sunlit fraction. Figure 14 shows a scatterplot of the logging road (pixel 5) and deciduous trees (pixel 6). Both plots show large residuals from the regression lines compared with that for the cedar canopy in Figure 11; thus, the change here depends on both reflectance and shadow, further complicated by the seasonal changes in greenness and vegetation structure, respectively. Residuals were distributed below (above) the regression line in spring (autumn) because the vegetation activity was delayed in response to the annual solar cycle.
Although we focused on corrections for the conifer canopy, we also provide results for some deciduous vegetation types. Pixel-based regression assumes a small seasonal change in the cedar tree. Thus, the results are unrealistic for deciduous vegetation, such as the grass on the logging road or deciduous trees. In the case of the logging road, the trends in the regression gain and offset differed in the north-south and east-west directions ( Figure 10). This was due to the relative position of the sun. The north-south direction was less affected by shadows because there was no canopy on the south side; thus, the seasonal change was mainly due to the vegetation. However, the east-west direction had canopies on the south side that were likely responsible for the large seasonal change in reflectance. It should be noted that the RMSR was larger in both directions of the logging road and for other deciduous vegetation ( Figure 12); that is, the seasonal changes in the reflectance were not clearly linear with respect to the sunlit fraction. Figure 14 shows a scatterplot of the logging road (pixel 5) and deciduous trees (pixel 6). Both plots show large residuals from the regression lines compared with that for the cedar canopy in Figure 11; thus, the change here depends on both reflectance and shadow, further complicated by the seasonal changes in greenness and vegetation structure, respectively. Residuals were distributed below (above) the regression line in spring (autumn) because the vegetation activity was delayed in response to the annual solar cycle.  Figure 15 shows the time series of observed and corrected NIR reflectance values of the same two pixels as shown in Figure 14. Although the reflectance values were increased by the large extrapolation from the lower sunlit fraction, seasonal change remained due to the larger regression residuals compared with those of the cedar canopy. However, the typical seasonal change in the reflectance of vegetation was nearly eliminated, as it was synchronized with the sunlit fraction.   Figure 15 shows the time series of observed and corrected NIR reflectance values of the same two pixels as shown in Figure 14. Although the reflectance values were increased by the large extrapolation from the lower sunlit fraction, seasonal change remained due to the larger regression residuals compared with those of the cedar canopy. However, the typical seasonal change in the reflectance of vegetation was nearly eliminated, as it was synchronized with the sunlit fraction.
Although we focused on corrections for the conifer canopy, we also provide results for some deciduous vegetation types. Pixel-based regression assumes a small seasonal change in the cedar tree. Thus, the results are unrealistic for deciduous vegetation, such as the grass on the logging road or deciduous trees. In the case of the logging road, the trends in the regression gain and offset differed in the north-south and east-west directions ( Figure 10). This was due to the relative position of the sun. The north-south direction was less affected by shadows because there was no canopy on the south side; thus, the seasonal change was mainly due to the vegetation. However, the east-west direction had canopies on the south side that were likely responsible for the large seasonal change in reflectance. It should be noted that the RMSR was larger in both directions of the logging road and for other deciduous vegetation ( Figure 12); that is, the seasonal changes in the reflectance were not clearly linear with respect to the sunlit fraction. Figure 14 shows a scatterplot of the logging road (pixel 5) and deciduous trees (pixel 6). Both plots show large residuals from the regression lines compared with that for the cedar canopy in Figure 11; thus, the change here depends on both reflectance and shadow, further complicated by the seasonal changes in greenness and vegetation structure, respectively. Residuals were distributed below (above) the regression line in spring (autumn) because the vegetation activity was delayed in response to the annual solar cycle.  Figure 15 shows the time series of observed and corrected NIR reflectance values of the same two pixels as shown in Figure 14. Although the reflectance values were increased by the large extrapolation from the lower sunlit fraction, seasonal change remained due to the larger regression residuals compared with those of the cedar canopy. However, the typical seasonal change in the reflectance of vegetation was nearly eliminated, as it was synchronized with the sunlit fraction.

SCS+C Correction
In SCS+C correction, the seasonal change in the corrected reflectance was similar to that for the observed reflectance (Figure 8), except for an increase in the average reflectance and a decrease in the reflectance variation during the winter season. However, corrected images showed a much higher local contrast, even in the same cedar tree plot (Figure 7). Figure 16 shows scatterplots of the observed and corrected NIR (band 8) reflectance values in summer and winter; note that the horizontal axis is the local solar incidence angle, not the sunlit fraction.

SCS+C Correction
In SCS+C correction, the seasonal change in the corrected reflectance was similar to that for the observed reflectance (Figure 8), except for an increase in the average reflectance and a decrease in the reflectance variation during the winter season. However, corrected images showed a much higher local contrast, even in the same cedar tree plot (Figure 7). Figure 16 shows scatterplots of the observed and corrected NIR (band 8) reflectance values in summer and winter; note that the horizontal axis is the local solar incidence angle, not the sunlit fraction. In summer, C and cos in Equation (3) were larger (1.50 and 0.94, respectively), and cos and cos varied within a narrower range; therefore, the correction term varied only slightly around a value of 1. In winter, C and cos values were smaller (0.12 and 0.51, respectively), and cos and cos varied more widely, resulting in a larger range for the correction term (0.7-4.5). The correction term in Equation (3) also increased exponentially as the incidence angle increased (i.e., as cos decreased). The rotation-like correction of the SCS+C method did not remove the shadow but rather normalized the terrain-induced variation in reflectance in the scene. Hence, the average reflectance in winter did not increase to the same level as that in summer. In summer, C and cos θ in Equation (3) were larger (1.50 and 0.94, respectively), and cos i and cos α varied within a narrower range; therefore, the correction term varied only slightly around a value of 1. In winter, C and cos θ values were smaller (0.12 and 0.51, respectively), and cosi and cos α varied more widely, resulting in a larger range for the correction term (0.7-4.5). The correction term in Equation (3) also increased exponentially as the incidence angle increased (i.e., as cos i decreased). The rotation-like correction of the SCS+C method did not remove the shadow but rather normalized the terrain-induced variation in reflectance in the scene. Hence, the average reflectance in winter did not increase to the same level as that in summer. Figure 17 shows seasonal profiles of the Level-1C (uncorrected) and Level-2A products of Sentinel-2 in the red and NIR bands, together with sample images, although some data for our research site were not available. After simultaneous corrections for the atmosphere and topography were applied in the Level-2A product, the seasonal profile was flattened in band 8. The image contrast was higher, particularly in winter (10th to 90th percentiles > 0.3), likely due to DEM-based topographic correction in winter. Atmospheric and topographic corrections under regular operation of a global satellite product can be quite challenging but are extremely valuable for users in terms of improving data quality.

Comparison with the Level-2A Product
Remote Sens. 2020, 12, x FOR PEER REVIEW 16 of 21 Figure 17 shows seasonal profiles of the Level-1C (uncorrected) and Level-2A products of Sentinel-2 in the red and NIR bands, together with sample images, although some data for our research site were not available. After simultaneous corrections for the atmosphere and topography were applied in the Level-2A product, the seasonal profile was flattened in band 8. The image contrast was higher, particularly in winter (10th to 90th percentiles > 0.3), likely due to DEM-based topographic correction in winter. Atmospheric and topographic corrections under regular operation of a global satellite product can be quite challenging but are extremely valuable for users in terms of improving data quality.

Reflection by Canopy
In the current methods, incident radiation was defined by the sunlit fraction, which was derived as a binary decision regarding an illuminated or shaded point in the cloud, as shown in Figure 18. However, if an illuminated canopy reflected light, the light would reach the lower canopy. In summer, because the angle between the sunlight and the forest stand is acute, incident radiation reaches deeper into the canopy. As a result, canopy reflection potentially reduces the difference in

Reflection by Canopy
In the current methods, incident radiation was defined by the sunlit fraction, which was derived as a binary decision regarding an illuminated or shaded point in the cloud, as shown in Figure 18. However, if an illuminated canopy reflected light, the light would reach the lower canopy. In summer, because the angle between the sunlight and the forest stand is acute, incident radiation reaches deeper into the canopy. As a result, canopy reflection potentially reduces the difference in brightness between illuminated and shaded areas. In winter, the contrast between illuminated and shaded areas would remain relatively high due to the larger angle between the sunlight and the canopy, as well as the reduced reflection by the canopy. This would result in larger regression gains in winter. Furthermore, this effect is more significant with higher canopy reflectance, as evidenced by the abovementioned sharp changes in bands 3 and 8 in summer (Figure 9c). Although diffuse solar radiation is a possible factor, it probably has limited impact, especially in the NIR band. In this study, the sunlit fraction was regarded as equivalent to the incident radiation in the extremely simplified case that the shaded point cloud was completely dark. To improve the representation, the horizontal axis of regression should be altered from a sunlit fraction to an equivalent scale to that of the incident radiation, because a shaded point should have some brightness due to the sunlight reflected by the surrounding point cloud. Thus, light reflection by the canopy should be considered to simulate non-zero brightness in a shaded canopy. One practical method for simulating reflection by the canopy would be to accumulate the number of illuminated point clouds from all surrounding directions. In addition, local topography should be considered in pixel-based correction, as shown in Figure  19. In summer, the acute angle between the sunlight and the forest stand reduces the relative difference in illumination conditions with respect to the slope direction; however, the difference would be expected to be more remarkable in winter, as there is a large difference in the illuminated area with respect to the slope direction when the binary decision is applied to the sunlit fraction. On the other hand, the contrast in brightness with slope direction would be enhanced in winter if canopy reflection was considered as an additional factor, as it is proportional to the illuminated area. The difference in the illuminated area with slope direction is larger in winter than in summer due to the obtuse angle between the sunlight and the canopy. The canopy reflection enhances this difference by increasing the brightness of the shaded canopy. Focusing on the seasonal difference, an anti-sunward slope shows a larger difference, which is enhanced by taking canopy reflection into consideration. To improve this correction method, the horizontal axis should be altered from a sunlit fraction to an incident radiation scale, similar to scene-based correction. After modification, the range of the horizontal axis would become wider (narrower) for anti-sunward (sunward) slopes due to the effect of canopy reflection. This improvement is the most important point for future research studies in this area. In addition, local topography should be considered in pixel-based correction, as shown in Figure 19. In summer, the acute angle between the sunlight and the forest stand reduces the relative difference in illumination conditions with respect to the slope direction; however, the difference would be expected to be more remarkable in winter, as there is a large difference in the illuminated area with respect to the slope direction when the binary decision is applied to the sunlit fraction. On the other hand, the contrast in brightness with slope direction would be enhanced in winter if canopy reflection was considered as an additional factor, as it is proportional to the illuminated area. The difference in the illuminated area with slope direction is larger in winter than in summer due to the obtuse angle between the sunlight and the canopy. The canopy reflection enhances this difference by increasing the brightness of the shaded canopy. Focusing on the seasonal difference, an anti-sunward slope shows a larger difference, which is enhanced by taking canopy reflection into consideration. To improve this correction method, the horizontal axis should be altered from a sunlit fraction to an incident radiation scale, similar to scene-based correction. After modification, the range of the horizontal axis would become wider (narrower) for anti-sunward (sunward) slopes due to the effect of canopy reflection. This improvement is the most important point for future research studies in this area. slope shows a larger difference, which is enhanced by taking canopy reflection into consideration. To improve this correction method, the horizontal axis should be altered from a sunlit fraction to an incident radiation scale, similar to scene-based correction. After modification, the range of the horizontal axis would become wider (narrower) for anti-sunward (sunward) slopes due to the effect of canopy reflection. This improvement is the most important point for future research studies in this area. This study used two sensors (Sentinel 2A and 2B) simultaneously. The calibration accuracies of both Sentinel 2A and 2B are always better than 5%, with 1-2% bias between two sensors (2A measures higher reflectance than 2B) [43][44][45]. Due to this high calibration accuracy, there was no practical fluctuation in the seasonal reflectance variation (e.g., Figure 8). Although we used data collected in three separate years (one image from 2017, five from 2018, and 14 from 2019), no obvious differences between years were detected due to stable climate conditions and a lack of land use changes during the study period.

Comparisons with the SCS+C Method and Level-2A Product
Unlike conventional DEM-based methods [3][4][5], the correction methods applied in the current study used a 3D point cloud to model terrain and canopy structures and shadows. DEM-based methods express the reduction in incident solar radiation in terms of the increases in the angle of the slope and solar direction; therefore, shadows are not corrected explicitly. As a result, the SCS+C This study used two sensors (Sentinel 2A and 2B) simultaneously. The calibration accuracies of both Sentinel 2A and 2B are always better than 5%, with 1-2% bias between two sensors (2A measures higher reflectance than 2B) [43][44][45]. Due to this high calibration accuracy, there was no practical fluctuation in the seasonal reflectance variation (e.g., Figure 8). Although we used data collected in three separate years (one image from 2017, five from 2018, and 14 from 2019), no obvious differences between years were detected due to stable climate conditions and a lack of land use changes during the study period.

Comparisons with the SCS+C Method and Level-2A Product
Unlike conventional DEM-based methods [3][4][5], the correction methods applied in the current study used a 3D point cloud to model terrain and canopy structures and shadows. DEM-based methods express the reduction in incident solar radiation in terms of the increases in the angle of the slope and solar direction; therefore, shadows are not corrected explicitly. As a result, the SCS+C method reduced reflectance variation due to topography and increased NIR-band reflectance by only about 0.06 in winter ( Figure 8). In the Sentinel-2 Level-2A product, NIR reflectance increased to 0.25-0.30 from the uncorrected Level-1C product in winter (0.10) through simultaneous atmosphere and topography correction (Figure 17), which was the same level achieved by pixel-based correction in the present study. However, local variation in reflectance due to topography was drastically enhanced in winter, perhaps due to the spatial resolution of the DEM used in Level-2A processing. Further analyses are required to confirm this hypothesis.
Canopy-induced shadow can occur despite a sunward slope or flat plane and generates seasonal changes in remotely sensed reflectance, even if surface conditions remain unchanged. We successfully used the canopy structure to simulate shadows, demonstrating that the proposed method can be applied to correct flat terrain in regions where a DEM cannot represent shadows. For example, the influence of the canopy shadow in an Amazon tropical forest [1] can be corrected using the pixel-based method, if the structure and reflectance remain stable for long enough to obtain sufficient time series data. Scene-based correction could also be applied because the solar elevation is much higher than that in high-latitude zones.

General Discussion
Quantitative evaluation of corrected reflectance is difficult in shadow or terrain correction due to the lack of a true reference. Other related research has used qualitative evaluations such as visual comparisons, correlation with the slope angle, spectral stability, and so on [6][7][8]. In this study, image appearance and seasonal profiles were adopted to evaluate the multi-temporal correction. These evaluations clearly showed that differences in observed reflectance values with slope direction can be greatly reduced in both the scene-based and pixel-based methods. In addition, the time series of the corrected reflectance values showed large seasonal discrepancy relative to the scene-based method. Existing shadow/terrain corrections have rarely been evaluated with respect to seasonal characteristics.
The seasonal evaluation in this study highlights the shortcomings of the methods examined and suggests how to proceed in terms of improving future shadow-correction schemes.
The pixel-based method was capable of suppressing excess overcorrection in winter, in contrast to the scene-based method. This is because the wide seasonal range in the sunlit fraction prevents large extrapolation in winter. However, scene-based and pixel-based corrections should be applied as complementary correction approaches. An advantage of the scene-based method is the convenience of single-scene applicability. This method should be applied in regions and seasons with high solar elevation such as at low or mid-latitudes in summer to avoid excessive extrapolation to a sunlit fraction of 1.0. Pixel-based correction is suitable for evergreen vegetation where both reflectance and land-surface structure are stable for sufficiently long periods. It is important to emphasize that deciduous vegetation should be carefully considered, given the structural and spectral changes that are synchronized to seasonal solar movement. Regression-based correction is probably not feasible without appropriate information, such as seasonal changes in the spectral reflectance of the target. In addition, we used the December 3D point cloud to calculate the sunlit fraction for summer, although deciduous vegetation showed a large seasonal change in structure. This choice also generated uncertainty; however, the synchronization of seasonal changes in reflectance and structure (i.e., shadows) in deciduous vegetation is a challenging issue that is beyond the scope of the current study.
Potential drawbacks of the SfM-based point cloud approach include the computational cost associated with a large number of photo images and restricted application over wide areas. UAV flights may be prohibited in some cases for legal or safety reasons. Airborne LiDAR may be an effective alternative, although the point density of that approach is generally much lower than that using a UAV/SfM-derived 3D point cloud. The density of the point cloud is a significant issue in the current approach because it has a large influence on the simulation of the sunlit fraction. If the point density is low, the probability of intersection between the line-of-sight vector and the point cloud will be reduced drastically. In such cases, it is probably advisable to switch from the point cloud approach to a canopy surface model. The effects of point cloud density on correction accuracy are unclear and should be investigated in a future study by reducing the density of the point cloud.

Conclusions
In this study, we investigated the feasibility of using an SfM-derived 3D point cloud to correct shadow effects on reflectance in a conifer canopy. Scene-based correction resulted in overcorrection in winter due to a lower sunlit fraction and greater variation in observed reflectance values relative to the change in the sunlit fraction. In summer, corrected reflectance values were almost the same as the observed reflectance values due to weaker correlation with the sunlit fraction. Pixel-based correction was applied based on the assumption that seasonal changes in reflectance values are much smaller than the shadow effect. Corrected images and time series reflectance data showed small seasonal changes in the cedar canopy and lower variation due to slope direction. These results show the benefit of using the sunlit fraction of incident radiation in dealing with the reflection of incident radiation by the canopy and the propagation of reflected light deeper into the canopy layer.
A novel finding of this research is the favorable performance of shadow correction in flat terrain due to detailed canopy structure modeling using a point cloud. The proposed method corrected canopy self-shadowing even on flat ground, unlike conventional DEM-based correction.
Shadow-corrected reflectance is essential for accurate extraction of land cover types or bio-physical attributes of vegetation using optical remote sensing. This study provides a basis for future utilization of 3D point clouds for shadow correction, with specific issues highlighted here to be addressed as next steps in further developing this approach. Additionally, the application of our proposed method over wide areas using airborne LiDAR is a potential next step in future studies.