Background on the Temporal Frequency of Surface Temperature Data

Landsat series between the 1980s and 2020 include Landsat 4 and 5 Thematic Mapper, Landsat 7 Enhanced Thematic Mapper Plus, and Landsat 8 Thermal Infrared Sensor. Each of these Landsat satellites has a 16-day revisit period, and the collected data also contain gaps because of clouds and shadows. The long-time record from Landsat U.S. ARD indicated complex patterns of Landsat record availability in space and time because of the geographic variations of Landsat overpass coverage, variable cloud and shadow obscuration at the time of overpass, past Landsat sensor acquisition and data quality issues, and the changing relative orientation of the fixed ARD tiles with respect to the Landsat orbit paths (Yan and Roy, 2020). All records from Landsat Collection 1 U.S. ARD ST in June 2018 in the eastern part of the Atlanta, Georgia, area (https://earthexplorer.usgs.gov/) are shown in figure 1A–F. These records were from Landsat 7 and 8 in different paths/rows. Atlanta is in the lower southeastern part of the tile. Landsat ST images on August 7 and 16 from Landsat 7 are shown in figure 1A and B, Landsat ST images on August 8, 15, 24, and 31 from Landsat 8 are shown in figure 1C–F, respectively. Only one cloud-free image covers the Atlanta area from Landsat 7 with the scan line corrector-off artifact. As a result, the image contains many data gaps. If all clear observations are used, they only represent Landsat ST on a specific day for the limited areas. The Landsat ST are not able to represent monthly means or cover entire urban areas. These gaps within time series can obscure parts of the landscape daily or seasonally. For many monitoring activities, such as UHI, these daily or seasonal data are critical for evaluating environmental response to climate variation.

Seasonal and Annual Variations of Land Surface Temperature Data

LST is highly related to incoming solar radiation and surface physical conditions. The high LST usually coincides with areas of high air temperature because higher LST can lead to higher air temperature through surface heat fluxes (Good and others, 2017; Xian and others, 2022). Air observational networks are, however, sparse and often nonrepresentative in complex anthropogenic landscapes. In the UHI study, for example, temperature differences between urban and rural areas are required; however, not every city has pairs of rural and urban weather stations, and it is difficult to rely on in situ data alone to obtain information about the UHI. Therefore, remotely sensed LST data are increasingly used as a convenient and accessible means to characterize UHIs or surface UHIs. However, the data gaps within Landsat ARD reduce the utility of these data sources to represent close to true seasonal or annual landscape conditions for modeling and monitoring environmental phenomena including surface UHIs. Magnitudes of Landsat ST and MODIS LST of a cloud-free pixel in 2016 from the same weather station site are shown in figure 2. Daily mean temperature data collected by the weather station (National Centers for Environmental Information, 2022) on the same day that Landsat cloud-free ST data were available also are shown in figure 2. Air temperatures that are observed above the ground are relatively larger than remotely sensed LST during cold months when insolation is low and atmospheric fluxes are large. During warm months, remotely sensed LSTs are relatively larger than surface air temperatures. Comparing with MODIS LST and climate observations, Landsat ST records, including Landsat 7 and 8, do not have enough clear records in every month for monthly means estimated from these records to represent the true thermal condition and its variation. The annual mean averaged from monthly means, therefore, will not be close to the true condition.

Gap-Filling Approach

The gap-filling method is based on a data-driven spatiotemporal gap-filling approach, which fills the missing pixel value using good pixels from the Worldwide Reference System-2 path overlap region. The satellite orbits are offset to allow 8-day repeat coverage of any Landsat scene area for one Landsat sensor. If more than two Landsat satellites are in operation, the revisit frequency will be larger than a single satellite. The gap-filling approach is based on the assumption that the orbit overlap region of the satellite’s remotely sensed imagery often maintains higher data density than the surrounding nonoverlap area (Zhou and others, 2020). Usually, a Landsat ARD tile contains two to three overlap regions. Many of these data in the overlap region can have similar spectral features with pixels that have the same land cover type within an ARD tile. Pixels in an ARD tile are usually from multiple Landsat paths/rows on different dates. For a given pixel time series in an ARD tile, good quality pixels can be modeled as a linear combination of a set of pixel time series from the orbit overlap region. Therefore, we can model the time-series seasonality for any given pixel using training pixels in the same month or season as the overlap region in the same tile.

Generally, the data gap is defined as absent pixels because of the swath coverage or cloud-contaminated observations determined from the Landsat pixel Quality Assessment band. The gap-filling approach used in this study contains the following steps: (1) training data selection, (2) training data refinement, (3) training sample selection, (4) linear regression model creation, and (5) time-series prediction.

Step 1. Potential training data selection.—This step identifies if data are missing in time-series images in an ARD tile. Pixels in the orbit overlap regions are then identified and extracted using the thresholds based on peak values of histograms of clear observations in overlap and nonoverlap regions.

Step 2. Reprocessing data in overlap regions.—Because time-series images in overlap regions may contain cloud-contaminated pixels, these data need to be processed to exclude these contaminated pixels before they can serve as training samples. A Harmonic Adaptive Penalty Operator (HAPO) algorithm (Zhou and others, 2022) is implemented to replace these contaminated pixels. The HAPO method first uses a seven-parameter linear harmonic model (eq. 1) to replace these pixels.

After fitting, HAPO adjusts parameters in the harmonic regression listed in equation 1 by removing the linear trend (c₀+c₁t) first. Then, HAPO normalizes the time series to ensure that the algorithm is robust to time series with different magnitudes. Two measures, the sum of squared residuals and a penalty function, are minimized in the calculation of regression model coefficients. The sum of squared residuals, which is defined in equation 2, minimizes the total difference between observations and model predictions. The penalty function, which is defined in equation 3, estimates the variance of the model curve.

Step 3. Clustering the training data.—After the potential training pixels are produced in step 2, the training data sources are clustered based on the Euclidean distance to meet time-series similarity so that the training samples can be stratified based on their similarity to the target pixel. The Kmeans analysis method is used to calculate the averages of time series in each cluster center, and the training data are clustered into 20 groups based on these means.

Step 4. Training data selection.—The training data are randomly stratified from multiple clusters so that some mixed pixel time series can have combined seasonality from different landscapes. The Euclidean distances from cluster centers to the target pixel are calculated and the inverse distances are used as the weights to randomly sample a total of 100 pixels from the training data. The stratified sampling strategy emphasizes training data with similar seasonality but is not restricted by clusters; therefore, the approach is robust when the target pixel has a similar distance to multiple cluster centers, which could be the case for mixed pixels.

Step 5. The full time-series prediction.—After training samples are selected, the samples are used to create linear regression models to predict time series of any target pixel. The least absolute shrinkage and selection operator regression (eq. 6) is used to build regression models to predicted values at gaps based on clear observations from the training pixels.

The workflow of the gap-filling procedures is shown in figure 3. The steps highlighted in the box with the dashed line are repeated 10 times for each target pixel. The median and standard deviation (STD) values from each time step are produced in the final prediction. The median value is used for gap filling, and the STD value is used as an indicator of uncertainty of the model prediction. Generally, the gap-filling procedures are only applied to these containment pixels and no-data pixels. Each ARD tile has a minimum of four Landsat ST tiles every month if records from only one Landsat sensor are available. In most of these ARD tiles, only a part of a tile has the Landsat ST record, as shown in figure 1. After gap filling, each ARD tile will be filled; therefore, at least four records can be used to calculate monthly means in every month. The maximum magnitude of a pixel in a month is used to represent the maximum Landsat ST for the month. The seasonal mean Landsat ST values are calculated from three monthly means, and seasonal maximum is averaged from the maximum values of 3 months.

Data Validation

After the gap-filling algorithm is used to produce filled Landsat ST magnitudes in an ARD tile, other ST records are used to evaluate filled Landsat ST data. Here, we use the Atlanta area gap-filling result to illustrate the validation result of filled data. The independent reference datasets for validation are from three major sources: the MODIS LST collection, surface air temperature records observed by the same weather station used in figure 2 (National Centers for Environmental Information, 2022), and gridded climate data from Daymet (Thornton and others, 2021). The MODIS record is used because it collected imagery consistently every 8 days and the data can be used as a reference to compare with the temporal frequency of filling data. Also, climate data are selected in the days that Landsat ST records are available. About 300 sampling pixels are randomly selected in the study tile. The climate records are used to evaluate temporal variation patterns between climate data and Landsat ST data. The comparison between gridded climate data and remotely sensed LST can generate several statistical parameters, including root mean square error (RMSE) and coefficient of determination (r²) between the reference and filled pixel values, to determine the connections between these data.

Gap-Filled Landsat Surface Temperature Data

The U.S. ARD tile (h24v13, where h denotes horizontal and v denotes vertical) that covers the southeastern side of the Atlanta area was chosen as the study tile. Most of Atlanta is included within the tile. The gap-filled procedures first used all available images in 1 month to produce an image to represent observed Landsat ST for that month. After gap filling, all derived clear pixels from different dates were averaged to generate the monthly Landsat ST imagery. These monthly Landsat ST data were averaged to produce seasonal and annual Landsat ST means to represent observed Landsat ST means in a specific season and year. The seasonal mean in this study is defined as follows: January–March is season 1 (SN1), April–June is season 2 (SN2), July–September is season 3 (SN3), and October–December is season 4 (SN4). The gap-filled Landsat ST in the same month and year as shown in figure 1 and the related uncertainty map are shown in figure 4A and B. The gap-filled algorithm filled every pixel that is either contaminated by cloud/shadow or has no data with extrapolated magnitudes in a tile from regression models. After gap filling, the tile was filled with estimated and original clear pixels from six dates; therefore, the monthly mean was calculated from the Landsat STs of these six dates to represent the monthly status of Landsat ST. The uncertainty map indicates that uncertainties are reduced for most of the urban area pixels. Higher uncertainties are observed in some nonurban pixels in the eastern and northwestern parts of the Atlanta metropolitan area. These areas have relatively large cloud cover or are without records (fig. 1). Most of the data used for gap filling are from the central part of the tile.

The monthly means of gap-filled Landsat ST were further averaged to have seasonal means. The Landsat ST seasonal means from gap-filled images were compared with Landsat ST seasonal means estimated directly from available clear pixels without gap filling. Landsat ST seasonal means in 2016 without gap filling are shown in figure 5A–D. The images cover most of the Atlanta urban area. Landsat ST in SN1 shows relatively large magnitudes that are unreasonably high in Landsat ST (fig. 5A). Landsat STs in SN2 and SN4 contain many pixels that are potentially affected by clouds/shadows, resulting in relatively small magnitudes of Landsat ST in most urban areas (fig. 5B, D). The Landsat ST seasonal means indicate urban warming in most urban areas, but they also consist of unreasonably low Landsat STs in some parts of the area. A clear seamline also is observed in SN1, SN2, and SN4. SN3 shows relatively large Landsat ST magnitudes in most urban areas (fig. 5C). The gap-filled Landsat STs fixed unreasonably high seasonal magnitudes in SN1 (fig. 6A) and SN4 (fig. 6D) and low seasonal magnitudes in SN2 (fig. 6B). The gap-filled Landsat ST in SN3 indicates that most pixels in the urban area have larger magnitudes than many surrounding nonurban areas (fig. 6C). The Landsat ST seasonal means do not contain pixels having relatively small magnitudes in the urban area, resulting in clear UHI spatial distribution patterns in all four seasons.

The uncertainty maps of gap-filled pixels in these four seasons are shown in figure 7A–D for SN1–SN4. Generally, the areas that have more data have fewer uncertainties and STD values are relatively small. The strips in the middle of the tile indicate relatively low uncertainties because the area represents the overlap area of multiple Landsat paths/rows. Several areas in the northwestern corner of the tile have relatively high uncertainty pixels. The uncertainties are relatively low in most of the urban areas.

Comparisons of Gap-Filled Data with Raw Landsat Surface Temperature Data

The annual variation of single pixel gap-filled Landsat ST at the same location used in figure 2 was compared with nongap-filled Landsat ST, MODIS LST, and climate records in 2016 (fig. 8). The gap-filled Landsat ST added more records in warm months. These increasing records improved the temporal density of Landsat ST, resulting in better estimates of monthly and seasonal means in these warm seasons. In cold months, the gap-filled Landsat STs are slightly elevated by including more cloud/shadow-free pixels that are warmer than cloud/shadow-contaminated pixels. The gap-filled Landsat STs have a temporal variation pattern that is similar to the pattern from MODIS and weather station records.

Validation of Gap-Filled Landsat Surface Temperature Data Using Moderate Resolution Imaging Spectroradiometer Land Surface Temperature and Gridded Climate Data

We used randomly stratified samples over urban and agricultural lands to compare gap-filled Landsat ST and Daymet daily maximum air temperature (number of samples [n] =86), gap-filled Landsat ST and MODIS LST (n=44), and MODIS LST and Daymet maximum air temperature records (n=44). Statistically correlated outcomes including correlation of determination (r²), standard error (SE), and RMSE are provided in table 1. The gap-filling Landsat ST samples over urban land have relatively higher correlations and lower SE and RMSE values than Daymet records and MODIS data over agricultural land. However, MODIS LSTs have similar correlations and SEs with air temperature but higher RMSE values in the urban area. Scatter plots of gap-filling Landsat ST versus air temperature (fig. 9A) and gap-filling Landsat ST versus MODIS LST (fig. 9B) in 2016 from randomly selected samples are shown in figure 9. The gap-filling Landsat ST and air temperature are highly correlated and have an r² value of 0.88, but Landsat ST values are greater than air temperature. The MODIS LST also has a relatively high and significant correlation at a less than 1-percent level (probability [or P] less than 0.01) with Landsat ST and has an r² value of 0.92.