Abstract The type of land use and land cover plays a decisive role in land surface temperature (LST). As cities are composed of varied covers, including vegetation, built-up areas, buildings, roads and areas without vegetation, understanding LST patterns in complex urban spaces is becoming increasingly important. The present study investigated the relationship between LST and albedo, NDVI, NDWI, NDBI and NDBaI in the period between 1994 and 2017. Images from Thematic Mapper (TM) and Thermal Infrared Sensor (TIRS) onboard the Landsat 5 and 8 satellites, respectively, were used in the study. The images were processed, resampled (spatial resolution of 120 m) in the environment of the QGIS 3.0 software and, finally, centroids were extracted resulting in a total of 1252 points. A classical regression (CR) model was applied to the variables, followed by spatial autoregressive (SARM) and spatial error (SEM) models, and the results were compared using accuracy indices. The results showed that the highest correlation coefficient was found between albedo and NDBaI (r = 0.88). The relationship between albedo and LST (r = 0.7) was also positive and significant at р < 0.05. The global Moran's I index showed spatial dependence and non-stationarity of the LST (I = 0.44). The SEM presented the best accuracy metrics (AIC = 3307.15 and R2 = 0.65) for the metropolitan region of Belém, explaining considerably more variations in the relationship between explanatory factors and LST when compared to conventional CR models.