This paper explores the potential of downscaling Land Surface Temperature, LST, based on land features multi-interaction with a spatial regression multi-modelling. The Radiative Transfer Equation first helped to create an LST15 m layer over Landsat-OLI/TIRS. Next, a bilinear assessment of LST is conducted over elevation and hillshade, so to adjust shadow/brightness. Then, interactions are modelled on a feature-to-feature linear basis between spectral indices, SI’s, representing vegetation, built-up, soil, water and shadow. A multilinear regression model is further built between combined pairs of interactions and LST15 m. The first principal Component, PC1, of all subtractions of each pair of interactions from others is stacked with individual SI’s, to build another multi-regression model around LST15 m. Each of the three models is individually subtracted from LST15 m, normalized, [0–1], and their sum serves as the residuals layer. The downscaling step uses coefficients of the interactions model with PC1 over the corresponding Sentinel2-MSI 10 m SI’s, and adds back the gaussian-kernel of residuals. The Normalized Urban-High Spatial Resolution-Land Surface Temperature, NU-HSR-LST10 m, is the final product, that sharpens hot/cold spots, with a highly spread of values among land features. As supporting results, directions of relations with vegetation and built-up were improved, while unexpected relations were alternatively revealed (water) or reversed (soil, shadow); determination coefficients, R2, shows a strong correlation of NU-HSR-LST10 m to LST30 m (R2:[0.7304–0.9844]), even stronger with a closest model (R2:[0.85–0.99]); a variance analysis between NU-HSR-LST10 m and LST30 m is quasi-insignificant between [0.0002–0.00297]; and a root mean square error computed in a war-disturbed urban context, was lower for NU-HSR-LST10 m, [0.057–0.096], than for LST30 m,[0.106–0.151], as stability in dynamics depiction. Finally, the machine learning algorithm of random forest based on different seeds achieved overall accuracy between [0.92–1]. From these results, the downscaling process is efficient in better distinguishing contributions per land feature in diverse urban environments, while more cross-validation based on meteorological stations is still needed.