High spatiotemporal land surface temperature (LST) datasets are increasingly needed in a variety of fields such as ecology, hydrology, meteorology, epidemiology, and energy systems. Moderate Resolution Imaging Spectroradiometer (MODIS) daily LST is one of such high spatiotemporal datasets that are widely used. But, it has a large amount of missing values primarily because of clouds, shadows, and other atmospheric conditions. Gapfilling the missing values is an important approach to create seamless high spatiotemporal LST datasets. However, current gapfilling methods have limitations in terms of accuracy and efficiency to assemble the data over large areas (e.g., national and continental levels). In this study, we developed a 3-step hybrid method by integrating daily merging (gapfilling missing values at one overpass using values at the other three overpasses each day), spatiotemporal gapfilling (estimating missing values based on values of their spatial and temporal neighbors), and temporal interpolation (gapfilling missing values based on values of their neighboring days), to create a seamless high spatiotemporal LST dataset using the four daily LST observations from the two MODIS instruments on Terra and Aqua satellites. We applied this method in urban and surrounding areas in the conterminous U.S. in 2010. The evaluation of the gapfilled LST product indicates its root mean squared error (RMSE) to be 3.3 K for mid-daytime (1:30 pm) and 2.7 K for mid-nighttime (1:30 am) observations. The method can be easily extended to other years and regions and is also applicable to other satellite products for large areas. This seamless daily (mid-daytime and mid-nighttime) LST product with 1 km spatial resolution is of great value for studying urban climate (e.g., quantifying surface urban heat island intensity, creating seamless high spatiotemporal air temperature dataset) and the related impacts on people (e.g., health and mortality), ecosystems (e.g., phenology), and energy systems (e.g., building energy use).