Satellite-based remote sensing technologies are utilized extensively to investigate urban thermal environments under rapid urban expansion. Current MODIS data is, however, unable to adequately represent the spatially detailed information because of its relatively coarser spatial resolution, while Landsat data can’t explore temporally the refined analysis due to the low temporal resolution. In order to resolve this situation, we used MODIS and Landsat data to generate “Landsat-like” data by using the flexible spatiotemporal data fusion method (FSDAF), and then studied spatiotemporal variation of land surface temperature (LST) and its driving factors. The results showed that 1) The estimated “Landsat-like” data have high precision; 2) By comparing 2013 and 2016 datasets, LST increases ranging from 1.8°C to 4°C were measurable in areas where the impervious surface area (ISA) increased, while LST decreases ranging from -3.52°C to -0.70°C were detected in areas where ISA decreased; 3) LST has a strongly negative relationship with the Normalized Difference Vegetation Index (NDVI), and a strongly positive relationship with Normalized Difference Built Index (NDBI) in summer; and 4) LST is well correlated with Building density (BD), in a complex conic mode, and LST may increase by 0.460°C to 0.786°C when BD increases by 0.1. Our findings can provide information useful for mitigating undesirable thermal conditions and for long-term urban thermal environmental management.