Automated Cloud Shadow Detection from Satellite Orthoimages with Un-corrected Cloud Relief Displacements

1. Introduction

The utility of optical satellites is increasing due to their ability to periodically observe extensive regions and to enable effective Earth observation missions [1]. However, despite these significant advantages, satellite imagery, particularly in the field of remote sensing, could present various challenges. One such challenge encountered in this field is shadows [2]. These shadows could be caused by different sources, such as buildings, clouds, and terrain. Shadows from these objects could lead to misleading outcomes like classification or object detection because the presence of shadows leads to a significant loss of radiometric information in the shadowed areas [2,3]. This vulnerability that also affects unmanned aerial vehicle (UAV) imagery [4,5]. Given the significant impact shadows could have on the accuracy and reliability of satellite imagery, many researchers perceive shadows as noise and targets for removal or correction in detection operations [6].

Building upon the general challenges posed by shadows, cloud shadows also present significant issues in optical satellite imagery. Due to the annual average cloud cover of 66% [7], a significant amount of information is lost as result of clouds and the corresponding cloud shadows. As noted above, the presence of cloud and cloud shadows in optical satellite images could interfere with ground observations [8,9,10]. Therefore, detecting cloud and cloud shadow is crucial for ensuring the reliability of satellite imagery [11].

To detect cloud shadows, which are perceived as noise, Zhu and Woodcock [12] conducted cloud and cloud shadow detection for the Landsat and Sentinel-2 satellite. Within the scope of cloud shadows, this method involves estimating the projected location of cloud shadows from cloud masks detected a priori by considering the solar and satellite azimuth and zenith angles. Foga et al. [13] reported that this method is well performed and tested by the United States Geological Survey (USGS). In high resolution satellite imagery, Le Hégarat-Mascle and André [14] and Fisher [15] developed a cloud shadow detection method focusing on the geometric relationship between clouds and their shadows for SPOT satellite images. Similarly, Zhong et al. [16] also developed a cloud shadow detection method based on the geometric relationships for Chinese satellite images.

As shown in the above developments, it is well known that the direction of shadow in raw satellite images depends on the direction of sun illumination and sensor viewing direction. It is also well known that ortho-rectified images remove relief displacements due to oblique sensor viewing direction. This implies that the direction of shadow and the direction of sun illumination should coincide in ortho-rectified images. However, we observed that this implication may not be realized for objects such as cloud and cloud shadows.

The reason is that the height of clouds is not considered during the ortho-rectification process, leading to relief displacement concerning cloud height in orthoimages. This phenomenon was also mentioned in the study by Pailot-Bonnétat et al. [17], however, without explicit experiments. This paper will highlight the necessity of the direction of sensor viewing as well as that of sun illumination for determining the direction of shadow of objects with un-corrected relief displacements such as clouds. We will show that relief displacement of clouds still exists in orthoimages and point out the importance of considering sensor geometry for cloud shadow detection from orthoimages. We will then propose an automated approach for cloud shadow detection.

In pursuit of cloud shadow detection, we utilized the OCM (Object-oriented Cloud and Cloud-shadow Matching) method previously researched by Zhong et al. [16] to generate a cloud map first. Candidate regions of cloud shadows were then projected from the cloud map along the direction of shadow and actual cloud regions were detected. We employed for experiments RapidEye satellite images taken at nadir and oblique viewing angles. To examine relief displacements of clouds in orthoimages, geometric correction and ortho-rectification were applied to the data. The direction of cloud shadows in orthoimages were estimated by two cases: the first case considered the direction of sun illumination only and the second case the direction of sun illumination and sensor viewing direction. The estimated direction of cloud shadows was compared with true direction manually measured. In the orthoimages made from nadir images, the angles from the two cases were similar to the true angle. However, in the orthoimages made from high-oblique images, there was a maximum difference of 21.3° between the angle from the first case and the true angle. The difference between the angle from the second case and the true angle was less than 4.0°. We performed cloud shadow detection using the shadow angles from the two cases. Accuracy results showed that the shadow detection considering the angle from the second case improved the average f1-score by 0.17 and increased the average detection rate by 7.7%, compared to the results from the first case. By considering both solar and sensor geometries, higher accuracy in cloud shadow detection was achieved. This result support the estimation of shadow direction from clouds with un-corrected relief displacements and cloud shadow detection proposed in this paper.

2. Materials and Methods

The cloud shadow detection method proposed in this study consists of three main steps: data preprocessing, estimation of the search range for cloud shadow detection, and detection of cloud shadows, as illustrated in Figure 1. In the preprocessing step, the pixel values of a satellite image are converted to Top-of-Atmosphere (TOA) values and cloud pixels are detected by the OCM method [16]. the search range estimation step, metadata at the time of satellite image acquisition, including the satellite and solar zenith and azimuth angles, are used to calculate the direction vector from clouds to cloud shadows. Subsequently, the search range where cloud shadows may be projected was calculated by presumed minimum and maximum cloud heights. In the cloud shadow detection step, the process involves shifting cloud pixels from the minimum to maximum range incrementally along the direction of cloud shadows. Cloud shadows are then detected by selecting the shift value that show the best spectral characteristics of cloud shadows.

2.1. Methods

2.1.1. Data Preprocessing

The pixel values (Digital Number, DN) in a satellite image contain distortions due to variations in solar angle resulting from temporal differences, changes in the distance between the Sun and Earth, and spectral differences [18]. Distortions result in variations in pixel values between images, making it challenging to accurately determine the physical quantities of the surface [19]. Therefore, a process is necessary to convert pixel values into reflectance values corrected for distortions. This study employed the formulas (1) and (2) to convert pixel values into Top-of-Atmospheric reflectance (TOA reflectance) values from the Product Specification provided by Planet Labs [20]. In Equation (1), $RAD\left(i\right)$ represents the radiance of the $i$-th band, where $i$ denotes the number of bands, and ScaleFactor denotes the coefficient for radiance conversion. In Equation (2), ${REF}_{TOA}\left(i\right)$ represents the Top-of-Atmospheric reflectance of the $i$-th band, $SunDist$ is the distance between the Earth and the Sun at the time of image acquisition, $EAI\left(i\right)$ denotes Exo-Atmospheric Irradiance, and $SolarZenith$ represents the solar zenith angle. The cloud map is generated using a cloud map with low confidence (CML) of the OCM method. The CML represents the operation of loosening the thresholds of the VNIR bands to sufficiently detect clouds [16].

$$RAD\left(i\right)=DN\left(i\right)\times ScaleFactor,$$

(1)

$${REF}_{TOA}\left(i\right)={\displaystyle \frac{RAD\left(i\right)\times \pi \times SunDist}{EAI\left(i\right)\times Cos\left(SolarZenith\right)}},$$

(2)

2.1.2. Calculation of Cloud Shadow Direction

The geometric principle behind shadow formation is that sunlight is obstructed by objects, resulting in the casting of shadows opposite to the solar direction. Shadows can be detected by utilizing the solar azimuth angle and the height of the object. However, as mentioned in the introduction, considering the geometry of the sensor is essential for cloud shadow detection. This is illustrated in Figure 2 and Figure 4. Figure 2 (a) illustrates the projected locations of cloud shadows based on the height of the clouds. ${\mathrm{C}}_0$ represents the true ground location of the clouds. ${\mathrm{C}}_1$ and ${\mathrm{C}}_2$ represent clouds locations at different heights ${∆\mathrm{h}}_1$ and ${∆\mathrm{h}}_2$. ${\mathrm{C}\mathrm{S}}_1$ and ${\mathrm{C}\mathrm{S}}_2$ denote the projected cloud shadows from ${\mathrm{C}}_1$ and ${\mathrm{C}}_2$ due to the oblique sun direction.

Figure 2 (b) illustrates the positions of clouds and cloud shadows projected in an orthoimage. If ortho-rectification is performed ideally by considering the height of clouds, all clouds would be positioned as at ${\mathrm{C}}_{\mathrm{o}}$ in Figure 2 (b). The cloud shadows for clouds with various heights would be projected along the line from the sun to ${\mathrm{C}}_{\mathrm{o}}$. However, in reality, the displacement due to cloud height will not be corrected during the ortho-rectification process. As a result, clouds depicted in the orthoimage are positioned along the line from the sensor to ${\mathrm{C}}_{\mathrm{o}}$on ${\mathrm{C}}_1$ to ${\mathrm{C}}_2$ according to their heights. With an increase in cloud height, relief displacement occurs in the orthoimages, posing a challenge in estimating cloud shadow positions solely based on the geometry of the sun [17].

Figure 3 illustrates clouds and their shadows in the original images before orthorectification (a) and in the orthoimage (b). The cloud shadow pixels in (b) have been ortho-rectified based on the height information from the corresponding surface heights. As a result, it is challenging to discern significant differences between the cloud and cloud shadow positions in (a) and (b). This indicates that the cloud relief displacement in the orthoimages has not been effectively corrected.

Figure 4 illustrates the cloud relief displacement caused by the satellite's viewing angle. Figure 4 (a) and (b) illustrate the positions of cloud relief displacement in vertical and oblique images, respectively. Here, ${\mathrm{C}}_{\mathrm{o}}$ represents the ground point directly beneath the cloud, and ${\mathrm{C}}_1$and ${\mathrm{C}}_2$ indicate clouds at different heights PC, ${\mathrm{P}\mathrm{C}}^1$, and ${\mathrm{P}\mathrm{C}}^2$ represent the positions of the projected ground point below the clouds and the positions of the clouds based on their heights. As the satellite captures images at more oblique angle, ${\mathrm{P}\mathrm{C}}^2$ in (b) exhibits a more pronounced relief displacement compared to ${\mathrm{P}\mathrm{C}}^2$ in (a). In summary, relief displacement in cloud height persists in orthoimages. As the satellite captures images at higher angles, the variation in cloud height becomes more pronounced. Therefore, it is crucial to consider the satellite's azimuth angle, which allows for the incorporation of cloud height variations.

To detect clouds that still exhibit relief displacement in orthoimages, it is necessary to determine the direction from clouds to their shadows. The calculation method for the cloud-to-shadow direction vector ($\overrightarrow{C}$) is illustrated in Figure 5. Since the relief displacement of clouds is not corrected in orthoimages, the clouds are projected onto the position of the Projected cloud in Figure 5. When applying the solar direction vector ($\overrightarrow{S}$), Projected

cloud does not align with the direction of the cloud shadow. Therefore, it is necessary to consider the sensor direction vector ($\overrightarrow{V}$). Both vectors ($\overrightarrow{S}$ and $\overrightarrow{V}$) can be calculated using the metadata at the time of image acquisition. $\overrightarrow{C}$ can be obtained through the difference operation of $\overrightarrow{S}$ and $\overrightarrow{V}$, as shown in Equation (3). In Equation (4), ${\theta }_S$, ${\phi }_S$, ${\theta }_V$​, and ${\phi }_V$ represent, respectively, the solar zenith and azimuth angles, and the sensor zenith and azimuth angles.

$$\overrightarrow{ C }=\overrightarrow{ S }-\overrightarrow{ V },$$

(3)

$$\overrightarrow{ C }=\pi +{\tan }^{-1}\left({\displaystyle \frac{(\mathrm{t}\mathrm{a}\mathrm{n}({\theta }_S)\times \mathrm{s}\mathrm{i}\mathrm{n}({\phi }_S\left)\right)-(\mathrm{t}\mathrm{a}\mathrm{n}({\theta }_V)\times \mathrm{s}\mathrm{i}\mathrm{n}({\phi }_V\left)\right)}{(\mathrm{t}\mathrm{a}\mathrm{n}({\theta }_S)\times \mathrm{c}\mathrm{o}\mathrm{s}({\phi }_S\left)\right)-(\mathrm{t}\mathrm{a}\mathrm{n}({\theta }_V)\times \mathrm{c}\mathrm{o}\mathrm{s}({\phi }_V\left)\right)}}\right),$$

(4)

2.1.3. Calculation of the Cloud Shadow Range

Once the cloud-to-shadow direction is determined, cloud shadow locations can be determined using the height of clouds [16]. However, in optical satellites, it is not possible to determine the height of the cloud. It is necessary to make assumptions about the minimum and maximum heights of the cloud to estimate the search range for cloud shadow detection. In mid-latitude regions of the Northern Hemisphere, excluding tropical areas, the maximum height of clouds is known to be 12 km [16,21]. The minimum height varies depending on the research; however, based on previous studies [22], we assumed a minimum cloud height of 200 m, as illustrated in Figure 6. The minimum-maximum search range ($D_{\left(m\right)}$) for cloud shadows is calculated using the minimum-maximum height of clouds ($h_{min}$ and $h_{max}$) and the cloud-to-cloud shadow direction angle ($\overrightarrow{C}$), as shown in Equation (5).

$$D_{\left(m\right)}=h\times ‖\overrightarrow{ C }‖,$$

(5)

2.1.4. Cloud Shadow Detection

The proposed cloud shadow detection method in this study involves shifting the clouds using the search range for cloud shadows and the cloud shadow direction angle. Subsequently, infrared band-based statistical analysis is conducted to identify the point with the minimum value, indicating the detection of cloud shadows.

Clouds have varying heights for each object; hence, cloud shadow detection needs to be performed for each cloud object. In this study, we employ the outline information of clouds within the cloud mask to individually extract cloud objects. The clouds extracted for each object are shifted fbased on the Ground Sample Distance (GSD) in the cloud shadow direction angle, ranging from the minimum to the maximum detectable distance. To calculate the shifted position of clouds, the ground coordinates of the shifted clouds are computed using Equations (6) and (7) as shown in Figure 7. These ground coordinates are then transformed into image coordinates.

Clouds have varying heights for each object; hence, cloud shadow detection needs to be performed for each cloud object. In this study, we employ the outline information of clouds within the cloud mask to individually extract cloud objects. The clouds extracted for each object are shifted based on the Ground Sample Distance (GSD) in the cloud shadow direction angle, ranging from the minimum to the maximum detectable distance. To calculate the shifted position of clouds, the ground coordinates of the shifted clouds are computed using Equations (6) and (7) as shown in Figure 7. These ground coordinates are then transformed into image coordinates.

$$\left[\begin{array}{c}{X}_{geo}\\ Y_{geo}\end{array}\right]=\left[\begin{array}{c}{GT}_0+\left({GT}_1\times {Boundar{y}_{TL}}_x\right)+D_{\left(m\right)}\times Cos\left(\overrightarrow{ C }\right)\\ {GT}_3+\left({GT}_5\times {Boundar{y}_{TL}}_y\right)+D_{\left(m\right)}\times Sin\left(\overrightarrow{ C }\right)\end{array}\right],$$

(6)

$$\left[\begin{array}{c}{X}_{Img}\\ Y_{Img}\end{array}\right]={\left[\begin{array}{cc}{\displaystyle \frac{(X_{geo}-{GT}_0)}{{GT}_1}}& {\displaystyle \frac{(Y_{geo}-{GT}_3)}{{GT}_5}}\end{array}\right]}^T,$$

(7)

Due to the unavailability of the true cloud height or the true distance of shadow projection for each cloud object, it is essential to incorporate a measure considering the spectral characteristics of cloud shadows. The formulation of the measure is as follows: First, cloud shadows typically exhibit low reflectance in the near-infrared (NIR) band. Second, to identify points with low reflectance values, it is imperative to exclude cloud areas with the cloud mask from operation. Third, the areas with the most affected by the shadow correspond to areas where the sum of NIR reflectance within the shifted cloud mask is at its minimum. Accordingly, we define the measure in Equation (8), and its explanation is provided below. $x\mathrm{r}\mathrm{e}\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{s}$ pixels, $A$ the region of a shifted cloud object, $B$ the cloud mask map, $h\left(x\right)$ the statistics of the NIR band as defined in Equation (9), $f\left(x\right)$ the region of cloud objects satisfying the condition, and $g\left(x\right)$ the objective function to satisfy the minimum value among the cloud object regions that meet the condition. In Equation (9), ${\mu }_{NIR}$ and ${\sigma }_{NIR}$ represent the mean and standard deviation of the NIR band, respectively. $\alpha$ denotes the Z-score, and a value of 1.96 was applied for $\alpha$.

$$g\left(x\right)=\left\{\begin{array}{c}if x \in A\setminus B, h\left(x\right)<0.17, D_{\left(m\right)}\le D_{max}, f\left(x\right)\\ if x \in A\setminus B, h\left(x\right)\ge 0.17, D_{\left(m\right)}>D_{max}, \ne Cloud shadow\\ otherwise, next loop\end{array}\right.,$$

(8)

$$h\left(x\right)=({\mu }_{NIR}+\alpha \times {\sigma }_{NIR}),$$

(9)

$$Potential Cloud Shadow Region\left(PCSR\right)=argmin\left[g\left(x\right)\right],$$

(10)

$$Potential Cloud Shadow=({NIR}_{TOA}<0.17 and {\displaystyle \frac{{NIR}_{TOA}}{{RED}_{TOA}}}>1.0) \in PCSR,$$

(11)

We detect cloud shadows in the following procedures. If the distance of the shifted cloud object is less than the maximum detectable distance ($D_{max}$), the iterative process is repeated by moving one pixel along the clod shadow direction. Through iterative operations, the region with the lowest statistics in the near-infrared band is identified and confirmed as the cloud shadow region. This region is formally defined as the Potential Cloud Shadow Region, as described by Equation (10). Finally, we assume that one cloud object corresponds to one cloud shadow object. Within the confirmed area, we detect Potential Cloud Shadow (PCS) pixels using Equation (11). Figure 8 (a) represents the state before noise removal, while (b) represents after noise removal. In PCS pixels, we extract an object with the largest outer contour, considering the rest as noise.

2.1. Materials and Validation Method

In this study, cloud shadow detection experiments were conducted using RapidEye images at 5m spatial resolution. As the images were provided in Level 1B, which is uncorrected for geometric distortion, geometric correction and ortho-rectification were applied based on the method developed in our previous studies [23,24,25]. To examine the relief displacement of clouds based on sensor geometry, this study collected images captured by sensors with different viewing geometry. Figure 9 shows images used for the experiments and Table 1 summarizes the metadata of the images used in the experiments. All three images were taken in September over the Dongducheon-si and Gaesong-si, whose latitude difference is less than 0° 7′ 35″ degree. Therefore, the sun illumination angles for the three images were very similar.

Image (a) was captured obliquely tilting the satellite to the west, with a sensor zenith angle of 16.3° and an azimuth angle of 281.3°. Image (b) was captured near vertically tilting the satellite to the east, with a sensor zenith angle of 3.8° and an azimuth angle of 99.8°. Image (c) was captured obliquely tilting the satellite to the east, with a sensor zenith angle of 17.1° and an azimuth angle of 98.8°. Note that the second image was captured with near nadir viewing and that the other two images was captured with significant viewing angles in opposite direction to each other. Images (a) ~ (c) corresponded to cloud cover percentages of 1.33%, 0.58%, and 12.66%, respectively. The cloud map was generated using the OCM method.

For verification of the method proposed here, reference cloud shadow directions were measured through visual interpretation of the angles from cloud to cloud shadow. Similarly, reference cloud shadow regions were extracted manually by visual interpretation. The accuracy assessment involved qualitative and quantitative analyses by comparing the reference data with the results from our method. Qualitative analysis was performed by visually interpreting the reference data. The accuracy assessment involved quantitative analysis using a confusion matrix, which includes True Positive (TP), True Negative (TN), False Positive (FP), and False Negative (FN). Precision, recall, and f1-score, calculated from the confusion matrix results (Equations 12, 13, and 14), represent the level of accuracy, with values closer to 1 indicating higher accuracy.

$$Precision={\displaystyle \frac{TP}{TP+FP}},$$

(12)

$$Recall={\displaystyle \frac{TP}{TP+FN}},$$

(13)

$$\mathrm{F}1 score=2\times {\displaystyle \frac{Precision\times Recall}{Precision+Recall}},$$

(14)

3. Results

3.1. Verification of Cloud Shadow Direction

In this section, verification of the direction angle for cloud shadow was conducted. Verification data was collected by visually inspecting 10 direction angles from cloud to cloud shadow in each image. The collection criteria involved connecting the center of cloud to the center of cloud shadow or connecting the edge of cloud to the corresponding edge of cloud shadow, depending on the shape of cloud and cloud shadows. Verification was performed by comparing the direction vector $\left(\overrightarrow{C_1}\right)$ of cloud shadow obtained using only the solar direction vector and the direction vector $\left(\overrightarrow{C_2}\right)$ obtained using both the solar direction vector and the sensor direction vector against the verification data ($\overrightarrow{T}$).

Figure 10 shows the 10 verification data extracted from Scene-1, which was obtained by tilting the satellite to the west by 16.3°. Figure 11 highlights the comparison of cloud shadow detection results for in Scene-1. The left images represent the NIR-R-G composite images indicating the two direction angles of cloud shadow, $\overrightarrow{C_1} \mathrm{a}\mathrm{n}\mathrm{d} \overrightarrow{C_2}$ and the verification data $\overrightarrow{T}$. The right image indicates the bounding boxes for cloud shadow shifted from the bound box of cloud along the two cloud shadow directions. As shown in Figure 11, $\overrightarrow{C_2}$ correctly points towards the center of the cloud shadow and the bounding box shifted along $\overrightarrow{C_2}$covers the entire cloud shadow region.

Table 2 shows the angle values of $\overrightarrow{C_1},\overrightarrow{C_2,}$ and $\overrightarrow{T}$, the angle differences between $\overrightarrow{C_1}$ and $\overrightarrow{T}$, and $\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{a}\mathrm{n}\mathrm{g}\mathrm{l}\mathrm{e}\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}\mathrm{e}\mathrm{r}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{e}\mathrm{s}\mathrm{b}\mathrm{e}\mathrm{t}\mathrm{w}\mathrm{e}\mathrm{e}\mathrm{n}\overrightarrow{C_2}$ and $\overrightarrow{T}$ for 10 verification data. The absolute mean angle difference between $\overrightarrow{C_1}$ and $\overrightarrow{T}$ was about 17.4° and that between $\overrightarrow{C_2}$ and $\overrightarrow{T}$ was about 3.2°. Figure 11 and Table 2 supports the determination of cloud shadow directions using the sun and sensor directions, as proposed in this paper.

Figure 12 shows the 10 verification data extracted from Scene-2, which was obtained by nearly vertical viewing angle. Figure 13 illustrates the comparison of cloud shadow detection results in Scene-2. For Scene-2, the differences between $\overrightarrow{C_1} \mathrm{a}\mathrm{n}\mathrm{d} \overrightarrow{C_2}$ were not as large as those in Scene-1. For near vertical images, the relief displacements are not severe and hence the shadow direction angles using the sun direction are somewhat similar to those using the sun and sensor directions. As before, Table 3 $\overrightarrow{C_1}$, $\overrightarrow{C_2}$, $\overrightarrow{T}$, $\overrightarrow{C_1}-\overrightarrow{T}$ and $\overrightarrow{C_2}-\overrightarrow{T}$ for Scene-2. The absolute mean angle difference between $\overrightarrow{C_1}$ and $\overrightarrow{T}$ was about 1.5° and that between $\overrightarrow{C_2}$ and $\overrightarrow{T}$ was about 2.3°. Figure 12 and Table 3 may indicate that sun directions can used directly as cloud shadow directions for near vertical images. They also support the proposed determination of cloud shadow directions using the sun and sensor directions for precise angle determination for near vertical images.

Figure 14 shows the 10 verification data extracted from Scene-3, which was obtained by tilting to the east by 17.1°. Figure 15 illustrates the comparison of cloud shadow detection results in Scene-3. For Scene-3, the differences between $\overrightarrow{C_1} \mathrm{a}\mathrm{n}\mathrm{d} \overrightarrow{C_2}$ were large as those in Scene-1. It is also notable that the location of shadow direction vectors, $\overrightarrow{C_2}$, and the corresponding shadow bounding boxes compared to the locations of $\overrightarrow{C_1}$ and the corresponding shadow bounding boxes in Scene-3 are in opposite to those in Scene-1. This is because Scene-3 was tilted to the east whereas Scene-1 to the west and hence relief displacements occurred in opposite directions. The proposed determination of cloud shadow directions can handle un-corrected cloud relief displacements in various tilt directions.

Table 4 shows $\overrightarrow{C_1},\overrightarrow{C_2,}$ $\overrightarrow{T}$, $\overrightarrow{C_1}-\overrightarrow{T}$, and $\overrightarrow{C_2}-\overrightarrow{T}$ for Scene-3. The absolute mean angle difference between $\overrightarrow{C_1}$ and $\overrightarrow{T}$ was about 21.3° and that between $\overrightarrow{C_2}$ and $\overrightarrow{T}$ was about 3.7°. The magnitude of the absolute mean angle difference between $\overrightarrow{C_1}$ and $\overrightarrow{T}$ for Scene-3 was larger than that for Scene-1. This is because the magnitude of sensor’s tilt angle for Scene-3 was larger than for Scene-1. The results and observations in Figure 14 and Table 4 strongly support that the proposed determination of cloud shadow directions using the sun and sensor directions should be valid.

3.2. Cloud Shadow Detection Accuracy

In this section, we conducted an accuracy analysis of cloud shadow detection results using $\overrightarrow{C_1}$ and $\overrightarrow{C_2}$. Figure 16 shows the intermediate process of cloud shadow detection using $\overrightarrow{C_1}$ and $\overrightarrow{C_2}$ for Scene-3. In this Figure 16, (a) represents the NIR-R-G composite image. (b) to (e) illustrate the location of a shifted cloud object using $\overrightarrow{C_1}$ and $h\left(x\right)$ as described by Equation (9), while (f) to (i) illustrate the same for $\overrightarrow{C_2}$. The NIR located in the top right corner of the figure represents the statistics as described in Equation (9), which are equivalent to $h\left(x\right)$. Orange pixels represent a shifted cloud object, while red pixels indicate PCSR as described in Equation (10). yellow directional vectors represent $\overrightarrow{C_1}$, and green directional vectors $\overrightarrow{C_2}$.

The distance for cloud shadow detection in C6 of Scene-3 was confirmed to be 989 m as seen in (h) using $\overrightarrow{C_2}$. However, it was observed that when using $\overrightarrow{C_1}$, the corresponding cloud shadow was not detected at distances similar to those of $\overrightarrow{C_2}$, and was skipped. For the same reasons mentioned in Section 4.1, a precise cloud shadow directional angle was not set. These results further support the validity of cloud shadow direction determination using the direction of the sun and sensor in the orthoimage.

Figure 17 illustrates the post-processing process following the detection of PCSR. Figure 17 (a) represents the reference data for clouds and cloud shadows, while (b) denotes before post-processing, and (c) post-processing. The red pixels in (b) and (c) represent final cloud shadow pixels, while yellow pixels denote remaining noise pixels such as terrain shadows or water bodies. PCS was detected within the boundary of the shifted cloud using Equation (11), such as in C6 of Figure 15, which is indicated by the green bounding box. Subsequently, the PCS of the object with the largest outer contour of the detected pixels, as shown in Figure 17 (C), is considered the final cloud shadow.

The results of cloud shadow detection are presented in Figure 18, Figure 19, and Figure 20 for Scene-1, 2, and 3, respectively. In these figures, (a) to (c) represent, in order, the reference cloud and cloud shadows, the detection results using $\overrightarrow{C_1}$, and the detection results using $\overrightarrow{C_2}$. White pixels represent clouds detected separately by the OCM method [16] and black pixels represent cloud shadows. One can observe that the shadow detection results using $\overrightarrow{C_1}$ could not find cloud shadows at all, in particular for Scene-1 and Scene-3. This is because the wrong cloud shadow directions could not guide the shifted cloud shadow bounding boxes to the location close enough to the real cloud shadow.

On the other hand, the cloud shadow detection using $\overrightarrow{C_2}$ has successfully identified the cloud shadows corresponding to the clouds, but there was small shape discrepancy in some results compared to reference data. This discrepancy could be attributed to the variation in the shadow's shape caused by the terrain relief. In the future research, there is a need to reduce this discrepancy through post-processing methods such as watershed algorithm.

Quantitative analysis results for the three scenes using $\overrightarrow{C_1}$ and $\overrightarrow{C_2}$ are summarized in Table 5. They indicate that the metrics values are consistently higher when using$\overrightarrow{C_2}$ for Scene-1 and -3. In the case of Scene-1, the detection rate was 97.5% for $\overrightarrow{C_2}$, while it was 86.2% for $\overrightarrow{C_1}$. In the case of Scene-3, the detection rate for $\overrightarrow{C_2}$ was 94.0%, while it was 80.8% for $\overrightarrow{C_1}$. For Scene-2, the shadow detection results from $\overrightarrow{C_1}$ and $\overrightarrow{C_2}$ were very similar to each other. This is because $\overrightarrow{C_1}$ and $\overrightarrow{C_2}$ were very similar and the location of shifted bounding boxes were also very similar.

The average f1-score for cases using $\overrightarrow{C_1}$ and $\overrightarrow{C_2}$ are 0.51 and 0.68, respectively. The detection accuracy of $\overrightarrow{C_2}$ improved by 0.17 compared to $\overrightarrow{C_1}$. The average detection rates for $\overrightarrow{C_1}$ and $\overrightarrow{C_2}$ are calculated as 88.8% and 96.5%, respectively. This result confirms that, even if the relief displacement of clouds is not corrected in the orthoimages, the proposed method allows for more accurate cloud shadow detection.

4. Discussion and Conclusions

This study applies a geometric method to detect cloud shadows satellite images. It is well-known that the direction of shadow in raw satellite images is determined by the direction of sensor and sun illumination. However, shadow directions in orthoimages has not been studied explicitly. This study proposed and verified that the direction of cloud shadow in orthoimages is also determined by both sensor and sun directions. This is because relief displacements due to cloud height above ground are not corrected through ortho-rectification processes. The findings in this paper should be applicable to other objects in ortho- or unmanned aerial vehicle (UAV) images with uncorrected relief displacements, such as high-rise buildings. This can be a new interesting research topic.

Experiment results support that cloud shadow directions cannot be determined solely by the sun illumination vector, in particular, in orthoimages generated from oblique images. The inclusion of sensor geometry improved the accuracy of the estimation of cloud shadow directions. This improvement contributed to significant improvements for detecting cloud shadows. By shifting a bounding box of a cloud region through the cloud shadow direction, we could locate the bounding box on cloud shadows. Cloud shadow regions were extracted precisely by identifying dark pixels within the box. In the future, there is a need to enhance accuracy of cloud shadow detection through post-processing methods such as watershed algorithm. The outcomes of this study are expected to be utilized as precise cloud shadow detection and cloud shadow correction.