An Improved Cloud Gap-Filling Method for Longwave Infrared Land Surface Temperatures through Introducing Passive Microwave Techniques

Abstract

Satellite-derived land surface temperature (LST) data are most commonly observed in the longwave infrared (LWIR) spectral region. However, such data suffer frequent gaps in coverage caused by cloud cover. Filling these ‘cloud gaps’ usually relies on statistical re-constructions using proximal clear sky LST pixels, whilst this is often a poor surrogate for shadowed LSTs insulated under cloud. Another solution is to rely on passive microwave (PM) LST data that are largely unimpeded by cloud cover impacts, the quality of which, however, is limited by the very coarse spatial resolution typical of PM signals. Here, we combine aspects of these two approaches to fill cloud gaps in the LWIR-derived LST record, using Kenya (East Africa) as our study area. The proposed “cloud gap-filling” approach increases the coverage of daily Aqua MODIS LST data over Kenya from <50% to >90%. Evaluations were made against the in situ and SEVIRI-derived LST data respectively, revealing root mean square errors (RMSEs) of 2.6 K and 3.6 K for the proposed method by mid-day, compared with RMSEs of 4.3 K and 6.7 K for the conventional proximal-pixel-based statistical re-construction method. We also find that such accuracy improvements become increasingly apparent when the total cloud cover residence time increases in the morning-to-noon time frame. At mid-night, cloud gap-filling performance is also better for the proposed method, though the RMSE improvement is far smaller (<0.3 K) than in the mid-day period. The results indicate that our proposed two-step cloud gap-filling method can improve upon performances achieved by conventional methods for cloud gap-filling and has the potential to be scaled up to provide data at continental or global scales as it does not rely on locality-specific knowledge or datasets.

1. Introduction

Land surface temperature (LST) is classified as a Global Climate Observing System Essential Climate Variable (GCOS-ECV) [1]. It is widely used in both research and commercial applications, with its key domains of relevance including agriculture [2], urban landscape management [3], disaster risk analysis [4], and investigations on heat flux and hydrological features across the globe [5,6]. However, satellite LST data derived from brightness temperatures (BT) recorded in the longwave infrared (LWIR) spectral region have the key limitation that widespread cloud can completely obscure land surface from the sensors’ view [7,8,9]. The resulting spatiotemporal ‘cloud gaps’ pose a challenge for applications reliant on regular and routine LST observations, both in terms of environmental models and derived data products [10,11]. The cloud gap problem may be further exacerbated when the sensor providing BT data is mounted on a polar orbiting satellite. In this scenario only a few (usually 1–2) views of an area in the tropical low latitudes and the temperate middle latitudes can be provided in every diurnal cycle, and these views can often be compromised by high (>50%) levels of cloud cover. This is a critical problem for both the monitoring of LST and its downstream users in applications, such as drought monitoring [12], urban island heat characterization [13], and fire detection [14,15].

1.1. Prior Cloud Gap-Filling Research

The filling of ‘cloud gaps’ in LWIR-derived LST data has been the subject of much research, with methods focusing primarily on data interpolation approaches based on temporal techniques such as the Savitzky–Golay filter [16,17], single spectrum analysis [18], and other similar methodologies [19]. The key assumption in such studies is a high continuity between daily LST dynamics at the scale of a remotely sensed pixel. However, daily continuity breaks down when a region suffers significantly from cloud shadowing, storm cells, differential urban heating, short-term precipitation events, and other environmental phenomena [20,21]. Therefore, there is an urgent need to enable a cloud gap-filling method that can account for these sub-daily changes in LST to minimize their negative effect on the intra-day continuity that temporal interpolation relies on. Geo-statistical interpolation is a variant of the above interpolation approaches for LST cloud gaps, in which areas of no-data are filled using functions built from the spatial distribution of the available clear-sky observations obtained at, or close to, the same local time at which the cloud gaps are observed [22,23,24]. However, this approach has been proven unreliable when the study area is affected by larger portions of contiguous data gaps [25] and is thus not suitable for use in areas suffering frequent cloud.

Another type of approach exploited to fill cloud gaps is based on LST variation and distribution patterns in the spatio-temporal space, sometimes referred to as “spatio-temporal data fusion” (STDF) methods. An STDF framework is usually constructed using clear-sky LST observations of spatially neighbouring pixels made at proximal or near-proximal dates. The methods deployed include linear models [26,27], Markov-based models [28], and those based on ratios between daily LST observations and the corresponding eight-day averages [29]. Improved performance of these methods, and thus more stable cloud gap LST estimates, have been obtained through the introduction of additional data beyond LST alone [25,30,31]. However, a flaw of this approach is highlighted by Song et al. [30] in that gap-filled pixels of MODIS LST data based on STDF methods represent an estimate of the LST that would be found under clear sky conditions, even if the gaps are caused by clouds. This is because the STDF approaches rely on relationships built from ‘clear sky’ LST observations. Therefore, a ‘clear sky bias’ within the STDF cloud gap-filled product can be introduced because the amount of outgoing longwave and incoming shortwave radiation in clouded areas is affected by the cloud cover in a manner not accounted for by the clear sky STDF relationship.

The use of physically based land surface energy balance (LSEB) modelling [31,32,33] and numerical weather prediction (NWP) modelling [34] are theoretically not affected by a clear-sky bias. However, since they require accurate local soil properties to work effectively, the LSEB models have large uncertainty on the model parameterization process. A recent investigation [35] has attempted to overcome this flaw by fusing LSEB theory with machine-learning-based methods. Despite the good accuracies (RMSE < 2.5 K) achieved in terms of results, their inputs still require spatially and temporally contiguous shortwave radiation products that are difficult to obtain globally with a high accuracy and at a high resolution from geo-stationary satellites [36]. In addition, validation of energy balance models of surface temperature has shown that there is a tendency for such models to suffer a cold bias [37,38]. The NWP modelling technique is much more sophisticated in describing the land surface and atmospheric conditions but has been shown to struggle in accuracy with regards to surface temperatures over a variety of landcovers [39,40]. Overall, physical models as they currently exist should be further improved to act as a suitable answer on their own to the filling of LST cloud gaps, especially on large (global) scales.

The final methodological approach available to tackle the cloud gap problem is to exploit passive microwave (PM) data. The advantage of PM-based approaches is the ability of microwave radiation to penetrate clouds, which allows the spaceborne sensor to record surface-emitted PM radiation and support the production of ‘all-weather’ LST data [41]. However, the ground pixel footprints of satellite-based PM instruments suitable for LST estimation are typically coarser than 10 × 10 km². This is a far worse spatial resolution than is provided by LWIR sensors, such as the Moderate Spatial Resolution Imaging Spectroradiometer (MODIS) and the Visible infrared Imaging Radiometer (VIIRS), and insufficient for many applications of LST data given the strong spatial heterogeneity typical of land surface temperatures. Furthermore, PM-derived LST estimates represent the temperature at some depth in the soil, typically 2 mm to 1 cm, which is different to the surface skin temperature estimates typically provided by LWIR-based LST methods [42,43,44].

Due to the known limitations of PM-derived LST data, numerous investigations have been performed with respect to how PM- and LWIR-derived LST data records might be combined to produce an ‘all-weather’ daily LST record at a spatial resolution the same or similar to that of the original LWIR-derived LST record.

Table 1. Summary of LWIR and Passive Microwave data fusion studies related to the production of ‘all-weather’ daily LST information at a higher spatial resolution than PM data alone can produce.

Study	Main Method Description	Spatial Scale	Limitations
1. Kou et al. [45]	A Bayesian Maximum Entropy (BME) blending approach to merge PM and LWIR data, achieving a Root Mean Square Error (RMSE) in LST of between 2.3–4.5 K.	A relatively small 100 × 100 km2 region	Only used on night-time data over a small region, with its universality requiring more validation.
2. Duan et al. [43]	An empirical model based on a digital elevation model (DEM) and clear sky LWIR-derived LST at neighboring pixels to downscale PM-derived data for achieving LST of cloudy LWIR pixels.	China	Downscaling of PM LST only relies on DEM, may be theoretically less effective in areas where the topographical variation is not important (e.g., low altitude plains).
3. Sun et al. [46]	A downscaling method for PM LST using NDVI and DEM data applied for gap-filling of LWIR LST at a finer resolution.	China	(1) Penetration depth difference between PM and LWIR not considered. (2) The method was applied over a large area but only validated at limited sites. (3) PM data were downscaled to 5 km but the feasibility of the method at a finer resolution (e.g., 1 km) requires further investigation.
4. Zhang et al. [47]; Zhang et al. [48]; Zhang et al. [49];	A temporal component decomposition developed for merging observations that achieved a 1-km all-weather daily LST data.	Northeastern China; the Tibetan Plateau	Both microwave data and reanalysis data are required as essential inputs to reconstruct the real LST under cloud. But this increases the complexity of the method and risks increased uncertainty from data inputs. Therefore, it should be cautiously discussed when applied in other regions.
5. Long et al. [50]	A data fusion method used to merge LWIR observations and PM-like coarse-resolution reanalysis datasets, based on correlations between images taken a limited time apart.	3 plots of 80 × 80 km2 in China	(1) This time-interpolation-like fusion method has strict requirement on the availability of its input datasets at neighboring dates [49]. (2) The method only addresses situations of temporally discontinuous pixel loss across relatively small study regions.
6. Yoo et al. [51]; Shwetha and Kumar [52]	Machine-learning based models to fuse LST at different spatial scales.	South Korea, about 10,000 km2; Cauvery river basin in India, about 80,000 km2	The physics behind machine learning models remains unclear, and it is difficult to justify the global universality of these models when the relationships they derive cannot be explicitly formalized.

summarises key studies in this field from the last five years, almost all of which rely on the environmental characteristics of the particular region studied. Most of the approaches did not address the potential universality of their method beyond the particular study region, which is likely a limitation to producing the all-weather daily LST at continental/global scales.

1.2. Research Objective

Based on the research described above, both conventional STDF approaches and PM-LWIR data fusion approaches have a mixture of merits and limitations. Therefore, a combination of these has potential for compensating for the ‘clear sky bias’ present in the STDF approaches and for removing the regional restriction present in many existing PM-LWIR fusion methods. In this study, we build on past efforts by developing a novel two-step framework for better estimation of LST in cloud covered areas at a relatively high spatial resolution, based on a combination of LWIR-derived LSTs and PM LST data. We applied the approach at a national scale in tropical Africa, and through our efforts aim to further enhance the universality of an approach to generate all-weather, fine spatial resolution (~1 km) LST data. We also aim to investigate the effect of cloud insulation on the gap filling of LST, answering whether and how PM observations might mitigate biases in the cloud gap filled LST data generated by the conventional STDF methodology. Our work may benefit future efforts to produce globally quasi-full coverage daily LST data at resolutions similar to those provided by current LWIR sensors.

2. Study Area and Datasets

2.1. Study Area and In Situ LST Ground Observations

The study area for this work covers Kenya (East Africa) with a total area of around 582,646 km². As illustrated in Figure 1, the far southwestern part of the country is characterized by highly vegetated cropland and forest, with complex topography at higher elevations (>1000 m). In contrast, the rest of the region is primarily dominated by flat plains and covered with sparser vegetation types (grass or shrub savannas). Most of Kenya is dominated by Arid (BSh, BWh) and Tropical (Af, Am, Aw) Köppen–Geiger climate zones with some Temperate (Cfb, Csb, Cwb) climate zones found in the south-west [53]. Snow and ice rarely occur except at high altitude. The local climate experiences a bi-modal rainfall pattern, with a season of ‘short rains’ typically between October and December and ‘long rains’ from April to late May or June. Other months are identified as dry seasons. Given the wide range of precipitation, topography and climate types, methodologies for improving LST cloud gap-filling developed in Kenya are likely applicable to many other areas globally, especially in the (sub) Tropics and the Temperate zones.

The ILRI Kapiti Research Station is located in south-west Kenya (−1.6083°, 37.1327°), as shown in Figure 1a. This station is run by the International Livestock Research Institute (ILRI) and provides observations of LST and land surface energy balance, as well as CO₂ and H₂O fluxes. The observation sites are instrumented using a series of mast mounted (Heitronics) LWIR radiometers to deliver upward (53° view zenith angle) and downward (15° view zenith angle) looking LWIR brightness temperatures. They are then processed to deliver logged LST estimates every 5-min according to the LST in situ measurement principles described in Gottsche et al. [54]. There are four sites in all (Masts 1–4 in Figure 1a) which can provide LWIR radiometer measurements suitable for the evaluation of satellite-based LST products since September 2018. The research station has a fenced 15 km × 20 km area of homogenous savannah, which can further be broken down into 94% grassland-soil complexes and 6% tree canopy cover.

We employ the in situ LST data on Masts 1–4 of ILRI Kapiti Research Station as our primary benchmark (validation) dataset, focusing on the period from October 2018 to February 2019. A brief description of the derivation and accuracy evaluation of our benchmark LST data is provided in the Supplementary Material (Figure S1). The certainty (range of possible error) of the ground observations ranges between 0.79 and 0.95 K per observation throughout the whole time series, across the four masts. The range of possible error consists of radiometric measurement uncertainty, surface emissivity uncertainty, and the uncertainty contribution of the logging equipment as driven by environmental temperature. This time period, encompassing both the short-rains and subsequent dry period, is sufficient in capturing the major climatic and vegetative cycles of an entire year (Figure 1b), by which land surface emissivity and the thermodynamic temperature are strongly affected [55]. In particular, the selected period includes the very rapid green-up (typically only a few days) when the rains start, the rainy sub-period accompanied by increasing vegetation canopy cover, and the sub-period of landscape drying when the rain ends. This phenological pattern is representative of the study area as a whole, as can be seen from the NDVI time series of five typical LWIR pixels in different locations and different land covers across Kenya in Figure 1b. Notice that the locations of the five pixels, i.e., Pixels A (at 3.63° N, 40.35° E), B (4.16° N, 35.23° E), C (0.78° S, 35.22° E), D (0.76° S, 40.84° E), and E (0.93° N, 39.37° E), have been marked in Figure 1c. They are respectively covered by different land use types of shrubland, grassland, cropland, forest, and urban area.

2.2. Satellite Datasets

Several satellite datasets were utilized for our study. The enhanced Aqua MODIS 1-km daily LST product (MYD21A1D/N.v061) was selected as the cloud gap-filling target. In addition, we also employed the MODIS MCD43A4 (v061) 500-m daily “Bidirectional Reflectance Distribution Function (BRDF)” adjusted reflectance dataset and the 90-m DEM data generated by the NASA Shuttle Radar Topography Mission (SRTM). The PM data employed are the Level-1R BT data collected by the Advanced Microwave Scanning Radiometer-2 (AMSR-2) operated onboard the Global Change Observation Mission1-Water (GCOM-W1) satellite. The polar orbiting GCOM-W1 was launched by the Japan Aerospace Exploration Agency (JAXA) in May 2012 and shares the same orbital track as Aqua, with a very closing equator-crossing time (about 1:30 p.m./a.m.). The similarity of AMSR-2 and Aqua MODIS observation times is beneficial for blending their observations to fill the gaps of MODIS LWIR LST under cloud.

In addition to the ~5-month in situ data record coming from the ILRI Kapiti Research Station, we also employed the LST data from the geostationary Meteosat platform’s Spinning Enhanced Visible and Infra-Red Imager (SEVIRI) [56] as a second benchmark validation dataset. SEVIRI-derived LST data covering an entire climatic year (from October 2018 to September 2019) were used across the whole Kenya. The SEVIRI LST product is provided every 15 min by EUMETSAT’s Land Surface Analysis Satellite Application Facility (LSA-SAF). Trigo et al. [56] describe the adapted split window algorithm for producing LST data from SEVIRI-measured LWIR brightness temperatures (BTs). Gottsche et al. [54] report that the product has an overall RMSE of 3.2 K compared to LWIR-radiometer derived LSTs collected at a West African validation station located in a tiger bush savanna biome. With a pixel size of around 4 km in Kenya, the primary role of Meteosat LST data was to aid evaluation of both the standard MODIS LST data and the cloud gap-filled MODIS LST data. The latter is possible because instances exist where MODIS pixels classed as cloudy (and thus with no MODIS LST data) are matched by Meteosat LST data that do have an LST observation classed as cloud free within a time lag of up to around 8 min, due to either cloud movements, differences in the cloud masking procedures, and/or due to the different viewing geometries of the two sensors (SEVIRI = ~40–42° view zenith angle (VZA) over Kenya, MODIS VZA = 0–65° depending on overpass). This can be observed in the contrast show in Figure 4c,e in Section 4.

3. Methodology

The flowchart of our cloud gap-filling approach for MODIS LST data is schematically shown in Figure 2. It generally comprises two steps, namely (1) a conventional realization of the STDF methodology and (2) a PM-driven bias adjustment process to calibrate the data coming out of step (1), accounting for the fact that these data should represent cloud-covered rather than clear-sky LSTs. Therefore, there are three datasets tested in the validation step: (i) the clear sky original MODIS product (MODIS_Clear), (ii) MODIS with gaps filled by step (1) (MODIS_STDF), and (iii) MODIS_STDF bias corrected for under-cloud pixels with passive microwave (AMSR-2) observations (MODIS_PMBC).

3.1. Step 1: STDF Methodology

3.1.1. Methodology Description

The 1 km MODIS LST data were subset to Kenya and pixels classed as cloud contaminated were removed using the “quality control (QC)” field contained within the MODIS LST product. In addition, all pixels flagged as ‘cloud-edge’ or low quality were also removed. A 500-m daily NDVI dataset for Kenya was calculated using the red-band and infrared-band reflectance data from the MODIS MCD43A4 product, which has been corrected for angular reflectance effects based on the BRDF. The resulting estimates were resampled using a spatial averaging procedure to the same 1-km resolution as the MODIS LST product. The spatial resampling process was also applied to the SRTM DEM dataset subset over Kenya. For each date (represented as t₁) in the study period (from October 2018 to September 2019), we used the STDF methodology version detailed in Song et al. [30] to fill cloudy pixels that have no LST data record within the MODIS daily LST product. The Song et al. [30] method was chosen as it has been applied to a large area in South China and obtained relatively good validation performance (with RMSEs between 1.5 and 3.5 K). A summary of the STDF method is given here, whilst the full details can be found in the work of Song et al. [30].

The STDF method first builds a transfer function based on the above-described datasets (LST, NDVI and DEM) as follows:

$$LS{T}^{\ast }{}_{t_1}=a\times LS{T}^{\ast }{}_{t_0}+b\times NDV{I}^{\ast }{}_{t_1}+c\times DE{M}^{\ast }+d$$

(1)

where the superscript “*” indicates that this variable has been normalized to the range 0 to 1.0, based on the maximum and minimum values of that variable found across Kenya during our study period. Parameters a, b, c, and d are coefficients fitted using all pixels that have a clear-sky LST estimate on a specific date t₁ (LST^*_t1) and at least one closely preceding date, t₀ (LST^*_t0) (see below for further detail on how this was calculated). NDVI^*_t1 indicates the corresponding (normalized) NDVI on the t₁ date. The NDVI data are available for almost every pixel continuously because the MCD43A4 daily product is generated using 16-day composited reflectance data, and therefore has good resistance to data loss from cloud contamination. Notice that day-time and night-time LST data must be separately treated in this transfer function.

After deriving the coefficients of a, b, c, and d, the Equation (1) transfer function was used to fill all cloudy MODIS LST pixels on the t₁ date. For any t₁ date included in the study period, the t₀ date was iterated among all neighbouring dates of t₁ meeting the condition |t₀ − t₁| <= 15 days (from the nearest date to the furthest date). An average of the estimated LST values for t₀ was then taken where a cloud gap pixel was filled more than once (based on the iterative t₀ dates), and the iteration was stopped when the fraction of pixels with effective LST values on t₁ was equal to or exceeded 0.9. In the original STDF method of Song et al. [30], the remaining small fraction (<0.1) of cloud gaps was then filled with LST data generated by a conventional inverse distance weighted (IDW) interpolation method [57]. However, the IDW step is not done here in order to avoid bringing in extra interpolation uncertainty that may negatively influence the validation step. The final modification to the original Song et al. [30] method is the use of iteration over all possible dates within a 30-day window (rather than only selecting one t₀ date that has a high fraction of clear-sky pixels) to increase the amount of data available to the procedure and therefore improve the applicability to the entire area of Kenya.

3.1.2. Method Sensitivity Analysis

In this section, we conducted sensitivity analysis for the STDF method expressed through Equation (1), in order to explore the existing uncertainty sources of the method which may influence performance of the subsequent bias adjustment process. Considering the linearity of Equation (1) and the fact that there are no interactive components between different input variables, the sensitivity of the output variable LST_t1 to each input variable X (denoting LST_t0, NDVI_t1, or DEM) can be expressed as the 1st order partial derivative to X, i.e., ($\frac{dLST\left(t_1\right)}{dX}$). It should be noted that these input variables are not normalized between 0–1. The sensitivity was calculated for all iterations conducted during Section 3.1.1. The statistical summary of this sensitivity analysis of the different iterations are reported in

Table 2. Statistical results for sensitivity analysis of the STDF method during the selected study period. The upper and lower threshold of the indicator are respectively defined as the 90th and 10th percentile of all result samples rather than the maximum and minimum values, in order to avoid outliers.

Sensitivity of LST(t1) on the Input Variables X, i.e., $(\frac{\mathit{d}\mathit{L}\mathit{S}\mathit{T}\left({\mathit{t}}_1\right)}{\mathit{d}\mathit{X}})$	Observation Time	Average	Standard Deviation	Upper Threshold (90th Percentile)	Lower Threshold (10th Percentile)
NDVI (K/0.1)	Day-time	−1.11	0.89	2.88	−4.28
	Night-time	−0.06	0.17	0.76	−0.93
DEM (K/100 m)	Daytime	−0.42	0.35	1.21	−6.30
	Night-time	−0.36	0.42	1.43	−5.69
LST(t0) (K/K)	Daytime	0.50	0.84	1.72	−0.80
	Night-time	0.59	0.25	1.49	−0.71

. From Table 2, it is apparent that the difference between daytime and night-time for sensitivity of output LST_t1 against LST_t0 and DEM is low. An LST_t0 difference of 1 K can cause a difference on average of 0.5–0.6 K for the (output) gap filled LST. In comparison an average decrease of about 0.3 K or 0.4 K was observed when elevation rises by 100 m. For NDVI, sensitivity is more important in the daytime compared with night-time, since an NDVI increase of 0.1 can result a decrease >1 K in the gap-filled outcome, while it has an almost negligible effect on the night-time outcomes.

3.2. Step 2: Cloud Shadowing Bias Adjustment with PM-Derived LST Data

As the STDF methodology is based on cloud-free LST observation data from MODIS but is attempting to estimate LST under cloud, biases may be introduced. Identifying and removing such bias is the focus of Step 2 and is based on the AMSR-2 PM-derived LST data. AMSR-2 provides both horizontally (H-) and vertically (V-) polarized BT observations in seven different frequencies [58]. LST data have been estimated by different combinations of these data [59,60,61,62,63], including at 6.9 GHz, 10.7 GHz, 18.7 GHz, 23.8 GHz, 36.5 GHz, and 89 GHz. The linear LST retrieval algorithm based on the 36.5 GHz V-polarized BT [61,64] is among the most widely used, though Song et al. [65] found it inferior to that based on the 18.7 GHz and 23.8 GHz bands [63,66], especially for pixels with higher fractions of water surface (including standing water and dynamic water signals like high levels of soil moisture). To minimize the influence from increased soil moisture on LST retrieval in the rainy season of our study area, we employ the AMSR-2 V- and H-polarized 18.7 GHz and 23.8 GHz BTs to retrieve LST using the methodology of Jones et al. [63]. The resulting PM-derived LST estimates have a spatial resolution of 25 km after the AMSR-2 Level1R BT data were resampled according to the global Equal Area Scalable Earth (EASE) Grid [67] projection system.

The PM BTs at 18.7 GHz and 23.8 GHz have a theoretical temperature sensing depth of around 1 cm [68], as compared to the MODIS LWIR-derived LST skin surface depth of ≤1 mm [69]. To deal with this we first averaged the MODIS 1-km LST data of Kenya to the 25-km EASE Grid and compared those cells that have more than 600 samples (i.e., >95%) of clear-sky MODIS pixels to the AMSR-2 retrieved LST for the same region. As seen in Figure 3, both daytime (ascending) and night-time (descending) MODIS and AMSR-2 LSTs show strong correlations but depart from the 1:1 line. This result agrees with the observations of Song et al. [65], who demonstrated that, apart from in desert and snow/ice covered areas, the LST derived from these two different sensors have a strong linear correlation at the global scale. The AMSR-2 observations have a negative bias of around −0.7 K in the daytime and a positive bias of around 2.5 K in the night-time. This is consistent with the logic that compared to temperature of sub-surface soils, that of skin surface is usually higher at mid-noon but lower at midnight. Therefore, the bias in Figure 3 can be largely ascribed to the different sensing depths. To enable the AMSR-2 LST data to provide the best estimate of what MODIS would have measured, we fit an independent linear equation to the data of day-time and night-time modes respectively, as shown by the orange text in Figure 3. These achieved high coefficients of determination (R² ≥ 0.9), with RMSEs between the AMSR-2 LST data adjusted via these equations and the true MODIS LST (RMSE_unbias) of less than 1.5 K. The small value of RMSE_unbias proves a high degree of consistency between PM (AMSR-2)-derived LST and MODIS LWIR LST, thereby confirming the feasibility of the penetration depth bias adjustment process.

Next, we inverted the linear relations in Figure 3 and thus transformed the AMSR-2 based LST data (LST_AMSR-2) into the corresponding best estimate of MODIS-derived LST at the 25-km pixel scale (LST_{MODIS_25km}):

$$LS{T}_{ MODIS\_25km}=k_0\times LS{T}_{AMSR-2}+m_0$$

(2)

where $k_0$ and $m_0$ are the linear relationship parameters found in Figure 3. This removes the penetration depth bias. Data at 1-km resolution can then be further calibrated to remove the cloud shadowing bias based on the relationship:

$$\frac{{\displaystyle \sum _{i=1}^{N1}LS{T}_{clear\_sky\_i}}+{\displaystyle \sum _{j=1}^{N2}(LS{T}_{cloud\_gap\_j}+\Delta LST)}}{N1+N2}=LS{T}_{25\_km}$$

(3)

where the number of clear-sky MODIS LST pixels (LST_{clear_sky}) within a 25-km pixel is N1, and the number of cloud gap-filled pixels (LST_{cloud_gap}) is N2. We then assume that all cloud gap-filled MODIS pixels within the 25-km AMSR-2 pixel have the same shadowing bias (ΔLST). Therefore, ΔLST for all cloud gap-filled pixels within a certain 25-km grid can be calculated via rearrangement of Equation (3) to:

$$\Delta LST=\frac{LS{T}_{25\_km}\times (N1+N2)-{\displaystyle \sum _{i=1}^{N1}LS{T}_{clear\_sky\_i}}-{\displaystyle \sum _{j=1}^{N2}LS{T}_{clear\_gap\_j}}}{N2}, \mathrm{if} \frac{{\displaystyle \sum _{j=1}^{N2}\Delta LST}}{N1+N2}>RMS{E}_{unbias}$$

(4)

Equation (4) is constrained with the premise of $\frac{{\sum }_{j=1}^{N2}\Delta LST}{N1+N2}\le RMS{E}_{unbias}$. RMSE_unbias denotes the fitting errors illustrated in Figure 3 and can be identified as the threshold value for random (non-systematic) error between LST retrievals made at different spatial resolutions. When $\frac{{\sum }_{j=1}^{N2}\Delta LST}{N1+N2}$ lies below this threshold the LST error is mainly driven by this random error. The influence of this random error on ΔLST can be exaggerated if (N1 + N2) >> N2. To mitigate such influences, we use (N1 + N2) to substitute N2 as the new denominator in Equation (5) where $\frac{{\sum }_{j=1}^{N2}\Delta LST}{N1+N2}\le RMS{E}_{unbias}$. This means that when random error dominates the overall shadowing bias of the 25-km grid, the bias will be equally allocated to all (N1 + N2) MODIS pixels within the grid regardless of whether they are clear sky or cloudy, as in Equation (5):

$$\Delta LST=\frac{LS{T}_{25\_km}\times (N1+N2)-{\displaystyle \sum _{i=1}^{N1}LS{T}_{clear\_sky\_i}}-{\displaystyle \sum _{j=1}^{N2}LS{T}_{clear\_gap\_j}}}{N2+N1}, \mathrm{if} \frac{{\displaystyle \sum _{j=1}^{N2}\Delta LST}}{N1+N2}\le RMS{E}_{unbias}$$

(5)

Finally, the STDF cloud gap filled LST values are calibrated by the application of ΔLST, as found by Equation (4) or (5) depending on the amount of cloud cover within an AMSR-2 25-km grid cell, to give MODIS_PMBC.

4. Results

4.1. Overview of Cloud Gap-Filling

Our two-step cloud gap-filling processes improved the spatial coverage of the daily daytime and night-time MODIS_Clear data of Kenya from <50% to >90% for the cloud gap-filled images, as can be seen in Figure 4a. However, it is important to note that the re-visit cycle of AMSR-2 at low latitudes such as those of Kenya is normally greater than a day, so not all cloud gap-filled MODIS LST pixels were able to be bias-adjusted based on a near-simultaneous AMSR-2 observation. Figure 4b shows the proportion distribution of bias adjustment fraction (BAF) for gap-filled pixels of each date. It indicates that BAF on most (≥49%) of dates is between 0.25 and 0.5. Although not high, this rate is sufficient to provide a basis for comparison between the bias adjusted and non-bias adjusted cloud gap-filled LST data. From Figure 4c–h, we employed the day-time data of 22 October 2018 as an example to demonstrate the spatial inputs, outcomes, and corresponding reference data of our methodology. Spatial distribution maps of LST are shown respectively for the non-bias adjusted result in Figure 4f and for the bias adjusted result in Figure 4g. The map of their difference is demonstrated in Figure 4h, revealing that some gap-filled LST values in the middle and south of the country are lower when bias adjustment is included. A series of quantitative validation steps against the reference datasets are carried out in the following sections.

4.2. Validation against In Situ LST Ground Observations

We first compared MODIS_Clear observations from clear-sky pixels, MODIS_STDF from Step 1, and MODIS_PMBC from Step 2, to the contemporaneously collocated in situ LST data recorded at the four masts of the ILRI Kapiti Research Station site (Figure 5 and

Table 3. Statistics resulting from the LST comparison shown in Figure 5. ** indicates significance at p < 0.05.

	MODIS LST Data Type	RMSE (K)	Mean Bias (K)	r	N
Day	MODISClear	2.7	−1.6	0.96 **	198
	MODISSTDF	4.3	2.5	0.94 **	71
	MODISPMBC	2.6	0.2	0.97 **	71
Night	MODISClear	1.1	−0.4	0.93 **	115
	MODISSTDF	1.0	−0.6	0.93 **	119
	MODISPMBC	0.8	−0.2	0.94 **	119

). As the local observation times of MODIS and AMSR-2 are reported to have a variation between ±15 min from 13:30 for different dates and different observation swaths, the corresponding in situ data were determined by averaging all in situ LST samples from 13:15 to 13:45 local time (UTC+3) in the day-time, and all in situ LST samples from 01:15 to 01:45 local time at night-time. Moreover, in situ LSTs from Mast 2 and Mast 3 were averaged rather than being used independently, because these masts are in the same 1 km MODIS pixel. The use of multiple masts across and within different pixels (along with the geometric illumination model upscaling method applied, see Section Supplementary Figure S1) means that surface heterogeneity is well accounted for in the landcover under observation. We also examined the gap-filling results against the in situ LST data from Mast 1, Masts 2 and 3 (averaged), and Mast 4 as a time series (Figure 6 and Figure 7). For the time series demonstration, each blue point (MODIS_PMBC) is accompanied by a corresponding green point (MODIS_STDF). Bias adjustments herein are made only on dates where LST estimates are available from both MODIS and AMSR-2. Bias-adjusted outcomes on dates when AMSR-2 data were not available are not shown.

At night, all three datasets show a high degree of agreements with the in situ LST data, with absolute values of mean bias less than 1 K, and RMSEs below 1.5 K. The day-time performance of our cloud gap-filling approach appears poorer than that of night-time, this is likely due to the stronger spatial and temporal variations of LST at noon than at mid-night. However, MODIS_PMBC agrees significantly better with the in situ LST data record in the daytime than does MODIS_STDF. The former has an RMSE of 2.6 K compared to 4.3 K for the non-bias adjusted data, and this is similar to the 2.7 K RMSE shown by MODIS_Clear. Based on all analyses above, the MODIS_STDF values tend to overestimate the in-situ data in general, whilst MODIS_PMBC values are in closer agreement in the daytime. At night, most of the cloud gap-filled satellite estimates are close to the in-situ values, regardless of whether they have been bias-adjusted or not.

4.3. Validation against SEVIRI Geostationary LST

Our cloud gap-filled MODIS LST data record was also validated against the LSA-SAF Meteosat SEVIRI LST dataset, based on the mean LST of all MODIS pixels within a SEVIRI pixel. This was carried out for the full twelve months of the year. Similar to the evaluation employed against the in situ LST data, SEVIRI data were compared in three independent scenarios. These scenarios are: (1) comparison against only pixels containing clear-sky MODIS_clear observations; (2) against only pixels that have undergone cloud gap-filling up to Step 1 (MODIS_STDF), and (3) against only pixels that have undergone cloud gap-filling up to Step 2 (MODIS_PMBC) (Figure 8). The temporal means of the SEVIRI-derived LST obtained from the local time of 1:15 p.m./a.m. to 1:45 p.m./a.m. were used for comparison against the spatial means of MODIS_Clear estimates within SEVIRI pixels produced at the corresponding time. To ensure a consistent comparison, only SEVIRI pixels for which ≥14 spatially coincident clear-sky MODIS observations were available (i.e., either all or almost all MODIS pixels within the SEVIRI pixel are identified as cloud free) were considered. The same strategy was also applied to the evaluation schema of MODIS_STDF and MODIS_PMBC. From Figure 8, it is clear that daytime results differ more between clear-sky observations and cloud gap-filled estimates than do the night-time results. In the daytime, the best correspondence between MODIS and SEVIRI LSTs is achieved under clear-sky conditions, with a negligible bias and an RMSE of 3.2 K. The performance of the bias-adjusted cloud gap-filled result is slightly worse, with an RMSE of 3.6 K. However, this is an improvement upon the performance of the non-bias adjusted cloud gap-filled data, where an RMSE of 6.7 K and a positive LST bias of about 2.6 K were observed. At night, the performance of MODIS_STDF and MODIS_PMBC correction is similar, with RMSE no higher than 3.0 K.

Based on the yearly Kenya-wide comparison in Figure 8, further detailed analyses were conducted with respect to independent land cover types and daily timescales, as shown in Figure 9 and Figure 10 respectively. The ESA CCI 300-m land cover map in 2018 (see Figure 1c) was exploited to provide land cover type information for the entirety of Kenya. The areal fraction of each land cover type in Figure 1c (except for water bodies) was calculated for all SEVIRI-viewed pixels (~4 km) within the country. A SEVIRI 4-km pixel containing any given land cover type with area fraction lower than 0.7 is identified as a ‘mixed pixel’ and thus eliminated, while the remaining ‘pure’ pixels were employed to evaluate the potential influence of land cover variation on the gap-filled results, as shown in Figure 9. The similar performance among different land covers was consistent with the overall RMSE results reflected in Figure 8. This suggests that the bias-adjusted cloud gap-filling method is effective over a wide range of vegetation covers in obtaining LST estimates with competitive accuracy against MODIS clear-sky LST observations. The poorest performance of MODIS_PMBC is found over bare soil/rock during the day (RMSE ≈ 4.2 K) and over sparse vegetation at night (RMSE ≈ 4.0 K). However, this is acceptable considering that the difference between MODIS_pmbc and MODIS_clear RMSE is insignificant (<1.5 K).

The daily time series of difference between RMSE for gap-filled LST (MODIS_STDF and MODIS_PMBC) and RMSE for clear sky LST (MODIS_Clear) are demonstrated in Figure 10. The fraction of clear sky LST pixels over the study area for each day are also shown using the red bar on the right-hand y-axis. From the contrast between blue (MODIS_PMBC) and orange (MODIS_STDF) lines, we can see that the relative performance between the bias adjusted gap-filling algorithm against the non-bias adjusted algorithm is stable over different seasons of the year. The performances of the non-bias adjusted vs. bias adjusted algorithms are not significantly influenced by the fraction of clear sky LST pixels that are needed to build the gap-filling model.

4.4. Influence of Cloud Duration on PM-Based Calibration

As the PM-based bias adjustment was found to have more of an effect in the daytime than at night-time, we further investigated the impact of cloud duration during the day and night on bias correction performance. Firstly, data from the four selected dates highlighted in Figure 6 were employed to calculate an evaluation metric, Δbias. This metric is defined as the difference between the absolute bias of the MODIS_STDF (|ΔLST_STDF|) and that of MODIS_PMBC (|ΔLST_PMBC|). The Δbias metric (|ΔLST_STDF| − |ΔLST_PMBC|) was computed over the pixels in which ground validation masts are located. It can be used as an indicator of the accuracy improvement found through applying the PM-based LST bias-adjustment approach of Step 2. A time series indicating the presence or absence of morning cloud was then generated (1: cloud, 0: clear sky) from 7:30 a.m. to 1:30 p.m. local time for each of the four dates selected. Cloud state was provided by data from the sky-pointing radiometer deployed on the masts at ILRI Kapiti Research Station (see Supplementary Material Section S2, Figure S2). The dates selected for this were chosen at random from the days within four equal divisi)ons of the available time period. Based on this time series, a daily cloud duration fraction (CDF) was calculated, defined as the number of stable (occluded sky for a minimum of 15 min prior to the given minute) cloud minutes from 7:00 a.m. to 1:30 p.m., divided by the total number of minutes in that period. The results are displayed in Figure 11 and suggest that there is a relationship to be explored between Δbias and CDF, especially considering the low Δbias of 1.3 K and a small CDF of 0.14 on 2nd November 2018 as compared to the far higher Δbias (4.8 K) and correspondingly high CDF of 0.85 of the morning of 15 January 2019. Results in Figure 11 suggest that the performance of our bias adjustment methodology in the daytime is influenced by the length of the cloud period experienced during the morning prior to the early afternoon MODIS overpass whose cloudy observations were filled by our proposed method.

To corroborate the relationship found in Figure 11, a further analysis was carried out using the Meteosat SEVIRI LST data, as this data contains cloud masking information that is available over the entirety of Kenya at 15-min intervals. This provides a more extensive test than with the in-situ data alone. A CDF based on the SEVIRI cloud mask was calculated for both daytime and night-time periods, following the same method as detailed above. For the day-time data, the CDF period was slightly extended to between 7:00 a.m. and 2:00 p.m. in order to account for differences in the observation time of the satellite sensors. Similarly, for the night-time data the period between 7:00 p.m. of the previous date and 2:00 a.m.of the observation date was used. Due to the temporal resolution of SEVIRI (15 min) as compared to that of the in-situ data (5 min), only 28 unique CDF values can be obtained in each seven-hour period.

The outcome (Figure 12) is that in the daytime, the mean Δbias rises as the CDF increases from 0 to 1, while the standard deviations are close for all CDF-based groups. This phenomenon, together with the results in Figure 11, indicate that better day-time performance of our bias adjustment strategy is achieved under conditions of higher CDF ranges in the morning, especially when CDF exceeds 0.5. Essentially, longer periods of cloud cover prior to the afternoon Aqua MODIS overpass tend to provide improved performance of the bias adjustment approach. This is logical because, up to a certain limit, the longer a pixel is covered by cloud, the longer it is shadowed from direct sunlight and therefore the larger the temperature difference will be relative to clear-sky conditions. By contrast, at night there is only a very slight increase in Δbias with increasing CDF, probably indicating that the change in LST is more influenced by the presence/absence of solar radiation in the day time rather than by the presence/absence of downwelling atmospheric radiation from cloud layers at night.

5. Discussion

5.1. Improvement and Universality of the Gap-Filling Methodology

Overall, our evaluation of the cloud gap-filled MODIS LST data against the benchmark data, based on both the in situ and Meteosat SEVIRI LST data records, provides a consistent narrative. For day-time Aqua MODIS overpasses gap-filled LST data using the conventional STDF methodology shows a positive bias. This supports our supposition that, if used alone, the STDF approach overestimates LST due to simulating what the LST would be under clear-sky type conditions, whereas in fact the location is or has recently been cloudy. Our results also reveal the effectiveness of using PM observations for calibration of the gap-filled LWIR-derived LSTs, especially during the daytime.

The implementation of this improved two-step cloud gap-filling framework on LST only depends on satellite remote sensing datasets. Compared to the existing methods (Table 1), it is easier to be implemented and is independent from all auxiliary datasets (e.g., reanalysis data) other than satellite remote sensing inputs. Our approach requires much less dependence on the particular environmental characteristics of the study area than some prior works (Table 1), especially when we take into account Figure 10 which shows that the accuracy of this method is not strongly influenced by the varied daily area fraction of cloudy pixels. The performance is stable under periods of continuous cloudy/rainy weather, as can be seen from the time series (e.g., between late November and early December 2018) in Figure 7. Taken together, these tests indicate the greater universality of the proposed method beyond its current study area, as compared to other existing methods.

It is noticeable that we find a relatively inconsequential cloud insulation effect on LST bias correction at night (Figure 12b), as opposed to the effect seen in the morning-to-noon window. The difference between daytime and night-time observations has rarely been discussed in previous studies related to cloud gap-filling of an LWIR-derived LST dataset, with some authors ignoring the problem by only tackling night-time cloud gap-filling [45]. Other works may have failed to find significant day/night differences, because they calibrated against air temperatures rather than the more physically accurate surface temperatures [30,70]. Our findings suggest that the difference between bias-adjusted and non-bias adjusted gap-filled LST, driven by cloud insulation, is less sensitive to ground long-wave outgoing radiation and downwelling atmospheric radiation at night, as compared with the shadowing of solar short-wave incoming radiation during the day.

From Section 2.1 we can see that Kenya is representative of a variety of different geographical settings due to the varied land covers, climate zones and topographical conditions found across the country. Based on the stable evaluation of results across the different land covers in Figure 9, we suggest that the improved two-step cloud gap-filling framework should be effective in the vast majority of vegetation covers ranging from sparsely vegetated soils and grassland to forest found in low- and middle-latitude regions, where most of the world’s human population resides. However, our current method is likely to be unsuitable for use in several land cover types, namely desert and snow/ice covered areas. This is because the temperature difference between microwave observations and LWIR skin surface observations is very high in such landcovers, whilst a conventional linear calibration (as applied in Equation (2)) insufficiently accounts for the depth-induced temperature differences encountered in these areas [65].

5.2. Uncertainty and Limitations of the Current Study

It is important to stress that neither the benchmark datasets (field radiometer derived LSTs and the Meteosat LST product) or the MODIS_Clear data employed in this study represent the ‘absolute truth’ in terms of surface temperature observations, and we do not consider them as such. Each sensor is subject to its own inherent measurement errors. Moreover, measurements observed by the different sensors may differ on account of the heterogenous spatial scales at which observations are made and aggregated over sensor footprint geometry, as well as changes in environmental conditions occurring between non-temporally contemporaneous measurements. As such, we focus here upon the relative difference in performance metrics between the gap filled data (MODIS_STDF and MODIS_PMBC) and MODIS_Clear, rather than the absolute performance of these gap filled products vs. the SEVIRI and field radiometer benchmark LST datasets themselves. For example, using the Meteosat LST product as a benchmark, a relatively large absolute daytime RMSE (3.6 K) was obtained for pixels gap-filled using the MODIS_PMBC dataset. This value is however much closer to the absolute daytime RMSE value (3.2 K) obtained from the ‘raw’ MODIS_Clear dataset than the RMSE value (6.7 K) obtained using the MODIS_STDF dataset, highlighting the effectiveness of our novel cloud gap-filling two-step framework. While the absolute RMSE found between the MODIS_Clear and the benchmark datasets during the daytime was observed to be high (>3 K), exploration of the causes of this discrepancy is beyond the scope of the current study. Furthermore, the validation effort of Gottsche et al. [54] found an average RMSE of 3.2 K for Meteosat SEVIRI data over a similar biome, indicating RMSE values of this magnitude are currently inherent in LWIR LST products in semi-arid savannah regions.

A limitation of the current study is the relatively small number of field-based LST observations available to validate our cloud gap-filling methodology. Radiometer measurements made at the ILRI Kapiti Research Station are temporally limited, spanning a ~5-month period. Fortunately, however, the period is sufficient to capture one entire climatic and vegetation cycle of local land surface (Figure 1b). In the spatial domain, the validation sites only coincide with ground footprints of three distinct 1-km MODIS pixels that are in close geographic proximity to one another relative to the size of the study area. Nevertheless, the establishment and maintenance of even a single site represents a substantial achievement. Indeed, only one other satellite validation site capable of producing similar LST validation data exists anywhere on the African continent and it is not currently operational [54]. In addition, the limited availability of ground validation data has been compensated for to some extent by the intercomparison of the MODIS-derived products against the independent Meteosat SEVIRI LST product [56] at the national scale through an entire climatic year. The intercomparison between the MODIS derived products and the ground/Meteosat-SEVIRI validation data provided consistent results. Prior studies on cloud gap-filling have primarily relied on validation strategies that artificially remove data to create “vacant sub-regions” of the imagery to replicate the presence of cloud gaps but for which the ‘true values’ of the data are actually known [26,71]. However, this strategy is unsuited to evaluating cloud gap-filling of LST data records, since the pixels for which the “true values” of LST are known have actually been observed under clear sky conditions and have thus not experienced the radiation interception effects of clouds that would be the case for ‘true’ cloud-covered LST pixels.

The sensitivity analysis results reported in Table 2 indicate that a theoretical error as large as 100 m in the DEM only leads to a difference of 0.3–0.4 K in LST on average. As the vertical error of the SRTM DEM is reported to be <16 m [72], the uncertainty of elevation data can be effectively neglected. On the other hand, the influence generated by uncertainty in the NDVI data can be much more significant. This is because the generation of daily-scale NDVI data from the 16-day composite MCD43A4 product relies on an assumption that vegetation change is not significant during the 16-day window. Unfortunately, vegetation changes rapidly during the onset of the rainy seasons in Kenya, which may result in a not-insignificant estimation error in NDVI estimates at certain dates. The resultant gap filled LST data is sensitive to such errors, especially during the daytime (as described in Section 3.1.2). Consequently, the differing performance between MODIS_STDF and MODIS_PMBC in the daytime cannot only be attributed to the effects of cloud as performance can also be impacted by NDVI uncertainty, particularly during the onset of the rainy season when rapid vegetation changes occur. However, the impact of cloud effects can still be identified as the primary driver of the performance differences between MODIS_STDF and MODIS_PMBC, based on the positive correlation observed between the CDF and Δbias of the day-time data in Figure 11 and Figure 12. Furthermore, the impact of the bias correction did not get significantly smaller at the rainy periods (typically between April-June and between October-December) in Figure 10.

Finally, despite the potential bias due to clouds as well as the uncertainty of NDVI input not being very significant for the night-time observations in this study, more validation experiments are required in other study areas to test the wider applicability of this finding. Currently, we recommend applying the PM bias correction process to both daytime and night-time data in future studies, in the event that larger night-time errors can arise from factors that were not present in the current study.

5.3. Future Work

Several points remain where further improvements to our approach might be investigated:

(1) Biases calculated at the coarse (25-km) resolution of AMSR-2 are linearly allocated to each cloud gap pixel at the fine (1-km) resolution in the current study, based on an assumption that cloud geometry (e.g., cloud top height, cloud thickness, etc.) is homogenous within the AMSR-2 pixel. However, a bias-adjustment approach that considers the influence of cloud geometry differences at the AMSR-2 sub-pixel scale may offer further performance improvements to the current two-step LST cloud gap-filling framework.

(2) If applying the methodology at continental to global scales, the derived relations between AMSR-2 and MODIS LSTs (Equation (2)) should be re-explored because they are likely dependent on seasonal features and climate types that will vary at the continental/global scales.

(3) As the re-visit cycle of AMSR-2 is greater than one day at low latitudes, better spatial coverage may be achieved by the fusion of AMSR-2 data and data from other PM radiometers.

(4) Cloud gap-filling of Terra MODIS LST data (MOD11A1) could be similarly conducted using data from a PM radiometer similar to AMSR-2, but with local time of overpasses around 10:30 a.m./p.m. One candidate is the microwave radiation imager onboard the Fengyun-3C satellite [73].

6. Conclusions

In this study, we have presented an effective two-step framework for improving the cloud gap-filling performance of the LWIR-based daily land surface temperature (LST) data delivered by Aqua MODIS. The primary innovation of the framework lies in the coupling of a conventional “spatial-temporal data fusion (STDF)” methodology to a cloud shadowing bias-adjustment process, based on passive microwave (PM) LST data derived from the AMSR-2 instrument. We have evaluated the resulting STDF-based cloud gap-filled and PM-bias adjusted LST data against both in situ LST data and the geostationary SEVIRI-derived LST data record. In the daytime, the STDF-based outcomes show an RMSE of 4.3 K against the in-situ data and 6.7 K against the SEVIRI LST data. The RMSEs of the PM bias-adjusted outcomes are 2.6 K and 3.6 K respectively. At night, the performance of bias-adjusted outcome and that of STDF-based outcome are very similar.

Overall, we conclude that, compared to the STDF approach alone, our two-step framework provides a means to improve the performance of cloud gap-filling on LWIR-derived LST through fusion with PM data, although the degree of improvement is far more significant for daytime compared to night-time observations. Finally, we have shown that the day-time accuracy improvements of our two-step approach to cloud gap-filling are increasingly apparent (relative to conventional STDF methods) under the condition of increased cloud cover residence time in the morning-to-noon timeframe. The proposed method has been shown to be stable in terms of RMSE change through time, even under spatially and temporally continuous pixel loss conditions, across a wide range of land cover types. In the future, this method could be applied to the next generation of sub-1 km all-weather LST products at the continental or global scale.