Makassar City is the capital of South Sulawesi Province as shown in Figure 1, located between 119º4ʹ29.038ʺE - 119º32ʹ35.781ʺE and 4º58ʹ30.052ʺS - 5º14ʹ0.146ʺS, with an area of 177.18 km2 (Makassar City Government, 2025). This city has elevations ranging from 0 to 25 meters with topographic characteristics of mostly flat to undulating, a gentle slope (0 – 5%), and a soil type dominated by Tropaquepts (67.04%) (Makassar Environmental Agency, 2024). In 2024, Makassar City had a population of 1,477, 861. The city’s mean annual air temperature is 27.75°C, and maximums often exceed 34°C during the dry season (April – October). The city’s mean relative humidity is lowest during the dry season, ranging from 69.94% to 80.34% (BPS-Statistics of Makassar Municipality, 2025).

The observational analysis period covers 2013 to 2024, using a multi-scale approach that starts at the city level to identify general trends, then progresses to the district level to examine spatial variations and comparisons between districts. Administrative units are adjusted according to changes in district boundaries. From 2013 to 2016, the analysis consisted of 14 districts. Then, from 2017 to 2024, the analysis covers 15 districts, including the Sangkarrang Islands district, which was officially formed in 2015 and separated from the Ujung Tanah district under Makassar City Regulation No. 3 of 2015 (Makassar City Government, 2015), and has been active since 2017.

This study uses data integration obtained from various sources, summarized in Table 1. Landsat 8 OLI/TIRS imagery serves as data for annual LST estimation and the calculation of spectral indices (2013 – 2024) during the dry season (April – October) to minimize the influence of cloud interference. The LST processing results were cross-validated using MODIS LST (MOD11A1) and correlated with air temperature data from the Indonesian Agency for Meteorological, Climatological, and Geophysics (BMKG) observation station and TerraClimate imagery. Meanwhile, for LCZ classification and validation, building morphology data from Google Open Buildings and very high-resolution imagery from aerial photographs, land use maps, and Esri World Imagery Wayback were utilized.

The data processing workflow for this study begins with LST estimation from Landsat 8 imagery and LCZ classification. These results are then integrated to calculate SUHII, followed by spatial pattern analysis. This entire process is illustrated in Figure 2.

The initial phase included a literature review and data processing, as detailed in Table 1 and Figure 2. For data processing, we used two platforms: ArcGIS Pro and Google Earth Engine (GEE). First, we use ArcGIS Pro to project datasets into the WGS 84/UTM 50 S coordinate system to ensure geometric consistency in area calculations. In parallel, the second platform, GEE, is used to derive spectral indices, including Normalized Difference Vegetation Index (NDVI), Normalized Difference Built-up Index (NDBI), Modified Normalized Difference Water Index (MNDWI), and Bare Soil Index (BSI), which represent vegetation, built-up areas, water bodies, and bare soil, respectively. These indices are prepared to calculate parameters used in the LCZ classification. Formulas and reference sources for these spectral indices are presented in Table 2.

To establish a uniform spatial framework, this study used a 100 m x 100 m grid clipped to the administrative boundaries of Makassar City. The grid size was matched based on the spatial detail of LCZ analysis for city-level assessment, with scale homogeneity for the environment, computational effectiveness, compatibility with supporting data for validation, and accordance with global standard protocols, World Urban Database and Access Portal Tools (WUDAPT) for LCZ analysis (Chen et al., 2023; Huang et al., 2023; Wang et al., 2024). This grid has served as the main spatial analysis unit, as container LCZ parameters and other variables such as LST and spectral indices were aggregated using Zonal Statistics (Mean). In parallel with this grid analysis, resolution harmonization is performed separately for cross-validation. Landsat LST data were temporarily downscaled to match the original resolution of the reference data (1000 m for MODIS, 4,638.3 m for TerraClimate) using a bilinear interpolation technique to enable direct statistical comparisons between pixels.

LST was derived from Landsat 8 TIRS Collection 2 Level 2 (Tier-1) imagery at Path 114/Row 64 during the dry season of 2013 – 2024 using GEE. The Single Channel Algorithm (SCA) method (Jiménez-Muñoz & Sobrino, 2003) was applied based on the availability of thermal and atmospheric thermal profile data, robustness, and compatibility with other data archives. To maximize data availability while guaranteeing quality, the image collection was filtered to include full scenes during the dry season. Scenes were filtered using a cloud cover threshold of ≤35%. Statistically, although scenes have a total cloud cover of up to 35%, the clouds are located outside the main study area (e.g., over Mount Lompobattang in Gowa Regency). For this reason, the QA_PIXEL band with bitwise operations was used to mask clouds and cloud shadows, providing data integrity in the study area despite potential cloudiness in other parts of the scene. Data processing involves extraction and calibration to obtain the original physical values of the intermediate bands by applying a predetermined scale factor (USGS, 2023). Specifically, Surface Thermal Radiance () from the ST_TRAD band, and both Upwelling () and Downwelling () atmospheric radiance from the ST_URAD and ST_DRAD bands, were calibrated using a factor of 0.001. Meanwhile, Atmospheric Transmittance () from ST_ATRAN band and Surface Emissivity () from ST_EMIS band were multiplied by a scale factor of 0.0001.

The retrieval of the LST () was performed using the SCA method on the calibrated intermediate bands. This process includes converting the surface radiance (), calculating the atmospheric correction functions () from the additional atmospheric bands, and integrating all components, including surface emissivity (), into the SCA formula as shown in Equations 1 - 4. The sensor-specific constants used for Landsat 8 band 10 are , and (Jimenez-Munoz et al., 2009; Nugraha et al., 2024; Sekertekin & Bonafoni, 2020).

where and are parameters derived from , , and . The LST values obtained in Kelvin were then converted to degrees Celsius.

To obtain a representative value and reduce residual noise, the annual median per-pixel value was calculated from all available cleaned LST observations, yielding a time series of 12 annual median LST (2013 – 2024). Due to the lack of in-situ LST measurements, accuracy assessment was conducted through multi source cross validation against: 1) MODIS LST products, 2) air temperature data from local BMKG meteorological stations (specifically the Paotere Maritime Meteorological Station in Makassar City, Sultan Hasanuddin Meteorological Station in Maros Regency and South Sulawesi Climatology Station in Maros Regency), and 3) air temperature data from the TerraClimate imagery dataset. The statistical accuracy of the LST processing results was evaluated using metrics such as correlation coefficient (Pearson and Spearman), the coefficient of determination (R2), Root Mean Square Error (RMSE), Mean Absolute Error (MAE), bias, linear regression slope, and intercept in the full Landsat scene. These validations assessed the regional-scale thermal consistency of the derived Landsat imagery rather than verifying pixel-level accuracy. Once the large-area consistency was confirmed, the LST time series was clipped to the administrative boundaries of Makassar City.

To analyze SUHII more objectively, this study uses the LCZ classification framework. LCZ is a classification system developed by Stewart and Oke (2012) to classify urban and natural landscapes into 17 standard classes based on surface structure, land cover, and thermal properties. This classification system is designed for local climate studies applied at neighborhood scales (100 m to several kilometers). The 17 LCZ classes are divided into 10 built types (LCZ 1-10) and 7 land cover types (LCZ A-G), as illustrated in Table 3. Each class represents a different condition, characterized by combinations of surface structure, land cover, and anthropogenic activities.

Source: Adapted from Stewart & Oke (2012)

LCZ classes were constructed based on a set of quantitative morphological and physical parameters. This study used six parameters within the global standard ranges for each LCZ class, as presented in Table 4. These parameters are derived for 2014, 2019, and 2024 to analyze temporal changes in the LCZ. Calculations were performed on a 100 x 100 m grid, as established in Section 2.4.2, which serves as the basic spatial unit for organizing all LCZ parameters. The data is used by integrating the Google Open Buildings Temporal V1 2.5D dataset with spectral indices from Landsat 8 to overcome limitations of historical building data. The formulas and data sources used to calculate each parameter are summarized in Table 5.

Source: Adapted from Stewart & Oke (2012).

The first phase of this parameter extraction process involved creating temporally consistent binary building masks for 2014, 2019, and 2024. This was achieved by integrating Google Open Buildings Temporal V1 2.5D (4 m resolution) as the main data, with Landsat-derived spectral indices as supporting data. The two primary bands from the main data used were building height, which contains height values, and building presence, which indicates the probability of a building pixel. The threshold of the building presence band was (≥ 0.34). Concurrently, a built-up spectral signature was identified from the corresponding year’s Landsat imagery using NDBI (> 0.2). To align these datasets both spatially and temporally, the 4 m building probability mask and its associated building height band were resampled to 30 m. The resampled building mask was then combined with the Landsat-derived NDBI mask using a logical AND operation. This final step produced a 30 m binary building mask where pixels meeting both criteria, a high building probability and a spectral built-up signature, were assigned a value of 1 (building), while all other pixels were assigned a value of 0 (non-building).

The second phase of the analysis focused on extracting built-up geometry and surface fractions within the standardized 100 m grid. Using the previously established 30 m binary building mask, the Height of Roughness Elements (HRE) was calculated by aggregating the underlying building height values using Zonal Statistics (Mean). Then, the Building Surface Fraction (BSF) value, expressed as a percentage, was derived by dividing the total building area by the grid area, which was summed using Zonal Statistics, as shown in Table 5. In that formula, is the binary pixel (1 for building, 0 otherwise), is the area per pixel (900 m²), and is the total area of the 100 m grid (10,000 m2). To characterize the other surface types, the Pervious Fraction (PSF) was calculated from the Fraction Vegetational Cover (FVC) derived from NDVI, the Water Surface Fraction (WSF) from MNDWI, and the Soil Surface Fraction (SSF) from BSI. Each index was normalized to 0-1 using min-max scaling, summed per pixel, and the mean value within each grid was used as the PSF. Finally, to ensure consistent total surface area, the Impervious Surface Fraction (ISF) was calculated by subtracting BSF and PSF from the total grid area.

The final stage consists of calculating the surface’s aerodynamic and radiative properties. Terrain Roughness (TR) was quantified as the effective roughness () for each grid cell, which was determined using the Davenport classification (Stewart & Oke, 2012) by calculating the maximum value between building roughness ( and vegetation roughness (), as shown in Table 5. Vegetation roughness () was determined based on the NDVI threshold corresponding to land cover types, as detailed in Table 6. In parallel, Surface Albedo (SA) was calculated as the weighted sum of Landsat 8 bands (Blue, Red, NIR, SWIR1, SWIR2), as specified in Table 4, to measure the ratio of reflected to incoming solar radiation. The resulting 30 m albedo raster was aggregated to the 100 m grid using Zonal Statistics (Mean).

The LCZ classifications for 2014, 2019, and 2024 were performed using a three-phase workflow that integrated hybrid Fuzzy Logic with a deterministic GIS-based rules approach and spatial aggregation, executed through a custom Python script. This method was chosen to accommodate continuous, gradual transitions in standard LCZ parameters Table 4 that involve overlapping class boundaries. Fuzzy Logic was calculating the degree of membership (μ) ranging from 0 to 1 to identify the transition zones that are uncovered or poorly represented by hard-threshold classifications (C. Chen et al., 2023; Estacio et al., 2019; Fonte et al., 2019). The challenge of this method with multisource input data is that it leaves unclassified grid cells. For that reason, the process was further processed using deterministic GIS-based rule classification (Y. Chen & Hu, 2022; Mo et al., 2024). The classification system uses a hierarchical decision tree with expert-designed hard thresholds to ensure complete spatial coverage and consistency.

This phase was carried out to identify the grid cells that best matched the ideal morphological range for each LCZ class. The spatial distribution maps for the six parameters across the three study years (2014, 2019, 2024) are provided in Appendix A. These parameters for each grid were first converted into degrees of membership (μ) using trapezoidal fuzzy membership functions (Table 7). Three function types that were adopted from (Chen, C., et al., 2023) were applied: (1) Standard Trapezoidal for parameters with an optimal central range (e.g., BSF), (2) Left-Shoulder for lower-is-better criteria (e.g., PSF in compact built-up classes), and (3) Right-Shoulder for higher-is-better criteria (e.g., HRE for high-rise classes). To accommodate spatial uncertainty and local variation, tolerance margins were applied to the standard LCZ thresholds (Kara et al., 2025): a 5% margin for percentage-based parameters (BSF, ISF, PSF), proportional margins for continuous physical-based parameters (HRE, SA), and an absolute margin of 0.5 for discrete Terrain Roughness (TR) classes. The final fuzzy suitability score for each grid cell for LCZ was calculated by aggregating the six membership degrees using the MIN operator Equation 5, meaning that a grid is assigned to an LCZ class only if all six parameters meet the class criteria. Thus, the overall score equals the minimum membership value across all parameters. Grid cells whose final highest score across all classes was below a threshold of 0.001 were labeled “Unclassified” and proceeded to the next phase for rule-based classification.

Source: Modified from C. Chen et al. (2023)

This deterministic phase classifies all grids labeled “Unclassified” in Phase I. This method is implemented using a decision tree with a hard threshold adapted from the LCZ standard and several previous studies, as presented in Figure 3 (Chen & Hu, 2022; Cilek & Cilek, 2021; Mo et al., 2024).

Source: Adapted from Chen & Hu (2022), Mo et al. (2024), and Cilek & Cilek (2021).

The final phase performs refinement to ensure spatial homogeneity and remove scattered salt-and-pepper pixels. Small or unclassified grid cells were merged into the neighboring polygon with which they shared the longest common boundary. This process evaluated the eight adjacent grids (the Moore Neighborhood) around each target cell, effectively smoothing the LCZ map and producing coherent, spatially contiguous zones.

The accuracy of the LCZ classification for 2014, 2019, and 2024 was validated against the 2024 Makassar City Block Aerial Photo, the 2024 Makassar City Land Use Map, and the Esri World Imagery Wayback for historical years (2014 and 2019). The stratified random sampling method was implemented to ensure proportional representation of all 17 LCZ classes. The total of validation sampling points exceeded the minimum target of 500 for each year’s observation, resulting in 563, 560, and 558 points for 2014, 2019, and 2024, respectively. The classification results were evaluated using a confusion matrix and statistical metrics, including Overall Accuracy (OA), Producer’s Accuracy (PA), User’s Accuracy (UA), and the Kappa Coefficient. The results are considered valid and acceptable for analysis if the OA is 80% or higher, in agreement with urban land cover studies (Ao et al., 2025).

Based on the LCZ framework applied by Stewart & Oke (2012), this study defines SUHII as the difference in Land Surface Temperature (LST) between specific Local Climate Zone (LCZ) classes. LCZ D (Low Plants) was selected as the reference zone for SUHII calculation in accordance with many previous studies (Deng et al., 2024; Han et al., 2024; He et al., 2024; Wei & Sobrino, 2024; Xi et al., 2024) that used this class as a reference zone. The selection of the reference zone is based on two main considerations: LCZ D is the dominant natural land cover class in Makassar City, and LCZ D has consistently been recorded as having a lower mean LST value compared to built-up land cover classes (LCZ 1-6, LCZ 8, and LCZ 10) throughout the 2013-2024 study period. SUHII for other LCZ classes () was calculated using Equation 6.

where is the mean LST of a specific grid cell (°C) and is the mean LST of the reference LCZ D class (°C).

The SUHII calculation involved temporal synchronization between LST and LCZ data because the LCZ classification map spans 5 years (2014, 2019, and 2024), while the LCZ data is available annually (2013-2024). To obtain conditions that represent urban morphology and thermal response, annual LST data are grouped and analyzed based on the LCZ map closest in time. This selection of temporal representations is based on the gradual development of urban structure and land cover patterns over a multi-year period, whereas LST fluctuates dynamically each year. The complete analysis period is defined in Table 8. The mean LST value is extracted for each 100 m grid cell by LCZ class using Zonal Statistics, enabling SUHII calculations for each year in the study period. However, the synchronization approach used in this study Table 8 assumes relative morphological stability within each LCZ period. Deviations between the actual morphological features and the classified LCZ map can bias SUHII estimates for the boundary years (2013, 2016, 2017, 2021, and 2022).

To quantify the spatial distribution and patterns of SUHII clustering, 100 m-resolution SUHII data for 2024 were aggregated to 300 m resolution to reduce microscale effects and improve computational efficiency for city-level analysis. This study used a two-tier autocorrelation approach. First, the analysis used Global Spatial Autocorrelation (Global Moran’s Index) to determine whether SUHII values were clustered, dispersed, or randomly distributed. This method was calculated using Equation 7, yielding values ranging from -1 (dispersed) to +1 (clustered) (Esri, 2025b; Getis & Ord, 1992).

where represents the deviation of the SUHII value at location from the global mean,  is the spatial weight features, is the total number of features, and is the aggregate of all spatial weights. Spatial relationships were defined based on K-Nearest Neighbors (k=8) with row standardization. The significance level is confirmed if the p-value is <0.05; a positive z-score indicates clustering, and a negative score indicates dispersion (Esri, 2025b).

Secondly, Hot Spot Analysis (Getis-Ord Gi*) involved pinpointing specific geographic thermal clusters. This method identified where high values (Hot Spots) or low values (Cold Spots) cluster significantly at the local scale. The Getis-Ord Gi* was calculated using Equation 8 (Esri, 2025a; Getis & Ord, 1992; Ord & Getis, 1995).

Where is the attribute value for feature , is the mean, and is the standard deviation. Significant Hot or Cold Spots are identified where p < 0.05 with high positive or low negative z-score, respectively. To ensure identified local clusters were valid, A False Discovery Rate (FDR) correction was applied to account for spatial dependency and multiple testing (Esri, 2025a).

The Landsat 8 LST processing results for Path 114 Row 64 were validated before being clipped to the administrative boundaries of Makassar City. This evaluation process consisted of cross-validation using multi-source reference data to ensure the LST's validity and consistency with actual field conditions. First, Landsat LST was validated against the MODIS LST product to assess consistency in LST values across a broader area.

As shown in Figure 4, the cross-validation indicates a very strong relationship between the two datasets in all years from 2013 to 2024. The Pearson correlation (r) consistently exceeded 0.98 annually. The coefficient of determination (R2) ranges from 0.91 to 0.99 (2023), indicating that Landsat data explains 91% to 99% of the variance in the MODIS reference data. Despite a persistent negative bias (Landsat cooler than MODIS) averaging between -0.58°C and -2.27°C, the overall precision is high, with RMSE values remaining low. The highest precision was achieved in 2023, with an RMSE of 1.72°C. However, this high correlation is affected by the aggregation of Landsat LST data (30 m) into coarser MODIS pixels (1,000 m), known as the Modifiable Areal Unit Problem (MAUP) (Briz-Redón, 2022). This spatial effect reduced variance in the comparison data. Therefore, these validation results are more indicative of regional-scale thermal consistency (≥ 1km) than absolute pixel-level accuracy at 30 m LST analysis scale.

Given the limited number of ground-based air temperature stations in the full Landsat scene (only three BMKG stations: Paotere Maritim Meteorological Station, Sultan Hasanuddin Station, and South Sulawesi Climatology Station), the validation strategy was augmented with a spatially continuous dataset. Figure 5(a) illustrates the relationship between Landsat LST and in-situ air temperature from these stations. As expected, due to fundamental physical differences between radiant skin temperature (LST) and ambient air temperature, the linear correlations are weak (r from -0.06 to 0.28). However, the agreement in magnitude is reasonable, with RMSE values ranging from 1.25°C at Paotere Maritime Meteorological Station to 2.34°C at South Sulawesi Climatology Station.

To ensure spatial and overall thermal consistency, TerraClimate air temperature data were also used as a validation reference. Prior to use, the reliability of this air temperature dataset was evaluated against the data of three BMKG stations Figure 5(b). The results showed a strong and significant association, particularly at the South Sulawesi Climatology Station (r = 0.81, p = 0.001). The annual mean temperature from TerraClimate captured annual trends but showed a warm bias. For instance, in 2019, Sultan Hasanuddin Meteorological Station recorded the mean air temperature of 27.42°C, while TerraClimate showed 28.27°C. The total bias across all stations was 0.57°C. Therefore, a bias correction was applied by subtracting 0.57°C from the TerraClimate data to align it with the local station data.

The relationship between Landsat LST and the bias-corrected TerraClimate air temperature was then assessed annually as an additional consistency check Figure 6. This comparison highlights the complex and scale-dependent relationship between two variables. Over the 2013-2024 period, Pearson correlation value (r) fluctuated between 0.44 and 0.71, indicating a general positive co-variation. The year 2018 showed the strongest linear relationship (r = 0.71, RMSE = 1.93°C).

Moreover, Figure 6 showed that several years (2014, 2017, and 2022) exhibited negative R2 values. This statistical result revealed significant physical and scale differences between the two datasets. Landsat LST captured spatial heterogeneity in surface temperature, while TerraClimate data represented homogeneous air temperatures at different vertical levels. The model is systematic, but it’s non-linear. The difference between the two datasets suggests that a simple linear regression model may perform worse than using the mean, resulting in a negative R2. While these results are tangled, the consistent positive correlation and reliability with MODIS LST and in-site stations confirm that Landsat-derived LST is statistically valid for SUHII analysis.

The spatial distribution of LST in Makassar City from 2013 to 2024 showed a fluctuating trend for the thermal environment, with the mean ranging from 25.80°C in 2022 to 32.28°C in 2023. During this 12-year observation period, 2023 was the hottest year with a maximum temperature of 40.15°C. This finding is supported by a recent study that reported a peak LST of 40°C in Makassar City in 2023 (Jumadi et al., 2025). The overall quantitative LST results, including minimum, mean, and maximum values, are summarized in Figure 7.

Based on Figure 7, during the observation period, the minimum LST ranged from 15.08°C to 24.59°C, while the maximum ranged from 31.77°C to 40.15°C. The mean annual LST has ranged from 25.80°C to 32.38°C, with 2022 and 2023 recording the lowest and highest LST values, respectively. These results indicated that certain years, such as 2015 and 2023, were influenced by climatic phenomena, whereas other years showed consistent intra-annual variability.

The spatial pattern shows persistent hotspots in the urban core.This study found that the mean LST ranged from 25.80°C to 32.28°C, with a significant increase of 6.48°C from 2022 to 2023. This is consistent with the findings of Jumadi et al. (2025), who reported an average LST range of 26°C to 35°C, with a 2023 increase of 8°C in Makassar City. 2023 was significantly affected by the El Niño phenomenon, which led to prolonged drought and elevated surface temperatures across the region (BMKG, 2025; NOAA, 2023).

Furthermore, despite fluctuating mean LST, Figure 7 also shows that the dynamic of maximum LST remains consistently high, frequently exceeding 34°C. Although the graph visualizes minimum, mean, and maximum values, the spatial pattern of LST is also found in the Standard Deviation (SD) values. SD has become an indicator to reveal thermal heterogeneity. Based on the above statistical records, the highest spatial variation occurred in 2017 (SD = 2.52) and 2013 (SD = 2.47), denoting a highly fragmented thermal pattern in which LST ranges vary markedly across land cover zones. In contrast, 2024 recorded the lowest variation (SD = 1.72), showing a more uniform LST distribution throughout the year.

The spatial distribution of LST is visualized in a multitemporal map in Figure 8. High temperatures (indicated by orange to red colors) are consistently concentrated in the city center, which is characterized by densely built-up areas. Conversely, lower temperatures (shown in green to blue) are observed in the eastern suburbs and northern coastal areas, where vegetation and water bodies provide a cooling effect on the surrounding environment. The minimum LST recorded was 15.08 °C in 2020. Given that the LST data used in this study were a composite of the dry-season median, these anomalous values reflect local cooling effects rather than transient sensor errors. This condition is associated with pixels located above water bodies or other natural areas, such as dense vegetation with high humidity, which maintain significantly lower temperatures even during dry months.

The accuracy of LCZ classifications for the years 2014, 2019, and 2024 was assessed using a Stratified Random Sampling approach. A total of 563 (2014), 560 (2019), and 558 (2024) sample points were validated against high-quality reference data. This included the 2024 Makassar City Aerial Photography, the Makassar City Land Use Map, and high-resolution historical imagery from Esri World Imagery Wayback for 2014 and 2019. The validation metrics are summarized in Table 9.

The results in Table 9 revealed highly consistent classification accuracy across the observation period. In 2014 and 2019, Overall Accuracy (OA) reached 94%, with Kappa coefficients of 0.92 and 0.93, respectively. Those values indicated nearly identical agreement between the classified maps and the reference data, which were assumed to represent actual ground conditions. Although the accuracy in 2024 (OA = 90%, Kappa = 0.88) falls within the range of results meeting the accepted classification threshold for land cover studies, namely 0.80 (80%) (Kotharkar et al., 2025; Liang et al., 2023).

Furthermore, the table shows different patterns in the evaluation of class-specific performance over the observation period. Certain classes appear to have very high accuracy, such as LCZ 1 (Compact high-rise), LCZ 4 (Open high-rise), LCZ 7 (Lightweight low-rise), and LCZ C (Bush, scrub), consistently achieve a Producer Accuracy (PA) of 100% with User Accuracy (UA) of 80% or higher, indicating that these different shapes are correctly identified. In contrast, some classes appear challenging due to fluctuating results, such as LCZ A (Dense trees), which show a stable but low PA of 60% across all years, indicating the classification model is less effective at correctly identifying this class despite its high UA (90%). This lower PA implies a 40% omission error, with dense trees likely misclassified into adjacent categories such as LCZ B (Scattered trees), potentially underestimating LCZ A’s total spatial extent and cooling contribution to the city. Similarly, LCZ 10 (Heavy industry) showed a maximum PA of 100%, but it varied and was relatively lower (60-80%), indicating that although all industrial areas are covered, they were sometimes identified as other building types. The stable PA values across almost every class from 2014 to 2024 indicated a robust and repeatable classification methodology, although UA fluctuated in some classes.

Analysis of the specific accuracy of each LCZ class was conducted to demonstrate the model’s reliability performance. LCZ 3 (Compact low-rise) is the dominant LCZ class in Makassar City, with User Accuracy (UA) ranging from 0.92 to 0.97 and a consistently high Producer Accuracy (PA) of 0.96 across all years. These results indicate that the class is reliably identified on maps and correctly classified in the field. Other building classes, such as LCZ 6 (Open low-rise) and LCZ 8 (Large low-rise), are also identified with high precision (UA > 0.83), with PA values at 0.84 and 0.94, respectively. Most importantly, LCZ D (Low plants), which serves as the reference zone for SUHII calculations, shows very high UA, with overall accuracy above 90%, reaching 0.99 in 2014, 0.97 in 2019, and 0.93 in 2024, while its PA remains stable at 0.82 across all observation years.

The figure illustrates the visual characteristics of the 17 identified LCZ classes, with each frame showing a specific urban built-up type or natural land cover type. The statistical validation results are visually presented in Figure 9, showing the 2024 LCZ classification overlaid on the corresponding aerial photographs. That figure represented the spatial delineation of each LCZ class in a 100-grid, showing the classified zones based on the observed urban morphological features in high-resolution imagery. These results provide visual evidence of the high statistical accuracy reported in Table 9. The high accuracy reflected in the PA and UA metrics for both built-up and natural land cover classes demonstrate the model's reliability in the LCZ classification, which accurately represents the physical urban morphology of Makassar City and effectively accommodates LST differences in the calculation and analysis of SUHII.

The validated multi-temporal LCZ maps for 2014, 2019, and 2024 are presented in Figure 10. The figure shows the dynamic spatial changes of Makassar City over the decade of observation. The figure shows that the city core remains densely built. While the peripheral area has faced expansion to impervious surfaces. The shift in urban structure from natural land cover to built-up land is clearly visible at several locations.

From the quantitative analysis of the overall area of the LCZ classification results, a clear and dominant pattern emerged: a systematic expansion of built-up areas (LCZ 1–LCZ 10) within the total study area of 17,728 ha. As detailed in Table 10, Makassar City’s urbanization is characterized by the spread and horizontal development of low to medium-rise buildings. The most obvious growth occurred in LCZ 6 (Open low-rise buildings), with an expansion of 948 ha, making this class represent the largest absolute increase of all classes. This is complemented by a steady expansion in LCZ 3 (Compact low-rise, +225 ha), LCZ 8 (Large low-rise, +268 ha), and LCZ 5 (Open midrise, +169 ha).

In contrast, natural land cover and water-permeable surfaces experienced a massive decline, as vegetation area used for agricultural purposes, such as rice paddies, was systematically converted into built-up areas. This conversion occurred in LCZ D (Low plants), decreasing the area by 870 ha as it transitioned to built-up zones, primarily in LCZ 6 (Open low-rise), LCZ 3 (Compact low-rise), and LCZ 8 (Large low-rise). Water bodies (LCZ G) also fell off significantly by 638 ha due to coastal reclamation and infrastructure development. Other natural classes, including LCZ B (Scattered trees, -93 ha) and LCZ F (Bare soil or sand, -270 ha), also experienced a reduction in area, further confirming the broader transformation.

The systematic land cover transition from natural areas, such as vegetation and water bodies, to built-up and impervious surfaces provides the basis for understanding LST changes in Makassar City. The replacement of surfaces with high evapotranspiration capacity with urban materials that absorb and store solar radiation forms, strengthens, and expands the SUHII area. Therefore, mapping this shift in urban morphology serves as a key indicator in monitoring and controlling the increase in SUHII.

The spatial distribution analysis continued with a quantitative description of the causes of morphological changes in Makassar City between 2014 and 2024. This section used the integration of multitemporal area trends, net changes, and detailed transition matrices to explain the systematic land conversion underlying LCZ changes. By directly identifying and explaining the “from-to” dynamics of LCZ changes, this analysis provided valid evidence needed to link physical urban growth with increasing LST in Makassar City.

Regarding temporal evolution, Figure 11 presents the multitemporal area changes for each LCZ, showing that the directional trend is consistent throughout the observation decade. This visualization reveals the continued growth of built-up areas, particularly LCZ 3 (Compact low-rise), which remains the dominant zone, covering over 6,000 ha and correlating with the continued decline of almost all-natural land cover. This demonstrated an overall pattern of expansion of built-up land classes across the study area. Furthermore, the results also show a clear shift towards a more fragmented landscape as natural land cover classes are systematically replaced by various urban morphologies ranging from open to dense.

Based on the trend shown in Figure 11, the analysis continued by examining the net area change for each LCZ class during the observation period shown in Figure 12. Figure 12 shows that each built-up land class recorded a net increase, with LCZ 6 (Open low-rise) experiencing the largest expansion (+948 ha), followed by LCZ 8 (Large low-rise, +269 ha) and LCZ 3 (Compact low-rise, +225 ha). In contrast, almost all natural land cover classes experienced substantial decreases, with the largest in LCZ D (Low plants, -870 ha) and LCZ G (Water, -638 ha). From these results, an inverse relationship is observed, with each class experiencing the greatest increase and decrease, indicating a direct, large-scale land conversion process.

The diagram illustrates the consistent directional shift in built-up over a decade.

To further understand the mechanism of LCZ change, the analysis is elaborated through the LCZ transition matrix in Figure 13. That figure provides definitive evidence by tracing the exact origin and destination of the converted land. The matrix shows that the dominant pathway of urbanization is the direct conversion of vegetated land into low-rise residential areas. Specifically, 670 ha of LCZ D (Low plants) were converted to LCZ 6 (Open low-rise), making it the largest contributor to land expansion for residential areas. This change was also supported by the conversion of LCZ D (Low plants) into 120 ha of LCZ 3 (Compact low-rise) and 84 ha of LCZ 8 (Large low-rise).

Another significant secondary pathway is the conversion of water bodies. This transition originated from LCZ G (Water), encompassing 140 ha, which was converted to LCZ E (Bare rock or paved) and 135 ha to LCZ D (Low plants). Related to this, 128 ha of water were converted to LCZ 6 (Open low-rise) and 91 ha to LCZ 8 (Low-rise large buildings), indicating a transition of coastal and inland water to accommodate land for the expansion of residential and industrial areas. These transformations represent a massive, structural decrease in the hydrological surface and in natural permeability, altering the energy balance of Makassar City.

Furthermore, Figure 13 also revealed the asymmetry and irreversibility of the transition during this observation period. It can be seen that more than 1,650 ha of natural land cover (LCZs B, D, F, and G) were converted to built-up land, while less than 30 ha of built-up land were converted back to natural area. This indicates that when natural land is converted to built-up land, it is highly unlikely to revert to natural land cover, as the land continues to be developed and used to meet the needs of the urban community. Moreover, the matrix also identified a stable urban core. This stability is evidenced by the large diagonal values, such as 5,886 ha of LCZ 3, which remained unchanged throughout the observation period. These results support the theory that urban expansion primarily occurs at the periphery through changes in surrounding natural lands.

Finally, based on the analysis above, Makassar City is undergoing a transition in urban morphology through a dual-lane conversion, characterized by expansion toward the urban fringe and the massive reclamation of water bodies in coastal and inland areas. This change involves replacing permeable surfaces with impervious materials. The danger of impervious materials was that they absorbed and stored heat, thereby altering the city’s surface energy balance. Consequently, the urban surface becomes more prone to heat retention and to the effects of climate change. These research findings also reveal the detailed, specific causes of the increase in LST and the intensification of SUHII through the transition pathway.

Further analysis was conducted to evaluate the relationship between LST and urban morphology, as represented by the LCZ. This analysis was based on previously obtained information indicating that the LCZ classification showed a massive transition from natural land cover to built-up land. Therefore, the relationship between LST and LCZ during the 2014 to 2024 observation period is important for systematically explaining the extent to which urban morphology influences LST values.

The annual LST distribution for each of the 17 LCZ classes is represented by Figure 14. Each of the twelve panels (a-l) represents each year of study. The thermal hierarchy is shown consistently. The boxplots in that figure showed LST values statistically across LCZ classes. They can be clearly distinguished into groups based on built-up areas, industrial zones, and natural land cover. This spatial structure serves as the basis for analyzing the SUHII pattern.

Statistical analysis of this multi-year dataset revealed stable LST ranking values. Densely built-up classes consistently record the highest LST. For instance, in the hottest year (2023), LCZ 3 (Compact low-rise) and LCZ 2 (Compact mid-rise) reached median LSTs of 33.30°C and 32.72°C, respectively, while in the other years, median LSTs were around 31°C. Moreover, LCZ 10 (Heavy industry) exhibited high LST, with a median reaching 34.96°C in 2023 and a maximum exceeding 38°C. This is consistent with its characteristics as a large-scale anthropogenic heat source. Another industrial area classified as LCZ 8 (Large low-rise) also exhibited high LST, but still lower than LCZ 10 (Heavy industry).

In contrast, natural land cover classes such as LCZ A (Dense trees), LCZ C (Bush, scrub), and LCZ G (Water) are consistently considered as urban cooling zones, with median LST often 5-6°C lower than the highest built-up zones (e.g., 27.35°C for LCZ C in 2023). Meanwhile, the selected class, LCZ D (Low plants), showed interannual variability in LST of ±1.5–2.2°C. The LST in LCZ D (Low plants) ranges from 24.37°C (2022) to 32.17°C (2023), with a median of 26.37°C. This 7.8°C fluctuation is not limited to LCZ D; 16 other LCZ classes also experienced a drastic increase in LST in 2023, with a mean of +6.45°C relative to 2022. Some classes show larger or comparable fluctuations (e.g., LCZ 10: +9.28°C, LCZ F: +7.72°C, LCZ 8: +7.14°C, and LCZ 9: +7.07°C). This indicated that the large LST change in 2022 – 2023 represents an externally driven anomaly (effect of the Super El Niño) rather than an inherent thermal characteristic. In addition, LCZ E (Bare rock or paved) and LCZ F (Bare soil or sand) also showed LSTs comparable to those of built-up zones in some observation years. This revealed that open impervious surfaces, although included in the natural land cover, contribute to heat, rather than act as cooling zones.

The temporal analysis in Figure 14 also confirmed the LCZ-based thermal hierarchy against interannual variability. Although the global El Niño phenomenon affecting regional climate increased median LST and compressed the Interquartile Range (IQR) in years such as 2015 and 2023, the widespread increase in LST maintained the relative ranking of the classes. This stability suggests that LST differences are related to the biophysical properties of each LCZ class. Moreover, the lower mean LST of LCZ D compared to the built-up LCZ class is relatively consistent, and its spatial representation as the dominant natural land cover in the study area further confirms its suitability as a reference zone for calculating SUHII. The systematic increase in LST in the built-up class is due to the city's morphological transition, as explained in Section 3.2.3.

Surface Urban Heat Island Intensity (SUHII) was calculated spatially and temporally in Makassar City. The SUHII was obtained by subtracting the mean LST of each pixel from the mean LST of LCZ D (Low plants), which served as the reference zone in this study. The calculation results provided information on how warm or cool a particular urban morphology is relative to the surrounding natural land cover.

To facilitate interpretation of the physical meaning of these calculations, the SUHII values are interpreted as thermal anomalies relative to the reference zone. A positive SUHII (> 0°C) represents the Surface Urban Heat Island (SUHI) effect, where the observed pixel is warmer than the reference zone. A positive SUHII is caused by zones that use impervious materials that store heat, such as LCZ 3 (Compact low-rise) and LCZ 10 (Heavy industry). Conversely, a negative SUHII (< 0°C) indicates the Surface Urban Cold Island (SUCI) effect, where the pixel is cooler than the reference zone. The negative SUHII was typically found in natural land cover. Finally, a neutral SUHII (around 0°C) represents a surface that is thermally equivalent to the reference zone, which serves as the surface energy balance point.

In spatial terms, the annual SUHII distribution map is presented in Figure 15, with a symmetrical scale of -9°C to 9°C to facilitate the representation of positive SUHII values as orange to brown. These positive SUHII values are concentrated in densely populated urban centers and gradually expand into the suburbs. This spatial expansion is related to the morphological transition documented in Section 3.2.3, in which natural land cover is systematically replaced by built-up areas that store and absorb heat. Conversely, negative SUHII values are shown in dark grey to black. These SUHII values are associated with water bodies and vegetated areas that provide a thermal counterbalance to the expansion of positive SUHII.

Regarding temporal dynamics, Figure 16 illustrates the interannual trends and annual variability of SUHII over the 12-year study period. Maximum intensity values were highest in 2017 (8.51°C) and 2021 (8.05°C), indicating greater thermal imbalance in those years. The SUHII distribution also shows that the period from 2017 to 2022 was distinguished by a consistent increase in SUHII, with a mean peak of 1.46°C occurring in 2018. A significant change in SUHII values occurred in 2023, when the mean SUHII dropped sharply to 0.01°C, indicating that in that year the mean LST of built-up was almost the same as the reference zone. However, a recovery was observed in 2024 to 0.87°C. Despite that variability, the stable, positive maximum values indicate that SUHII occurs with high intensity every year. Meanwhile, the negative minimum value throughout the observation period (reaching -7.43°C in 2014) indicates that the LST in other natural land cover classes, such as LCZ G (Water) and LCZ A (Dense tree), is lower than that in the reference zone, LCZ D (Low plants).

The sharp decrease in the mean SUHII to 0.01°C in 2023 was caused by the occurrence of a super El Niño. This condition caused absolute temperatures across the LCZ to reach their highest intensity. For example, LCZ 3 (Compact low-rise) reached a median of 33.30°C. Then, severe drought also encompassed the reference zone, LCZ D (Low plants). The loss of evaporative cooling in this reference zone caused a sharp increase in LST. Therefore, the LST between built-up and natural land cover was almost equal, resulting in a mean SUHII approaching zero. However, the maximum SUHII persistence remained high at 7.63°C in 2023. This indicates that, despite a narrowing of the mean gradient, intense anthropogenic hotspots remained highly active under these conditions.

Ultimately, based on spatial and temporal analysis, it was found that the intensification of SUHII in Makassar City is a direct consequence of its continuously evolving urban morphology, and in certain years, such as 2015 and 2023, SUHII values were also significantly influenced by the El Niño phenomenon. The transition from natural land cover to built-up, driven by the increasingly widespread use of impervious materials, led to the mean SUHII value remaining continuously positive and to the area expanding. This condition prompted a warning to increase vigilance for heat hazards that were not only focused on the urban core but also a phenomenon that developed and expanded its impact as urban physical growth expanded.

This section continues the analysis by grouping SUHII by LCZ classification across three distinct time periods (2014, 2019, and 2024). The analysis aims to uncover how specific morphological configurations and material properties contribute to SUHII. A detailed temporal segmentation of these time periods enables an examination of specific land conversions and their impact on the surface energy balance. Ultimately, the results reveal a hierarchy of thermal contributors and sinks that have evolved alongside Makassar City's rapid urban expansion.

During this early period, the spatial distribution of SUHII is visualized in Figure 17. That figure showed a massive positive SUHII expansion to the east and south of the city. Thus, positive SUHII is no longer confined to the city core but has expanded significantly into the suburbs. This shift described the conversion of natural land cover to impervious surfaces and to built-up areas.

The statistical data in Figure 18 show that LCZ 3 (Compact low-rise) had the highest SUHII in 2013, with a mean of 3.10°C. This high intensity is related to the finding that LCZ 3 is the most dominant built-up class in Makassar City, with an area of 5,887 ha. Positive SUHII in this zone is caused by using zinc or concrete roofs combined with low-density building structures with airtight surface materials. In this zone, the density of built-up areas, besides storing and absorbing heat, also limits surface air circulation due to a lack of vegetation.

Furthermore, LCZ 1 (Compact high-rise) and LCZ 2 (Compact mid-rise) showed stable, positive SUHII values of around 2.00°C. This condition reflects the high heat storage capacity of high-rise building materials. Meanwhile, morphologies with lower density, such as LCZ 9 (Sparse buildings) and LCZ 8 (Large low-rise), showed the highest intensity at 0.60°C in 2013. In addition, LCZ E (Bare rock or paved) was also found to produce positive SUHII values above 0.70°C consistently. Although LCZ is included in the natural land cover class, the expansion of impervious surfaces, such as asphalt road networks and concrete parking areas, tends to strengthen SUHII rather than reduce it. Moreover, LCZ F (Bare soil or sand) also began to shift from negative to positive SUHII in 2016 due to seasonal moisture variability. This pattern emphasizes that the focus of SUHII prevention and management should not be solely on built-up areas, but also on non-residential impervious surfaces and unmanaged open land. These two zones are the potential sources of thermal degradation in Makassar City.

Regarding the natural cooling mechanism, LCZ G (Water) and LCZ A (Dense trees) remain consistent heat sinks with negative mean SUHII values of -2.55°C and -2.31°C, respectively. However, the lowest negative SUHII was observed in LCZ C (Bush, scrub) in 2014, reaching -5.66°C. The fluctuation occurred mainly due to the high evapotranspiration in this LCZ C class and the influence of sea breezes, consistent with the finding that LCZ C in Makassar City is in an island area.

Regarding the mid-study expansion, the spatial distribution of SUHII in Figure 19 showed that thermal concentration shifted towards industrial zones and high-density built-up areas. The largest positive SUHII identified in Figure 20 was experienced by LCZ 10 (Heavy industry), which reached a mean peak of 3.91°C in 2018. This value represents the largest positive SUHII within the LCZ class across all observation periods. This spike was caused by heat absorption from the metal roofs, extensive impermeable concrete surfaces, and the heat emitted by anthropogenic activities in that large-scale industrial zone.

Meanwhile, in LCZ 1 (Compact high-rise) and LCZ 3 (Compact low-rise), they consistently showed the highest mean positive SUHII values of 2.77°C and 3.06°C, respectively. These results indicate that the dense, high-rise building development phase effectively and intensively traps and stores solar radiation, depending on the materials used in the buildings. In addition, a transition in SUHII values was observed in LCZ F (Bare soil or sand), which in the previous LCZ period had a neutral or negative value, but shifted to a positive value, reaching 1.39°C during this period. This shift was found in land clearing without vegetation cover, thus turning this area into an active heat contributor.

Although LCZ A (Dense trees) and LCZ G (Water) continue to serve as cooling zones, their combined area appears to be declining, a serious concern for urban resilience. In 2019, LCZ A (Dense trees) had a negative SUHII of -3.44°C, acting as a buffer against the massive warming of the Makassar mainland. The cooling capacity of LCZ A is under pressure, as evidenced by a significant reduction in area from 301 ha in 2014 to 267 ha in 2019. This reduction in natural absorption requires concrete intervention from the local government through targeted city planning and policies to protect the remaining green corridors.

In the final assessment phase, an analysis was conducted using the 2024 LCZ period to evaluate contemporary SUHII patterns in Makassar City. As shown in Figure 19 (f-h), positive SUHII is increasing and expanding. This condition is also influenced by the increasingly unstable climate from 2022 to 2024. The statistical details presented in Figure 21 indicate that 2022 and 2024 remain the positive SUHII hierarchy, with the built-up class having a mean SUHII value above 1.49°C, except for LCZ 7 (Lightweight low-rise). LCZ 7 in this study was found on islands, not on the mainland, so the SUHII in this class is relatively lower than in other built-up classes because it receives a cooling effect from wind-sea. However, in 2023, there appears to be thermal equalization. This is identified as being influenced by the super El Niño, causing the mean SUHII to decrease to 0.01°C. This decrease in value is not due to increased urban cooling, but the reference zone, LCZ D (Low plants), recorded an increase in LST to 32.28°C. This LST value narrows the difference between the built-up zone and natural land cover.

Despite a decrease in the mean SUHII value during 2023, LCZ 10 (Heavy industry) continues to dominate with a mean positive SUHII value of 2.87°C and 2.75°C in 2024. These results reinforce the finding that heavy and large-scale industrial materials in Makassar City absorb heat. This zone is certainly very dangerous if not controlled, especially when global climate change affects the region, increasing the LST of the reference zone. A similar thing happened to LCZ 3 (Compact low-rise), which experienced a decrease in the mean SUHII to 0.94°C in 2023, but increased again to 1.85°C in 2024. This trend shows that densely populated residential areas remain a major contributor to the heat that city residents face.

In contrast, natural land cover, particularly LCZ G (Water) and LCZ C (Bush, scrub), continued to record negative SUHII values. In 2023, the LCZ C reached a mean SUHII of 4.96°C, becoming the highest cooling that occurred on islands influenced by the surrounding sea. This revealed that a dominant wind-sea-cooling effect surrounds the islands. Meanwhile, on the mainland, water bodies (LCZ G) serve as the main heat sink, with a mean SUHII of -4.14°C in 2023. Moreover, the negative SUHII of LCZ A (Dense trees) in 2023 at -3.74°C also confirms that dense vegetation maintains a superior cooling pattern rather than low-density vegetation. These findings demonstrate that preserving maritime assets and natural infrastructure on the mainland is essential to preventing total thermal saturation, especially under intense climate anomaly conditions.

This section examines the spatial and temporal variations of SUHII across all districts in Makassar City from 2013 to 2024. The analysis was conducted to identify the district administrative areas with the greatest impact on positive SUHII. Ultimately, these results illustrate the development trends of thermal characteristics over a 12-year period.

The distribution of the mean SUHII per district in Makassar City is presented in Figure 22. That figure shows that the positive SUHII is found in the urban core. The districts included Makassar, Bontoala, Mamajang, and Rappocini, consistently recorded mean positive SUHIIs, often exceeding 2.00°C. This strongly confirms the intensification of the SUHI phenomenon. The highest positive SUHII value was found in the Mamajang district in 2013, with a mean of 4.39°C. Districts in the urban center are characterized by high population density and a dense building morphology, with low- and mid-rise buildings, which are identified as the main thermal contributors in Section 3.5. In Contrast, districts on the urban periphery, including Tallo and Tamalanrea, show much lower mean SUHII, even negative in some years. This is because these areas have an extensive cooling zone. While the Sangkarrang Islands generally have a negative SUHII, in certain years, such as 2021 and 2022, the value is above 0. However, the SUHII in this area still appears to be under control, as it returned to negative values in 2023 and 2024.

Districts with high positive SUHII values, such as Bontoala, Makassar, Mamajang, and Rappocini, had the highest mean SUHII values throughout the observation year (Figure 23). Mean SUHII in those districts exceeding 2.00°C indicates thermal resilience to annual volatility. Conversely, an increase in positive SUHII was observed in Manggala and Tamalate after 2020, while peripheral expansion zones such as Biringkanaya and Tamalanrea showed a gradual increase, reflecting rapid industrialization and land conversion. Moreover, other districts, particularly Ujung Pandang and Panakukkang, exhibited high variability in mean SUHII due to their high sensitivity to changes in surface properties. Finally, Tallo and Tamalanrea consistently recorded the lowest values, often falling into the negative range due to their greater proportions of vegetated land, water bodies, and impervious surfaces. Therefore, although Tamalanrea experienced an increase in positive SUHII, the cooling zone in the region is still functioning well. It should also be noted that Tamalanrea is the second-largest district in Makassar, and it is divided into several land zones.

The impact of the 2023 El Niño was clearly visible as a regional event affecting SUHII values. This year marked a reset point, with mean SUHII in almost all sub-districts dropping sharply to near zero. For example, in the Makassar district, the mean SUHII value dropped from 3.96°C in 2022 to only 0.78°C in 2023. This trend indicates that this phenomenon lowered SUHII values as the reference zone experienced increased intensity. However, a direct recovery was observed in 2024, so the 2023 SUHII value was heavily influenced by climate effects and not the permanent SUHII value in Makassar City.

The SUHII spatial configuration analysis was performed at a 300 m resolution. The original 100 SUHII data was aggregated to a 300 m scale. This was performed to reduce microscale noise and optimize computational efficiency for evaluating city-scale patterns. This aggregation process resulted in an area difference of approximately 2 ha (0.01%), increasing the total study area from 17,728 ha to 17,730 ha. This increase was stated not to significantly affect the results and validity of the analysis because the percentage difference was very small. The SUHII spatial analysis was conducted in two stages: a global evaluation of spatial autocorrelation, followed by the identification of significant cluster locations.

Global spatial patterns were evaluated using the Global Moran’s Index. This analysis aimed to statistically confirm the spatial clustering pattern in the SUHII processing results for Makassar City. As illustrated in Figure 24, representative graphs for 2013 and 2024 show that the results fall within the positive tail of the normal distribution. This finding indicates that the clustering pattern of the SUHII distribution in Makassar City is structurally organized and not random.

As summarized in Table 11, Moran’s Index results show high values between 2013 and 2024, ranging from 0.70 to 0.79. With a p-value of 0.0000 and a z-score consistently above 65. Therefore, the clustering pattern is a result of chance, with a probability of less than 1%. The results indicate a clustered, strong, non-random spatial structure of SUHII in Makassar City, with distribution throughout the city. It can be seen that, throughout the observation period, SUHII in Makassar City showed a continuous clustering pattern. This means that areas with similar heat intensity were geographically clustered throughout the study period.

After the global spatial pattern was statistically confirmed through the Global Moran’s Index, a Getis-Ord Gi* analysis was conducted to map the geographic evolution of SUHII clusters during the study observation period. The results of this hotspot analysis are visualized in a spatial significance map in Figure 25. It can be seen that significant positive SUHII (Hot Spot) with 99% confidence is consistently found in densely populated urban centers, including Makassar, Bontoala, Mamajang, and Rappocini districts, as well as industrial clusters in the eastern expansion zone. Meanwhile, a significant negative SUHII (Cold Spot) is observed in the northern part, including the Tallo, Tamalanrea, and Biringkanaya districts. The Hot Spot and Cold Spot cluster patterns are geographically resilient and are associated with specific LCZ classes. Therefore, the spatial distribution of SUHII in Makassar City shows high stability despite small annual fluctuations.

Significant Cold Spots are mostly located in water bodies, primarily in the Tamalanrea and Tallo districts. These blue-green infrastructure locations are important to preserve as they serve as a resource to mitigate the increase in Hot Spots. While the identified Hot Spots in the urban center remain geographically stable, Figure 25 also reveals annual fluctuations in suburban areas, which are not significant zones. These transition areas are found to be sensitive to rapid land-use changes and regional climate events, such as El Niño in 2023. These results indicated that both the urban core and the suburban area require proactive land management policies to prevent further hotspot expansion.

While the spatial map in Figure 25 distinguishes among different levels of statistical confidence (90%, 95%, and 99%), the area percentage in Figure 26 combines these categories into Hot Spot and Cold Spot categories. These spatial significance levels of Makassar City’s annual SUHII are classified into Hot Spot, Cold Spot, and Not Significant zones. The integration provides a clearer pattern of the overall thermal trend. The Not Significant zone accounts for the majority of the urban surface, ranging from 66% to 72% in most years.

The findings show that the coverage of the Hot Spot area decreased from 28.10% (4,986 ha) in 2023 to 20.60 (3,654 ha) in 2024. Likewise, the Cold Spot area also showed a decrease from 21.40% (2013) to 17.90% (2024). Meanwhile, the Not Significant zone expanded by 61.50% of the total urban area in 2024, as shown in Figure 26. This indicates that the shift in urban morphology has become more homogeneous in its thermal characteristics.

This shift in the significance pattern indicates continued thermal homogenization across Makassar City. While the decrease in Hot Spot area initially appears to be a positive improvement, it suggests that SUHII is becoming more diffuse and widespread across the entire area. In this regard, the decrease in Cold Spot coverage is an important warning. It reveals the loss of the cooling zone as a protective barrier for the city due to the massive expansion of built-up land. Furthermore, the substantial increase in the Not Significant zone indicates that Makassar City is transitioning into “sea of heat”, with thermal contrast disappearing and the reference zone becoming more uniform.

The relationship between urban morphology and the SUHII Hot Spot is further analyzed through visualization in Figure 27. The figures show the annual percentage of Hot Spot coverage for each LCZ. These results reveal a consistent and systematic pattern. In the built-up class, LCZ 3 (Compact low-rise) and LCZ 2 (Compact mid-rise) show the highest and most stable percentages of Hot Spot areas across all years. In addition, LCZ 10 (Heavy industry) and LCZ 1 (Compact high-rise) also show Hot Spot dominance in some years, reaching 100% hotspot saturation. This very high concentration indicates that industrial activity and high-density urban structures are the main physical drivers of the SUHII Hot Spot clustering in Makassar City.

In contrast, natural land cover LCZ classes such as LCZ A (Dense trees), LCZ D (Low plants), and LCZ G (Water) consistently showed very low or near-zero SUHII Hot Spot coverage percentages. This suggests the need to preserve these natural land cover classes as they provide urban cooling. The stability of these patterns revealed that despite the rapid urban land-use transition, their thermal characteristics remain structurally tied to their physical morphological properties. Therefore, urban planning efforts should prioritize preserving these natural coolers while targeting densely built-up areas for rapid, targeted heat mitigation.

In addition, the spatial analysis continued with a Thermal Exposure Classification analysis for 15 districts in Makassar City, as shown in Figure 28. This classification categorizes exposures as Very High (> 50%), High (25 – 50%), Moderate (0 – 25%), and Low (0%). Exposure class is defined as the percentage of Hot Spot pixels within each district's total administrative area. Based on Figure 28, SUHII Hot Spots were found in districts in the urban center, especially Bontoala and Makassar. These districts remained in the Very High category from 2013 to 2022, with hotspot saturation exceeding 50%. This indicates that the dense urban structure acts as a SUHII Hot Spot cluster, largely unaffected by annual climate variations, but rather by the urban morphology. Similarly, Rappocini and Mamajang districts showed high thermal exposure, often at the Very High exposure level, throughout the study period. By 2024, these districts had the highest exposure.

Although suburban districts such as Biringkanaya and Tamalanrea are mostly classified as having moderate exposure, they exhibit high sensitivity to climate change and development. This is evident in the sharp increase to Very High or High exposure levels in 2023. Meanwhile, the Sangkarrang Islands had low exposure to the SUHII Hot Spot throughout the observation period. The location of the islands surrounding the ocean acts as a heat sink, preventing the formation of the Hot Spot SUHII. This cooling effect ensures that the islands are stable and unaffected by mainland SUHII.

The findings of this study reveal four interconnected outcomes. First, annual LST analysis confirms a general upward thermal trend across Makassar City between 2013 and 2024, punctuated by an anomalous spike in 2023 attributable to the super El Niño event. Second, the 5-year interval LCZ classification documents a clear morphological shift from natural land cover toward compact and open low-rise built-up classes over the observation period. Third, LCZ D (Low plants) functioned as a reliable thermal reference zone throughout most of the study period, though its buffering capacity was critically compromised during 2023. Fourth, spatial analysis reveals SUHII distribution is transitioning from a concentrated hotspot structure toward diffuse thermal homogenization.

During the observation period (2013-2024), Makassar City underwent an evolution, marked by the transition from natural land cover to impervious surfaces and built-up areas. In 2024, LCZ 3 (Compact low-rise) become the largest built-up area, covering 6,112 ha, followed by LCZ 6 (Open low-rise), at 3,865 ha. These increases in built-up areas were in line with a decrease in natural land cover of more than 1,650 ha, particularly in LCZ D (Low plants) and LCZ (Water). Despite significant decreases, LCZ D and LCZ G retained the largest areas in 2024, at 2,060 ha and 1,924 ha. To quantitatively evaluate this trajectory, a point-biserial correlation analysis was performed to examine the relationship between urban surface types (coded as built-up versus natural land cover) and mean SUHII across the observation periods (N = 51). The analysis revealed a strong, highly significant positive correlation (), providing empirical evidence that the systemic conversion from natural to impervious classes directly drives the city’s thermal intensification. However, that massive conversion altered the urban surface energy balance. Research by Liong et al. (2021) confirmed that changes in land use and land cover are driving factors for the increase in LST in Makassar City. This is further supported by Arifin et al. (2023), who identified that massive urban expansion resulted in a significant increase in total built-up area compared to previous decades. The conversion of permeable and evapotranspirable zones into heat-absorbing and heat-storing materials was the main factor causing the intensification and expansion of positive SUHII in Makassar City.

In parallel, a stable thermal hierarchy throughout the study period (2013-2024) confirms the LCZ framework as a robust indicator of SUHII distribution. LCZ 3 (Compact low-rise) and LCZ 10 (Heavy industry) consistently function as major heat emitters. LCZ 3 contributes to positive SUHII due to its dense material configuration. These materials absorb heat, such as in zinc and concrete roofs. Moreover, the lack of vegetation prevents efficient heat dissipation. Consequently, this area is warmer than the reference zone. This is in line with the findings of Ardiyansyah et al. (2021), who identify LCZ 3 as the largest built-up LCZ class in Makassar City and noted that LCZ 3 and LCZ 10 encompass more than 80% of the total built-up category.

On the other hand, natural land cover such as LCZ A (Dense trees), LCZ C (Bush, scrub), and LCZ G (Water) still function as cooling zones, providing a cooling gradient of 5-6°C relative to built-up areas. For example, in the hottest year of 2023, mean LST reached 35.14°C for LCZ 10 (Heavy industry), 33.34°C for LCZ 1 (Compact high-rise), and 33.21°C for LCZ 3 (Compact low-rise). Natural land cover is cooler with 27.35 for LCZ C, 28.14°C for LCZ G, and 28.53 for LCZ A. That cooling effect is provided through vegetation shading, evapotranspiration, and water thermal characteristics (Deng et al., 2026; Faridatul et al., 2024; Isinkaralar et al., 2025). Moreover, these zones are affected by the connectivity of their location. Increasing connectivity and maintaining large and continuous areas of vegetation and water can maximize their cooling benefits (Zhang et al., 2026).

Although LCZ A consistently functions as a primary cooling zone due to its high evapotranspiration rates, its spatial extent may be systematically underestimated given the 40% omission error documented in Section 3.2.1, where dense tree canopy was likely misclassified into adjacent categories such as LCZ B (Scattered trees). This misclassification implies that the mapped distribution of dense trees underrepresents the true extent of the city’s green cooling infrastructure. Consequently, the total cooling capacity contributed by LCZ A in Makassar City is likely greater than what is captured in this analysis. However, since the UA for LCZ A remains high (90%), the LST values extracted from successfully identified patches still reliably reflect the true thermal behavior of dense tree zones, ensuring that thermal characterization of LCZ A remains valid despite its spatial underestimation.

However, the stability of the thermal hierarchy was disrupted by the 2023 super El Niño event, which led to near-complete thermal equalization between built-up and natural land cover, bringing the mean SUHII to 0.01°C. During this period, mean LST across Makassar City ranged from 25.80°C to 32.28°C, with an increase of 6.48°C from 2022 to 2023. This result is consistent with Jumadi et al. (2025), who reported a mean LST range of 26°C to 35°C for Makassar City with a spike of 8°C in 2023. This disruption was driven by a disproportionate heating of LCZ D (32.28°C), the reference zone, which reduced its thermal contrast with built-up areas such as LCZ 3 (33.22°C) to only 0.94°C. Notably, this response was not observed during the 2015 El Niño event, despite stronger dry-season (April-October) oceanic forcing, with Oceanic Niño Index (ONI) values reaching 1.0–2.5 in 2015 compared to 0.6–1.8 in 2023 (NOAA, 2026). The asymmetric LST response of LCZ D between these two El Niño events suggests that oceanic forcing intensity alone does not determine the thermal behavior of the vegetated reference zone. Instead, the 2023 collapse is more likely attributable to compound climate conditions, including a strong positive Indian Ocean Dipole (IOD) that amplified rainfall deficits during the dry season, an abrupt soil moisture deficit resulting from the transition from the preceding triple-dip La Nina (2020 – 2022) to El Nino conditions which suppressed evapotranspiration in LCZ D, and the progressive reduction of LCZ D spatial extent by approximately 870 ha due to urbanization which weakened its thermal buffering capacity. These findings suggest that the effect of El Niño on SUHII reference zone stability is mediated by compound climate forcing and concurrent LST conditions, limiting the direct generalizability of the 2023 response to future El Niño events of similar or greater oceanic forcing magnitude.

Furthermore, the spatial configuration analysis shows that SUHII in Makassar City is structurally clustered, with a trend warranting attention that points toward thermal homogenization. Although the Hot Spot area (90-99% confidence) decreased from 28.10% in 2013 (4,986 ha) to 20% (3,654 ha) in 2024, the expansion of the Not Significant zone to 61.5% indicates a shift toward thermal homogenization or a uniform “sea of heat”. This occurs with fewer urban cooling zones in different locations. This finding is supported by the observation that high-density impervious surfaces produce higher temperatures. It also contributes to thermal uniformity across the urban surface (Ardiyansyah et al., 2021; Xu et al., 2025). This process towards homogenization reveals an increase in the temperature of the reference zone to a more uniform level. That is very dangerous because it weakens thermal resilience and increases residents’ heat exposure. Research by Ali et al. (2019) supports this by showing that even water bodies in Makassar City can lower LST, yet the large-scale expansion of built-up areas exacerbates SUHII. In addition, Jumadi et al. (2025) emphasize that built-up density contributes to LST across the city.

The physical transformations documented in this study have direct consequences for human comfort and public health in Makassar City. Heat-retaining materials, such as concrete, asphalt, and corrugated iron roofs, that are dominated by LCZ 3 (Compact low-rise), actively contribute to thermal discomfort. High-density areas with a lack of vegetation trap solar radiation and re-emit long-wavelength heat, exposing residents to high levels of radiation (de Almeida et al., 2025; W. Zhang et al., 2024). That exacerbates physiological heat stress (da Silva Espinoza et al., 2023). It is because the body’s primary cooling mechanism, sweating, becomes less effective when surrounded by hot, dry, impervious surfaces that radiate heat back to the skin. Previous studies have confirmed that materials such as asphalt increase heat stress, while grass and high-albedo surfaces reduce it (Irmak et al., 2017; Tyagi & Danish, 2025). Consequently, thermal discomfort reduces quality of life, reduces outdoor activity, and increases reliance on high energy consumption for intensive cooling. Importantly, this burden falls unequally on vulnerable groups, including children, the elderly, and those with pre-existing health conditions (Aghazadeh et al., 2025; Fong et al., 2023).

Finally, the analysis of SUHII exposure at the district level indicates that districts in the urban center, such as Rappocini and Mamajang, experienced Very High (>50%) SUHII Hot Spot exposure in 2024. In contrast, peripheral districts such as Biringkanaya and Tamalanrea become sensitive areas. They were characterized by frequent shifts between Moderate (0-25%) and High exposure (25-50%) levels due to rapid industrialization and residential expansion. The patterns were confirmed to significantly increase LST (Yang, Bai, et al., 2023; Zou et al., 2024). Those effects have been shown to pose poor risks to urban residents and urban sustainability (G. Wang et al., 2025; Yang, Bai, et al., 2023). Although absolute SUHII values at fine spatial scales should be interpreted with caution due to unquantified uncertainty from the lack of in-situ LST validation, the macro-scale thermal pattern across the district provides a baseline for a mitigation framework for Makassar City. The main focus in urban districts is structural repairs and the implementation of cooling roofs as soon as possible. Meanwhile, peripheral areas should focus on protecting remaining natural infrastructure to prevent further expansion of hot spots.

A further methodological limitation relates to the temporal synchronization between LST and LCZ data. As shown in Figure 13, substantial land conversion occurred over the 2014 – 2024 period, yet LCZ maps were available only for 2014, 2019, and 2024. Consequently, annual LST data from the boundary years (2013, 2016, 2017, 2021, and 2021) were paired with LCZ maps that may not fully reflect actual LST conditions at the time of LST acquisition. For years when LST preceded LCZ maps (2013, 2017, and 2022), the assigned LCZ reflected conditions of greater urbanization than existed at that time, potentially leading to an overestimation of SUHII. Conversely, in years when the LST aligns with the LCZ map (2016 and 2021), urban development may have exceeded the map's coverage, potentially leading to an underestimation of SUHII. Future studies would benefit from LCZ mapping at a higher frequency, such as a biennial or triennial classification interval, to mitigate these morphological differences. While creating an annual LCZ is technically feasible, it is still constrained by the limited availability of high-resolution spatial data needed to parameterize the classification for a specific study area.