A cloud detection method based on a time series of MODIS surface reflectance images

Abstract

The Moderate Resolution Imaging Spectroradiometer (MODIS)-Terra surface reflectance product (MOD09A1), with bands 1 to 7, is a gridded, eight-day composite product derived from the MODIS-Terra top of atmosphere reflectance swaths. It performs cloud detection and corrects for the effects of atmospheric gases and aerosols. The cloud mask (CM) algorithms for MODIS are based on empirical thresholds on spectral reflectance and brightness temperature. Since the spatial resolution of the thermal band is 1000 m, while that of MOD09A1 is 500 m, many undetected and false clouds are observed in MOD09A1. These errors always result in temporal and spatial inconsistencies in higher-level products. In this paper, a cloud detection algorithm (TSCD) based on a MOD09A1 time series is introduced. Time series cloud detection (TSCD) algorithm is based on the relative stability of ground reflectance and the sudden variations in reflectance that result from cloud cover. The algorithm first searches the clear-sky reference data, and then discriminates clouded and unclouded pixels by detecting a sudden change of reflectance in the blue wavelength and spectral correlation coefficient at the pixel level. Compared with cloud cover assessments obtained from MODIS' original CM, TSCD provides similar or better discrimination in most situations when the land surface changes slowly.

Keywords:

• MODIS
• surface reflectance
• time series
• cloud detection
• BRDF

1. Introduction

Due to the effects of cloud cover, cloud detection is very important in processing satellite remote sensing data. Surface reflectance products are the basis for characterizing the global land surface, and they are often positively biased by undetected clouds. Due to the moderate resolution of The Moderate Resolution Imaging Spectroradiometer (MODIS), discrimination between cloud cover and clear pixels is very difficult.

The current cloud detection methods based on top of atmosphere (TOA) reflectance can be grouped into three categories of methods based on spatial information, spectral information, and temporal information. Methods in the latter two categories are commonly used, since spatial information is difficult to use at full resolution.

Currently operational cloud mask (CM) algorithms for low orbiting sensors, including the Advanced Very High Resolution Radiometer (AVHRR) Clouds from AVHRR (CLAVR) (McClain 1993) and MODIS algorithms (Ackerman et al. 1998, 2006), use threshold-based spectral reflectance and brightness temperature, and perform processing pixel by pixel.

Cloud reflectance is usually higher than ground reflectance in the visible and near infrared (VNIR) spectral region. Thresholds on reflectance in the blue band are suitable for detecting low clouds, but they may fail when the Earth's surface is bright (Bréon and Colzy 1999). Shortwave infrared (SWIR) bands are often used to discriminate snow from clouds: these targets have similar reflectance ranges in the VNIR, but in SWIR snow reflectance is much lower than that of clouds (Dozier 1989). Among the SWIR bands, the 1380-nm band is located in a very strong water vapor absorption band, where only the upper layers of the atmosphere are visible and the background is completely black. This band has been successfully used on MODIS products to detect high clouds (Gao, Goetz, and Wiscombe 1993). Thermal infrared bands are used to detect clouds colder than the Earth surface, which corresponds to almost all types of clouds except thin or low clouds (Saunders and Kriebel 1988; Ackerman et al. 1998). Since then, there have been many improved methods of this kind (Frey and Ackerman 2008; Patrick, Qing, and SunMack 2008; Luo et al. 2008; Liu et al. 2007).

These methods even work well with the thin, cirrus clouds, which are very difficult to detect. However, they use a generic land type classification without a priori knowledge about specific surface reflectance and therefore have to deal with large uncertainties caused by wide natural variability of both land surface and clouds.

On the other hand, the time series data coming from the satellites with a high revisit frequency have provided possibilities to detect clouds using criteria based on time series consistency.

Rossow and Garder (1993) proposed a method that builds a clear-sky composite map from previous measurements, inferring the CM for every pixel by comparing a current measurement with the clear-sky reference value. The uncertainty of the reference value, caused by natural variability and sensor noise, is directly calculated from the measurements.

Another method assumed that the presence of a cloud is likely if a low correlation in the red band is observed at the local scale between two successive images of the same zone (Lyapustin et al. 2008). But it also indicates that a few pixels may be cloudy even in high correlation and may be clear even in low correlation.

Hagolle et al. (2008, 2010) used the blue band reflectance variations in detecting cloud and clear sky on the pixel level. They proposed a threshold related to the number of days between the current day and the most recent clear sky. The day would be labeled as cloudy if the difference of the blue band reflectance between the two dates was greater than the threshold. This method was successfully applied to images from high-resolution sensors, such as FORMOSAT-2, VENµS, Landsat (5 and 7), and Satellite Pour l'Observation de la Terre (SPOT4), whose viewing angles are constant with time.

These time series methods provide a new way to detect clouds. However, modifications are necessary to be applied to moderate resolution images with large and varying viewing angles, especially for MOD09A1 images.

Derived from these methods, a time series cloud detection (TSCD) algorithm is introduced in this paper. The new algorithm will be adopted as a pre-processing step in the “Generation and Application of Global Products of Essential Land Variables” project of China in order to improve the surface reflectance (Liang et al. 2013).

In the following sections, the datasets, new cloud detection method, and experiment results will be described in detail.

2. Data-sets used in the study and problems

MOD09A1 is a surface reflectance product of MODIS-Terra composed of the best observations of Terra during an 8-day period, determined by the overall pixel quality and observational coverage. The product is an estimate of the surface spectral reflectance for each band as it would have been measured at ground level. It corrects for the effects of atmospheric gases and aerosols.

Table 1 lists the important components of MOD09A1. MOD09A1 is a surface reflectance product including seven bands with 500-m resolution. First, it is atmospherically corrected. Second, it is gridded into a map projection (Wolfe, Roy, and Vermote 1998). There is one additional data-set of importance named the 1-km State QA, flagging pixels as land, ocean, or containing cloud or snow. The original cloud mask (MODCM) of each pixel in the State QA includes three parts: cloud state, cirrus, and the internal cloud algorithm flag (Vermote, Kotchenova, and Ray 2008).

Table 1. Important components of MOD09A1.

Band	Bandwidth (nm)	Original resolution (m)	Description
1	620–670	250	Resampled to 500-m resolution
2	841–876	250	Resampled to 500-m resolution
3	459–479	500
4	545–565	500
5	1230–1250	500
6	1628–1652	500
7	2105–2155	500
State		1000

Although MODCM values are provided, there are still some problems in the data flags. First, the resolutions of the image bands are 500 m, but that of MODCM is 1000 m, so there must be some inconsistencies. Second, the algorithm of MODCM is constructed to be cloud-conservative, minimizing false cloud detections but missing clouds that resemble clear conditions (Ackerman et al. 2006). As a result, the processing of the remaining clouds is left to the user.

In the “Generation and Application of Global Products of Essential Land Variables” project, the surface reflectance data are one of the inputs of algorithms of the products including leaf area index (LAI), emissivity, and photo-synthetically active radiation (PAR) . Remaining clouds will result in temporal and spatial inconsistencies of those products. Therefore, it is necessary to refine the cloud detection of MOD09A1 data.

There are three main cases in the MOD09A1 data

1. Clear pixels are labeled as clear. The pixel is atmospherically corrected and provides correct surface reflectance data.

2. Clear pixels are labeled as cloudy. The pixel is atmospherically corrected. It has a correct surface reflectance, though its CM is wrong.

3. Cloudy pixels. The pixel contains atmospherically corrected cloud reflectance. Though it is little different from the TOA cloud reflectance, it is not very different from TOA cloud properties, which still can be used to tell clouds from clear sky.

The quality assurance information has two possible cases

	The pixel was cloudy or other abnormalities, such as polar night data or bad data. The pixel keeps the original value because atmospheric correction is meaningless.
	A clear pixel is labeled as cloudy and remains the original TOA reflectance. We found that this was almost nonexistent after a wide survey.

Since case (4) ought to be discarded or labeled and case (5) seldom exists, we can still make use of the real surface reflectance and atmospherically corrected cloud reflectance to identify cloud and clear sky after cases (4) and (5) are discarded.

Furthermore, since errors in the original cloud masks or abnormalities in surface reflectance data from the preceding processing steps influence subsequent applications, a pre-processing step in the “Generation and Application of Global Products of Essential Land Variables” project is necessary in order to improve the quality of MOD09A1.

MOD09A1 products are gridded and registered with an accuracy better than 0.2 pixels. Therefore, it is feasible to detect clouds using time series methods.

In this study, we selected land data from 2000 to 2010 in a study area including Asia, South America, and North America. The land covers include grasses, forests, deserts, and croplands. Waters and high-latitude snow data were ignored in our study.

Most data were standard MOD09A1 products, and some came from MODIS Land Product Subsets (http://daac.ornl.gov/MODIS/modis.shtml), but with the same processing.

3. Methodology

Generally, the surface reflectance of a land pixel varies with time very slowly. A sharp change in reflectance or spectral correlation could possibly be due to cloud cover or an abnormality.

Hagolle et al. (2008, 2010) calculated the difference in the blue band between the current day and a day with clear sky to detect cloud and clear sky on the pixel level. This method has been successfully applied to high-resolution images of FORMOSAT-2, ENµS, Landsat, and SPOT4 sensors with consistent viewing angles. In this study, the idea of Hagolle et al. (2008, 2010) was adopted for MOD09A1 images, which cover the 4592155-nm spectral ranges, but to apply this methodology to MOD09A1, the bidirectional reflectance distribution function (BRDF) effect had to be considered. Further, the temporal correlation of the spectrum could also be used as a criterion for ground continuity detection, owing to the seven bands of MOD09A1.

There are two main criteria to detect MOD09A1 clouds in TSCD: the threshold for the reflectance increase in band 3 (blue for MOD09A1), and the threshold for the spectral correlation coefficient variation. The reason for selecting the blue band is that the reflectance of land pixels in the blue band is always lower than the other bands and its variation is also slower. On the contrary, the cloud reflectance in the blue band is usually very high. The spectral correlation coefficient ensures consistency in the spectral domain.

In order to compute the variations and to detect clouds for a pixel on date D, a cloud-free reference pixel is needed. For each date D, the clear reference pixel contains the most recent cloud-free reflectance obtained in the time series before date D. A pixel is labeled as clear for the time series criterion if:

(1)

(2)

where ρ_b3(D) is the surface reflectance of the pixel in the blue band at date D; D_r is the date of the most recent clear data before D; D–D_r is the number of days between them; and corr(D, D_r) is the normalized spectral correlation coefficient, where θ_D,Dr is the spectral angle between ρ(D) and ρ(D_r). K₁, L₁, K₂, and L₂ are coefficients whose threshold values depend on D–D_r. K₁ and K₂, which characterize the change of the spectra and correlation coefficient with time, are related to the land type and season alteration. L₁ and L₂ represent the effects of BRDF and noise.

For original data from a routine day, the solar zenith and view zenith vary greatly with satellite orbits. The BRDF has a great effect on the surface reflectance, especially for data with a high resolution. Due to the optimizing selection of MOD09A1, the effects of solar zenith and view zenith variation are weakened. Therefore, the selection of L₁ was simplified in our project as follows:

(3)

The 500-m resolution of MOD09A1 results in mixed pixels. Since D is close to D_r, the selection of K₁ was also simplified in our project as follows:

(4)

corr is an assistant criteria to estimate clouds. We found that the influence of BRDF on corr was smaller than 0.01. At the same time, D was close to D_r. Therefore, the threshold of corr was selected as 0.96.

If these criteria are used, attention should be paid to the following points:

3.1. Label snow pixels including those under clouds

Clouds are apparently similar to snow. Our objective in the project was to discriminate clouds and non-snow objects, so we first needed to label the snow pixels regardless of cloud cover. For MOD09A1 data, the reflectance of snow and atmospherically corrected clouds has the following features:

• a high-reflectance value in the VNIR region (bands 1, 2, 3, and 4), and

• the snow reflectance in the SWIR region (bands 5 and 6) is much lower than that of clouds.

As a result, discriminating clouds from snow with band 3 was noneffective, as was corr. However, our objective was to discriminate clouds from non-snow objects, so our scheme on snow processing was simple. Based on the MOD09A1 snow mask, we used information on the spatial and temporal continuity of snow and followed these criteria to refine the snow mask:

(5)

These criteria do well for pure snow, but not for mixed snow. Spatial and temporal knowledge of the two targets can be used to distinguish them. For example, snow will not melt until spring in high-latitude areas, so we can obtain the transition of pure snow turning to mixed snow according to temporal continuity.

3.2. Find clear sky

There are two kinds of clear sky: the first clear sky in the time series, and all subsequent clear skies.

The conditions for finding the first clear sky are rigorous. Data of two successive dates (D₁ and D₂, where each is the best among eight days) should satisfy the following criteria in order to be selected as clear sky:

• the data difference in band 3 is smaller than 0.025;

• the average value in band 3, , is smaller than 0.35;

• there is no band 3 value smaller than within date [D₁−30, D₂+30];

• there are more band 3 values located in within date [D₂, D₂+30];

• the correlation coefficient of spectra is smaller than 0.96.

Once the first clear sky is determined, the subsequent clear sky is easier to obtain. In order to save computing time, a subsequent pixel will be labeled as clear if Equations (1) and (2) are satisfied. Then the reference clear sky will be substituted with the new clear sky. Searching for subsequent clear sky will continue until a new clear sky could not be found for more than 80 days. At this time, the object may have changed suddenly because of land cover change or because of cropping, for instance, so the idea is to restart the search for a new first clear sky since the last clear sky.

Since the criteria for the first clear sky are strict, some clear skies may be skipped. So we search backward from the first clear sky just like searching for subsequent clear skies.

3.3. Processing high-reflectance objects

For MOD09A1, the band 3 values of most objects are smaller than 0.3, and few objects exceed 0.3. Sand, for example, may sometimes exceed 0.35. However, clear sky with a band 3 value smaller than 0.35 can always be found. Once found, the skipped clear skies with a band 3 value higher than 0.35 can be recruited by searching backward.

Based on the continuity of MODIS band 3 values and the spectral correlation coefficient, we can identify clouds and clear sky pixels through a time series of data. To summarize, the general flowchart of the TSCD algorithm is given in Figure 1.

1. As the first step, find and discard snow and abnormal data such as noise bands and uncorrected data.

2. Beginning from day D₀, find the first clear sky pair at days D₁ and D₂.

3. Search backward for clear skies from days D₁ to D_0.

4. Beginning at D₂, search the subsequent clear sky and substitute the new clear sky for the reference clear sky until a new clear sky cannot be found for more than 80 days.

5. If the interval exceeds 80 days, search again for a new first clear sky pair since the last clear sky again. Let the last clear sky date be the new D₀, and repeat steps (2) through (4).

6. Finally, the TSCD CM is output.

Figure 1. The general flowchart of TSCD algorithm.

4. Results and discussion

The TSCD algorithms have been applied in the “Generation and Application of Global Products of Essential Land Variables” project. A large dataset from 2000 to 2010 was used to test our algorithms. The dataset includes about 100 tiles located in middle- and low-latitude areas. There are more than 400 images for each tile. The size of each image is 2400 by 2400.

There is no reliable independent CM source at a given hour at the given position. In order to assess the algorithm performance, a visual comparison of time series cloud detection cloud masking (TSCDCM) against the true-color MODIS Red, Green, Blue (RGB) image as well as the MODIS CM (MODCM) was conducted. The agreement between TSCDCM and MODCM was obvious in most areas, though sometimes, TSCDCM seemed better than MODCM. Since thresholds in the MODCM method are empirically adjusted, thin clouds are often labeled as clear, which leads to a fragmentized CM in continuous clouds.

For sake of conciseness, only part of the results is presented here. First, a subset of the time series data was used to test the algorithm. Then, several typical tile data, such as h27v06 in 2007, h11v04 in 2005, h13v09 in 2008, and h22v07 in 2008, were processed as a test dataset.

Two time series datasets of Changbaishan, a Jilin Province subset, and one time series subset of Beijing, from March 2001 to the end of 2004 were selected. Figures 2–4 show the cloud detection results. The x-coordinate represents about 180 dates (at 8 days each), while the y-coordinate represents the band 3 values, MODCM, and TSCDCM. The ‘+’ at value 0.1 represents MODCM, and the ‘x’ at value 0.05 represents TSCDCM. TSCDCM generally agreed with MODCM, and it was sometimes even better than MODCM. The differences between them are represented by ‘A’, 'B’, and 'C’ in the figures. ‘A’ represents the summer clouds whose MODIS flag is clear; ‘B’ represents pixels with a cloudy MODCM flag but a clear TSCDCM flag (these pixels are almost clear judging from the comparison with pixels of nearby times in the current year or same time in other years); and ‘C’ represents pixels that may be clouds and for which the MODCM flag is cloudy and the TSCDCM flag is clear. Their influences on surface reflectance are less than 0.02, which is equivalent with the noise level.

Figure 2. One time series (March of 2001 to the end of 2004) of the Changbaishan subset, x-coordinate represents about 180 dates (each per 8 days), y-coordinate represents the band 3 values, MODCM and TSCDCM. ‘+’ at value 0.1 represents MODCM, and ‘x’ at value 0.05 represents TSCDCM.

Figure 3. Time series (March of 2001 to the end of 2004) of the Changbaishan subset.

Figure 4. Time series (March of 2001 to the end of 2004) of the Beijing subset.

Figures 5–10 are visual comparisons of TSCDCM with the true-color MODIS RGB image as well as the MODCM. Figures 5 and 6 are time series of cloud detection results for one year. The left column is the RGB images, the middle column is the TSCDCM, and the right column is the MODCM. White represents clouds. Figures 7 and 8 show the details of days 129 and 193. In Figure 7, TSCDCM might underestimate the cloud cover in the top-right corner where haze and thin clouds exist. In other places, TSCDCM seems better than MODCM because visibly continuous clouds appear fragmentized in MODCM.

Figure 5. A subset of h27v06 from 009 day to 361 day for 2007, size 200×200.

Figure 6. A subset of h11v04 form 009 day to 361 day for 2005, size 200×200.

Figure 7. A subset of h27v06 for date 129 in 2007, size 200×200.

Figure 8. A subset of h27v06 for date 193 in 2007, size 200×200.

Figure 9. Two subsets of h22v07 for date 145 (left) and 193 (right) in 2008 size 200×200, the two RGB images have a different normalization for helping visual distinction between snow and Clouds.

Figure 10. A subset of h13v09 for date 201 in 2008, size 400×400.

The result may be incorrect on snow because of simplified treatment. However, our project aimed to label the snow regardless of cloud cover on it, so we used the MODIS snow mask and sr_NDSI criteria together with spatial and temporal continuity information to settle this problem.

The algorithm does not work well above inland water surfaces because of sun glint, turbidity, and foam. Besides this, thin clouds and high aerosol optical thicknesses might be confused – some clouds might be too thin to be detected. The threshold of the band 3 value <0.35 might also result in some overestimated clouds, which was described in section 3.

The x-coordinate represents about 180 dates (eight days each); the y-coordinate represents the band 3 values, MODCM, and TSCDCM; the ‘+’ at value 0.1 represents MODCM; and ‘x’ at value 0.05 represents TSCDCM.

5. Conclusion

A TSCD method was developed for MOD09A1 products. The TSCD method was successfully applied to a project in China titled “Generation and Application of Global Products of Essential Land Variables.” Compared to heritage CM algorithms, which are empirical, the TSCD method is self-adaptive to practical objects on the basis of their changing rules.

The MODCM served as a reference in the comparison with TSCDCM in this study. TSCDCM was mostly consistent with MODCM except on some occasions. Generally, MODCM is always more conservative in order to avoid false cloud detection, but it results in clouds being undetected. In some cases, TSCDCM classifies thick aerosol plumes as clouds, but, as heavy aerosol plumes cannot be corrected with good accuracy, the confusion of aerosol plumes and clouds was not a problem in this project.

The TSCD method has provided good performance in most cases, especially in low-latitude forests and middle-latitude grasses and crops. Sometimes, TSCDCM was better than MODIS CM. On the other hand, the TSCD method might overestimate clouds on water, bright objects, and other objects that change suddenly, and it might underestimate thin clouds, which tend to be confused with noises and uncorrected aerosol plumes.

Improvement of the TSCD algorithm is still in progress. For example, uncorrected, thick (opaque) aerosol plumes are often confused with clouds. The coefficients in formulas (1) and (2) could be processed in much more detail in future versions according to land cover, season, zone, and so on. It is certain that appropriate parameters exist for each instance.

Currently, the method is being applied as a post-processing procedure of MOD09A1. At the same time, its principles are also suitable for other multi-spectral data as well as TOA reflectance, as long as the thresholds are reconsidered.

Funding

This work was supported by the “Generation and Application of Global Products of Essential Land Variables” project of China [grant number 2009AA122100] under the “State Program for High-tech Research and Development” (863 Program).

Additional information

Funding

Funding: This work was supported by the “Generation and Application of Global Products of Essential Land Variables” project of China [grant number 2009AA122100] under the “State Program for High-tech Research and Development” (863 Program).