Cloudy Region Drought Index (CRDI) Based on Long-Time-Series Cloud Optical Thickness (COT) and Vegetation Conditions Index (VCI): A Case Study in Guangdong, South Eastern China

: Widespread and long-lasting drought disasters can aggravate environmental degradation. They can lead to signiﬁcant economic losses and even a ﬀ ect social stability. The existing drought index mostly chose arid and semi-arid regions as study areas, because cloudy weather in humid and semi-humid regions hindered the satellite in its attempts to obtain the surface reﬂectivity. In order to solve this problem, a cloudy region drought index (CRDI) is proposed to estimate the drought of the clouded pixels. Due to the cumulative e ﬀ ect of drought, the antecedent drought index (ADI) has a certain impact on the calculation of the current drought. Furthermore, cloud is the only source of natural precipitation, and it also a ﬀ ects the evaporation and emission process on the ground. Therefore, based on the remote sensing drought index, ADI and cloud optical thickness (COT) are used to estimate the drought of pixels with missing data due to cloud occlusion. In this paper, a case study of the cloudy Guangdong, which is located in a humid area, is presented. First, we calculated the CRDI using Moderate Resolution Imaging Spectroradiometer (MODIS) data from 2003 to 2017, and then discussed the e ﬀ ect of CRDI with the data from 2016 as examples. Through the analysis of the parameters of regression equation, ﬁlling e ﬃ ciency, rationality of the estimated value, the continuity of CRDI and the rationality of CRDI spatial distribution results, it is concluded that CRDI can e ﬀ ectively estimate the drought severity of the cloud-covered pixels, and more comprehensive drought data can be obtained by using CRDI. The successful application of CRDI in Guangdong shows it is robust and ﬂexible, suggesting high e ﬃ ciency and great potential for further utilization.


Introduction
Drought, as a serious natural disaster, has a serious impact on social economy and people's lives [1]. It is estimated that the economic losses caused by drought around the world are far more than for other meteorological disasters [2].
There are many ways to monitor droughts. For a single sample, we can obtain soil moisture data through a soil moisture monitoring station [3][4][5][6]. Meteorological station data can also be used to monitor drought by calculating drought index, such as from the Palmer drought severity index (PDSI), the standardized precipitation index (SPI), the surface water supply index (SWSI) and so on [7][8][9][10][11]. We can also observe the growth of single plants to determine whether they are suffering from

.2. Data Source
The Moderate Resolution Imaging Spectroradiometer (MODIS) is a sensor mounted on Aqua nd Terra. It provides observation data of the earth's surface every 1-2 days, and the minimum spatial esolution of the observation data is 250 m. MODIS data have moderate resolution, long coverage ime and a short revisit period and have great potential for large-scale drought monitoring. The main ata used in this paper are MODIS surface reflectance data (MOD09 and MYD09) with 250 m esolution in 8 days, COT data (MOD08 and MYD08) with 1 degree resolution in 8 days, and land se/land cover (LULC) data (MCD12Q1) with 500 m resolution.
In this paper, the MODIS surface reflectance and COT products with a period of 8 days from 003 to 2017 are used. We decided to use 8-day data since droughts are generally persistent events ith varying durations, and most of the drought events can be monitored using 8-day data. For COT ata, the best value of different satellite synchronization data was selected. For the surface reflectance ata, after selecting the best value, the basic drought index VCI should be calculated by using the urface reflectance values of red and near infrared bands (NIR). In addition, land use/land cover data elected in this paper are the international geosphere biosphere programme (IGBP) products that can efine the surface classification of ecosystems. The research object of this paper is the area with egetation cover. In the subsequent use of land cover data, it is necessary to integrate the original and cover data of the IGBP classification scheme (Table 1). For the convenience of analysis and alculation, we resample the low-resolution COT and land use data to 250 m, so as to establish the ne-to-one correspondence relationship between different data sets.

Category
Proportion (%) Category Corresponding to IGBP In this paper, the MODIS surface reflectance and COT products with a period of 8 days from 2003 to 2017 are used. We decided to use 8-day data since droughts are generally persistent events with varying durations, and most of the drought events can be monitored using 8-day data. For COT data, the best value of different satellite synchronization data was selected. For the surface reflectance data, after selecting the best value, the basic drought index VCI should be calculated by using the surface reflectance values of red and near infrared bands (NIR). In addition, land use/land cover data selected in this paper are the international geosphere biosphere programme (IGBP) products that can define the surface classification of ecosystems. The research object of this paper is the area with vegetation cover. In the subsequent use of land cover data, it is necessary to integrate the original land cover data of the IGBP classification scheme (Table 1). For the convenience of analysis and calculation, we resample the low-resolution COT and land use data to 250 m, so as to establish the one-to-one correspondence relationship between different data sets.

Methodology
A new method monitor drought in a cloudy area using ADI and COT was proposed. It included the following procedures: calculation of VCI and ADI, discussion on the correlation between VCI and ADI, establishment of regression equation and construction of a CRDI model.

Calculate VCI and ADI
NDVI is the most widely used vegetation index, which is related to vegetation coverage and mainly reflects the physiological and ecological characteristics of vegetation and also reflects the background impact of the plant canopy, such as soil, ground tide, snow, dead leaves, coarseness, etc. [37].
NDVI was first proposed by Rouse [38] based on the following equation: where NIR stands for the reflectivity of near infrared band, and R stands for the reflectivity of red band. The range of NDVI is between −1 and 1, the negative value indicates that the ground is covered by clouds, water, snow, etc.; 0 indicates that there is rock or bare soil, etc.; the positive value indicates that there is vegetation coverage, and its NDVI increases with the increase of vegetation coverage. VCI is a remote sensing index based on NDVI. It is mainly used to reflect the relative growth of vegetation. VCI was first proposed by Kogan [39] and based on the following equation: where, i ranges from 1 to 46 and is represented for the current period in a year. NDVI i represented the NDVI value of period i in a specific year. NDVI max and NDVI min represented the maximum and minimum NDVI values of period i of all years and reflect the best and worst conditions of vegetation growth, respectively. Furthermore, the difference represents the habitat of vegetation, such as soil moisture, nutrition and sunshine. The smaller the difference between NDVI i and NDVI min , the worse the vegetation growth in the current period [39]. The range of VCI is between 0 and 100. The smaller VCI is, the worse the vegetation condition is, and the more serious the drought is [40,41]. Due to the cumulative effect of drought, the antecedent drought has a certain impact on the current drought [30]. Antecedent drought used in this paper is obtained from VCI based on the following equation: where j represented the current period. The minimum value of j is 2, and the maximum value of j is the length of our data. We used 15 years of data, and there are 46 issues per year, so the length of our data is 690.

Correlation between VCI and ADI
The relationship between VCI and ADI of forest, grassland and agricultural land ( Figure 2) was discussed by using the data of 2016. The points in the scatter diagram are randomly selected. It can be seen that the scattering points are basically concentrated in the upper right corner, and VCI is positively correlated with its ADI. However, for the scatter diagram of the forest, some scatter points are concentrated in the upper left and lower right. This is due to the area and species of woodland in Guangdong. Some of the points concentrated in the upper left and lower right may be related to different forest types. However, most of the points are still concentrated in the upper right corner, and there is a positive correlation on the whole. For three different land cover types, the aggregation of points in the scatter diagram are similar. However, whether the same equation can be used for different land cover types needs further investigation. Remote Sens. 2020, 12, x FOR PEER REVIEW 5 of 19 In order to explore the differences of scatter plot of different land cover types and different periods, the scatter plots of random points of the 19th, 33rd and 43rd period were selected from the 46 periods in 2016. Furthermore, period 19, 33 and 43 correspond to late-April, mid-September and early-December in 2016, respectively. The points in the scatter plot tend to be concentrated in the upper right corner, but it is obvious that the aggregation of different land cover types and different periods is different ( Figure 3). Therefore, regression equations should be constructed for different land cover types and different periods.

Estimation Method of Cloudy Pixel
In this paper, ADI and COT are used to estimate the drought value of the pixel in a cloudy region. Using ADI and COT to estimate the drought value of pixels in cloudy regions mainly includes the following steps. Taking the agricultural region as an example, firstly, the agricultural region in current period is extracted, and the regression equation is established according to the COT, ADI and In order to explore the differences of scatter plot of different land cover types and different periods, the scatter plots of random points of the 19th, 33rd and 43rd period were selected from the 46 periods in 2016. Furthermore, period 19, 33 and 43 correspond to late-April, mid-September and early-December in 2016, respectively. The points in the scatter plot tend to be concentrated in the upper right corner, but it is obvious that the aggregation of different land cover types and different periods is different ( Figure 3). Therefore, regression equations should be constructed for different land cover types and different periods. In order to explore the differences of scatter plot of different land cover types and different periods, the scatter plots of random points of the 19th, 33rd and 43rd period were selected from the 46 periods in 2016. Furthermore, period 19, 33 and 43 correspond to late-April, mid-September and early-December in 2016, respectively. The points in the scatter plot tend to be concentrated in the upper right corner, but it is obvious that the aggregation of different land cover types and different periods is different ( Figure 3). Therefore, regression equations should be constructed for different land cover types and different periods.

Estimation Method of Cloudy Pixel
In this paper, ADI and COT are used to estimate the drought value of the pixel in a cloudy region. Using ADI and COT to estimate the drought value of pixels in cloudy regions mainly includes the following steps. Taking the agricultural region as an example, firstly, the agricultural region in current period is extracted, and the regression equation is established according to the COT, ADI and

Estimation Method of Cloudy Pixel
In this paper, ADI and COT are used to estimate the drought value of the pixel in a cloudy region. Using ADI and COT to estimate the drought value of pixels in cloudy regions mainly includes the following steps. Taking the agricultural region as an example, firstly, the agricultural region in current period is extracted, and the regression equation is established according to the COT, ADI and VCI Remote Sens. 2020, 12, 3641 6 of 19 corresponding to the cloudless pixels, and then the parameters in the regression equation are obtained. Finally, according to these parameters, the drought value of cloudy pixels in agricultural region is estimated and recorded as DI. The equation for estimating DI is as follows: where, a is the constant term, b is the coefficient of ADI, c is the coefficient of COT and τ is the residual of the regression equation. Finally, it is necessary to judge whether the estimated value DI is reasonable. VCI ranges from 0 to 100, so the boundary of DI should also be within 0 and 100. When the value of DI exceeds the boundary but does not exceed 5%, DI is assigned a boundary value. When the estimated DI value exceeds 5% of the boundary value, the current pixel is considered to be inestimable. After determining the regression equation of each period and different land cover types, we can further obtain the CRDI of each pixel according to the flow chart ( Figure 4). Furthermore, the specific equation of CRDI is as follows: CRDI = VCI , cloudless region DI(COT, ADI), cloudy region (5) the basic remote sensing drought index VCI is still used in cloudless pixels; for the cloudy pixels which are unable to obtain drought data, the estimated drought value was used. Among them, DI(COT, ADI) is shown in Equation (4).
obtained. Finally, according to these parameters, the drought value of cloudy pixe region is estimated and recorded as DI. The equation for estimating DI is as follow where, a is the constant term, b is the coefficient of ADI, c is the coefficient of COT an of the regression equation. Finally, it is necessary to judge whether the estim reasonable. VCI ranges from 0 to 100, so the boundary of DI should also be within the value of DI exceeds the boundary but does not exceed 5%, DI is assigned a boun the estimated DI value exceeds 5% of the boundary value, the current pixel is inestimable.
After determining the regression equation of each period and different land co further obtain the CRDI of each pixel according to the flow chart ( Figure 4). Further equation of CRDI is as follows: CRDI = VCI , cloudless region DI(COT, ADI), cloudy region the basic remote sensing drought index VCI is still used in cloudless pixels; for which are unable to obtain drought data, the estimated drought value was use DI(COT, ADI) is shown in Equation 4.
According to the regression equation, parameters corresponding to each perio type can be obtained, including constant term, coefficient of COT, coefficient of A Therefore, it is possible to determine the form of specific regression equations for land cover type ( Figure 4). The current estimated DI can also be used as ADI estimation.

Results
Taking the data of 2016 as an example, the parameters of regression equation, validity of estimated value, the continuity of CRDI and the rationality of CRDI sp results are analyzed. According to the regression equation, parameters corresponding to each period and land cover type can be obtained, including constant term, coefficient of COT, coefficient of ADI and residual. Therefore, it is possible to determine the form of specific regression equations for each period and land cover type (Figure 4). The current estimated DI can also be used as ADI for the next DI estimation.

Results
Taking the data of 2016 as an example, the parameters of regression equation, the efficiency and validity of estimated value, the continuity of CRDI and the rationality of CRDI spatial distribution results are analyzed.

Analysis of VCI Efficiency
According to the statistics of total VCI data, it is found that most of pixels are missing 10-70 periods of 690 periods (Figure 5a). Among them, most of the pixels lost 15-35 periods. There were few pixels with less than 10 or more than 70 missing periods. There are fewer missing data in eastern Guangdong and Leizhou Peninsula, which are basically within 20 periods, while there are more missing data in northern Guangdong, ranging from 40 to 70 periods.

Analysis of VCI Efficiency
According to the statistics of total VCI data, it is found that most of pixels are missing 10-70 periods of 690 periods (Figure 5a). Among them, most of the pixels lost 15-35 periods. There were few pixels with less than 10 or more than 70 missing periods. There are fewer missing data in eastern Guangdong and Leizhou Peninsula, which are basically within 20 periods, while there are more missing data in northern Guangdong, ranging from 40 to 70 periods.
The loss rate of each period was calculated. Periods with a loss rate greater than 70% are not available in arid land analysis; the quality of the data with a loss rate of less than 2% is better, with a total of 503 periods (Figure 5b). The periods with a loss rate of more than 2% have a greater impact on the overall interpretation of the drought situation, there are 187 periods of such data, accounting for more than 27% of the total periods. The main problem to be solved in this paper is the estimation of missing pixels in such periods.

Analysis of Parameters
The parameters in Equation 4 can be obtained by the regression equation. The statistical values of parameters of forest, forest grass mixture and agricultural land in 2016 have certain regularity ( Table 2). Among the three land cover types, the average, median, variance and standard deviation of the P value and mean residual value are small. The absolute residual is around 10. It shows that 10% of them cannot be explained by the parameters in this paper. The average and median of ADI coefficients are similar, and the variance and standard deviation are very small. The ADI coefficient has little difference in different land cover types and different periods. The average and median values of COT coefficients were not significantly different, but the variance and standard deviation of COT coefficients were larger, which indicated that the coefficient of COT varied greatly among different land cover types or periods. Compared with other types of forest, the average and median values of constant and mean residual of forest are quite different. The regression parameters of forest grass mixture and agricultural are better, because the vegetation growth in agricultural land is relatively regular. The vegetation coverage rate of forest grass mixed type is relatively low, and the difference within the vegetation condition is small. The loss rate of each period was calculated. Periods with a loss rate greater than 70% are not available in arid land analysis; the quality of the data with a loss rate of less than 2% is better, with a total of 503 periods ( Figure 5b). The periods with a loss rate of more than 2% have a greater impact on the overall interpretation of the drought situation, there are 187 periods of such data, accounting for more than 27% of the total periods. The main problem to be solved in this paper is the estimation of missing pixels in such periods.

Analysis of Parameters
The parameters in Equation (4) can be obtained by the regression equation. The statistical values of parameters of forest, forest grass mixture and agricultural land in 2016 have certain regularity ( Table 2). Among the three land cover types, the average, median, variance and standard deviation of the P value and mean residual value are small. The absolute residual is around 10. It shows that 10% of them cannot be explained by the parameters in this paper. The average and median of ADI coefficients are similar, and the variance and standard deviation are very small. The ADI coefficient has little difference in different land cover types and different periods. The average and median values of COT coefficients were not significantly different, but the variance and standard deviation of COT coefficients were larger, which indicated that the coefficient of COT varied greatly among different land cover types or periods. Compared with other types of forest, the average and median values of constant and mean residual of forest are quite different. The regression parameters of forest grass mixture and agricultural are better, because the vegetation growth in agricultural land is relatively regular. The vegetation coverage rate of forest grass mixed type is relatively low, and the difference within the vegetation condition is small.
The histogram of parameters of the regression equation for each period of 2016 and different land cover type has certain regularity ( Figure 6). Most of the parameters conform to normal distribution. Among them, the distribution of the constant term and absolute residual of the forest is different from other land cover types and other parameters, and they are not normal distributions. This may be due to the different sensitivity of different forest types to drought. Compared with the other two land cover types, the deviation of the regression estimation value of forest land may be relatively large.  The histogram of parameters of the regression equation for each period of 2016 and different land cover type has certain regularity ( Figure 6). Most of the parameters conform to normal distribution. Among them, the distribution of the constant term and absolute residual of the forest is different from other land cover types and other parameters, and they are not normal distributions. This may be due to the different sensitivity of different forest types to drought. Compared with the other two land cover types, the deviation of the regression estimation value of forest land may be relatively large.

Effectiveness Analysis of Filling
According to the filling situation of CRDI to the original VCI in 2016, the filling rate of each period was calculated (Table 3). Among them, the average filling rate of total data, forest, forest grass mixed and agricultural were 98.0%, 99.1%, 97.5% and 99.5%, respectively. The filling rate of each period is different, but the difference is not significant.
There are 31 periods with a filling rate greater than 98%, and 10 periods between 95% and 98% of total data. There are 41, 28 and 42 periods with a filling rate greater than 98%, and 2, 12 and 3 periods between 95% and 98% of forest, forest grass mixed and agriculture, respectively. The filling efficiency of these periods is quite good. The filling rate of less than 80% is poor, which was not found in 2016. The CRDI model can fill in the missing pixels of the traditional remote sensing drought index.

Effectiveness Analysis of Filling
According to the filling situation of CRDI to the original VCI in 2016, the filling rate of each period was calculated (Table 3). Among them, the average filling rate of total data, forest, forest grass mixed and agricultural were 98.0%, 99.1%, 97.5% and 99.5%, respectively. The filling rate of each period is different, but the difference is not significant.
There are 31 periods with a filling rate greater than 98%, and 10 periods between 95% and 98% of total data. There are 41, 28 and 42 periods with a filling rate greater than 98%, and 2, 12 and 3 periods between 95% and 98% of forest, forest grass mixed and agriculture, respectively. The filling efficiency of these periods is quite good. The filling rate of less than 80% is poor, which was not found in 2016. The CRDI model can fill in the missing pixels of the traditional remote sensing drought index.

Rational Analysis of the Estimated Value
It is not comprehensive to evaluate the CRDI model from the perspective of filling efficiency, and whether the distribution of estimated values is reasonable should be considered. We took the CRDI in 2016 as the current CRDI, VCI in 2016 as the original VCI, and the CRDI in 2015 as the previous CRDI ( Figure 7). The histogram distribution of different land cover types was analyzed and the rationality of the estimated value of CRDI is discussed through the histogram. The histogram distribution of CRDI in 2016 for different land cover type is similar to that of the original VCI and the previous CRDI.

Rational Analysis of the Estimated Value
It is not comprehensive to evaluate the CRDI model from the perspective of filling efficiency, and whether the distribution of estimated values is reasonable should be considered. We took the CRDI in 2016 as the current CRDI, VCI in 2016 as the original VCI, and the CRDI in 2015 as the previous CRDI (Figure 7). The histogram distribution of different land cover types was analyzed and the rationality of the estimated value of CRDI is discussed through the histogram. The histogram distribution of CRDI in 2016 for different land cover type is similar to that of the original VCI and the previous CRDI.
In 2016, the total number of filled pixels accounted for 4.4% of the effective pixels, and the estimated CRDI values were mainly distributed between 40 and 80. The number of filled pixels in forest type accounted for 5.4% of the total number of forest pixels, and the estimated values in the CRDI model mainly ranged from 40 to 90. In the forest grass mixed type and agricultural type, the number of filled pixels accounted for 4.3% and 2.8% of the total number of forest grass mixed and agricultural areas respectively, and the estimated values in the CRDI model were mainly distributed between 40 and 80. The CRDI in 2016 did not change the overall distribution of the original VCI, and the range of estimated values in the CRDI model was reasonable.

Continuity Analysis of CRDI
The filling efficiency of the CRDI model is quite good, and the range distribution of the estimated value is also reasonable. In order to further explore the effect of the CRDI model, it is necessary to investigate the spatial and temporal continuity of CRDI. In the study area, the CRDI values of the pixels passed by the section line and the time series curve of the sample points were extracted. On In 2016, the total number of filled pixels accounted for 4.4% of the effective pixels, and the estimated CRDI values were mainly distributed between 40 and 80. The number of filled pixels in forest type accounted for 5.4% of the total number of forest pixels, and the estimated values in the CRDI model mainly ranged from 40 to 90. In the forest grass mixed type and agricultural type, the number of filled pixels accounted for 4.3% and 2.8% of the total number of forest grass mixed and agricultural areas respectively, and the estimated values in the CRDI model were mainly distributed between 40 and 80. The CRDI in 2016 did not change the overall distribution of the original VCI, and the range of estimated values in the CRDI model was reasonable.

Continuity Analysis of CRDI
The filling efficiency of the CRDI model is quite good, and the range distribution of the estimated value is also reasonable. In order to further explore the effect of the CRDI model, it is necessary to investigate the spatial and temporal continuity of CRDI. In the study area, the CRDI values of the pixels passed by the section line and the time series curve of the sample points were extracted. On the one hand, the trends of the curve itself were investigated, and on the other hand, the CRDI values of the current, previous and the same period previously were compared.
In the field of remote sensing, it is usually necessary to reconstruct the vegetation index curve, especially the time series curve. This is because some noises in remote sensing data will make the data fluctuation unneglectable [42]. There are many methods to reconstruct remote sensing data. We choose the WS (Whittaker smoother) method which was introduced by Eilers P.H. in 2003 [43].WS has the advantages of fast and convenient, simple parameters and it is easy to implement [44][45][46]. It has been widely used in the analysis of time series curve in the remote sensing field.
We selected the 36th period (early September) of 2016 as the sample data. The loss rate of VCI in this period is at a medium level. There is no significant difference between the original VCI and the previous VCI. Therefore, the estimated value of this period should be basically consistent with the corresponding previous and contemporaneous values. The reference value of this period of data is good.
The selected section line should pass through the missing pixels of the original VCI as much as possible and should also include some pixels that are not missing in the original VCI. Line AA' and BB' are located in the north and east of Guangdong respectively (Figure 8). The width of the section line is defined as five pixels, and the average value is taken as the value of the ordinate. When three or more pixels are missing, it is considered that the current pixel is missing. The sample points should be evenly distributed in Guangdong and represent different land use types. Point C is located in Leizhou Peninsula, point D is located in the west of Guangdong, and point E is located in the central part of Guangdong. The width of a single point is defined as 3 × 3 pixels, and the average value is taken as the value of ordinate. When there are five or more missing pixels, it is considered that the current pixel is missing. The first period data is taken as the starting point, and the abscissa is the number of periods of the current data.
Remote Sens. 2020, 12, x FOR PEER REVIEW 10 of 19 the one hand, the trends of the curve itself were investigated, and on the other hand, the CRDI values of the current, previous and the same period previously were compared.
In the field of remote sensing, it is usually necessary to reconstruct the vegetation index curve, especially the time series curve. This is because some noises in remote sensing data will make the data fluctuation unneglectable [42]. There are many methods to reconstruct remote sensing data. We choose the WS (Whittaker smoother) method which was introduced by Eilers P. H. in 2003 [43].WS has the advantages of fast and convenient, simple parameters and it is easy to implement [44][45][46]. It has been widely used in the analysis of time series curve in the remote sensing field.
We selected the 36th period (early September) of 2016 as the sample data. The loss rate of VCI in this period is at a medium level. There is no significant difference between the original VCI and the previous VCI. Therefore, the estimated value of this period should be basically consistent with the corresponding previous and contemporaneous values. The reference value of this period of data is good.
The selected section line should pass through the missing pixels of the original VCI as much as possible and should also include some pixels that are not missing in the original VCI. Line AA' and BB' are located in the north and east of Guangdong respectively (Figure 8). The width of the section line is defined as five pixels, and the average value is taken as the value of the ordinate. When three or more pixels are missing, it is considered that the current pixel is missing. The sample points should be evenly distributed in Guangdong and represent different land use types. Point C is located in Leizhou Peninsula, point D is located in the west of Guangdong, and point E is located in the central part of Guangdong. The width of a single point is defined as 3 × 3 pixels, and the average value is taken as the value of ordinate. When there are five or more missing pixels, it is considered that the current pixel is missing. The first period data is taken as the starting point, and the abscissa is the number of periods of the current data. Among all the pixels passing through line AA', the missing pixels account for 68.1%. The trend in the current CRDI is similar to that of the previous CRDI, but the position of peak valley is slightly offset (Figure 9a). Line AA' was divided into three parts: A1, A2 and A3. The difference between the Among all the pixels passing through line AA', the missing pixels account for 68.1%. The trend in the current CRDI is similar to that of the previous CRDI, but the position of peak valley is slightly offset (Figure 9a). Line AA' was divided into three parts: A1, A2 and A3. The difference between the current and previous CRDI of A1 is quite small, and the trend is similar. The CRDI value of A2 in the previous is higher than that in the current period, but the trend is still very similar, especially in the missing pixels of the original VCI data. Instead, the pixels where the data are not missing are different. The current CRDI of A3 is similar to that of the previous CRDI, but the trend of them is slightly different. In the area with more missing data in the middle section of A3, the details of the previous CRDI are more and the fluctuation is more obvious, while the current CRDI lacks some details here, and the trend is relatively stable. slightly different. In the area with more missing data in the middle section of A3, the details of the previous CRDI are more and the fluctuation is more obvious, while the current CRDI lacks some details here, and the trend is relatively stable.
Among all the pixels passed by line BB', 59.1% were missing pixels. The trend of the current CRDI is similar to that of the previous CRDI, the peak valley position is slightly shifted (Figure 9b). However, in some areas, there is no obvious law of the two curves. Line BB' was divided into four parts: B1, B2, B3 and B4. There was no significant difference between the current CRDI and previous CRDI of B1. The trend of the front and middle segments of B1 is similar. The previous CRDI of the tail segment has more details, and the curve has some fluctuations. As a whole, the current CRDI of B2 is higher than that of the previous CRDI. It can be seen that the floating range of the two curves is small. CRDI fluctuated more in the previous period, but it was more stable in the current period. Furthermore, the current CRDI is also less detailed than the previous CRDI. The value of current CRDI of B3 is close to that of the previous CRDI, and their trends are also very similar. The current CRDI of B4 was also close to the previous CRDI. There are some differences between the two curves in B4, but the trends of them are similar in the missing pixels of the middle and rear segment.
In general, the current CRDI and previous CRDI values of the two lines are different in different regions. However, the trend of the two lines is similar on the whole.
(a) (b) Figure 9. Comparison of CRDI between the 36th and 35th periods of the section line AA' (a) and BB'(b). The grey parts refer to the missing data of the original 36th VCI. According to the trend of the 36th and 35th periods, the CRDI can be divided into several sections separated by purple lines, which were labeled as A1, A2 and A3 and B1, B2, B3 and B4, respectively. The blue line refers to CRDI of the 36th period, and the red line refers to CRDI of the 35th period.
In 2016, point C has 3 of the 46 periods are missing, and the loss rate is 6.5% (Figure 10a). The CRDI values of point C in the current period and the same period vary from 50 to 80, both of which are small. The trend of current CRDI in the first three periods is opposite to that in the same period. However, there is little difference between the estimated value of the 2nd and 3rd period and the value of CRDI in the same period. Furthermore, the trends of the 2nd and 3rd periods of current CRDI are consistent with the 1st and 4th period. The estimated value of the 41st period is basically consistent with the data of the same period and has good continuity with the data before and after. The estimation results of the three periods are all good. The 2nd and 3rd periods of point C are continuous missing pixels, and compared with the 2nd period, the estimated value of the 3rd period has a certain deviation.
In 2016, point D has 3 of the 46 periods missing, with the loss rate of 6.5% (Figure 10b). The range of CRDI of point D in the current period is small, basically between ~70 and ~95, while the range of CRDI in the same period is large, basically between ~25 and ~80. Although the fluctuation ranges of the two curves is quite different, the trend is a little similar. There is a certain gap between the estimated values of the 2nd, 29th and 34th periods and the data of the same period, but the consistency is better with the earlier and later periods. The trend of the current CRDI before and after The grey parts refer to the missing data of the original 36th VCI. According to the trend of the 36th and 35th periods, the CRDI can be divided into several sections separated by purple lines, which were labeled as A1, A2 and A3 and B1, B2, B3 and B4, respectively. The blue line refers to CRDI of the 36th period, and the red line refers to CRDI of the 35th period.
Among all the pixels passed by line BB', 59.1% were missing pixels. The trend of the current CRDI is similar to that of the previous CRDI, the peak valley position is slightly shifted (Figure 9b). However, in some areas, there is no obvious law of the two curves. Line BB' was divided into four parts: B1, B2, B3 and B4. There was no significant difference between the current CRDI and previous CRDI of B1. The trend of the front and middle segments of B1 is similar. The previous CRDI of the tail segment has more details, and the curve has some fluctuations. As a whole, the current CRDI of B2 is higher than that of the previous CRDI. It can be seen that the floating range of the two curves is small. CRDI fluctuated more in the previous period, but it was more stable in the current period. Furthermore, the current CRDI is also less detailed than the previous CRDI. The value of current CRDI of B3 is close to that of the previous CRDI, and their trends are also very similar. The current CRDI of B4 was also close to the previous CRDI. There are some differences between the two curves in B4, but the trends of them are similar in the missing pixels of the middle and rear segment.
In general, the current CRDI and previous CRDI values of the two lines are different in different regions. However, the trend of the two lines is similar on the whole.
In 2016, point C has 3 of the 46 periods are missing, and the loss rate is 6.5% (Figure 10a). The CRDI values of point C in the current period and the same period vary from 50 to 80, both of which are small. The trend of current CRDI in the first three periods is opposite to that in the same period. However, there is little difference between the estimated value of the 2nd and 3rd period and the value of CRDI in the same period. Furthermore, the trends of the 2nd and 3rd periods of current CRDI are consistent with the 1st and 4th period. The estimated value of the 41st period is basically consistent with the data of the same period and has good continuity with the data before and after. The estimation results of the three periods are all good. The 2nd and 3rd periods of point C are continuous missing pixels, and compared with the 2nd period, the estimated value of the 3rd period has a certain deviation.
In 2016, point D has 3 of the 46 periods missing, with the loss rate of 6.5% (Figure 10b). The range of CRDI of point D in the current period is small, basically between~70 and~95, while the range of CRDI in the same period is large, basically between~25 and~80. Although the fluctuation ranges of the two curves is quite different, the trend is a little similar. There is a certain gap between the estimated values of the 2nd, 29th and 34th periods and the data of the same period, but the consistency is better with the earlier and later periods. The trend of the current CRDI before and after the 34th periods is similar to that of the same period, the estimation result of the 34th period is better. It should be noted that CRDI in 2015 showed a significant trough between the 17th and 28th periods. It can be inferred that there may be drought near point D in this period, but there is no obvious drought in 2016.
In 2016, point E has 4 of the 46 periods missing, and the loss rate is 8.7 % (Figure 10c). The value of CRDI in the current period and that in the same period fluctuated greatly, basically between 40 and 90. Compared with the same period, the estimated values of the 1st, 21st, 29th and 45th period have large gaps. However, the trend of the current period is consistent with that of the previous period in the 1st, 21st and 45th periods. These three periods have good estimation effect.
Remote Sens. 2020, 12, x FOR PEER REVIEW 12 of 19 the 34th periods is similar to that of the same period, the estimation result of the 34th period is better. It should be noted that CRDI in 2015 showed a significant trough between the 17th and 28th periods. It can be inferred that there may be drought near point D in this period, but there is no obvious drought in 2016.
In 2016, point E has 4 of the 46 periods missing, and the loss rate is 8.7 % (Figure 10c). The value of CRDI in the current period and that in the same period fluctuated greatly, basically between 40 and 90. Compared with the same period, the estimated values of the 1st, 21st, 29th and 45th period have large gaps. However, the trend of the current period is consistent with that of the previous period in the 1st, 21st and 45th periods. These three periods have good estimation effect. From the perspective of spatial continuity, the estimation effect of CRDI is quite good. There is a good consistency between the estimated data of the three points and the trend of the previous and later data. From the perspective of time continuity, the estimation effect of CRDI is also good.

Comparison of Spatial Distribution Between CRDI and VCI
Take the 3rd (mid-January) and 15th (early April) period of 2016 as an example; in the enlarged area where VCI deficiency is concentrated, the distribution of the current CRDI, the original VCI and the previous CRDI are compared.
The value of VCI in the 3rd period of 2016 is relatively high, basically distributed above 80. The missing pixels of this period are concentrated in Leizhou Peninsula and eastern Guangdong, and the missing pixels account for 3.33% of the effective pixels. However, the loss rate of CRDI is 0.04%, most of the missing pixels can be effectively filled by CRDI model (Figure 11). The estimated values are basically distributed between 70 and 90.
In the 3rd period, the missing pixels in the pink region are mainly forest type. In the brown region, the missing pixels were mainly those of an agricultural type. In the blue region, the missing pixels were mainly those of the forest grass mixed type. The filling effects of these three regions of the estimated value are quite good. The spatial distribution of CRDI after filling are also good. There is a certain difference between the current CRDI and the corresponding previous CRDI in the second region, because the original VCI is different from the previous CRDI itself. There is a small area of CRDI in the middle of the graph with low value in the previous CRDI in the third region, but this information is lost in the current CRDI. The previous data can only be used as reference, not exactly the same as the current data. From the perspective of spatial continuity, the estimation effect of CRDI is quite good. There is a good consistency between the estimated data of the three points and the trend of the previous and later data. From the perspective of time continuity, the estimation effect of CRDI is also good.

Comparison of Spatial Distribution Between CRDI and VCI
Take the 3rd (mid-January) and 15th (early April) period of 2016 as an example; in the enlarged area where VCI deficiency is concentrated, the distribution of the current CRDI, the original VCI and the previous CRDI are compared.
The value of VCI in the 3rd period of 2016 is relatively high, basically distributed above 80. The missing pixels of this period are concentrated in Leizhou Peninsula and eastern Guangdong, and the missing pixels account for 3.33% of the effective pixels. However, the loss rate of CRDI is 0.04%, most of the missing pixels can be effectively filled by CRDI model (Figure 11). The estimated values are basically distributed between 70 and 90.
In the 3rd period, the missing pixels in the pink region are mainly forest type. In the brown region, the missing pixels were mainly those of an agricultural type. In the blue region, the missing pixels were mainly those of the forest grass mixed type. The filling effects of these three regions of the estimated value are quite good. The spatial distribution of CRDI after filling are also good. There is a certain difference between the current CRDI and the corresponding previous CRDI in the second region, because the original VCI is different from the previous CRDI itself. There is a small area of CRDI in the middle of the graph with low value in the previous CRDI in the third region, but this information is lost in the current CRDI. The previous data can only be used as reference, not exactly the same as the current data. Remote Sens. 2020, 12, x FOR PEER REVIEW 13 of 19 Figure 11. Comparison of spatial distribution of original VCI, current CRDI and previous CRDI of the 3rd period of 2016.
The high value of VCI in the 15th period of 2016 was more than 80. Furthermore, there are also some low value distributions below 20, but not so much. The distribution of missing pixels is relatively scattered, accounting for 3.27% of the effective pixels. However, the loss rate of CRDI is 0.01%, most of the missing pixels can be effectively filled by CRDI model. The estimated values are basically distributed between 70 and 90.
In the 15th period, most of the missing pixels are forest or forest grass mixed type. The filling effect of these three regions of the estimated value are quite good (Figure 12). The spatial distributions of CRDI after filling are also good. There are some differences between the estimated data of the second and third regions and the corresponding previous CRDI. This is because some of the missing pixels of the current VCI are also missing in the previous period. In other words, some of the missing pixels in the 15th period were also missing in the 14th period, and these pixels were continuously missing. The high value of VCI in the 15th period of 2016 was more than 80. Furthermore, there are also some low value distributions below 20, but not so much. The distribution of missing pixels is relatively scattered, accounting for 3.27% of the effective pixels. However, the loss rate of CRDI is 0.01%, most of the missing pixels can be effectively filled by CRDI model. The estimated values are basically distributed between 70 and 90.
In the 15th period, most of the missing pixels are forest or forest grass mixed type. The filling effect of these three regions of the estimated value are quite good ( Figure 12). The spatial distributions of CRDI after filling are also good. There are some differences between the estimated data of the second and third regions and the corresponding previous CRDI. This is because some of the missing pixels of the current VCI are also missing in the previous period. In other words, some of the missing pixels in the 15th period were also missing in the 14th period, and these pixels were continuously missing. Remote Sens. 2020, 12, x FOR PEER REVIEW 14 of 19 Figure 12. Comparison of spatial distribution of original VCI, current CRDI and previous CRDI of the 15th period of 2016.

The Estimation of Extremum Value
Compared with the histogram of current CRDI, original VCI and the same period CRDI ( Figure  7), for the CRDI of 2016, the number of estimated pixels distributed between 40 and 80 for total and different land cover types was more than the original data and the previous data. The amount of pixel distribution between 80 and 90, even between 90 and 100, is relatively small. According to the

The Estimation of Extremum Value
Compared with the histogram of current CRDI, original VCI and the same period CRDI (Figure 7), for the CRDI of 2016, the number of estimated pixels distributed between 40 and 80 for total and different land cover types was more than the original data and the previous data. The amount of pixel distribution between 80 and 90, even between 90 and 100, is relatively small. According to the distribution trend of current VCI and CRDI in previous, more estimated pixels should be distributed between~80 and~100, but the number of estimated pixels in this region is less, especially in the range of 90~100. This shows that the actual estimates of many pixels that should be between 80 and 100 are smaller. In the range of 80~100, especially 90~100, the CRDI estimation results are not as accurate as other intervals.
CRDI is not so accurate for the estimation of extreme value. Usually, the same set of parameters is used to estimate the same land cover type of the same period. The data used to calculate the parameters of regression equation are distributed between 0 and 100, and the parameters are inclusive, so the parameters usually have some deviation for the estimation of extreme value.

The Influence of Continuous Loss of CRDI
For the accuracy of CRDI, the first problem of continuous missing pixels is that if the pixel cannot be estimated in the previous period, the pixel cannot be estimated in the current period ( Figure 11). The second problem is that the residual of CRDI will accumulate. Due to the lack of correction of the real VCI value in the continuous missing pixels, the estimation residual will become larger and larger (Figures 11 and 12).
The number of continuous missing pixels accounted for 0.47% of the total effective pixels of the 12th, 13th and 14th periods of 2016. The estimated values of continuous missing pixels and the histogram corresponding to the data before and after showed certain regularity (Figure 13a). Before and after the continuous loss of pixels, the VCI values were basically distributed between 40 and 70, and the distribution was relatively uniform. The estimated values are also distributed in this interval, but more concentrated between 50 and 60. of 90~100. This shows that the actual estimates of many pixels that should be between 80 and 100 are smaller. In the range of 80~100, especially 90~100, the CRDI estimation results are not as accurate as other intervals.
CRDI is not so accurate for the estimation of extreme value. Usually, the same set of parameters is used to estimate the same land cover type of the same period. The data used to calculate the parameters of regression equation are distributed between 0 and 100, and the parameters are inclusive, so the parameters usually have some deviation for the estimation of extreme value.

The Influence of Continuous Loss of CRDI
For the accuracy of CRDI, the first problem of continuous missing pixels is that if the pixel cannot be estimated in the previous period, the pixel cannot be estimated in the current period ( Figure 11). The second problem is that the residual of CRDI will accumulate. Due to the lack of correction of the real VCI value in the continuous missing pixels, the estimation residual will become larger and larger (Figures 11 and 12).
The number of continuous missing pixels accounted for 0.47% of the total effective pixels of the 12th, 13th and 14th periods of 2016. The estimated values of continuous missing pixels and the histogram corresponding to the data before and after showed certain regularity (Figure 13a). Before and after the continuous loss of pixels, the VCI values were basically distributed between 40 and 70, and the distribution was relatively uniform. The estimated values are also distributed in this interval, but more concentrated between 50 and 60.
Among them, the distribution of the first period estimates is relatively average, which is more in line with the transition trend from the previous period to the next period ( Figure 13b). The histogram of the estimated values of continuous missing pixels between 46 and 74 has certain regularity. Estimate values for the second and third periods are more concentrated than for the previous period. Compared with the first period, the second period is more concentrated between 56 and 58. Furthermore, compared with the second period, the third period is more concentrated between 56 and 60. When there are continuous missing pixels, the cumulative residual will be larger and larger, and the estimated value will be more and more concentrated on a certain value, and the difference between the missing data pixels will be smaller and smaller.
(a) (b) Figure 13. Distribution of estimated values of continuous missing pixels: (a) compare the continuous missing data with the data before and after (b) comparison between consecutive missing data.

A Large Number of Data Missing Periods
If the number of missing pixels is too large, the accuracy of CRDI will be affected. The number of reference pixels used to calculate the model parameters is relatively small and concentrated. Therefore, the deviation of the calculated parameters is relatively large (more than 50%), and the Among them, the distribution of the first period estimates is relatively average, which is more in line with the transition trend from the previous period to the next period (Figure 13b). The histogram of the estimated values of continuous missing pixels between 46 and 74 has certain regularity. Estimate values for the second and third periods are more concentrated than for the previous period. Compared with the first period, the second period is more concentrated between 56 and 58. Furthermore, compared with the second period, the third period is more concentrated between 56 and 60. When there are continuous missing pixels, the cumulative residual will be larger and larger, and the estimated value will be more and more concentrated on a certain value, and the difference between the missing data pixels will be smaller and smaller.

A Large Number of Data Missing Periods
If the number of missing pixels is too large, the accuracy of CRDI will be affected. The number of reference pixels used to calculate the model parameters is relatively small and concentrated. Therefore, the deviation of the calculated parameters is relatively large (more than 50%), and the estimated value also has a certain deviation. There are 11 periods with the loss rate more than 50%, accounting for about 1.6% of the total 690 periods.
68.15% of the valid pixels were missing in the 10th period of 2016, and the missing pixels were very concentrated. There is a large deviation in the estimated data of this period. However, the estimated data still have some continuity (Figure 14). In the 10th period of 2016, almost all the pixels in northern Guangdong and eastern Guangdong are missing. Due to the small number of original VCI pixels used to calculate parameters in this period, ADI plays a smaller role while COT plays a more important role in the estimation of missing pixels in this period. As a result, there are obvious edges the CRDI period below, which are the edges of the COT period. Even so, the estimated values still follow a certain distribution law, and even the left and right ends of the COT pixel edge still have a certain continuity.
Remote Sens. 2020, 12, x FOR PEER REVIEW 16 of 19 estimated value also has a certain deviation. There are 11 periods with the loss rate more than 50%, accounting for about 1.6% of the total 690 periods. 68.15% of the valid pixels were missing in the 10th period of 2016, and the missing pixels were very concentrated. There is a large deviation in the estimated data of this period. However, the estimated data still have some continuity (Figure 14). In the 10th period of 2016, almost all the pixels in northern Guangdong and eastern Guangdong are missing. Due to the small number of original VCI pixels used to calculate parameters in this period, ADI plays a smaller role while COT plays a more important role in the estimation of missing pixels in this period. As a result, there are obvious edges the CRDI period below, which are the edges of the COT period. Even so, the estimated values still follow a certain distribution law, and even the left and right ends of the COT pixel edge still have a certain continuity.

Evaluation of CRDI and Future Research Directions
The main advantage of CRDI is that it can solve the problem of large blank of drought data caused by cloud cover in a traditional remote sensing model. CRDI has a good application prospect and strong applicability. In this paper, VCI, which is more advantageous in vegetation area, is selected as the remote sensing drought index. When applying the CRDI model, other drought indexes, such as CWSI, SPI or PDSI, can be selected according to the specific situation of the study area.
During the discussion of the results, we found that the CRDI model is not so accurate for the estimation of the extreme value, the continuous loss pixels and the large number of data for the missing periods. In order to improve the accuracy of the CRDI model, the following aspects can be considered for future research: 1. The accuracy of basic data needs to be improved. The spatial resolution of VCI and COT is 250 m and 1 degree, respectively, COT should be resampled to 250 m. In the follow-up study, higher resolution data can be considered. The accuracy of the estimation results depends on the precision of land use/land cover data. In order to improve the accuracy, more precise land use/land cover data can be considered. 2. The effect of the CRDI model in small area may be better. In calculating the parameters of the regression equation, all the pixels without missing data are used, and then the values of all missing pixels are estimated with the obtained parameters. Therefore, the intra class differences of different land cover type in the study area will affect the results of the model. The smaller the study area is, the smaller the intra class differences of different land cover type may be. Therefore, the parameters obtained are more representative of the study area. Therefore, the algorithm of the model can be improved. For example, an n × n window can be opened with the target pixel as the center, and the parameters of the model can be calculated according to the data in the window, and then the drought value of the target pixel can be estimated. However, this

Evaluation of CRDI and Future Research Directions
The main advantage of CRDI is that it can solve the problem of large blank of drought data caused by cloud cover in a traditional remote sensing model. CRDI has a good application prospect and strong applicability. In this paper, VCI, which is more advantageous in vegetation area, is selected as the remote sensing drought index. When applying the CRDI model, other drought indexes, such as CWSI, SPI or PDSI, can be selected according to the specific situation of the study area.
During the discussion of the results, we found that the CRDI model is not so accurate for the estimation of the extreme value, the continuous loss pixels and the large number of data for the missing periods. In order to improve the accuracy of the CRDI model, the following aspects can be considered for future research: 1.
The accuracy of basic data needs to be improved. The spatial resolution of VCI and COT is 250 m and 1 degree, respectively, COT should be resampled to 250 m. In the follow-up study, higher resolution data can be considered. The accuracy of the estimation results depends on the precision of land use/land cover data. In order to improve the accuracy, more precise land use/land cover data can be considered.

2.
The effect of the CRDI model in small area may be better. In calculating the parameters of the regression equation, all the pixels without missing data are used, and then the values of all missing pixels are estimated with the obtained parameters. Therefore, the intra class differences of different land cover type in the study area will affect the results of the model. The smaller the study area is, the smaller the intra class differences of different land cover type may be. Therefore, the parameters obtained are more representative of the study area. Therefore, the algorithm of the model can be improved. For example, an n × n window can be opened with the target pixel as the center, and the parameters of the model can be calculated according to the data in the window, and then the drought value of the target pixel can be estimated. However, this improvement has higher requirements for the running efficiency of the algorithm, and the running efficiency of the algorithm should also be considered. Furthermore, how to determine the size of window is also a difficult problem.

Conclusions
In this paper, a drought retrieval model CRDI based on COT and ADI is established for the pixels whose drought value cannot be obtained due to cloud cover during the monitoring period.
The drought values of the missing pixels of different land cover types were estimated for each period. The effect of the CRDI model was evaluated from the filling efficiency of the model, the rationality of the estimated value and the temporal and spatial continuity of CRDI. Furthermore, the estimations of large missing and continuous missing periods are analyzed. The main conclusions are as follows: 1.
The filling efficiency of CRDI is high. CRDI can fill most of the missing data. The average filling efficiency of total data, forest, forest grass mixed and agricultural was as high as 98.0%, 99.1%, 97.5% and 99.5%. Even for the most inefficient periods, 80% were filled.

2.
The estimated value of CRDI is reasonable. The range distribution of CRDI and its corresponding original VCI is similar, and the estimated value is also within a reasonable range. However, the estimated value is more concentrated than the corresponding previous and the same period data, and the estimation of extreme value is not as accurate as others.

3.
The continuity of CRDI is quite good. The continuity of CRDI data is analyzed from the perspective of space and time by comparison section line and sample points. The trend of the section line between the current CRDI and the previous CRDI is very similar, and it is reasonable for CRDI to fill the missing areas. The trend of the CRDI curve of sample points is similar to that of the same period, and the continuity of estimated data is better than that of previous and next curve. Considering both time and space, the CRDI model has good continuity. In addition, by comparing the spatial distribution of CRDI and VCI, it is found that the estimated values have good continuity with the surrounding drought values.
Author Contributions: J.Y. proposed the research idea. W.L. and Y.W. determine the specific research methods. W.L. carried out the RS data processing and analysis and wrote the manuscript. All authors have read and agreed to the published version of the manuscript.