The recently released Landsat Analysis Ready Data (ARD) over the United States provides the opportunity to investigate landscape dynamics using dense time series observations at 30-m resolution. However, the dataset often contains data gaps (or missing data) because of cloud contamination or data acquisition strategy. We present a new algorithm that focuses on data gap filling using clear observations from orbit overlap regions. Multiple linear regression models were established for each pixel time series to estimate stable predictions and uncertainties. The model's training data came from stratified random samples based on the time series similarity between the pixel and data from the overlap regions. The algorithm was evaluated using four tiles (5,000 × 5,000 30-m pixels for each tile) from 2018 land surface temperature data (LST) in Atlanta, Georgia. The accuracy was assessed using 1,000 randomly masked pixels and daily air temperature from eight ground stations. Both assessments showed the r2 value above 0.75, except two stations with mixed Landsat pixels. We also compared our results with the eMODIS LST product in terms of annual mean temperature. The two maps showed a similar spatial pattern at the region level, but our results showed more spatial detail in the urban area that matched the pattern of impervious surface. We also applied the method on ARD surface reflectance bands at Fairbanks, Alaska, to illustrate its improvements in surface reflectance products and in land change modeling. This approach can also be applied to other datasets, vegetation indexes, or spectral reflectance bands of other sensors.