Decadal Urban Heat Island Intensification in a Rapidly Growing Tropical Secondary City: Integrating Land Use Dynamics and Spatial Autocorrelation Analysis in Kendari, Indonesia

2.1 Study area

The research location is at coordinates 3°58.540' South Latitude – 122°31.036' East Longitude and 3°58.095' South Latitude – 122°31.923' East Longitude. Administratively, Kendari City is located in Southeast Sulawesi Province, covering eleven subdistricts with an area of 267.37 km² (Figure 1). All spatial analyses were performed using ArcGIS 10.8 and QGIS 3.22. This location was chosen because it is a coastal city with a tropical climate that is undergoing rapid urbanization accompanied by significant land cover change.

1.png

Figure 1. Research location

2.2 Data and methods

In this study, data collection was conducted through a combination of field measurements and satellite image analysis. Using a single satellite platform (Landsat 8 OLI/TIRS) for all observation years improves radiometric consistency and minimises methodological bias in inter-temporal LULC and LST comparison (Table 1). Landsat 8 Collection 2 Level-2 Surface Reflectance products were utilised in this study because they provide atmospherically corrected surface reflectance and improved radiometric calibration compared to the previous Collection 1 datasets. The use of Collection 2 Level-2 products enhances the reliability and inter-temporal consistency of LULC classification, NDVI calculation, emissivity estimation, and LST retrieval for multi-year UHI analysis.

All Landsat scenes were selected during the dry season (September–October) to minimize seasonal variability in vegetation moisture, precipitation, and atmospheric conditions. Images with cloud cover below 15% were prioritized to reduce cloud contamination and ensure inter-annual comparability of thermal conditions.

2.3 Land use classification

The land use classification was done using the supervised Maximum Likelihood Classification (MLC) technique on multi-temporal Landsat 8 OLI (Operational Land Imager) images for the years 2014, 2019, and 2024, employing the same methodology, band combinations, and spatial resolution to ensure methodological consistency across the years of observation. The study identified six land use categories: settlement (residential, commercial, industrial); forests (natural forests and plantation forests); gardens/agricultural land (mixed agricultural land); rice fields (paddy fields); water bodies (rivers, lakes, coastal areas); and vacant land (open land, barren land and construction sites).

It had an overall accuracy of 86.67% (2014), 87.33% (2019) and 91.67% (2024) and Kappa coefficients of 0.83, 0.84 and 0.89, respectively, using 300 validation points.

To evaluate changes in land use, three indicators were computed for each land cover class for each year: absolute area change (ha), relative percentage area change, and proportional land share change. The general trends in urban growth and landscape change for the study area between 2014 and 2024 were identified.

2.4 Land surface temperature retrieval

Land surface temperature (LST) is important in environmental monitoring and studying land surface-atmosphere interactions. It is useful in climate system modeling, studies on the UHI effect, and ecological modeling studies [29, 30]. Numerous techniques and methods for LST analysis exist. The retrieval of the LST requires the use of satellite remote sensing. The use of thermal infrared information enables the availability of the LST values from different land covers using algorithms such as the split-window method, the single-channel method, and other methods, taking into account the effect of atmospheric interference and emissivity [31-33]. Steps to derive LST from Landsat imagery include converting Landsat DN (Digital Number) to radiance, converting band radiance to BT (Brightness Temperature), and retrieving the Land Surface Emissivity (LSE). The LST retrieval process was performed using Landsat 8 Thermal Infrared Sensor (TIRS) Band 10. Only Band 10 was utilized because Band 11 is affected by larger calibration uncertainty associated with stray-light contamination, as documented in USGS technical recommendations and previous thermal remote sensing studies. The use of a single thermal band approach improves the reliability and consistency of inter-temporal LST estimation. The Digital Number (DN) values from Band 10 were converted into Top of Atmosphere (TOA) spectral radiance using the following equation:

$L \lambda=M L \times Q c a l+A L-O i$    （1）

where, Lλ = TOA spectral radiance (W·m⁻²·sr⁻¹·µm⁻¹); ML = band-specific multiplicative rescaling factor (RADIANCE_MULT_BAND_10 from image metadata); Qcal = quantized calibrated pixel value (DN); and AL = band-specific additive rescaling factor (RADIANCE_ADD_BAND_10 from image metadata).

TOA spectral radiance was converted to at-sensor brightness temperature (BT) in degrees Celsius using the thermal calibration constants provided in the image metadata:

$B T=\frac{K_2}{\ln \left(\frac{K_1}{L_\lambda}+1\right)}-273.15$      （2）

where, BT = at-sensor brightness temperature (℃); K₁ = band-specific thermal conversion constant (K1_CONSTANT_BAND_10 = 774.8853 W·m⁻²·sr⁻¹·µm⁻¹); and K₂ = band-specific thermal conversion constant (K2_CONSTANT_BAND_10 = 1321.0789 K).

LSE (ε) was estimated following the three-stage NDVI Threshold Method (NDVI_THM) [34].

First, the Normalized Difference Vegetation Index (NDVI) was computed from the atmospherically corrected surface reflectance values of OLI Band 5 (NIR) and Band 4 (Red):

$N D V I=\frac{N I R-R E D}{N I R+R E D}=\frac{B 5-B 4}{B 5+B 4}$       (3)

where, NDVI values range from −1 to +1, with values approaching +1 indicating dense photosynthetically active vegetation and values near 0 or negative indicating bare soil, impervious surfaces, or open water.

Second, the fractional vegetation cover (PV) was derived from NDVI using the Carlson and Ripley (1997) formulation:

$P V=\left(\frac{N D V I-N D V I_{\min}}{N D V I_{\max}-N D V I_{\min}}\right)^2$       (4)

where, NDVI_min and NDVI_max represent the minimum and maximum NDVI values observed within the study area for each respective acquisition year, used to normalize PV to the interval [0, 1]. For this study, NDVI_min and NDVI_max were determined empirically from each Landsat scene rather than using fixed global thresholds in order to account for inter-annual variability in phenological and atmospheric conditions.

Third, LSE (ε) was estimated as a linear function of PV following the parameterization [34]:

$986.0+\Lambda d \times 500.0=3$      (5)

The empirical intercept of 0.986 and the slope coefficient of 0.004 were derived from regression analysis that relates broadband emissivity to fractional vegetation cover across various land [34], representing land surface types, representing a weighted composite of bare soil emissivity (~0.986) and the incremental emissivity contribution of vegetation cover.

LST was calculated from the emissivity-corrected brightness temperature using the following equation:

$L S T=\frac{B T}{1+\left(\frac{\lambda \times B T}{\rho}\right) \times \ln (\varepsilon)}$          (6)

where, LST = land surface temperature (℃); BT = at-sensor brightness temperature (℃); λ = effective wavelength of emitted thermal radiance for Band 10 (10.895 µm); and ρ = h·c/σ_B = 1.4388 × 10⁻² m·K = 14,388 µm·K, in which h = Planck's constant (6.626 × 10⁻³⁴ J·s), c = speed of light in vacuum (2.998 × 10⁸ m·s⁻¹), and σ_B = Boltzmann constant (1.381 × 10⁻²³ J·K⁻¹).

2.5 Urban Heat Island Intensity classification

To enable a fair comparison of the spatial distribution of the UHI intensity of the three years under evaluation without being affected by differences in LST absolute values, the five UHI intensities examined were classified according to relative thresholds calculated as a function of both the annual mean (μ) and the annual standard deviation (σ) of the LST. The method employed is in line with Chen et al. [28], although it corrects for interannual differences in absolute LST values, thus making it possible to analyze the correspondence in the spatial distribution of the UHI levels between the 3 years under consideration. The classification thresholds and their interpretation are summarized in Table 2.

Table 2. Urban Heat Island (UHI) intensity classification thresholds and physical interpretations

It is important to note that because this classification employs year-specific relative thresholds, the boundaries between categories shift with each observation year in accordance with the annual LST distribution. Consequently, changes in the spatial extent of each UHI category between years reflect genuine redistribution of relative thermal intensity rather than changes in fixed absolute temperature ranges. The Natural Breaks (Jenks) algorithm was additionally applied as a complementary method for cartographic visualization purposes only, to produce maps that capture the natural clustering structure of the LST data distribution for each year.

2.6 Spatial autocorrelation analysis

The Moran index serves to detect the presence or absence of spatial autocorrelation in data globally. The value of this index ranges from −1 < I < 1 with a standardized spatial weight matrix. A negative value (−1 < I < 0) indicates negative spatial autocorrelation or a dispersive pattern, while a positive value (0 < I < 1) indicates positive spatial autocorrelation or a clustered pattern. If the value is close to zero, it means that no spatial pattern has been formed. To calculate spatial autocorrelation, one commonly used method is the Moran Index, which can be formulated through a mathematical equation as stated in previous studies [27, 28, 35]. The spatial unit consisted of 30 m × 30 m raster grid cells derived from Landsat-based LST pixels. A queen-contiguity spatial weights matrix was applied to define neighborhood relationships among adjacent grid cells in the Global Moran’s I analysis. Global Moran's I (I) was computed to quantify the overall degree of spatial autocorrelation in the LST distribution across the study area for each observation year:

$I=\frac{n}{\sum_i \sum_j w_{i j}} \times \frac{\sum_i \sum_j w_{i j}\left(x_i-\bar{x}\right)\left(x_j-\bar{x}\right)}{\sum_i\left(x_i-\bar{x}\right)^2}$       (7)

where, I = Global Moran's Index; n = total number of grid cells; wᵢⱼ = element of the queen-contiguity spatial weights matrix for cells i and j; xᵢ and xⱼ = mean LST values at cells i and j, respectively; and x̄ = global mean LST across all cells.

Moran's I values range from −1 to +1:

I > 0 → positive spatial autocorrelation (clustered pattern)

I < 0 → negative spatial autocorrelation (dispersed pattern)

I ≈ 0 → random spatial distribution

2.7 Spatial association analysis between Land Use/Land Cover and Urban Heat Island

To analyze the spatial relationship between the LULC types and UHI intensity classes, a post-classification cross-tabulation analysis was performed using the Tabulate Area tool in ArcGIS10.8. A cross-tabulation matrix was produced for each year of observation to assess the distribution of LULC classes across different UHI intensity classes. Temporal trend analysis (2014-2019-2024) was undertaken to identify systematic changes in the LULC-UHI spatial association over the decade-long study period.

Land cover and land use due to rapid urbanization (2014 to 2024) have had an important impact on UHI intensification. The rapid urban land use growth caused an increase in surface LST and the spatial distribution of thermal hotspots in densely populated subdistricts in Kendari City. These observations are highly consistent with the assumptions of UHI theory that urban warming occurs due to removal of vegetated surfaces and replacing them with impervious surfaces, resulting in decreased evapotranspiration, increased thermal mass and increased anthropogenic heat build-up [1, 2, 4]. The findings also indicate that urbanization-induced land-use changes are an important driver of developing UHI in Kendari, given the high LST values in densifying areas.

The increased magnitude of UHI in Kendari City can be explained using the classical theory of UHI proposed by Oke [1] and further elaborated by Voogt and Oke [2]. This theory involves the alteration of the energy balance of the urban micro-climate caused by changing land surface cover. In addition, when vegetation is replaced by built structures, it reduces evapotranspiration and the release of latent heat, which provides natural cooling in urban areas. For example, asphalt and concrete have a low albedo and high heat storage capacity, and thus they can absorb and retain solar energy during the day. Additionally, anthropogenic heat flux from urban activities such as the traffic intensity and energy consumption in buildings also increases urban temperature. The higher LST of the densely populated districts in Kendari is caused by the loss of the cooling function of vegetation and the strengthened heat storage and anthropogenic heat emissions within the canopy layer of the urban area.

Furthermore, the increase in Moran’s I indicates stronger spatial concentration of thermal hotspots over time, suggesting that high-temperature zones are becoming increasingly clustered within continuously urbanized areas. This finding is consistent with the urban thermal concentration hypothesis, whereby rapid urbanization promotes the spatial consolidation of heat accumulation as impervious surfaces become dominant within the urban landscape. This finding is consistent with the urban thermal concentration hypothesis, whereby active urbanization leads to a continuous chain of thermal hotspots as imperviousness becomes the dominant land cover in the local urban setting [1, 2]. Similar patterns of spatial thermal coalescence have been reported for rapidly urbanizing cities, such as Guangzhou, China [7], and Delhi, India [6], where built-up land expansion was associated with spatial coalescence of thermal hotspot areas and increased fragmentation of urban cooling areas. In Kendari, the loss of areas of vegetated and open-land patches appears to have weakened the spatial continuity of natural cooling systems and, in turn, contributed to an accumulation of heat in urban centers and development corridors.

Relative to other cities in Asia experiencing rapid urbanization, the increase in overall UHI in Kendari between 2014 and 2024 is relatively high, and relative to studies from Guangzhou [7], Delhi [6], and Harbin [22], the results suggest that secondary tropical cities may be particularly sensitive to a rapid thermal response under unstructured urbanization conditions. It has been suggested that the city's bay-enclosed landscape may limit ventilation by the atmosphere, thus trapping warm air masses within the urban canopy layer, as compared to open-plan urban centers [1, 2]. Second, the fragmented, market-driven urbanization of Kendari city, combined with limited application of climate-adaptive spatial planning and urban thermal mitigation strategies, may have contributed to the development of irregular urban morphology, low sky view factor, and meaningful longwave radiation trapping in highly built-up areas [19, 20]. Third, the limited extent and spatial fragmentation of green infrastructure restrict the effectiveness of evapotranspirative cooling and ecological buffering at the city scale. In line with previous studies, secondary cities with rapid urbanization, low topographical relief, weak spatial planning controls, and limited capacity for green infrastructure are particularly vulnerable to accelerated UHI intensification.

The findings of this study highlight the important ecological role of agricultural land, urban vegetation, and water bodies as natural thermal regulators within Kendari City’s urban landscape. The substantial reduction in rice fields and plantation areas during the 2014–2024 period indicates not only land-cover transformation but also the progressive loss of urban thermal buffering capacity. Rice fields, in particular, are good at cooling because standing water and thick vegetation support continuous evapotranspiration processes that keep surface temperatures lower than those of impervious urban materials [36]. However, the results suggest that the cooling effectiveness of these landscapes has become increasingly weakened as agricultural and vegetated areas are progressively enclosed and fragmented by expanding built-up development. This condition reflects the gradual degradation of the “urban cooling island” function, where formerly connected green and blue spaces become isolated within a warmer urban thermal matrix. The decline in the thermal performance of remaining forest and water-body areas further indicates that the effectiveness of urban ecological cooling systems depends not only on total area coverage, but also on spatial connectivity, landscape continuity, and protection from surrounding thermal pressure. These findings demonstrate that fragmented green spaces embedded within dense urban environments may no longer provide sufficient cooling benefits to offset the accelerating thermal accumulation associated with rapid urbanization.

In addition to biophysical features, UHI in Kendari City has an impact on public health, urban energy demand, hydrology, and urban sustainability. Rising surface air temperature in urban areas will increase the risk of exposure to heat stress and cause an increase in energy demand for cooling, surface runoff, and the risk of flooding due to the continuous formation of impervious areas. The study findings indicate that climate responsive spatial planning, urban green infrastructure and blue green corridor development should be integrated into future urbanization strategies. Furthermore, from a scientific perspective, the study contributes to the literature on UHI in tropical secondary cities through successfully integrating multi-temporal analysis of LULC change, satellite-based LST retrieval, and spatial autocorrelation within a single modeling framework. The findings therefore provide an empirical basis for supporting Sustainable Development Goal (SDG) 11 on sustainable cities and SDG 13 on climate action, while also offering spatially explicit evidence that can inform the revision of Kendari’s urban spatial planning and climate adaptation policies.