Algorithms for Moderate Resolution Imaging Spectroradiometer cloud-free image compositing

Abstract The presence of clouds is the biggest obstacle in the investigation of land cover, and many techniques have been developed to detect clouds. However, few indicators have been proposed for the detection of cloud-free conditions. To address this, we propose two indicators for use in compositing 8-day cloud-free images: the B 1 / 7 ratio is the ratio of the band 1 reflectance to the band 7 reflectance of Terra surface reflectance images (MOD09GA), and saturation refers to the color saturation of these images. Here, we describe the principles underlying these two indicators and analyze their characteristics for vegetation, water, urban, and nonvegetation pixels under cloud-free, cloud shadow, and cloudy conditions using MOD09GA from October 16 to 23, 2007, in North China (sample A) and using data published by the U.S. Geological Survey (USGS). We found that the B 1 / 7 ratio and saturation are suitable for extracting cloud-free pixels over land and water, respectively; therefore, we combined these two indicators to develop a single-unified model. In particular, our results demonstrate that the pixels exhibiting the lowest B 1 / 7 ratios should be adopted as cloud-free pixels over land when the value of B 1 / 7 for land surfaces is between 0 and 1, and the surface reflectance of Moderate Resolution Imaging Spectroradiometer band 1 is less than 0.3. Otherwise, the pixels exhibiting the greatest saturation values should be adopted. We used our model to composite cloud-free images for two additional regions in China: the Tarim basin (sample B) from October 2 to 9, 2012, and the coastal areas of southeastern China (sample C) from April 15 to 23, 2013. We compared the cloud-free images of these regions with the 8-day surface reflectance product (MOD09A1) with respect to clouds, cloud shadow, and cirrus clouds, and we found that our proposed cloud-free image compositing approach can accurately eliminate both clouds and cirrus clouds. Specifically, the percentage of residual cloud pixels in sample C was found to be less than that in MOD09A1. Moreover, in the cloud-free images obtained using our newly developed method, cloud-free pixels are typically associated with greater sensor zenith angles and smaller scatter angles than those in MOD09A1. However, our method retains some limitations. In particular, 9.68, 33.22, and 33.00% of cloud-shadow pixels remain in the cloud-free images for samples A, B, and C, respectively.

Resolution Imaging Spectroradiometer (MODIS); therefore, accurate cloud detection is considered to be a crucial component in eliminating the impact of clouds on remote-sensing products. 3,4 Three types of information can be used to detect clouds: spectral, texture, and multitemporal information. Spectral information was first utilized to distinguish clouds from the Earth surface using classification methods, [5][6][7][8] in which thresholds in the value of one band or the difference between two bands were considered. 3,4,9 Subsequently, texture information was used to distinguish clouds from the surface based on thresholds in heteroscedasticity. Such identification is possible because the reflectance of cloud is greater than that of the land surface, which results in greater variance along cloud edges, particularly over the oceans. [10][11][12] With increases in the availability of remote-sensing data, researchers have started to use multitemporal information to distinguish clouds from the Earth surface based on the differences between appropriate benchmarks and series of daily images; such methods are based on the assumption that the reflectance of the Earth surface is more stable than that of cloud-affected pixels, which exhibit frequent variations. [13][14][15][16] Many satellite-based instruments, ranging from the Advanced Very High Resolution Radiometer (AVHRR) to MODIS, employ these methods. For instance, the International Satellite Cloud Climatology Project (ISCCP) 17,18 and the AVHRR Processing Scheme Over Cloud Land and Ocean (APOLLO) cloud detection algorithm 19,20 employ only spectral information, and the Cloud Advanced Very High Resolution Radiometer (CLAVR) employs a series of spectral and spatial variability tests to detect clouds. 1,[21][22][23] Some researchers have also considered the use of temporal information for cloud detection. 4,9,[24][25][26] However, the computational power required for such methods increases with the complexity of the cloud detection algorithms adopted. Ultimately, obtaining composite cloud-free images based on the results of cloud detection has become a time-consuming process. This is evidenced particularly by the case of MOD09A1 published by the United States Geological Survey (USGS), in which only the pixels with the highest quality are considered; in this instance, cloud-cover pixels, cloud-shadow pixels, high-angle pixels, and high-aerosol pixels have been removed. 27 Obtaining data relating to land surface characteristics is one of the primary challenges in land cover research. To address this, some algorithms have been proposed, in which cloud-free pixels are selected based on quantitative physical indexes such as the maximum normalized difference vegetation index (NDVI) 28 and the lowest MODIS band 3 surface reflectance obtained during a specific time period. 29 Maximum NDVI has been used previously as an indicator in gridded vegetation indices (MOD13) for vegetation monitoring. However, the use of this index typically leads to the selection of the pixel with the greatest sensor angle. Conversely, the minimum reflectance value for MODIS band 3, which represents the blue band, tends to consider shadow pixels as cloud-free pixels. Ma et al. used min½tan −1 ðB1∕ðB2 − B1Þ (where B1 and B2 are the bands 1 and 2 references of MODIS, respectively) to select cloud-free pixels. 4 However, this model was shown to be ineffective in distinguishing between clouds, nonvegetated land, and shadows on vegetated land. The Canada Centre for Remote Sensing (CCRS) developed a new compositing technology that employs six criteria to obtain pixels of six different classes, in addition to mixed pixels. The CCRS researchers also compared max(NDVI), minðB1Þ, and max½maxðB2; B6Þ∕B3 with their own multicriteria scheme and reported some details regarding the effects of the bidirectional reflectance distribution function (BRDF); 9 however, many criteria and classification steps are required for the CCRS compositing technology, and this method is associated with considerable uncertainty.
The present study aims to obtain highly accurate composite cloud-free images by employing the fewest criteria possible. According to Kaufman et al. 30 and Levy et al., 31 the reference ratio B 1∕7 between reflectances of 0.66 and 2.1 μm (corresponding to bands 1 and 7 of MODIS, respectively) is a good indicator of cloud cover for land areas. Moreover, as clouds are typically white or gray color, their saturation is lower than that of water pixels in a hue-saturation-intensity (HSI) model, making saturation a good indicator of cloud-free regions for water areas. Therefore, in the present study, we propose B 1∕7 ratio and saturation as appropriate indicators for the selection of cloud-free pixels on land and water, respectively, and present a set of algorithms for compositing cloud-free images for China.

Underlying Principles of the B 1∕7 Ratio
Previous studies have shown that 2.1-μm-wavelength images are not affected considerably by atmospheric scattering; 30 thus, this wavelength can pass through all types of aerosol except dust. Moreover, for vegetated areas, the apparent reflectance at 2.1 μm is twice that at 0.66 μm when a given region is cloud free. Accordingly, the apparent reflectance at 2.1 μm is often used to estimate the surface reflectance at 0.66 μm. 30 These wavelengths (i.e., 0.66 and 2.1 μm) are the central wavelengths of MODIS bands 1 and 7, respectively, and are typically expressed as B 1 and B 7 . The reflectance ratio of these bands, B 1∕7 , for a given land surface has been shown to be positively correlated with scattering angle and the vegetation index derived from the MODIS near-and short-wave infrared bands. 31 where Θ is the scattering angle and defined as in Eq. (2). B NDVI SWIR 1∕7 is a piecewise function related to the vegetation index for shortwave infrared radiation [Eq. (3)], and NDVI SWIR is the vegetation index derived from MODIS near-and short-wave infrared bands [see Eq. (4)].
In Eq. (2), θ 0 , θ, φ 0 , and φ are the solar zenith angle, the sensor zenith angle, the azimuth of the sun, and the azimuth of the sensor, respectively. 8 > > < > > : In Eq. (4), ρ 5 and ρ 7 are the apparent reflectance of MODIS bands 5 and 7, respectively. It is clear from Eq. (3) that B NDVI SWIR 1∕7 ranges from 0.48 to 0.58 and is correlated positively with NDVI SWIR . Moreover, according to a previously described geometric optics model, 32 NDVI SWIR is proportional to the sensor zenith angle (i.e., NDVI SWIR increases with increasing sensor zenith angle for a given vegetation pixel), and B 1∕7 increases with NDVI SWIR .
As the Terra satellite is sun synchronous, θ 0 and φ 0 can be considered constants, such that θ and φ are the primary factors determining Θ. Observations collected by the Terra satellite can be separated into two groups based on (φ − φ 0 ): the backward hemisphere with (φ − φ 0 ) ∈ 0 to 90 deg or 270 to 360 deg and the forward hemisphere with (φ − φ 0 ) ∈ 90 to 270 deg. Values of both Θ and B 1∕7 are typically smaller for the forward hemisphere observations than for the backward hemisphere; therefore, the pixels of the forward hemisphere observations are typically selected first. To sum up, selecting the pixel with the smallest B 1∕7 for a given period to represent a cloudfree pixel is equivalent to choosing the pixel in which the sum of NDVI SWIR and Θ is the smallest.

Underlying Principles of Saturation
The HSI model scales color in hue, saturation, and intensity spaces. Hue characterizes the main color received by the observer. Conversely, saturation represents the degree to which the solid color is diluted by white, and intensity is a measure of the brightness of the color. Thus, saturation decreases as the amount of white in the mix increases 33 and can be defined by a saturation equation as follows: where S represents saturation, and R, G, and B represent red, green, and blue rays, respectively; all variables lie between 0 and 1. The minðR; G; BÞ function finds the variable with the lowest value. Pixels can be classified as cloud-free, cloud shadow, or cloud pixels depending on the locations of the sensor, sun, and clouds ( Fig. 1). The saturation of cloud-free, cloud shadow, and cloud pixels can be expressed as shown in Eqs. (6), (7), and (8), respectively.
where S a is the saturation of a cloud-free pixel, and r R , r G , and r B are the reflectances of MODIS bands 1, 4, and 3, respectively.
where S b represents the saturation of a cloud-shadow pixel, t is the transmissivity of the cloud, and t R , t G , and t B are the transmissivity values of cloud for MODIS bands 1, 4, and 3, respectively.
In Eq. (8), S c is the saturation of a cloud pixel, ρ is the reflectance of cloud, and ρ R , ρ G , and ρ B are the reflectance of the cloud for MODIS bands 1, 4 and 3, respectively.
The diameters of droplets and dust in clouds and fog in the atmosphere are typically in the range 5 to 100 μm, which is much greater than the wavelength range of visible light. Therefore, random scattering occurs when visible light passes through these particles. 34 The scattering coefficients are the same for the red, green, and blue bands; therefore, Eqs. (7) and (8) can be simplified to produce Eqs. (10) and (11), respectively. Fig. 1 A model map of cloud-free, cloud shadow, and cloud pixels.
When the cloud is so thick that t approaches 0, S 0 b also approaches 0; in such cases, the true value of S 0 b is 0 and S 0 c tends to 0, which is in accordance with our knowledge that clouds and fog are typically white or gray. Moreover, when thinner clouds or fewer aerosols are present in the atmosphere, more visible light will reach the Earth's surface; in such cases, S 0 b will be equal to S a , and S 0 c will be smaller than S a according to Eqs. (10) and (11).

Algorithm for Integrated Cloud-Free Image Synthesis
We developed a combined algorithm for cloud-free image compositing by combining B 1∕7 with saturation. Our model adopts the land surface reflectance images for a specific period as the candidate dataset (MOD09GA) and calculates B 1∕7 from the surface reflectance data of bands 1 and 7 and pixel saturation S from the surface reflectance data of bands 1, 4, and 3. Assuming a land surface mask (LSM) of 1, which means it is land pixel, B 1∕7 of (0, 1], and B 1 < 0.3, the pixel exhibiting the lowest value of B 1∕7 is selected as a cloud-free pixel; otherwise, the pixel exhibiting the greatest saturation is selected as a cloud-free pixel. The LSM is a land/ water flag and can be extracted from the 1-km State QA Descriptions (16-bit) layer in MOD09GA. Our newly developed model can be used to extract an index image, HC_time, in which the value of each pixel corresponds to the date of the image least contaminated by cloud, selected from a set of multitemporal remote-sensing images obtained within a certain period. Then, cloud-free images for the surface reflectance (LSREF_clear) can be built based on this index image. The compositing workflow is presented in Fig. 2.
3 Data Processing and Analysis

Data Processing
We employed three samples in the present study, hereafter referred to as samples A, B, and C (Fig. 3). We used sample A to assess the feasibility of distinguishing cloud-free pixels for different land cover types, considering a particular region in northern China (34°53′ to 43°26′ N, 113°49′ to 122°41′ E) that includes Hebei Province, Beijing, Tianjin, Shandong Province, and Bohai Bay. Sample A includes nine land cover classes, divided according to the Land Cover Type Yearly L3 Global 500-m product (MCD12Q1) (Fig. 4).
We used samples B and C to validate the algorithm for deserts, snowy regions, and tropical areas often covered by clouds. The region considered for sample B is located to the south of the Tarim basin and the north of the Tianshan Mountains; 74% and 10% of this sample can be considered deserts and snowy areas, respectively. Conversely, the study region for sample C is located in southern China, adjacent to the South China Sea, in an area often covered by cloud; thus, it was difficult to obtain cloud-free images for this region.
We employed three separate MODIS products, all of which can be downloaded freely using the Warehouse Inventory Search Tool (WIST): the Surface Reflectance 8-Day L3 Global 500 m (MOD09A1), the Surface Reflectance Daily L2G Global 1 km and 500 m (MOD09GA), and the MODIS Land Cover Type Yearly L3 Global 1 km SIN Grid product (MCD12Q1). MO09A1 provides bands 1 to 7 at a resolution of 500 m in an 8-day gridded level-3 product in sinusoidal projection. Each MOD09A1 pixel contains the best possible level-2 gridded (L2G) observation for a given 8-day period, selected based on high-observation coverage, low-view angle, the absence of clouds or cloud shadow, and aerosol loading. MOD09GA provides bands 1 to 7 in a daily-gridded L2G product in sinusoidal projection including 500-m reflectance values and 1-km observation and geolocation statistics. MCD12Q1 incorporates five different land cover classification schemes, derived using a supervised decision-tree classification method. We used land cover type 4 of the MODIS-derived net primary production (NPP) scheme and type 5. Further details of the data used are presented in Table 1.
All preprocessing of these data was performed using the MODIS Reprojection Tool (MRT). Sample A was used to conduct geographic transformations and to mosaic the four parts of the data. Finally, spatial subsets of the images were obtained using a 500 × 500 m mask for the study area. Samples B and C were transformed from the hierarchical data format (HDF) to the tagged image file format (TIFF) for which the projection was sinusoidal and the resolution was 463.31 × 463.31 m.

Analysis of B 1∕7 Ratio
To estimate the ability of B 1∕7 to identify cloud-free pixels according to the MODIS land cover type 4 for NPP, we reclassified the land cover into four groups: water, vegetation (including all    Fig. 5. Some fundamental conclusions can be deduced from the comparison of Figs. 5(a) and 5(b). B 1∕7 was found to be greater than 1 for cloud pixels and very high (i.e., greater than that under cloud-free conditions) for water areas under cloudy conditions. For urban and vegetation areas under cloud-free conditions, B 1∕7 was found to be less than and greater than 1, respectively. Conversely, under cloudy conditions, we found B 1∕7 to be around 1 and greater than the value under cloud-free conditions. Comparison of Figs. 5(c) and 5(d) demonstrates that B 1∕7 of shadow pixels for land areas is around 1 and is higher than B 1∕7 under cloud-free conditions. However, for water areas, we found B 1∕7 to be less than 0 (primarily because B 1 is less than 0), except under cloud-free conditions. In summary, our results demonstrate that B 1∕7 is typically higher for cloud pixels than for cloud-free pixels ( Table 2).  Note: For water pixels under cloud shadow conditions, we found B i to be −0.01; therefore, we did not calculate the slope for this relationship. We also omitted the calculation of the slope for cloud pixels in cloud shadow.
To investigate our results further, we plotted B 1∕7 values for cloud-free, cloud shadow, and cloud conditions for the different land cover types (Fig. 6) using the data presented in Table 2. We found B 1∕7 to be in the range 0 to 1 for land areas under cloud-free conditions; moreover, our results suggest that B 1∕7 is always lower under cloud-free conditions than under cloudy conditions. Comparison of shadow and cloud-free pixels demonstrates the existence of several pixels, in which B 1∕7 is lower for shadow pixels than for cloud-free pixels. However, no such pixels were observed over water areas. Therefore, in cases where a given pixel could be considered to be cloud-free, cloud shadow, or cloud depending on the observation time, we selected the period with the minimum value of B 1∕7 to represent the period with clear conditions.

Analysis of Saturation
To investigate whether the pixel exhibiting the maximum saturation can be considered a cloudfree pixel, we calculated saturation from MODIS land surface reflectance using Eq. (5) (Fig. 7). In particular, we calculated R, G, and B from bands 1, 4, and 3, respectively. In general, we found that the saturation of cloud pixels to be lower than that of cloud-free pixels for all types of land cover, although the saturation for cloud shadow pixels for water and nonvegetation areas was found to be lower than that for cloud-free pixels. Moreover, our results suggest that the saturation of shadow pixels increases with increasing percentage of vegetation cover [Figs. 7(a), 7(c), and 7(d)]. For land cover classes including vegetation, urban, and nonvegetation, the saturation of cloud-free pixels is typically greater than that of cloud and shadow pixels. Based on the above analysis, we suggest that the maximum saturation can act as a good cloud-free pixel index for water and nonvegetation areas.

Results and Discussion
The algorithm proposed here for compositing cloud-free images is similar to that of MOD09A1, in which both select cloud-free pixels from MODIS daily images and combine these pixels. We investigated samples A, B, and C using both the MOD09A1 method and our compositing technique both qualitatively and quantitatively through visual evaluation and analysis of indicators.

Visual Evaluation
Based on visual evaluation, our proposed algorithm performs better in the selection of cloud-free pixels over water (including marine, coastal, lake, reservoir, and river environments) than the MOD09A1 product (Fig. 8). For example, MOD09A1 erroneously defined a cloud pixel as water under the clear-sky conditions [red circles in Figs. 8(a1) and 8(a2)]. Thus, our algorithm performs better than MOD09A1 in such cloud regions, although our method is less spatially consistent than MOD09A1 [e.g., the rectangle for sample B in Figs. 8(b1) and 8(b2)]. Moreover, for sample C [Figs. 8(c1) and 8(c2)], our algorithm removed the cloud in the lower left of the image, whereas MOD09A1 did not. Finally, the ocean color in the cloud-free images obtained using our method is more spatially homogeneous than that obtained from the MOD09A1 product.

Indicator Evaluation
Residual cloud, cloud shadow, and mixed pixels can be found in the composite image, and quantitative evaluation of the relative proportions of these residual pixels is essential. The 1-km State QA Descriptions (16-bit) dataset of the MOD09GA product includes 11 indicators, including five quantitative indicators that are useful for the estimation of the quality of image composition: cloud state, internal cloud, cloud shadow, cirrus detected, and pixels adjacent to cloud.
In the original dataset of daily images, each pixel includes three classes: every day is cloudy; every day is noncloudy; there are at least one cloudy day and at most seven cloudy days throughout the 8-day period. For this third class, the accuracy of the cloud-free composited image has been evaluated previously according to Eq. (12). where i is a state of cloud-free, cloud, cloud shadow, cirrus, and so on, M i represents the pixel number of a given state remaining in the cloud-free image, and M clear oud i represents the number of pixels exhibiting the state in at least one day and at most seven days throughout the 8-day period. Based on Eq. (12), we conducted a comparative evaluation of the quality of the images produced using our compositing technique (Table 3). It is clear that our method performs well in removing cloud and cirrus pixels, but less well for cloud-shadow pixels. In the composite image, the internal cloud in the coastal areas of southeastern China (sample C) is lower than that in the MOD09A1 product, indicating that fewer residual cloud pixels remain in the composite image than in the MOD09A1 product. Moreover, slightly less high cirrus can be found in the composite image than in the MOD09A1 product, although there is slightly more medium height and small cirrus in the composited image. However, residual small cirrus is less than 5% in samples A and B, and that in sample C is 7.15%. More residual cloud shadow pixels are present in the composite image than in the MOD09A1 product with residual cloud shadow greater than 30% in the Northwest Territories (sample B) and the coastal areas of southeastern China (sample C). Table 4 presents average values for sensor zenith angle, scattering angle, and relative azimuth angle for samples A, B, and C for both the composite image and the MOD09A1 product. It is clear that, regardless of whether the B 1∕7 ratio or the maximum saturation method is adopted, the average sensor zenith angle is greater for the composite image than for the MOD09A1 product. This occurs because our method exhibits preference for nonnadir pixels, whereas the MOD09A1 product exhibits preference for nadir pixels. Furthermore, the average scatter angle for the composite image is less than that for the MOD09A1 product, and the average scatter angle obtained according to the maximum saturation method is less than that obtained using the B 1∕7 ratio method. This may have occurred because our method mistakes cloud-shadow pixels for cloud-free pixels, particularly when the maximum saturation method is adopted. The average relative azimuth angle for the composite image is approximately 90 to 270 deg, whereas that of the MOD09A1 product is less than 45 deg or greater than 315 deg. This likely occurs because our method preferentially selects the forward hemisphere pixels as cloud-free pixels, whereas the  MOD09A1 product preferentially selects backward hemisphere observations pixels as cloudfree pixels. Thus, it is clear that our method exhibits obvious differences from the MOD09A1 algorithm. The selection of the maximum saturation pixel as a cloud-free pixel is likely the key factor causing the differences that are apparent from Table 4. The atmosphere has a difference to the extinction coefficient of the red, green, and blue, the longer the light transmission path in the atmosphere, the greater the energy difference among the red, green, and blue, and then the greater the saturation. Therefore, pixels exhibiting large zenith angles are likely to be selected as cloud-free pixels, although the maximum saturation method is more likely to preferentially select such pixels than the B 1∕7 method. Thus, when backward (i.e., back to the sun) observations are used, the transmission distance for reflected sunlight is short, and the pixel saturation is low. Conversely, when forward (i.e., toward the sun) observations are used, the transmission distance for scattered sunlight is long, and the pixel saturation is high. Therefore, our method preferentially selects forward observation pixels as cloud-free pixels. The incident light of cloud-free pixels is primarily reflected sunlight, whereas that of shadow pixels is primarily light scattered by the surrounding objects. Therefore, the saturation of cloud-shadow pixel is greater than that of cloud-free pixel, and then shadow pixels can be more easily mistaken for cloud-free pixels.
To verify whether our results have been affected by land-use type, we selected two indicators of internal cloud and cloud shadow and used Eq. (12) to calculate P i for residual cloud and shadow for the composited cloud-free image for the water, vegetation, urban, snow, and nonvegetation land cover categories (Table 5). For this, we used land cover data from the MCD12Q1 type 5 dataset. Water bodies, nonvegetated areas, urban and built-up areas, and snow and ice were extracted directly from the dataset; the vegetation category was obtained by combining the evergreen needleleaf trees, evergreen broadleaf trees, deciduous needleleaf trees, and deciduous broadleaf trees classes.
It is clear that, in the cloud-free images from North China (sample A) and the Northwest Territories (sample B), the percentage of residual cloud for the vegetation, built-up, bare ground, and other land use types is small but still greater than corresponding percentages for the MOD09A1 product. In the cloud-free image for the Northwest Territories (sample B), the percentage of residual cloud pixels in snow areas was found to be 43.29% (27.3%) using our method (the MOD09A1 product). Thus, it is clear that the MOD09A1 product performs significantly better than our method in this respect. In the cloud-free image of the coastal areas of southeastern China (sample C), considerably more residual cloud pixels remain over water than in the MOD09A1 product. Moreover, in the Northwest Territories (sample B) and coastal areas of southeastern China (sample C), the number of residual cloud shadow pixels is relatively high with percentages greater than 20%. In the Northwest Territories (sample B), the composited image was produced primarily using the B 1∕7 method; here, P i of the vegetation were found to be significantly lower than those of water, indicating that the B 1∕7 method performs better than the saturation method in removing cloud shadow.

Conclusions and Future Work
Our newly developed integrated cloud-free image compositing algorithm requires only bands 1, 3, 4, and 7 of MODIS images, all of which are readily available, and the algorithm complexity is smaller than MOD09A1. The B 1∕7 method is suitable for extracting cloud-free pixels for land areas, whereas the saturation method is suitable for extracting cloud-free pixels for marine environments, coastal areas, and other water bodies. Our proposed method can eliminate cloud effects and produces slightly higher values for the percentage area ratio of residual cloud pixels and pixels adjacent to cloud in cloud-free images than those in the MOD09A1 product. Moreover, our proposed method performs well in terms of removing cirrus, particularly high cirrus. The calculation accuracy of our method depends on land cover type: the algorithm performs best in the removal of cloud and shadow for vegetation, followed by urban areas, bare land, and water, but performs relatively poorly for snow and ice. Conversely, the MOD09A1 product exhibits no such variation in performance with land cover type. Moreover, our algorithm preferentially selects cloud-free image pixels with larger zenith angles and smaller scattering angles and prefers forward observations to backward observations, which is in stark contrast to the MOD09A1 product. However, our method has some limitations in compositing cloudfree images. For example, the method can only be used to process multiday images obtained in daytime; it also produces more residual cloud pixels than the MOD09A1 product, although it performs similarly to the MOD09A1 product in terms of removing cloud and shadow in ice and snow areas. These limitations need further improvements.