A New Fully Gap-Free Time Series of Land Surface Temperature from MODIS LST Data

Abstract

Temperature time series with high spatial and temporal resolutions are important for several applications. The new MODIS Land Surface Temperature (LST) collection 6 provides numerous improvements compared to collection 5. However, being remotely sensed data in the thermal range, LST shows gaps in cloud-covered areas. We present a novel method to fully reconstruct MODIS daily LST products for central Europe at 1 km resolution and globally, at 3 arc-min. We combined temporal and spatial interpolation, using emissivity and elevation as covariates for the spatial interpolation. The reconstructed MODIS LST for central Europe was calibrated to air temperature data through linear models that yielded R² values around 0.8 and RMSE of 0.5 K. This new method proves to scale well for both local and global reconstruction. We show examples for the identification of extreme events to demonstrate the ability of these new LST products to capture and represent spatial and temporal details. A time series of global monthly average, minimum and maximum LST data and long-term averages is freely available for download.

1. Introduction

Remote sensing based time series are becoming increasingly available, and this tendency will continue to grow not only because of new Earth observation satellites being launched, but because of the availability of new methods to harmonize their data [1,2] and reconstruct incomplete records [3,4,5,6,7] along with the growing demand of different sectors for the monitoring of environment, analysis of trends and patterns, and forecasting. In this scenario, air temperature is as an essential climatic and ecological driver, one of the most important variables in climate research and global change. It controls many biological and physical processes between the hydrosphere, atmosphere and biosphere [8]. Therefore, the availability of spatially and temporally continuous temperature time series at high spatial and temporal resolution is crucial for a wide range of disciplines including hydrology, meteorology and ecology [9,10], as well as for numerous applications fields.

While meteorological stations can provide accurate air temperature measurements (near-surface temperature, usually measured about 2 m above ground), they only represent single locations with their immediate surrounding area of influence. Moreover, the distribution of meteorological stations is far from optimal for most areas of the world. This type of point data needs to be spatialized (usually performed through geostatistical approaches or thermodynamical modelling) to characterize regions. Regardless of the spatialization method used, unrepresentative smooth spatial patterns may occur due to the lack of dense data [11,12,13]. In the end, the obtained precision depends on the quality, representativeness and spatial distribution of the input network(s) of stations [12,13,14]. Besides, processing such data can be quite time-consuming regarding computation time (especially to obtain higher spatial resolution datasets), difficult to automate (i.e., data requires a lot of “curation”), and therefore, it is not easily updated with new incoming data [14,15,16].

Land surface temperature (LST) as recorded by remote sensing instruments on board of satellites, on the other hand, is intrinsically spatial and provides the spatial coverage that meteorological stations lack. As such, satellite LST is able to capture at greater detail local differences in temperature originating from varying meteorological conditions, environmental differences and/or active heat sources (e.g., urban areas, land cover classes) [15]. The most commonly used LST products are those from the MODerate resolution Imaging Spectroradiometer (MODIS) instruments on board of Terra and Aqua NASA satellites. The combination of these satellite data provides spatial estimates of LST at high temporal (four Earth coverages per 24 h; approximately 01:30, 10:30, 13:30, and 22:30, local solar time) and high spatial resolution (nominally 1 km) across the world [17,18]. Remotely sensed LST, however, as all remotely sensed data in the optical/thermal portions of the spectrum, suffers from cloud contamination, i.e., occurrence of gaps produced by cloud cover and/or presence of outliers in undetected cloudy pixels [5,19] which reduces the usage rate and hinders any subsequent interpretation [20]. The percentage of valid LST values varies in space and time depending on the region, but in any case, the gaps must be filled to render the data useful for other applications or research fields (e.g., for deriving climatic indices, modeling, etc.).

Previous LST related studies have attempted to obtain spatially and temporally continuous estimates of temperature for areas of varying sizes using different methods, time series of different lengths and different LST products. There have been two approaches aiming at two different outputs: (1) to use MODIS LST as one of the explanatory variables in statistical models to obtain gridded time series of air temperature [9,15,16,21,22]; and (2) to directly reconstruct LST products either with or without covariates [2,5,19,23,24]. In the first case, approaches go from linear regression models to more complex spatio-temporal regression-kriging interpolations, from regional (country scale) to global predictions and from only one to 10 years of output air temperature time series. In all of the first cases, these approaches use MODIS LST products (either daily or eight-day composites, both from Terra and Aqua satellites or only from Terra) to enhance the spatial detail attained with meteorological stations with the help of one or more predicting variables in a statistical model. In the second case, the reconstruction of satellite LST data, diverse approaches have been used. For example, a MODIS LST time series has been reconstructed and compressed using temporal Fourier transform by Scharlemann et al. [23] who processed MODIS LST at a full global extent for the years 2001 to 2005 and extracted parameters for the annual, bi-annual and tri-annual cycle through Fourier transform at a resolution of 1 km. Another temporal interpolation method replaced/estimated missing Aqua LST values with Terra LST values [24], however the resulting dataset still contained gaps. On the other hand, Neteler [19] used 3D spline interpolation to fully reconstruct LST (2000–2008, daily, at 200 m spatial resolution) for a small alpine region in northern Italy with an algorithm optimized for very complex topography. For this latter reason, the method could not be extended to larger areas, and motivated the study by Metz et al. [5], combining temporal averaging with multiple regression and spline interpolation to produce a seamless and gap-free daily reconstructed LST for Europe at 250 m spatial resolution for the period 2000–2013 (“EuroLST” dataset). For a summary of previous approaches to reconstruct LST time series, see Table 1 in Metz et al. [5].

All papers cited above have made use of LST MODIS collection 5. In 2017, NASA completed the reprocessing of the entire MODIS archive and made available a new collection 6 of MODIS products. This new version includes different kinds of fixes and improvements, including removal of cloud-contaminated LST from MODIS level 2 and 3 products; updated coefficients for the split-window algorithm; adjustments in the classification-based surface emissivity values, etc. [18]. On the software side, new methods became available in the last years for the gap-filling of time series. For example, the long standing free and open source software GRASS GIS [25,26], offers new enhanced modules and methods to efficiently handle and process huge amounts of geospatial time series data [27,28,29,30,31].

In this study, we present a new method to reconstruct MODIS LST collection 6 data that combines temporal interpolation to fill small gaps in time and spatial interpolation with covariates to fill the remaining gaps in space. We used a digital elevation model and one of the emissivity bands delivered along with MODIS LST products as covariates in the spatial interpolation step. The use of emissivity allows us to characterize the thermal properties of different land covers at the same temporal resolution at which we have LST records while providing complementary information to that contributed by elevation data, a well known predictor of temperature. To our knowledge, emissivity has not been used before in the reconstruction or interpolation of temperature surfaces. We applied the new method to the MOD11A1/MYD11A1 LST products for Central Europe (2003–2016) with a spatial resolution of 1 km, as well as to global LST products at 3 arc-min (nominally 5.6 km at the equator). Therefore, this new approach delivers seamless and gap-free LST time series at two different spatial resolutions with comparable satisfying quality and can easily be applied as-is to other parts of the world. Furthermore, it can be automated and continuously updated (with new incoming MODIS LST data) without much delay.

Having a long time series of daily observations allows for different aggregations and extraction of other (synthetic) variables. In this paper, we use the identification of extreme events to demonstrate how this new LST product covering a relatively short time span of 14 years is able to detect well-known anomalous events with high spatial and temporal detail.

2. Materials and Methods

2.1. Data

We used the MODIS LST products MOD11A1/MYD11A1 collection 6 acquired from the LP DAAC data pool (https://lpdaac.usgs.gov/dataset_discovery/modis/modis_products_table) starting with the year 2003, the first year fully covered by both products. These two products combined provide four LST values per day at a spatial resolution of nominally 1 km. We used the tiles h18v03, h18v04, h19v03, h19v04 covering central Europe (data provided in MODIS sinusoidal projection, the extents are 2223.9 km × 2223.9 km with a raster resolution of 926.6 m). The spatial coverage comprised a total of 5,760,000 pixels, of which 4,342,823 pixels were on land (75.4%). Additionally, we applied the same reconstruction method to the global MODIS LST products MOD11C1/MYD11C1 at a spatial resolution of 3 arc-min (approximately 5.6 km at the equator). As elevation model we used the Global Multi-resolution Terrain Elevation Data 2010 grids (GMTED2010; resolution 30 arc-s) available from the U.S. Geological Survey (https://lta.cr.usgs.gov/GMTED2010). Monthly ground meteo-station records of air temperature were obtained from the German Weather Service (DWD) for the year 2016 (the number of active DWD stations per month varied between 483 and 489, see

Table 1. Difference of the reconstruction over simulated gaps with the raw data in K. For reference, the difference of the reconstructed land surface temperature (LST) values without simulation (only outlier filtering) to the raw data is also provided.

Hole Size in Pixels	Simulation-Raw LST	Reconstructed-Raw LST
9	−1.17	−0.42
61	−1.91	−0.07
161	−2.44	−0.07
305	−2.57	0.03
509	−2.63	0.01
761	−2.61	0.03
1069	−2.65	0.02
1425	−2.73	0.01

in the Results, Section 3.5). Land cover information was obtained from the MODIS product MCD12Q1 collection 051 for the year 2013 (which is the latest year and collection available at the time of writing), using land cover type 4: MODIS-derived Net Primary Production (NPP) scheme (spatial resolution: 500 m).

2.2. LST Reconstruction Method

All pixels with an LST error >3 K were filtered out using the corresponding MODIS LST QA layers (quality assurance layers in the MODIS LST HDF files). Thus missing values were all pixels for which no LST value had been produced as well as all pixels with an LST error >3 K. We used this rather liberal threshold because preliminary tests showed that a lower threshold for LST error would discard also LST values that appear to be correct. Therefore we decided to keep potentially valid LST values at the expense of including more outliers. Outlier filtering is performed at a later stage.

Variation of daily surface temperature can be quite high. Any interpolation in time should therefore allow for fast fluctuations and should avoid smoothing of the time series in order to preserve extreme events. For temporal interpolation we used local weighted regression (LWR; GRASS GIS add-on r.series.lwr [29,30]) with polynomial order 2, considering the five nearest neighbors in time. Since the intra-day and inter-day fluctuations might be mixed up in the interpolation by the existence of any missing values, we treated each daily time step (01:30, 10:30, 13:30 and 22:30) as a separate time series. Intra-day fluctuations in LST can vary considerably among days, i.e., day-time and night-time temperatures can be nearly identical or quite different. Therefore, temporal interpolation of all four daily overpasses together over several days would require a priori assumptions on intra-day fluctuations for each day separately, which are impossible to determine if none or only one observation for a certain day is available. There exist models based on diurnal temperature cycle, but those require several temperature measures per day (see, for example, Liu et al. [32]), which are not available in the case of MODIS sensors. Gaps larger than seven days were not interpolated and left as missing values. We used thin plate spline interpolation (TPS; GRASS GIS add-on r.resamp.tps [29,30]) with elevation data from GMTED2010 and band 31 emissivity from the MODIS LST products as covariates to fill any remaining gaps in space. For details on TPS, see [33,34,35]. Here we used an adaptive window size determined by the spatial distribution of the nearest 150 points. Preliminary tests showed that considering emissivity increased the contrast between urban areas, vegetated areas, and water bodies. Emissivity also exhibited missing values due to cloud cover. However, since emissivity does not change as long as land cover remains constant, it can be easily reconstructed with temporal interpolation. We used LWR with polynomial order 1, using the nearest 7 neighbors in time and filling gaps up to 30 days long to obtain gap-free emissivity that was then used as covariate in the spatial interpolation of LST. Low outliers were defined as values larger than Q3 + 1.5 × (Q3 – Q1), with Q1 and Q3 being the first and third quartile of reference data [5,19]. As reference data we used the time steps of the same day, the last time step of the previous day and the first time step of the next day. The workflow of our MODIS LST reconstruction method is illustrated in Figure 1. Importantly, since valid pixels in the MODIS Land Surface Temperature products represent clear-sky conditions, our reconstruction method as well as most other studies attempting to fill missing values in LST, also represent clear-sky conditions.

To evaluate the robustness of the reconstruction method, we created gaps of different sizes in the raw LST data, reconstructed these gaps and compared the resulting reconstructed LST values to the raw LST values. These differences of the simulated gaps were then compared to the differences of the originally reconstructed LST values to the raw LST values.

2.3. Software

For mosaicking, we used the GDAL tool gdalbuildvrt (GDAL, http://www.gdal.org). All further processing was performed in GRASS GIS 7.2 [25,26]. The MODIS LST time series was processed on a Linux cluster computing system, managed by the “Son of Grid Engine” job scheduler (SGE, https://arc.liv.ac.uk/trac/SGE).

2.4. Assessment of LST Derived Indicators: Identification of Extreme Events

Considering climatological standards, our time series is quite short, spanning only 14 years so far. Therefore, we wanted to test if it is possible to identify extreme events with this short time series as reference.

Extreme events in climate science are usually defined by three criteria, i.e., rarity, intensity and severity. Definitions of “rare” vary, but an extreme weather event would normally be as rare or more rare than the 10th or the 90th percentile. Rarity can thus be assessed with the calculation of the corresponding quantile for the current observation according to a reference time series. A quantile of 0.5 means that the current observation is equal to the median, values larger than 0.5 mean larger than median and values smaller than 0.5 smaller than median. A value of 0.9 for example means that the current observation is equal to the 90% percentile of the reference time series. Intensity, in the context of extreme events, refers to their magnitude, i.e., events that have large magnitude deviations from the norm. Intensity can be assessed with standardized anomalies: the difference of the current observation to the long-term average, divided by the long-term standard deviation. These standardized anomalies are thus the number of standard deviations away from the average, with the sign indicating if the current observation is larger or smaller than the average. Finally, severity is related to the socio-economic losses that a rare and intense event may cause [36].

3. Results

The newly developed LST reconstruction method was applied to central Europe covered by four MOD11A1/MYD11A1 tiles with approximately 1 km resolution and furthermore globally applied to MOD11C1/MYD11C1 images with a resolution of 3 arc-min. Regarding central Europe, there are 20,456 time steps covering the years 2003 to 2016 (four records per day), corresponding to 384 GiB of selected raw data (i.e., LST, quality assurance, and emissivity layers).

The results consist of two parts: (1) an analysis of the MODIS LST products used here and of the LST reconstruction method; and (2) the identification of extreme events as an exemplary use case.

3.1. Missing Values

For the four tiles covering central Europe during 2003–2016, 61.4% of the temperature values were missing from the original MODIS LST collection 6 products (after filtering with the corresponding QA flags as described in the Methods Section 2.2). The smallest percentages of missing values occurred at high altitudes and along edges of water bodies. The water bodies themselves had a slightly higher percentage of missing values than their edges, probably due to fog over water bodies. The highest percentage of missing values appeared in low-lying inland areas away from oceans such as central Germany.

An inspection of months aggregated over all available years showed that most missing values appeared in the winter months November–February (71–74%), and least missing values appeared in the summer months June–September (49–55%). A comparison of missing values in the months January, April, July, and October (average over all years) revealed that not only the percentage of missing values but also their spatial distribution differs between months (Figure 2). In January (74.4% missing), most missing values appeared in lowland areas (up to 1000 m a.s.l.). In April (56.9% missing), most missing values were observed at higher altitudes. In July (49.8% missing), the lowest percentage of missing values appeared in southern Europe and, in October (62.3% missing), a high percentage of missing values appeared in north-eastern Europe and parts of France, Switzerland, Germany, and Austria. The seasonal differences in the amount of missing LST values is most probably caused by climate conditions of the current region (central Europe). For other parts of the world, different seasonal patterns in the amount of missing LST values are to be expected. To understand if missing values were equally common and equally distributed in different times of the day, the number of missing values was counted for each pixel and each time of the day over all processed days (Figure 3 with 01:30 and 13:30). High altitude areas such as the Alps showed few (<40%) missing values at night and average (50–60%) missing values during the day, while northern parts of central Europe had about 10% more missing values during the day than at night.

Temporal interpolation with LWR reconstructed 68% of all missing values. In summer, more missing values could be temporally reconstructed than in winter, e.g., averaged over all years, 85% of missing values were reconstructed in August and 48% in January.

An example for the temporal and spatial reconstruction of missing values is given in Figure 4. The proportion of pixels filled with temporal interpolation varies between time steps, depending on the coverage of nearby days.

3.2. Outliers

For the four tiles covering central Europe from 2003 to 2016, 5.9% of all surface temperature values were identified as outliers in the original MODIS LST products. Generally, the highest percentage of outliers appeared at high altitudes and in water bodies. Regarding monthly outlier counts aggregated over years, most outliers were detected in the months April–July (6.2–6.4%), and the lowest percentage of outliers was detected in the months November–February (5.2–5.6%).

Figure 5 shows a comparison of the percentage of outliers in the months January (5.2% outliers), April (6.4% outliers), July (6.4% outliers), and October (6.0% outliers) aggregated over the years 2003 to 2016. In January, most outliers appeared in northern and north-eastern Europe. A notable exception are water bodies where only few outliers were detected. In April and July, the highest percentage of outliers appeared at high altitudes and over water bodies. In October, outliers appeared at high altitudes and also in eastern and northern Europe, while water bodies did not show peculiar percentages of outliers.

To compare the four different overpass times per day in terms of outlier occurrence, the outlier counts were aggregated over all days and years for each overpass time of the day, and then divided by the corresponding number of valid LST pixels. The most obvious difference in the percentage of outliers was observed over water bodies where only few outliers at night and many outliers during the day occurred, e.g., in southern Sweden and at the border between Estonia and Russia (Figure 6).

3.3. Effect of Emissivity

Since emissivity differs among land cover types, particularly among urban areas, water bodies, and vegetated areas, we compared reconstructed LST data with land cover types from the MODIS product MCD12Q1. The effects of emissivity on LST reconstruction are clearly visible for urban areas and larger water bodies. In winter, both urban areas and water bodies appear hotter during the day due to the high specific heat capacity when compared to both vegetated areas and the result of LST reconstruction without using the emissivity layer (Figure 7).

3.4. Simulation of Missing Data

We evaluated the robustness of the reconstruction method by creating holes of different sizes in the raw LST data and then reconstructing these holes. The difference of the reconstruction over simulated gaps to the raw data was then compared to the difference of the original reconstruction (only outlier filtering, no gap filling) to the raw data (Table 1). As test data, we used an area without gaps in the Aqua daytime overpass at 13:30 for 23 June 2016. While the difference of the original reconstructed LST values to the raw data was always less than 0.5 K, the results over simulated gaps where between 1.2 and 2.7 K smaller than the raw LST values.

3.5. Calibration to Station Data

We compared time series of the new LST reconstruction considering emissivity (this paper), the corresponding raw data of the MODIS 11A1 products, the old “EuroLST” reconstruction not considering emissivity [5], and air temperature records from the German Weather Service (DWD) ground stations, using monthly averages for the year 2003 for the city of Bonn, Germany. Both reconstruction methods are closely following the raw LST data, but tend to be higher than the raw LST data due to the removal of low outliers. It appears that values of the new LST reconstruction are equal to or slightly lower than those of the old LST reconstruction. All LST time series are slightly above the values of ground station records in summer and slightly lower than ground station records in winter (Figure 8), indicating that for each time step, surface temperature needs to be individually calibrated to air temperature because the bias of surface temperature to air temperature is different for each time step.

We used monthly averages for 2016 as test data for calibrating reconstructed LST to air temperature recorded at ground stations. Air temperature values were obtained from DWD meteorological station records, which varied for 2016 from 483 to 489 stations among different months. Because the bias of surface to air temperature is different for each time step, surface temperature was calibrated to air temperature for each month separately with a linear model (

Table 2. Calibration of reconstructed monthly MODIS LST data for 2016 against air temperature data of meteorological stations from the German Weather Service (DWD) and elevation. Values in brackets are the results without elevation as additional explanatory variable. AIC: Akaike information criterion; RMSE: root mean square error.

Month	n Stations	R2	AIC	RMSE
Jan	486	0.85 (0.83)	−433.2 (−379.1)	0.64 (0.67)
Feb	487	0.86 (0.82)	−690.6 (−554.8)	0.49 (0.56)
Mar	487	0.90 (0.84)	−795.2 (−580.6)	0.44 (0.55)
Apr	484	0.85 (0.62)	−688.4 (−231.6)	0.49 (0.78)
May	484	0.91 (0.72)	−732.7 (−193.7)	0.47 (0.82)
Jun	483	0.89 (0.74)	−698.7 (−302.2)	0.48 (0.73)
Jul	484	0.86 (0.67)	−643.2 (−220.3)	0.51 (0.79)
Aug	486	0.84 (0.69)	−656.9 (−330.2)	0.51 (0.71)
Sep	489	0.86 (0.61)	−685.2 (−207.1)	0.49 (0.81)
Oct	485	0.90 (0.52)	−937.0 (−183.0)	0.38 (0.82)
Nov	489	0.83 (0.79)	−701.9 (−594.9)	0.48 (0.54)
Dec	486	0.71 (0.40)	−91.3 (258.7)	0.90 (1.30)

). The influence of elevation as additional variable in the linear models was also investigated. We used R², root mean square error (RMSE), and Akaike information criterion (AIC) to evaluate the goodness of the calibration. Calibration succeeded with R² in the range of 0.83 to 0.91 and RMSE in the range of 0.38 to 0.64, excluding December. The month of December presented by far the lowest R² (0.71), the highest RMSE (0.90), and the highest AIC (−91.3). The percentages of missing values and outliers were not unusually high for December 2016. The poor match between MODIS LST and the reconstruction method and air temperature recorded at ground stations for this month could be caused by a failure of the reconstruction method or by errors in the ground station data. Generally, elevation improved the calibration slightly. For a few months (April, September, October, and December), the improvements when including elevation were considerable (ΔR² > 0.2).

The new LST reconstruction method was applied to both central Europe using the MODIS products MOD11A1/MYD11A1 and globally using the MODIS products MOD11C1/MYD11C1. These MODIS LST products are created independent of each other. Comparing an intra-day time series for December 2016 revealed some overshoots in the global reconstruction (Figure 9b), where the global reconstruction had lower values at night and higher values during the day. These differences could well be caused by differences in the raw data (Figure 9a).

3.6. Extreme Events

Extreme events can be identified by comparing current observations to an adequate reference time series. The reconstructed LST time series spanning 14 years is relatively short compared to climatological standards. Nevertheless, extreme events can be identified using this time series as reference. Extreme events can be highly localized in both space and time (a small area affected in a short period of time). The high spatial and temporal resolution of the reconstructed MODIS LST time series allows to detect such localized events as illustrated with the example of the heat wave in Europe in August 2003. This heat wave is visible in both the global and the European reconstructed LST time series (Figure 10). Both time series were able to show the high intensity of the heat wave with a very similar spatial pattern despite different spatial resolutions (nominally 5.6 km vs. nominally 1 km). Notably, the global time series (Figure 10a) shows that this heat wave was restricted to western Europe, illustrating that the reference time series spanning 14 years is long enough for spatial discrimination.

Next, we compared average monthly temperatures for January and July 2016 with the time series for the years 2003–2016 with regard to intensity and rarity. For January 2016, temperatures in western and southern parts of our study area were considerably higher than the 2003–2016 average and median (Figure 11). These high temperatures were thus an intense and rare event. In July 2016, eastern parts of our study area were warmer than the 2003–2016 average and median. Similar analyses on extreme events can also be performed on shorter time scales such as weeks or days.

4. Discussion

Land surface temperature (LST) is a key variable in the physical processes of surface energy and water balance at different scales [37,38]. As such, many research fields and applications, namely evapotranspiration, climate change, hydrology, vegetation monitoring, urban climate, public health, among others, require gap-free time series of this variable as inputs to better understand the spatio-temporal changes of different study targets [39,40,41,42].

In this paper, we presented a new method to obtain a fully gap-free time series of gridded daily surface temperature from MODIS collection 6 LST products (tiled MOD11A1/MYD11A1 data at 1 km resolution and global MOD11C1/MYD11C1 data at 3 arc-min resolution). The spatial support has been enhanced by using elevation and emissivity as covariates. This resulted in plausible differences among distinct land cover types (e.g., cities and water bodies) and also between low and high lying areas. The former were evident when we compared reconstructed LST including and omitting emissivity as covariate in the spatial interpolation step (Figure 1). By using one of the emissivity bands delivered with the MODIS LST products, we were able to account for spatial differences and temporal changes in land cover, thus catching for example seasonal changes due to snow cover or loss of leaves in trees. This would not have been the case if we would have used constant LULC maps produced annually (i.e., MCD12Q1 of which the production ended with the year 2013).

The highest percentage of missing values in central Europe was observed in the month of January (average based on 2003–2016). While this is probably an expected result given the climatic conditions common of that month in central Europe, it means that short-term fluctuations in temperature are probably not well represented for gaps that are large both in time and in space, eliminating any possible unrecorded short-term fluctuations. Interpolation of gaps that are small in time or space can preserve short-term fluctuations. A 68% of all missing values present in MODIS LST products were interpolated in time through LWR, while the remaining gaps were filled by means of TPS in the spatial interpolation step. The fact that there were more missing values during day than at night might be explained by the MODIS cloud detection algorithm that uses thermal and reflective bands for day overpasses, but only thermal bands for cloud detection at night, i.e., clouds are more easily detected during the day than at night [43,44]. The difference in the percentages of missing values among water bodies and their borders might also be related to the MODIS cloud detection algorithm that uses different spectral thresholds for different types of covers [45]. Therefore, if the land cover layer used is not up-to-date or water bodies change their size, differences in cloud detection over water bodies and their borders might be expected. This was observed by Kotarba [43] while analyzing MODIS collection 5 cloud product for the period 2004–2009 over Poland. However, they detected more clouds over water bodies than over water body borders. According to the authors, the approach used by the MODIS project for operational production of level 3 cloud data (used afterwards for other level 3 products such as LST) tends to over-estimate cloud amount by 6.5% in Europe (annually). In MODIS collection 6 some of the artifacts (decreased and spatially biased data availability) have apparently been reduced but not eliminated [46].

A high percentage of outliers was observed over water bodies in April and July when outlier percentages were aggregated per month, and at 13:30 when outlier percentages were aggregated per overpass time. Considering that water bodies do not change their temperature quickly (i.e., within a day), it is possible that the applied outlier detection might have resulted in false positives for water bodies, particularly at mid-day (13:30) and in summer, when land surfaces are typically hotter than water bodies. On the other hand, water bodies are frequently covered by clouds or fog, but not that many missing values were detected over water bodies as compared to other surfaces. Therefore, the observed higher percentage of outliers over water bodies might be explained by undetected clouds (a task marked as difficult in MOD35 ATBD). These undetected clouds will have significantly lower surface temperature and thus are identified as outliers by our method. This might partially explain why the percentage of outliers detected was higher over water bodies, but it does not fully explain why outliers were more common in summer months in which water bodies should not be so frequently covered by clouds or fog [45]. However, if those cells were indeed covered by clouds but not detected [47], it is expected that they are filtered as outliers given the much lower temperature of clouds. Outlier detection remains a challenge and will be the focus of future improvements.

In the standard MOD/MYD11 LST and emissivity products, emissivities are assigned a priori based on land cover classification maps, i.e., the MCD12Q1 products. The latter product however was only generated until 2013. Therefore, some differences can be expected in places where LULC have changed to a drastically different cover in terms of emissivity. In the estimation of emissivity for MOD/MYD11 products, some other adjustments are made for the thermal infrared bidirectional reflectance distribution function (TIR BRDF), snow (from MOD10_L2 product), and green vs. senescent vegetation. This a priori approach is known to perform well for surfaces where emissivity can be correctly assigned based on the classification, i.e., it is best suited for land-cover types such as dense evergreen canopies, lake surfaces, snow, and most soils, all of which have stable emissivities error levels known to be within 0.01. However, it is significantly less reliable over arid and semi-arid regions [48]. Furthermore, errors in LST estimation may be introduced when the emissivity changes from day to night observation (e.g., due to condensation or dew), sudden natural changes (e.g., due to precipitation, wildfires), and from undetected nighttime cloud.

A LST time series reconstructed from remotely sensed records such as MODIS LST products, including the one we are presenting here, provides more spatial detail than gridded air temperature time series interpolated from ground station data (point observations), not only for central Europe with a relatively dense coverage, but more importantly also for other areas of the world with a lower density of ground stations. Commonly used climatic datasets with a comparable spatial resolution are WorldClim (http://worldclim.org/, [49]), and CHELSA (http://chelsa-climate.org/, [50]). These datasets, however, provide long-term statistics for past time periods such as the average for the years 1970–2000 for average monthly temperatures and not up-to-date time series with, e.g., monthly average temperature for the year 2016, not to mention daily temperature time series. Therefore, while highly used, most of the times they do not match the time frame of the phenomenon under study (see for example [51]). The long-term summaries we provide might be used to estimate the same temperature-based bioclimatic variables as those provided by WorldClim or CHELSA datasets.

We have demonstrated that our reconstructed LST time series is capable of identifying extreme events such as the well known heat wave that affected central Europe in 2003. LST might also be used to detect/monitor (surface) Urban Heat Islands (UHI) [52] either by estimating the difference in LST between urban areas and their rural surroundings (see Tomlinson et al. [53] and references therein) or by estimating the number of consecutive tropical nights from night-time LST [40]. Further use cases might include for example, the field of agriculture. Here, growing degree day (GDD) maps are relevant for the phenological assessment of budburst, flowering and crop maturing. GDD maps derived from reconstructed MODIS LST maps can account for local variations, i.e., they are likely to provide more detail. In addition, seasonal gradient maps (intra-annual short term trends) are of interest for agriculture. Furthermore, using MODIS LST time series, the shift of the upper forest limits may be assessed, likewise the suitability for the presence of tree (or other plants) species in terms of sensitivity to temperature. In the vast areas of permafrost, on the other hand, where almost no meteorological stations are available, MODIS LST data may be used to determine the duration of phases over 0 °C and their spatial extent [54]. Moreover, for applications that require higher spatial detail, the method that we present allows enhancing the spatial resolution and detail of the reconstructed LST with the use of a higher resolution elevation model as covariate in the TPS interpolation. However, this enhancement of the spatial resolution will affect only areas with diverse topography [5,19].

The new LST reconstruction method presented here improves over our previous LST reconstruction [5] in various ways. It requires less external data: previously climatic data and elevation, now only elevation. The temporal interpolation has been improved to better represent short-term fluctuations: previously weighted average, now local weighted regression (LWR) with order 2. We preferred LWR [55] over HANTS [56] or Fourier transforms (as used in [23] and others) because HANTS and Fourier transforms are only suitable to detect seasonal and long-term fluctuations with a regular pattern. LWR in turn is able to represent and preserve also irregular short-term fluctuations, particularly when used with polynomial order 2 or 3. Furthermore, the spatial interpolation has been improved by using TPS [33,34,35] instead of bspline interpolation with Tykhonov regularization [57]: TPS interpolation preserves existing values which can be optionally smoothed, whereas bspline always smoothes existing values. Moreover, the integration of covariates into TPS interpolation can further improve the results. Regarding spatial interpolation, the most important change is that no longer a fixed interpolation window is used, but instead an adaptive window, accounting for gaps of different size in the spatial domain. Using an adaptive window instead of a fixed window helps to reduce artifacts.

The entire workflow presented here was carried out with the open source geospatial suites GDAL and GRASS GIS along with newly published GRASS GIS add-on modules [29,30]. Therefore, the procedure can be replicated for other areas of interest and the reconstructed LST time series can be continuously updated with newly arriving data. The new LST reconstruction method scales well both locally (sub-continental) with a spatial resolution of nominally 1 km, and globally with a spatial resolution of 3 arc-min. Importantly, the reconstructed time series can be kept up-to-date with commodity hardware (processing one month takes about one day with an 7th generation Intel Core i5). The workflow could also be used to reconstruct other LST time series, for example the AVHRR LST time series from NOAA satellites or the Sentinel-3 SLSTR, which would extend the temporal coverage of the MODIS LST time series.

5. Conclusions

In this paper, we present a novel method to gap-fill LST time series at sub-continental and global scale. In total, more than 20,000 maps were reprocessed to obtain a gap-free daily time series of LST as well as aggregated products at 1 km and 3 arc-min resolution. We have shown the potential of our method and resulting LST time series spanning so far 14 years for the identification of extreme events. Moreover, such a long time series of daily observations is relevant for numerous applications in environmental monitoring, agriculture, urban planning and public health.

Finally, the use of free data and of open source suites guarantees the reproducibility of the proposed method along with the possibility to apply it to other areas and future LST time series. Global monthly reconstructed minimum, average, and maximum LST data (2003–2016) are available from https://www.mundialis.de/lst/ and Zenodo (http://0-doi-org.brum.beds.ac.uk/10.5281/zenodo.1115666) in GeoTIFF format.