Multitemporal Analysis of Urban Heat Island Dynamics in Response to Land Use/Land Cover (LULC) Changes in Bukidnon Province, Philippines (2017–2024)

1. Introduction

Urbanization is a defining characteristic of the modern era, significantly altering landscapes worldwide. By 2050, over 68% of the global population is expected to reside in urban areas, necessitating expanded infrastructure to support human activity. This expansion often results in land cover changes, such as converting vegetation and agricultural lands into built-up environments, leading to increased land surface temperatures (LST) and intensifying the Urban Heat Island (UHI) effect [1,2,3,4,5]. The UHI phenomenon, wherein urban areas experience higher temperatures than surrounding rural regions, has broad implications, including increased energy consumption, air pollution, and adverse health outcomes [6,7,8,9,10,11,12,13,14,15,16].

While extensive research has addressed UHI effects globally, limited analysis exists on the dynamics of UHI in provincial urbanizing regions such as Bukidnon. Previous studies often focus on metropolitan centers, overlooking larger and rapidly developing areas where land cover transformation and agricultural expansion play a critical role [17,18]. Bukidnon is Northern Mindanao’s most populous province, with over 1.5 million residents (Philippine Statistics Authority, 2020). Its evolving landscape, driven by urban expansion, poses challenges to local climate patterns, necessitating an in-depth analysis of how land cover changes contribute to LST variations and UHI intensification [19,20,21,22,23]. Bukidnon’s dual identity as an agricultural hub and an emerging urban center positions it as a key area for studying the intersection of urbanization, land cover change, and climate variability. Its significant contributions to national agricultural output contrast with the accelerating pace of urban infrastructure development, emphasizing the need for localized UHI assessments [24,25].

Previous studies have highlighted the role of remote sensing technologies, such as Landsat and Sentinel satellites, in monitoring land cover and temperature dynamics [26,27,28,29,30,31,32,33,34,35]. However, challenges such as accessing high-resolution LST data and errors in index-based LST estimation persist [36,37,38,39]. Addressing these challenges is crucial for ensuring that urban development aligns with the Sustainable Development Goals (SDGs), particularly SDG 11, which focuses on creating sustainable and resilient cities (United Nations, 2019).

Given the rapid urban transformation and declining forest cover, understanding UHI dynamics in Bukidnon is crucial for formulating sustainable urban strategies. This study addresses these gaps by conducting a multitemporal analysis, examining LST trends, and offering insights to guide climate-resilient urban planning. Specifically, this study aims to conduct a multitemporal analysis of land cover change and its impact on UHI dynamics in Bukidnon from 2017 to 2024. By leveraging satellite imagery and Geographic Information Systems (GIS), this research seeks to fill critical knowledge gaps, offering insights into the relationship between land cover changes and UHI effects [13,15,20,23,40,41,42,43]. The findings will provide essential data to guide urban planners and policymakers in developing strategies for mitigating UHI impacts and enhancing climate resilience in the region [44,45,46,47].

2. Study Area

The study was conducted in the Province of Bukidnon, Philippines. It is situated on Mindanao Island, the southern island of the Philippines; with a population of over 1.5 million, Bukidnon is the most populated province in Northern Mindanao (Philippine Statistics Authority, 2020). In 2023, it led the country in agricultural and fisheries production, contributing approximately 6.8% to the national total, with an output valued at ₱155.28 billion (Philstar, 2024). This significant contribution highlighted Bukidnon’s role as a major producer of crops such as rice, corn, sugarcane, and various fruits and vegetables. The province’s agricultural prominence is further highlighted by its substantial land area dedicated to farming. As of 2022, Bukidnon’s gross value added (GVA) in agriculture, forestry, and fishing was ₱125.4 billion, accounting for 7% of the sector’s national GVA (PSA, 2023). This extensive agricultural activity has led to land use changes, including the conversion of forestlands into agricultural areas.

The province’s agricultural landscape is characterized by a mix of traditional farming and large-scale plantations. In 2005, about 16.4% of Bukidnon’s land was devoted to agriculture, with corn cultivation occupying 7.3% (approximately 66,400 hectares), making it the second-largest corn-producing province in the Philippines (Environmental Science for Social Change, 2005). This extensive agricultural activity has implications for land use and environmental sustainability.

The conversion of land for agricultural purposes has raised concerns about deforestation and its environmental impacts. In 2010, Bukidnon had 448,000 hectares of natural forest, covering more than 50% of its land area. However, by 2018, the remaining natural forest was about 213,066 hectares, indicating significant deforestation and land conversion over the years (Global Forest Watch, 2018; Forest Foundation Philippines, 2024). This reduction in forest cover is attributed to the expansion of agricultural lands and other land use changes.

The province’s growing urban areas are likely experiencing significant land cover changes, contributing to the exacerbation of UHI effects. These changes can lead to heightened energy demands, increased air pollution, and adverse health outcomes. Understanding these impacts through precise, localized studies is essential for developing effective urban planning and mitigation strategies that align with sustainable development goals.

Figure 1. Satellite imagery showing the geographical boundary: (a) Province of Bukidnon, the study area, and (b) its location within the Philippines.

Figure 1. Satellite imagery showing the geographical boundary: (a) Province of Bukidnon, the study area, and (b) its location within the Philippines.

3. Materials and Methods

As shown in Figure 1, this study utilized Landsat and Sentinel imagery. Landsat 8, and Sentinel-2 and Sentinel-3 data were sourced from their respective platforms: USGS Earth Explorer (https://earthexplorer.usgs.gov/) and Copernicus Data Space (https://browser.dataspace.copernicus.eu/).

Land cover maps were generated using Landsat 8 and Sentinel-2 imagery. These datasets were mosaicked and corrected for the years 2017, 2020, and 2024. Meanwhile, Sentinel-3 data captured Land Surface Temperature (LST). Landsat imagery has a spatial resolution of 30m, while Sentinel-2 offers a 10m resolution. To ensure consistency across datasets, Sentinel-2 data were resampled to 30m to match Landsat imagery. Resampling Landsat to 10m was not performed, as it would not enhance the actual resolution, making the resampling of Sentinel-2 a more appropriate and reliable choice to maintain data integrity.

LST data were derived from Sentinel-3, which provides daily surface temperature measurements through the Sea and Land Surface Temperature Radiometer (SLSTR) sensor. This allowed for a continuous dataset throughout the year. To calculate the mean LST, daily LST data were acquired and averaged. Sentinel-3 images with distortions or sensor errors, such as negative temperatures, fragmented images, and no data images, were excluded to ensure accuracy and reliability. On average, between 15 and 30 LST observations per month were included in the analysis. The data were reprojected to the appropriate coordinate reference system to ensure compatibility across datasets. Using GIS software, the monthly average LST values were computed by applying standard mean formulas, resulting in an accurate representation of temperature trends across the study period.

Figure 2. Methodological Framework of the Study outlining the workflow and key outcomes of the research.

Figure 2. Methodological Framework of the Study outlining the workflow and key outcomes of the research.

3.1. Land Cover Classification

Multispectral images from Landsat and Sentinel satellites were obtained via the USGS Earth Explorer and Copernicus Browser platforms. These images underwent preprocessing steps to ensure accuracy and consistency, including atmospheric and geometric corrections. The land cover classification was performed using advanced algorithms within Geographic Information System (GIS) software, categorizing the landscape into four primary classes: Cropland, Woodland, Water Body, and Built-up Land. This classification methodology was adapted from [48] which focused on similar studies in Vietnam. While the Water Body class was included in the classification, it was not analyzed in detail due to its minimal area and limited influence on land surface temperature [48]. The classified results were mapped and visualized using GIS tools, generating detailed land cover maps that highlight spatial distributions and temporal trends. The classification accuracy was assessed and validated using a confusion matrix, comparing ground truth data and classified images within the GIS software, ensuring the reliability of the findings.

3.2. Land Surface Temperature Retrieval

The land surface temperature (LST) data were acquired from the Copernicus Browser, specifically utilizing Sentinel-3’s Sea and Land Surface Temperature Radiometer (SLSTR) sensor. A series of LST datasets were downloaded for the years 2017, 2020, and 2024. Approximately 30 LST images per month were available for download; however, not all could be included in the analysis due to distortions and potential sensor errors.

An example of such issues includes occurrences of negative temperature values (e.g., below 0 °C) within the Province of Bukidnon. As the Philippines is a tropical country, such negative temperature readings are unrealistic and were excluded from the analysis. On average, a minimum of 15 and a maximum of 30 valid LST datasets per month were processed.

To ensure compatibility and consistency, the LST data were reoriented by reprojecting them to the appropriate coordinate reference system. Using the GIS calculator, monthly average LST values were computed by applying the basic mean formula, resulting in a more refined and accurate representation of temperature trends.

3.3. Urban Heat Island Index Calculation

Urban Heat Island (UHI) is a phenomenon in which urban areas exhibit higher temperatures than their surrounding rural regions (NASA, 2024). Quantifying the effects of UHI can be challenging; however, using Land Surface Temperature (LST) data derived from Sentinel-3 can be simplified by calculating the difference between urban and rural temperatures. This difference, known as the Urban Heat Island Index (UHII), is straightforward for analyses involving two distinct areas. However, the process becomes significantly more complex when considering larger regions, such as provinces with diverse characteristics encompassing both urban and rural areas. To address this challenge, the researchers propose a technique for quantifying UHI values by adopting the formula introduced by [48] in their research and as presented by [49] and [50] expressed as follows:

UHII = (LST_i − LST_m)/SD,

(1)

where UHII is the Urban Heat Island (UHI) Index, LST_i and LST_m are the Land Surface Temperature (LST) of a particular area and the mean LST of the area, respectively, and SD is the computed Standard Deviation of LST.

Table 1. Urban Thermal Field Variance Threshold.

Urban Thermal Field Variance Index Threshold	Urban Heat Island Phenomenon	Ecological Conditions Evaluation
< 0	None	Excellent
0.000–0.005	Weak	Good
0.005–0.010	Middle	Normal
0.010–0.015	Strong	Bad
0.015–0.020	Stronger	Worse
>0.020	Strongest	Worst

Table 1. Urban Thermal Field Variance Threshold.

Urban Thermal Field Variance Index Threshold	Urban Heat Island Phenomenon	Ecological Conditions Evaluation
< 0	None	Excellent
0.000–0.005	Weak	Good
0.005–0.010	Middle	Normal
0.010–0.015	Strong	Bad
0.015–0.020	Stronger	Worse
>0.020	Strongest	Worst

3.3. Urban Thermal Field Variance Index Calculation

The Urban Thermal Field Variance Index (UTFVI) is a quantitative measure used to evaluate the impact of urban heat intensity on urban surfaces, calculated using the equation below [51]. While the Urban Heat Island Index (UHII) quantifies the temperature difference between distinct regions, it does not provide insight into how these temperature variations affect ecosystems. By calculating UTFVI, the ecological implications of urban heat intensity can be visualized, offering a clearer understanding of both ecological conditions and the Urban Heat Island (UHI) effect. The UTFVI is determined using the equation (2) adopted from [51], with the threshold levels categorized into six classes, as outlined in Table 1.

UTFVI = (LST_i − LST_m)/LST_i,

(2)

where UTFVI is the Urban Thermal Field Variance Threshold Index, LST_i and LST_m are the Land Surface Temperature (LST) of a particular area and the mean LST of the area, respectively.

4. Results

4.1. Land Cover Classification and Analysis

Figure 3. Land Cover Classes of the Province of Bukidnon for the year (a) 2017, (b) 2020, and (c) 2024.

Figure 3. Land Cover Classes of the Province of Bukidnon for the year (a) 2017, (b) 2020, and (c) 2024.

A supervised classification technique was employed using the support vector machine algorithm to analyze data from the years 2017, 2020, and 2024. Four distinct land cover classes were identified: Cropland, Woodland, Water Body, and Built-up Land. The descriptions of these classes are provided in Table 2.

Table 2. Land cover types.

Classes	Description
Cropland	Lands used for agriculture, horticulture, and gardens, including paddy fields and irrigated and dry farmland and vegetation. Lands covered with wetland plants and water.
Woodland	Lands covered with trees, including deciduous and coniferous woodlands and sparse woodland.
Water Body	Water bodies in the land area, including rivers, lakes, reservoirs, and fish ponds.
Built-up Land	Lands modified by human activities, including all kinds of habitation, industrial and mining areas, transportation facilities, and interior urban green zones and water bodies.

Table 2. Land cover types.

Classes	Description
Cropland	Lands used for agriculture, horticulture, and gardens, including paddy fields and irrigated and dry farmland and vegetation. Lands covered with wetland plants and water.
Woodland	Lands covered with trees, including deciduous and coniferous woodlands and sparse woodland.
Water Body	Water bodies in the land area, including rivers, lakes, reservoirs, and fish ponds.
Built-up Land	Lands modified by human activities, including all kinds of habitation, industrial and mining areas, transportation facilities, and interior urban green zones and water bodies.

The accuracy assessment of the land cover classification for 2017, 2020, and 2024 was evaluated using a confusion matrix based on 1,000 randomly selected points compared with the classified map. The results demonstrated acceptable classification performance across all years. The User’s Accuracy for the Cropland class was consistently high, ranging from 0.810 in 2017 to 0.859 in 2020, with a slight decrease to 0.815 in 2024. The Woodland class also performed well, achieving User’s Accuracy values between 0.777 in 2020 and 0.840 in 2017. In contrast, the Built-up class showed variability, with its lowest accuracy recorded in 2020 at 0.619 but recovering to 0.720 in 2024. The Waterbody class displayed significant variation, achieving perfect accuracy (1.000) in 2017 but declining to 0.700 in 2020 and slightly improving to 0.714 in 2024. The Producer’s Accuracy remained relatively consistent across all years, with values close to 0.82, indicating a reliable performance in correctly identifying ground truth pixels. The Kappa index, which measures classification agreement while accounting for chance, exceeded 0.61 for all years, with values of 0.655 (2017), 0.648 (2020), and 0.647 (2024). According to the threshold established by [52], these Kappa values indicate a substantial level of agreement between the classified maps and reference data.

From 2017 to 2024, significant shifts in land cover were observed, marked by an increase in built-up areas and a reduction in vegetative cover, particularly woodland. Built-up land expanded from 21,256 hectares in 2017 to 32,099 hectares by 2024 (see Table 4), indicating rapid urbanization in Bukidnon. Woodland areas experienced a decline, shrinking from 392,463 hectares in 2017 to 367,737 hectares in 2024 (refer to Figure 3). This reduction suggests a substantial conversion of natural landscapes into urban and agricultural areas.

Table 3. Accuracy assessment using confusion matrix for land cover classification (2017, 2020, and 2024).

		2017	2020	2024
User’s Accuracy	Cropland	0.810084	0.858716	0.815466
	Built-up	0.769231	0.619048	0.720000
	Woodand	0.840426	0.777518	0.823383
	Wate Body	1.000000	0.700000	0.714286
Producer’s Accuracy		0.821535	0.817547	0.815553
Kappa		0.654973	0.648250	0.647020

Table 3. Accuracy assessment using confusion matrix for land cover classification (2017, 2020, and 2024).

		2017	2020	2024
User’s Accuracy	Cropland	0.810084	0.858716	0.815466
	Built-up	0.769231	0.619048	0.720000
	Woodand	0.840426	0.777518	0.823383
	Wate Body	1.000000	0.700000	0.714286
Producer’s Accuracy		0.821535	0.817547	0.815553
Kappa		0.654973	0.648250	0.647020

Table 4. Land cover classifications and changes from 2017 to 2024 in hectares.

Land cover class	2017	2020	2024
Built-up Land	21,256	26,321	32,099
Woodland	392,463	381,339	367,737
Cropland	513,212	521,200	532,008
Waterbody	3,111	5,315	3,586

Table 4. Land cover classifications and changes from 2017 to 2024 in hectares.

Land cover class	2017	2020	2024
Built-up Land	21,256	26,321	32,099
Woodland	392,463	381,339	367,737
Cropland	513,212	521,200	532,008
Waterbody	3,111	5,315	3,586

4.2. Land Surface Temperature Analysis

Figure 4 illustrates the spatiotemporal variations in Land Surface Temperature (LST) across the study area for the years 2017 (a), 2020 (b), and 2024 (c). This multitemporal analysis highlights shifts in surface temperature, reflecting the influence of land cover changes and urban development over the period.

In 2017, LST values ranged from 15.9783 °C to 34.4996 °C, with higher temperatures concentrated in the southern and central regions. Areas with lower LST values, depicted in green, are primarily distributed in the northern and mid-central portions, likely corresponding to woodland or heavily vegetated and less urbanized zones.

By 2020, the LST range increased to 17.5329 °C–36.9591 °C, indicating a notable rise in land surface temperatures. This trend suggests intensified urban heat island (UHI) effects or reduced vegetation cover, with expanded red and orange zones across the map.

In contrast, the year 2024 reveals a reduction in maximum LST to 31.0354 °C, while the minimum LST decreases to 12.9585 °C. Despite the lower overall LST values, the spatial distribution of high-temperature zones remains widespread, highlighting persistent thermal hotspots likely driven by ongoing land use dynamics.

Figure 4. Land Surface Temperature (LST) of the Province of Bukidnon for the year (a) 2017, (b) 2020, and (c) 2024.

Figure 4. Land Surface Temperature (LST) of the Province of Bukidnon for the year (a) 2017, (b) 2020, and (c) 2024.

4.3. Urban Heat Island Index and Urban Thermal Field Variance Index Analysis

Figure 5. Computed Urban Heat Island (UHI) Index of the Province of Bukidnon for the year (a) 2017, (b) 2020, and (c) 2024.

Figure 5. Computed Urban Heat Island (UHI) Index of the Province of Bukidnon for the year (a) 2017, (b) 2020, and (c) 2024.

Figure 5 presents a comparative analysis of the Urban Heat Island Index (UHII) across three temporal maps: 2017 (a), 2020 (b), and 2024 (c). The UHII is spatially visualized, illustrating variations in surface temperature intensities across the study area.

In 2017, higher UHII values, indicated by red to orange hues, were predominantly concentrated in the central and southern regions, signifying elevated land surface temperatures. Conversely, the northern and select mid-central areas exhibit cooler UHII indices, marked by green shades. By 2020, there is a noticeable reduction in the extent of high UHII zones, suggesting potential improvements in vegetation cover or land-use changes. However, the 2024 projection reveals a resurgence of elevated UHII values, indicating a possible increase in urbanization or reduction in vegetative cover.

The UHII values reflect a subtle fluctuation over time. In 2017, the highest recorded UHII value was 2.16133, which slightly decreased to 2.00838 in 2020, before rising to 2.06475 in 2024. Conversely, the lowest UHII demonstrates a gradual increase from -3.53954 in 2017 to -2.9478 in 2020, with a slight decrease to -3.11711 by 2024.

Figure 6. Zonal distribution of Urban Thermal Field Variance (UTFVI) Index values across Barangays of the Province of Bukidnon for the year (a) 2017, (b) 2020, and (c) 2024.

Figure 6. Zonal distribution of Urban Thermal Field Variance (UTFVI) Index values across Barangays of the Province of Bukidnon for the year (a) 2017, (b) 2020, and (c) 2024.

Figure 6 illustrates the zonal distribution of the Urban Thermal Field Variance Index (UTFVI) across barangays for the years 2017 (a), 2020 (b), and 2024 (c). This visualization provides insight into the spatial and temporal shifts in ecological conditions and the intensity of the Urban Heat Island (UHI) phenomenon.

In 2017, a significant portion of the study area exhibited “Excellent” ecological conditions, represented by dark blue, indicating no UHI phenomenon. However, localized areas of weak to strong UHI intensity (depicted in purple to orange) appear along the periphery, and selected barangays corresponding to urbanized zones.

By 2020, there is a noticeable expansion of zones classified as “Normal” to “Bad” ecological conditions (yellow and light orange), reflecting an increase in middle to strong UHI phenomena. This trend suggests urban expansion and the corresponding rise in surface temperatures. The presence of “Worst” ecological conditions (represented by yellow) becomes more pronounced, highlighting critical areas experiencing the most severe thermal stress.

The 2024 illustration reveals a mixed pattern. While some regions show improvements, returning to “Excellent” conditions, areas with “Worst” ecological conditions persist in densely populated or highly urbanized sectors. The spread of “Weak” and “Middle” UHI zones suggests ongoing urbanization and land use change, necessitating targeted mitigation efforts. The shift in land surface temperatures observed in 2024, which subsequently affects the UTFVI index, indicates the presence of additional influencing factors. One notable factor could be the persistence of tropical cyclones in the Philippines during this period, leading to increased rainfall in the study area. This heightened precipitation likely enhances the cooling effect, contributing to the observed changes in surface temperature and ecological conditions.

4.3. Correlation between LandCover Classes and Land Surface Temperature (LST)

Pearson correlation analysis was employed to examine the relationship between Land Surface Temperature (LST) and land cover classes, incorporating population and elevation as additional variables. The Province of Bukidnon comprises a total of 464 barangays, providing 464 data points for each variable in the study. Population data were sourced from the Philippine Statistics Authority (PSA). Since PSA data were available only for 2015 and 2020, linear interpolation was performed at the barangay level to estimate population values for 2017 and 2024. Elevation data were derived from the Philippine Geoportal (https://www.geoportal.gov.ph/) and processed using zonal averaging techniques. For land cover classification, spatial normalization was conducted to express the data as percentages. Tables 5 to 7 present the correlation matrices for the years 2017 through 2024, respectively.

The correlation matrices for the years 2017, 2020, and 2024 reveal significant relationships between Land Surface Temperature (LST) and various variables. Elevation consistently shows a strong negative correlation with LST across all years (-0.817 in 2017, -0.775 in 2020, and -0.761 in 2024), indicating that higher elevations are associated with lower temperatures. Population exhibits a weak positive correlation with LST (0.165 in 2017, 0.152 in 2020, and 0.190 in 2024), suggesting a minor influence of population density on temperature. Built-up areas demonstrate a moderate positive correlation with LST, ranging from 0.197 to 0.283, highlighting the impact of urbanization on temperature increase. Cropland consistently shows the strongest positive correlation with LST among all land cover types, peaking in 2020 (0.667) compared to 2017 (0.521) and slightly declining in 2024 (0.551), highlighting the role of agricultural activities in temperature dynamics. Water bodies display a weaker positive correlation with LST, fluctuating across the years (0.233 in 2017, 0.109 in 2020, and 0.167 in 2024), reflecting a lesser influence on temperature patterns. While woodland demonstrates a strong negative correlation with LST, reflecting its cooling effect. This correlation strengthens from 2017 (-0.683) to 2020 (-0.794) before slightly weakening in 2024 (-0.697).

The results consistently reveal the variables that are highly correlated, with a notable observation that cropland and woodland consistently exhibit a strong negative relationship. This is evidenced by correlation values of -0.864 in 2017, -0.880 in 2020, and -0.798 in 2024. These findings indicate that areas with a higher proportion of cropland tend to have significantly less woodland, suggesting an ongoing land conversion process where cropland expansion occurs at the expense of woodland areas.

Table 5. Correlation Matrix for Year 2017.

	LST Average	Elevation	Population	Built-Up	Cropland	Water Body
Elevation	-0.817
Population	0.165	-0.008
Built-Up	0.283	-0.078	0.120
Cropland	0.521	-0.454	-0.006	-0.288
Water Body	0.233	-0.322	0.119	0.014	0.136
Woodland	-0.683	0.514	-0.063	-0.232	-0.864	-0.198

Table 5. Correlation Matrix for Year 2017.

	LST Average	Elevation	Population	Built-Up	Cropland	Water Body
Elevation	-0.817
Population	0.165	-0.008
Built-Up	0.283	-0.078	0.120
Cropland	0.521	-0.454	-0.006	-0.288
Water Body	0.233	-0.322	0.119	0.014	0.136
Woodland	-0.683	0.514	-0.063	-0.232	-0.864	-0.198

Table 6. Correlation Matrix for Year 2020.

	LST Average	Elevation	Population	Built-up	Cropland	Water Body
Elevation	-0.775
Population	0.152	0.018
Built-up	0.197	-0.101	0.109
Cropland	0.667	-0.498	0.060	-0.320
Water Body	0.109	-0.217	0.035	0.141	-0.012
Woodland	-0.794	0.582	-0.117	-0.163	-0.880	-0.134

Table 6. Correlation Matrix for Year 2020.

	LST Average	Elevation	Population	Built-up	Cropland	Water Body
Elevation	-0.775
Population	0.152	0.018
Built-up	0.197	-0.101	0.109
Cropland	0.667	-0.498	0.060	-0.320
Water Body	0.109	-0.217	0.035	0.141	-0.012
Woodland	-0.794	0.582	-0.117	-0.163	-0.880	-0.134

Table 7. Correlation Matrix for Year 2024.

	LST Average	Elevation	Population	Built-up	Cropland	Waterbody
Elevation	-0.761
Population	0.190	-0.009
Built-up	0.243	-0.083	0.165
Cropland	0.551	-0.433	0.049	-0.286
Waterbody	0.167	-0.223	0.113	0.167	-0.018
Woodland	-0.697	0.486	-0.156	-0.346	-0.798	-0.140

Table 7. Correlation Matrix for Year 2024.

	LST Average	Elevation	Population	Built-up	Cropland	Waterbody
Elevation	-0.761
Population	0.190	-0.009
Built-up	0.243	-0.083	0.165
Cropland	0.551	-0.433	0.049	-0.286
Waterbody	0.167	-0.223	0.113	0.167	-0.018
Woodland	-0.697	0.486	-0.156	-0.346	-0.798	-0.140

5. Discussions

The results of the study reveal substantial land cover changes in Bukidnon Province from 2017 to 2024, driven by urban expansion and agricultural development. The increase in built-up areas and the corresponding decline in woodland highlights the extent of urbanization, which is a significant contributing factor to the intensification of the Urban Heat Island (UHI) effect observed in the area. This aligns with global trends indicating that urban growth often significantly alters local land cover, leading to elevated land surface temperatures (LST) [6,7,8].

Land cover classification revealed that built-up land expanded by nearly 10,000 hectares between 2017 and 2024, while woodland areas experienced a considerable reduction. These shifts indicate a progressive conversion of natural and agricultural lands into urban environments. The reduction of vegetative cover, which serves as a natural cooling mechanism, has been directly linked to increased LST across Bukidnon. This pattern is consistent with existing literature demonstrating that replacing vegetation with impervious surfaces leads to higher heat retention and exacerbates UHI effects [11,13].

The analysis of LST data shows a clear upward trend in surface temperatures, particularly in urbanized and deforested areas. From 2017 to 2020, LST values increased notably, reflecting the peak of urban expansion during this period. Interestingly, by 2024, there was a slight decrease in maximum LST, which may be attributed to increased rainfall. Despite this decline, the spatial distribution of high-temperature zones remained widespread, indicating persistent thermal hotspots. When comparing the Land Surface Temperature (LST) across different land cover types, it is evident that woodland areas consistently maintain lower LST values. In contrast, croplands and built-up areas consistently exhibit higher LST, reflecting the influence of vegetation density and urbanization on surface temperature dynamics. The decline in observed LST in 2024 may be attributed to persistent rainfall associated with tropical cyclones. This is supported by PAGASA records, which indicate that 2024 experienced lower average atmospheric temperatures. The increased rainfall and cooler atmospheric conditions likely contributed to the reduction in surface temperatures. However, the trend in LST and corresponding land cover remains consistent across 2017, 2020, and 2024, with woodland areas consistently exhibiting lower temperatures than other land cover types. Conversely, the water bodies cannot be clearly analyzed in relation to LST due to their minimal coverage, making them nearly indistinguishable in the maps.

The Urban Heat Island Index (UHII) analysis further highlights the relationship between land cover changes and surface temperature dynamics. High UHII values were consistently recorded in urban centers and the cropland area, while cooler indices were observed in woodland and rural areas. This variation highlights the critical role of vegetation in mitigating the effects of urban heating [7,35]. The resurgence of elevated UHII values in 2024 suggests that urban expansion continues to exert pressure on local temperature patterns. The Urban Heat Island Intensity (UHII) follows a similar trend to LST, with higher LST values correlating with increased UHII. These elevated temperatures are consistently observed in cropland and built-up areas. The Urban Thermal Field Variance Index (UTFVI) provides additional insights into the ecological implications of rising temperatures. The expansion of areas labeled as having “Bad” and “Worst” ecological conditions suggests that the intensifying Urban Heat Island (UHI) effect contributes to environmental degradation [51]. This raises serious concerns about how these changes could affect biodiversity, agricultural productivity, and public health.

The Pearson correlation analysis highlights the complex interplay between Land Surface Temperature (LST), land cover classes, elevation, and population in Bukidnon Province. The findings offer valuable insights into how these variables influence thermal environments and underscore the implications of land use changes on local climate dynamics. The consistent strong negative correlation between elevation and LST across all years (-0.817 in 2017, -0.775 in 2020, and -0.761 in 2024) underscores the well-documented relationship where higher altitudes experience lower temperatures. This trend is expected given the adiabatic cooling effect of elevation, emphasizing the role of Bukidnon’s mountainous terrain in moderating LST. The weak but positive correlation between population and LST (0.165 in 2017, 0.152 in 2020, and 0.190 in 2024) suggests that population density alone may not be a primary driver of LST changes in the province but could indirectly contribute through urbanization and land use modifications. Built-up areas demonstrate a moderate positive correlation with LST, which ranges from 0.197 to 0.283 over the study years, reflecting the contribution of urbanization to the Urban Heat Island (UHI) effect. This result aligns with global findings where impervious surfaces such as concrete and asphalt exacerbate temperature increases. Similarly, cropland exhibits the strongest positive correlation with LST among land cover types, peaking in 2020 (0.667) and slightly declining in 2024 (0.551). This trend highlights the influence of agricultural activities on LST, possibly due to reduced vegetation cover and changes in soil moisture levels associated with farming practices.

Woodland consistently exhibits a strong negative correlation with LST (-0.683 in 2017, -0.794 in 2020, and -0.697 in 2024), reflecting its cooling effect. This finding underscores the critical role of forest cover in mitigating LST through evapotranspiration, shading, and other ecosystem services. The inverse relationship between cropland and woodland (-0.864 in 2017, -0.880 in 2020, and -0.798 in 2024) is particularly striking. It reveals an ongoing land conversion process where cropland expansion occurs at the expense of forested areas. This pattern contributes to higher LST and raises concerns about biodiversity loss, carbon sequestration, and watershed stability.

This study underscores the need for significant urban planning and the development of green infrastructure to help combat UHI effects. Increasing green spaces in cities, preserving existing forests, and encouraging sustainable farming practices can minimize the negative impacts of urbanization on local climates. Looking ahead, research should aim to improve predictive models that take into account the socioeconomic factors driving changes in land use and their long-term effects on urban heat dynamics [37,38], particularly in Bukidnon and similar areas.

6. Conclusions

This study highlights the significant impact of land cover changes on Land Surface Temperature (LST) dynamics in the Province of Bukidnon from 2017 to 2024, driven primarily by urban expansion and agricultural development. The findings reveal that built-up areas and cropland consistently exhibit higher LST values, while woodland demonstrates a strong cooling effect, emphasizing its critical role in regulating local temperatures. The consistent negative correlation between elevation and LST underscores the influence of Bukidnon’s topography in moderating surface temperatures. In contrast, the weak positive correlation with population suggests that urbanization, rather than population density itself, is a key driver of LST increases.

The reduction in woodland and cropland expansion highlights an ongoing land conversion process, which contributes to elevated LST and raises concerns about biodiversity loss, carbon sequestration, and watershed stability. Furthermore, the intensification of the Urban Heat Island (UHI) effect, as evidenced by the Urban Heat Island Index (UHII) and Urban Thermal Field Variance Index (UTFVI), suggests that urbanization and agricultural practices are placing significant pressure on the local climate and ecological balance.

These findings underline the urgent need for sustainable land management practices, such as preserving forested areas, expanding urban green spaces, and adopting sustainable agricultural techniques. Addressing these challenges requires proactive urban planning and policies to mitigate the adverse effects of land cover changes on local climates. Future research should focus on refining predictive models incorporating socioeconomic drivers of land use change and their implications for urban heat dynamics. This will be essential in developing effective strategies to balance urban and agricultural development with environmental sustainability in Bukidnon and similar regions.