Spatial Modelling of Land Surface Temperature Responses to Vegetation and Built up Intensity Across Oyo State, Nigeria: A Multi-Source Remote Sensing and Geospatial Statistical Approach

Full text

Text extracted from the article PDF for reading and search. Tables, figures and equations are in the PDF, which remains the authoritative version.

1. Introduction

Land surface temperature (LST) is an important indicator of the thermal condition of the Earth's surface and provides useful information for understanding spatial variations in surface heating. Satellite thermal infrared remote sensing has become an important approach for observing LST because it provides spatially continuous information over large areas and enables relationships between surface temperature and landscape characteristics to be examined [13]. LST has consequently been applied extensively in studies of urban thermal environments, vegetation conditions and land-use/land-cover (LULC) characteristics.

Different land-cover types have different physical and biophysical properties and therefore exhibit different thermal responses. Vegetated surfaces can influence surface temperature through processes such as evapotranspiration, whereas built-up and exposed surfaces can exhibit greater surface heating. Spectral indices provide useful measures for examining these relationships. The Normalized Difference Vegetation Index (NDVI) is commonly used to characterize vegetation, while the Normalized Difference Built-up Index (NDBI) is used to represent built-up characteristics. In Abuja, Nigeria, Adeyeri et al. found negative relationships between LST and vegetation-related indices and a strong positive relationship between LST and NDBI, demonstrating the importance of surface-cover characteristics in explaining thermal variation. [1] Research in southwestern Nigeria has similarly demonstrated the relationship between land-cover change and surface temperature. Salam et al. examined urban expansion, LULC dynamics and LST in Oyo town using multi-date Landsat imagery and reported increasing built-up areas alongside changes in vegetation and increasing surface temperature [12]. Although that study provides useful evidence for the Oyo environment, its focus was on urban expansion and temporal LULC dynamics within Oyo town rather than the spatially varying relationships between LST, vegetation intensity and built-up intensity across the entire state. A further consideration in LST analysis is that relationships between environmental variables may not be spatially uniform. Global regression models such as ordinary least squares (OLS) provide a single coefficient for the study area and may therefore conceal local differences in the strength or direction of relationships. Spatial dependence can also occur when nearby observations have similar thermal characteristics. Geographically weighted regression (GWR) provides a means of addressing spatial non-stationarity by estimating local relationships across geographic space. Eshetie, for example, compared OLS and GWR in Addis Ababa and found that GWR provided a substantially better representation of the spatially varying relationship between LST and explanatory variables [4].

Despite the increasing use of remote sensing for LST assessment, there is still a need for integrated analyses that combine LULC-based temperature comparison with continuous measures of vegetation and built-up intensity while explicitly examining spatial dependence and spatial heterogeneity. This is particularly relevant in Oyo State, where built-up areas, agricultural land, vegetation, bare surfaces and water bodies occur across a large and heterogeneous landscape. A state-wide assessment can therefore provide a broader understanding of surface thermal variation than analyses restricted to individual urban centres. This study investigates the spatial relationships between LST, vegetation intensity and built-up intensity across Oyo State, Nigeria, using multi-source remote sensing data and geospatial statistical techniques. Specifically, the study (i) classifies and assesses the accuracy of major LULC classes, (ii) examines differences in LST among the identified LULC classes, (iii) evaluates the relationships between LST, NDVI and NDBI using global statistical techniques, (iv) assesses spatial autocorrelation and residual patterns, and (v) applies GWR to determine how the relationships between LST and NDVI/NDBI vary across Oyo State.

2. Materials and Methods

2.1. Study Area

Oyo State is located in southwestern Nigeria and covers approximately 27,514 km². The state has a heterogeneous landscape consisting of built-up areas, cropland, dense vegetation, bare land and water bodies. The diversity of these land-cover types provides an appropriate setting for examining variations in land surface temperature and their relationships with vegetation and built-up intensity.

2.2. Data Acquisition and Pre-processing

Landsat 8 and Landsat 9 Collection 2 Level-2 imagery were used as the principal satellite data source. Surface reflectance bands were used for LULC classification and spectral-index derivation, while the surface-temperature product was used for LST estimation. January 2025 was selected for the thermal analysis, with 22 suitable scenes acquired between 1 and 26 January 2025, comprising 12 Landsat 8 and 10 Landsat 9 scenes [9,11].

For LULC classification, SR_B2–SR_B7, NDVI, NDBI and NDWI were combined to form a nine-feature dataset. The imagery was subjected to quality screening to reduce the effects of clouds, shadows and radiometric contamination before classification and analysis.

2.3. LULC Classification and Accuracy Assessment

Five LULC classes were considered: built-up, dense vegetation, cropland, bare land and water. A Random Trees classifier was applied using 100 trees, a maximum tree depth of 30 and up to 1,000 samples per class [2]. Classification accuracy was evaluated using 310 independent validation points, separate from the training samples.

2.4. Derivation of LST, NDVI and NDBI

LST was derived from the Landsat surface-temperature product and converted from Kelvin to degrees Celsius. NDVI was derived from the red and near-infrared bands to represent vegetation intensity, while NDBI was calculated from the near-infrared and shortwave-infrared bands to represent built-up intensity. The resulting datasets were clipped to the study area and prepared for statistical analysis.

2.5. Sampling and Statistical Analysis

A total of 750 equalized stratified-random points, consisting of 150 points from each LULC class, were used for statistical analysis. LST, NDVI and NDBI values were extracted at these locations.

The Kruskal–Wallis test was used to assess differences in LST among LULC classes because the data did not satisfy the assumptions required for conventional parametric comparison [8]. Significant differences were further examined using Dunn's post-hoc test with Holm correction [3,7]. Pearson and Spearman correlation coefficients were calculated to evaluate relationships between LST, NDVI and NDBI.

2.6. Global and Local Spatial Modelling

An ordinary least squares (OLS) model was used to estimate the global relationship between LST and the two explanatory variables, NDVI and NDBI. Global Moran's I was subsequently applied to evaluate spatial autocorrelation in both LST and OLS residuals [10]. Getis–Ord Gi* analysis was used to identify spatial clusters of positive and negative regression residuals [6] Geographically weighted regression (GWR) was then applied to examine whether the relationships between LST, NDVI and NDBI varied across space [5]. The model used a Gaussian distribution, bisquare weighting and 46 neighbouring observations. The neighbourhood size was selected based on the lowest AICc obtained during the search procedure. The analyses were treated as observational, and relationships between LST, NDVI and NDBI were interpreted as associations rather than definitive causal effects.

3. Results

3.1. Land Use/Land Cover Classification and Accuracy

The 2025 land use/land cover classification of Oyo State produced five classes: built-up, dense vegetation, cropland, bare land and water. The classification achieved an overall accuracy of 93.55% with a Kappa coefficient of 0.9067, indicating a high level of agreement between the classified image and the independent validation data. The producer's accuracy was highest for water (100%), followed by cropland (95.86%), bare land (93.85%), built-up (91.89%) and dense vegetation (86.79%). User's accuracy ranged from 89.47% for built-up to 100% for water. These results indicate that the classified LULC dataset was sufficiently reliable for subsequent analysis of LST differences across land-cover classes.

3.2. Spatial Distribution of Land Surface Temperature

The January 2025 LST analysis produced values ranging from 19.64°C to 49.56°C, with a mean temperature of 36.53°C and a standard deviation of 3.02°C. The resulting LST surface showed substantial spatial variation across Oyo State.

Higher LST values were generally associated with locations characterized by built-up and bare surfaces, whereas lower temperatures occurred mainly over water bodies and areas of dense vegetation. This spatial pattern provides an initial indication of the thermal differences associated with land-cover conditions across the state.

Fig 3: Land Surface Temperature of Oyo State, Nigeria, 2025

3.3. Land Surface Temperature Variation among LULC

Classes The equalized sample of 750 observations revealed clear differences in mean LST among the five LULC classes. Bare land recorded the highest mean temperature (37.98°C), followed by built-up areas (37.73°C) and cropland (36.50°C). Dense vegetation had a substantially lower mean temperature (33.43°C), while water recorded the lowest mean temperature (27.45°C).

Table 2: Descriptive statistics of LST by LULC class LULC Class n Mean (°C) SD (°C) Minimum (°C) Maximum (°C) Built-up 150 37.73 2.31 31.05 42.80 Dense vegetation 150 33.43 2.44 28.84 41.09 Cropland 150 36.50 2.91 28.35 42.67 Bare land 150 37.98 2.28 32.43 43.14 Water 150 27.45 1.95 25.04 35.41 The Kruskal–Wallis test indicated a statistically significant difference in LST among the five LULC classes, H(4) = 473.80, p < 0.001. The estimated epsilon-squared effect size was approximately 0.631, indicating a very large effect of LULC grouping on the observed distribution of LST. Dunn's post-hoc comparisons with Holm correction showed significant differences for most class pairs. The only non-significant comparison was between built-up and bare land (Holm-adjusted p = 0.581), despite their slightly different mean temperatures. This indicates that these two surface types exhibited statistically comparable thermal conditions within the sampled observations.

3.4. Relationships between LST, NDVI and NDBI

The correlation analysis showed different relationships between LST and the two continuous explanatory variables. NDBI exhibited a strong positive relationship with LST, with Pearson's r = 0.862 (p < 0.001) and Spearman's ρ = 0.837 (p < 0.001). A simple regression of LST on NDBI produced an R² of 0.742, indicating a strong association between built-up intensity and LST.

The overall Pearson correlation between NDVI and LST was r = 0.491 (p < 0.001). However, this overall positive relationship should be interpreted cautiously because the observations were distributed across substantially different LULC classes. Within most individual classes, NDVI showed a negative association with LST, indicating that the global relationship was influenced by differences among land-cover contexts.

When NDVI and NDBI were considered simultaneously in an OLS regression, the model explained 77.32% of the variation in LST (R² = 0.7732; adjusted R² = 0.7726). Both explanatory variables were statistically significant. The estimated coefficients were 5.340 for NDVI and 31.477 for NDBI, with both variables having robust p-values below 0.001.

3.5. Spatial Autocorrelation and OLS Diagnostics

Global Moran's I for the sampled LST values was 0.7728, with a z-score of 44.75 and p < 0.001, indicating strong positive spatial autocorrelation. Thus, locations with similar LST values tended to occur near one another rather than being randomly distributed.

The OLS residuals also displayed significant spatial autocorrelation, with Moran's I of 0.5381, z = 31.20, and p < 0.001. This indicates that although the OLS model explained a substantial proportion of LST variation, considerable spatial structure remained unexplained.

The Koenker statistic was significant (p < 0.001), indicating spatial non-stationarity and/or heteroscedasticity in the relationships. The Jarque–Bera statistic was not significant (p = 0.260), indicating no significant departure from normality of the OLS residuals.

3.6. Spatial Clustering of OLS Residuals

The Getis–Ord Gi* analysis identified spatial clusters of positive and negative OLS residuals across Oyo State. Positive residual clusters indicate locations where observed LST was higher than predicted by the global NDVI–NDBI model, whereas negative residual clusters indicate locations where observed LST was lower than predicted. The spatial distribution of these residual clusters, particularly the concentration of positive residuals in parts of southern Oyo State and negative residuals in northeastern and eastern areas, further demonstrates that the relationships between LST, vegetation and built-up intensity vary spatially.

3.7. Geographically Weighted Regression

GWR substantially improved the model compared with the global OLS approach. Using a bisquare weighting scheme and 46 neighbouring observations, the GWR model produced an R² of 0.9424, an adjusted R² of 0.9272, and an AICc of 2566.73. In comparison, the OLS model produced an R² of 0.7732 and an AICc of 3314.10.

The reduction in AICc of approximately 747 points indicates a substantially better model fit when spatially varying relationships were incorporated. The improvement also demonstrates that the relationships between LST, NDVI and NDBI were not adequately represented by a single global regression equation.

The GWR residuals had a Moran's I of −0.0239, z = −1.31, and p = 0.191, indicating that the residuals were not significantly spatially autocorrelated. This represents a major improvement over the OLS residuals and suggests that GWR successfully accounted for much of the spatial structure remaining in the global model.

3.8. Spatial Variation in GWR Coefficients

The GWR coefficients demonstrated considerable spatial variation in the relationships between LST and the explanatory variables. The NDVI coefficient ranged from −67.52 to 35.66, with a mean of −7.72, while the NDBI coefficient ranged from −38.38 to 47.82, with a mean of 12.43.

The local R² values ranged from −0.205 to 0.990, with a mean of 0.752. The wide range of local coefficients confirms that the strength and direction of the relationships varied across the study area rather than remaining spatially constant. The local significance assessment further demonstrated this spatial heterogeneity. For NDVI, 288 locations (38.40%) showed significant negative associations, 155 locations (20.67%) showed significant positive associations, while 307 locations (40.93%) were not significant based on the approximate ±1.96 t-statistic screening criterion. For NDBI, 445 locations (59.33%) showed significant positive associations, 42 locations (5.60%) showed significant negative associations and 263 locations (35.07%) were not significant.

4. Discussion

4.1. LULC and Land Surface Temperature

The results demonstrate clear thermal differences among the five LULC classes. Bare land and built-up areas recorded the highest mean LST, while water and dense vegetation were considerably cooler. The similar temperatures observed for bare land and built-up areas suggest that surfaces with limited vegetation and moisture can exhibit comparable thermal responses.

The lower temperature associated with dense vegetation is consistent with the cooling role of vegetation through shading and evapotranspiration. Thermal remote sensing studies have similarly established strong links between land-cover characteristics and surface thermal conditions [13,14].

4.2. Influence of Vegetation and Built-up Intensity

NDBI showed a substantially stronger positive association with LST than the overall NDVI relationship. This indicates that increasing built-up intensity was generally associated with higher surface temperatures across the study area. This finding agrees with previous work in Nigeria showing relationships between urban development, surface temperature and vegetation characteristics [1,12]. The overall NDVI relationship requires more careful interpretation. Although the global correlation was positive, the GWR results showed that the NDVI coefficient varied in both magnitude and direction across Oyo State. This suggests that vegetation does not have a uniform relationship with LST throughout the state and that local environmental and land-cover conditions influence the observed relationship.

4.3. Spatial Dependence and Model Performance

The significant Moran's I for both LST and OLS residuals confirms that the observations were spatially dependent. Consequently, a single global regression model was insufficient to fully represent the spatial structure of LST. GWR produced a much higher explanatory power and substantially lower AICc than OLS, while its residuals showed no significant spatial autocorrelation. This demonstrates the importance of accounting for spatial non-stationarity when examining relationships between LST, vegetation and built-up intensity. Similar advantages of spatially varying regression approaches have been reported in recent LST studies [4].

4.4. Implications

The findings indicate that thermal conditions across Oyo State are strongly associated with land-cover characteristics, particularly built-up intensity, but these relationships vary geographically. Areas experiencing greater built-up intensity may therefore require greater attention in land-use planning and urban environmental management, while maintaining vegetation can contribute to moderating surface thermal conditions.

Overall, the results support the use of multi-source remote sensing and spatial modelling for understanding spatial variations in LST and provide a more location-specific interpretation than a conventional global regression model.

5. Conclusion

This study assessed the spatial relationships between land surface temperature, vegetation and built-up intensity across Oyo State, Nigeria, using Landsat 8/9 imagery, LULC classification, statistical analysis and spatial regression. The LULC classification achieved a high overall accuracy of 93.55% with a Kappa coefficient of 0.9067. Clear differences in LST were observed among the five LULC classes. Bare land and built-up areas recorded the highest mean temperatures, while water and dense vegetation recorded the lowest. The Kruskal–Wallis test confirmed that these differences were statistically significant. NDBI showed a strong positive association with LST, while the NDVI relationship varied across land-cover contexts. Significant spatial autocorrelation was observed in LST and OLS residuals, demonstrating that spatial dependence could not be ignored. GWR substantially improved model performance, increasing R² from 0.7732 to 0.9424, while its residuals showed no significant spatial autocorrelation. The findings demonstrate that the relationship between LST, vegetation and built-up intensity is spatially variable across Oyo State. Therefore, spatially explicit approaches such as GWR provide a more informative framework for understanding thermal patterns than a single global regression model. The results can support land-use planning and environmental management aimed at reducing excessive surface heating in rapidly developing areas.