This study aims to quantify the impact of multiscale morphological features on urban heat heterogeneity by comparing linear (Multiscale Geographically Weighted Regression -MGWR-) and non-linear (eXtreme Gradient Boosting -XGBoost-) statistical models to decode global-local relationships and develop cooling strategies. The research followed a five-stage methodology: (1) hexagon-based sample selection, (2) urban heat extraction, (3) multiscale morphological measurements, (4) application of MGWR and SHAP-explained XGBoost, and (5) sensitivity analyses. Both models computed the model accuracy by employing three parameters, the coefficient of determination, the adjusted coefficient of determination, and the root mean square error, enabling the comparison of the linear-based and non-linear-based approaches. This framework addressed knowledge gaps related to model-specificity, statistical model performances, building granularity and urban network parameters impact, spatial heterogeneity in linear models, and interpretability challenges in machine learning outputs. XGBoost outperformed at modelling discrete spatial heterogeneity due to its tree-based algorithm, while MGWR effectively modelled continuous spatial heterogeneity. Both models consistently identified ground area ratio (GAR) and number of plots (UBP) as major contributors to urban heat. GAR and UBP showed a strong non-linear influence on urban heat. The non-linearity extracted by XGBoost initially showed an upward trend in temperatures, but a diminishing return at a higher value of GAR and UBP. Although building features showed a low individual impact, their negative correlation with urban heat suggested a cumulative cooling potential. The research underscores model-specific non-stationary relationships between urban morpho-space and localised urban heat, promoting a tailored examination of mitigation potentials over a 'one-size-fits-all' approach, guiding urban planners to optimise built-fabric for effective heat mitigation and improved urban thermal quality.