Land Use/Land Cover Changes and the Relationship with Land Surface Temperature Using Landsat and MODIS Imageries in Cameron Highlands, Malaysia

Mountainous regions are more sensitive to climatic condition changes and are susceptible to recent increases in temperature. Due to urbanization and land use/land cover (LULC) issues, Cameron Highlands has been impacted by rising land surface temperature (LST) variation. Thus, this study was carried out to explore the impact of the LULC change on LST in the Cameron Highlands from 2009 to 2019 using remote sensing images acquired from Landsat 7 ETM+, Landsat 8 Operational Land Imager (OLI/TIRS), and Moderate Resolution Imaging Spectroradiometer (MODIS) 11A Thermal sensors. A split-window algorithm was applied to Landsat 8 images (2013–2019) to derive the LST. Air temperature data of the study area were also obtained to cross-validate data sources. Based on the validation results, the accuracy of LULC and LST outputs were more than 94.6% and 80.0%, respectively. The results show that the current trend of urban growth continues at a rate of 0.16% per year, and the area experienced an LST increase of 2 ◦C between 2009 and 2019. This study is crucial for land planners and environmentalists to understand the impacts of LULC change on LST and to propose appropriate policy measures to control development in Cameron Highlands.


Introduction
Little is known about deforestation and land cover dynamics and their effects on land surface temperature (LST), particularly in highland areas. Cameron Highlands is situated in Malaysia at an elevation between 300 and 1800 m. The average temperature recorded at the main town, namely Tanah Rata, is 18 • C. The daily minimum temperature is 15 • C, with a deviation ranging from 5 to 7 • C. The temperature has steadily increased since the mid-1970s [1]. The warmest year was 1998 and 2002 was the second warmest. Recent deforestation has possibly contributed to the warming and increased LST. Previous studies using satellite observations indicate that tropical deforestation results in warmer and drier conditions at the local scale. LST is related to the surface energy balance and the integrated thermal conditions of the atmosphere within the planetary boundary layer [2]. The rise in LST could be a direct result of deforestation due to the expansion of urban zones and the overall global surface temperature rise caused by the reduction of the ozone layer.
Recent global warming events have impacted this microclimate and its habitats. However, the magnitude and spatial scale of this impact are unknown; hence, further study is needed to observe and address these issues. Despite their relatively small surface area, mountains are well understood to

Satellite (Landsat) images from 2009 to 2019 (Enhanced Thematic Mapper+ (ETM+) and
Operational Land Imager/Thermal Infrared Sensor (OLI/TIRS)) were collected from the United States Geological Survey (USGS) website. Landsat imagery was chosen because it is readily available and has long-term historical coverage. Dates in the first quarter of the year were chosen to maintain a constant seasonal condition, thus avoiding biases, which is essential for any change detection analysis. Images were also screened for minimum cloud and haze coverage. Hence, if the imagery in the first quarter did not satisfy these conditions, the second quarter for that particular year was chosen.
Landsat 7 was used for 2009 to 2012, whereas Landsat 8 was used for 2013 to 2019. MODIS imagery was used for 2009 to 2019 (Table 1).

Remote Sensing Image Acquisition
Satellite (Landsat) images from 2009 to 2019 (Enhanced Thematic Mapper+ (ETM+) and Operational Land Imager/Thermal Infrared Sensor (OLI/TIRS)) were collected from the United States Geological Survey (USGS) website. Landsat imagery was chosen because it is readily available and has long-term historical coverage. Dates in the first quarter of the year were chosen to maintain a constant seasonal condition, thus avoiding biases, which is essential for any change detection analysis. Images were also screened for minimum cloud and haze coverage. Hence, if the imagery in the first quarter did not satisfy these conditions, the second quarter for that particular year was chosen.
Landsat 7 was used for 2009 to 2012, whereas Landsat 8 was used for 2013 to 2019. MODIS imagery was used for 2009 to 2019 (Table 1).

Land Use/Land Cover Change Mapping
Prior to image classification, the presence of clouds was completely masked while haze was successfully removed using the MODTRAN codes available in the ATCOR module in ERDAS Imagine. An extra step was included for Landsat 7 images due to the error caused by the Scan Line Corrector (SLC) failure on board the satellite; the image chosen had to be composited with an image acquired from the previous week using ENVI version 5.3 [28]. Image selection was made based on similarities in the air temperature readings and climatic conditions. To ensure that the landcover was correctly assigned, Google Earth time series images are used as a cross-reference. From ground verification, it was concluded that the main five land cover classes ( Table 2) were classified. Numerous details and various land cover classes, primarily vegetation areas of specific varieties, were omitted due to the limitations of the Landsat 30 m resolution data. Table 2. Description of land use/cover categories.

Primary forests
Forests of native tree species, where there are no clearly visible indications of human activities and the ecological processes are not significantly disturbed.
Farm and mixed vegetable Comprises of agriculture land for orchards as well as vegetative crops.

Urban areas
Comprises an area such as settlement areas, roads, and industries.

Cleared lands
Comprises various types of old-growth forests that have become so severely degraded by logging and no longer resemble the spectral signatures of forests.
Water bodies Water features, including, rivers, lakes, and artificial ponds.

Classification of Images for LULCC
This study used the Object-Based Image Assessment (OBIA) method to classify the Land Use/Land Cover Change (LULCC). Trimble eCognition Developer Software was used for this purpose. The OBIA classification method used texture, shape, position, digital number, and a series of band combinations to identify the various details on an image. The entire work process was conducted similar to the workflow of Wang Xin et al. [29]. Different band combinations led to different segmentation patterns due to the type of wavelength that corresponds to the digital image number; hence, numerous combinations Land 2020, 9, 372 5 of 23 were tested. Based on visual interpretation, band combinations of 6,5,4 and 5,4,3 for Landsat 8 and, band combinations 5,4,3 and 4,3,2 for Landsat 7 provided the best discrimination ability between urban and vegetated areas due to the presence of NIR bands that were best to highlight vegetation or forested areas (see Figure 2) [30]. These band combinations were chosen due to the nature of the data image, whereby forest cover exceeds urban area. Training of the software was performed by assigning training areas from the segmented shapefiles onto associated land cover classes. Since the OBIA-produced image required additional checking, it was subjected to a toolbox program within ArcMap, explicitly designed for assessing the accuracy of OBIA images, namely RSOBIA by Oceanwise [31].
Land 2020, 9, x FOR PEER REVIEW 5 of 28 checking, it was subjected to a toolbox program within ArcMap, explicitly designed for assessing the accuracy of OBIA images, namely RSOBIA by Oceanwise [31]. To ground-truth the images, GPS points were obtained from a previous visit to the field site. The points were used for the 2019 images, for which Google Earth was used as additional verification. However, as points were not available for 2009 and 2014 images, only Google Earth was used for this purpose. ArcMap was then used for the final stages of map making and labelling. A summary of the method is shown in Figure 3.

Land Surface Temperature (LST) Calculation
In this study, we chose both Landsat thermal bands of Landsat 7 ETM+; namely the low gains band 6 (6L) and high gains band 6 (6H). Band 6L detects areas where surface brightness is high (nonvegetated areas), whereas band 6H detects areas where surface brightness is low (vegetated areas) [32]. A study by Nguemhe Fils et al. [33] had used only band 6H as their project was based on a high vegetation area where surface brightness was lower. However, the decision to use both bands for this project was applied, with both band values added and divided to get an average. This is done, due to the study area comprising of both vegetative and non-vegetative details. Similarly, this was done by Tavares et al. [32], Hesslerová et al. [34], Kumari et al. [35], and Walawender et al. [36].  can be distinguished by the slightly reddish tone with a mixture of light green) and (C) (cleared lands seen with a bare earth color, white-gray amidst primary forest around). Band combination 6,5,4 shown in image (D) (cleared lands has a light brown color in a lighter green tone of farms. Cleared lands are also identified as a lighter brown-white color tone when it has been cleared for weeks and left under the sun), (E) (urban areas in a purple tone surrounded by farms seen as light pink), (F) (water bodies-Ulu Jelai hydroelectric plant. Urban area towards the east of the dam in a light brown-purple color tone. The area is surrounded by primary forests seen as a darker shade of green; similarly, in (D)).
To ground-truth the images, GPS points were obtained from a previous visit to the field site. The points were used for the 2019 images, for which Google Earth was used as additional verification. However, as points were not available for 2009 and 2014 images, only Google Earth was used for this purpose. ArcMap was then used for the final stages of map making and labelling. A summary of the method is shown in Figure 3.

Land Surface Temperature (LST) Calculation
In this study, we chose both Landsat thermal bands of Landsat 7 ETM+; namely the low gains band 6 (6L) and high gains band 6 (6H). Band 6L detects areas where surface brightness is high (non-vegetated areas), whereas band 6H detects areas where surface brightness is low (vegetated areas) [32]. A study by Nguemhe Fils et al. [33] had used only band 6H as their project was based on a high vegetation area where surface brightness was lower. However, the decision to use both bands for this project was applied, with both band values added and divided to get an average. This is done, For Landsat 8 LST calculations, the chosen method is the Split-Window Algorithm (SWA), whereby both bands 10 and 11 are applied; the split-window algorithm was recommended over other methods [1,2,12,28,37,38]. Despite band 11 suffering from uncertainty of stray light, researchers demonstrated that the accuracy produced from the SWA is better than the use of single channel [20,39]. As for MODIS, daily LST's were acquired from the NASA Earth Data website. Data type TERRA 11A (Land Surface Temperature/Emissivity Daily L3 Global 1km) was chosen with the acquired data closely related to the dates of Landsat imagery. Then MODIS images were collected for each year of the study, subsequently processed in ArcMap software and converted into °C unit from K unit.

Conversion of Digital Numbers to Top-of-Atmosphere (TOA) Spectral Radiance
Digital numbers (DNs) store thermal data in the Landsat sensor and provide a manner of representing pixels that have not yet been calibrated and converted into radiance units. To calculate the LST, first, this process must be followed using an equation (for Landsat 8 Thermal Infrared Sensor (TIRS)) to convert DN to radiance: where is the TOA spectral radiance in units W/(m2.sr. m), is the band-specific multiplying rescale factor in the metadata, is the band-specific additive rescaling factor in the metadata, is the calibrated and quantized pixel values (DN), and represents the offsets supplied by USGS for the TIRS band calibration [15,34,40]. To convert the DN to radiance for Landsat 7 ETM+, the "gain and bias" method was used. The equations of this method are shown below [11,40]:  For Landsat 8 LST calculations, the chosen method is the Split-Window Algorithm (SWA), whereby both bands 10 and 11 are applied; the split-window algorithm was recommended over other methods [1,2,12,28,37,38]. Despite band 11 suffering from uncertainty of stray light, researchers demonstrated that the accuracy produced from the SWA is better than the use of single channel [20,39]. As for MODIS, daily LST's were acquired from the NASA Earth Data website. Data type TERRA 11A (Land Surface Temperature/Emissivity Daily L3 Global 1km) was chosen with the acquired data closely related to the dates of Landsat imagery. Then MODIS images were collected for each year of the study, subsequently processed in ArcMap software and converted into • C unit from K unit.

Conversion of Digital Numbers to Top-of-Atmosphere (TOA) Spectral Radiance
Digital numbers (DNs) store thermal data in the Landsat sensor and provide a manner of representing pixels that have not yet been calibrated and converted into radiance units. To calculate the LST, first, this process must be followed using an equation (for Landsat 8 Thermal Infrared Sensor (TIRS)) to convert DN to radiance: where L λ is the TOA spectral radiance in units W/(m2.sr. µm), M L is the band-specific multiplying rescale factor in the metadata, A L is the band-specific additive rescaling factor in the metadata, Q cal is the calibrated and quantized pixel values (DN), and O i represents the offsets supplied by USGS for the TIRS band calibration [15,34,40]. To convert the DN to radiance for Landsat 7 ETM+, the "gain and bias" method was used. The equations of this method are shown below [11,40]: where L λ is the spectral radiance at the sensor's aperture in W/(m2.sr. µm), gain is the rescaled gain of the "Level 1 product" header in W/(m2.sr. µm), bias is the rescaled bias in W/(m2.sr. µm), and QCAL is Land 2020, 9, 372 7 of 23 the quantized calibrated pixel value in DN. LMAX λ is the spectral radiance that is scaled to QCALMAX in W/(m2.sr. µm), QCALMIN is the minimum quantized calibrated pixel value, and QCALMAX is the maximum quantized calibrated pixel value that corresponds to LMAX λ in DN = 255 [11,40].

Conversion of Radiance to Brightness Temperature
After the digital numbers (DNs) are converted to reflection, the thermal band's spectral radiance must be converted to brightness temperature (BT). This is commonly accomplished using the thermal constants provided in the metadata file. The following equation is used in the tool's algorithm to convert reflectance to BT: where K 1 and K 2 are constant band-specific thermal conversions from the metadata. However, this result must be converted to Celsius. The radiant temperature is revised by adding absolute zero (approx. −273.15 • C).

Atmospheric Transmittance
Because water vapor content affects imagery in the form of clouds and haze, an estimation of Atmospheric Transmittance (AT) is needed to nullify that uncertainty. It is defined as an atmospheric capacity in which electromagnetic energy is transmitted through means of air by solar radiation, water vapor, and dust. As it fluctuates with various wavelengths through the various layers of the atmosphere, there is a relationship in which, for each small light scattering and atmospheric absorption, there is a higher air transmittance. To entirely derive AT, the calculation of water vapor content and mean atmospheric temperature must be previously calculated. In addition, estimations of relative humidity (RH) and the temperature of the upper air are needed. The following equations outline the process to obtain the equation required to calculate AT [18,41,42]: where e w (p,t) is the saturation vapor pressure of moist air; f (p) is the pressure fraction; e w (t) is the saturation vapor pressure in the pure phase; t is the ambient temperature; p is the atmospheric pressure; W is the atmospheric vapor pressure content (in g.cm −2 ); and P a is the standard atmospheric temperature of 1013.25 hPa. Hence, with data obtained from Met Malaysia, such as atmospheric pressure, air ambient temperature at a daily average, and relative humidity (RH), we were able to compute the necessary calculations with in situ data to yield an accurate estimation for the study site.

Emissivity Estimation Using the NDVI Technique
Landsat red and near-infrared bands were used for calculating the Normal Difference Vegetation Index (NDVI). The basic premise of calculating the NDVI is essential because the amount of vegetation present is an essential factor, and the NDVI can be used to understand the general vegetation condition. The calculation of the NDVI is essential because, subsequently, the proportion of the vegetation (PV) Land 2020, 9, 372 8 of 23 should be premeditated, which is closely related to the NDVI. In addition, emissivity (ε) should be calculated, which is related to the PV: where NIR = near-infrared band; Red = red band. Landsat 8 data used Bands 5 (infrared) and 4 (red), whereas Landsat 7 used Bands 4 and 3, to calculate the NDVI values. From the NDVI values obtained, the next step in acquiring Land Surface Emissivity (LSE), namely the calculation of the proportion of vegetation (P v ), was performed.

Calculating the Proportion of Vegetation (P v )
The proportion of vegetation (P v ) is calculated according to the equation below:

Land Surface Emissivity (LSE) Assessment
Land surface emissivity (LSE (ε)) is essential to estimate LST, because LSE is a proportionality factor that scales blackbody radiance (Planck's law) to predict emitted radiance and represents the efficiency of transmitting thermal energy across the surface into the atmosphere. Using a formula used by Ogunode and Akombelwa [18], Dutta [43] and Jimenez-Munoz et al. [44], LSE can be calculated using: where ε V and ε s are the vegetation and soil emissivity, respectively. The emissivity constants used for Band 10 are ε s = 0.971 and ε V = 0.987, and for Band 11 are ε s = 0.977 and ε V = 0.989 [45].

LSE Mean and Difference
To retrieve LST, the mean and difference of LSE must be obtained because they vary at varying distances. This allows or an estimation of LSE to utilize the NDVI values for the retrieval of soil and vegetation proportions-that is, minimum and maximum NDVI values. Both the LSE mean and difference can show a relationship between the TIR bands, as shown in the formula below:

Derivation of Atmospheric Vapor Pressure Content (W)
The obtained value of W from Equation (8) is used in the SW algorithm shown in Equation (14) for the LST retrieval process.

Coefficient
Value

Land Cover Change Analysis
The highest class recorded for 2009 was primary forests (Figure 4) (578.15 km 2 , 86.33% of total land area) followed by urban areas (50.10 km 2 , 7.48%) and farm and mixed vegetables (33.00 km 2 , 4.93%). The two minor land cover classes, namely cleared lands (7.84 km 2 ) and water bodies (0.60 km 2 ), were also recorded. The farm and mixed vegetables class includes tea plantations and other vegetation crops, such as corn, grapes, flowers, tomato, and strawberry farms. Within the farm and mixed vegetables land cover type, "other vegetation", which comprises uses that are not individually categorized, account for only 1/10 of the total area. This is due to the lack of growth in area of these crop types in Cameron Highlands. Due to the limitations of Landsat 30 m resolution, these finer details are not able to be classified; hence, the difference can only be seen when higher resolution data are used.   In 2014, image classification results ( Figure 5 and Table 4) show that the highest class was also primary forest (541.76 km 2 , accounting for 80.90% of the total land area and representing a decline of 36.39 km 2 compared to 2009), followed by urban areas (59.50 km 2 , 8.88%), and farm and mixed vegetables (57.03 km 2 , 8.52%). The remaining land cover classes were cleared lands (11.06 km 2 ) and water bodies (0.34 km 2 ). Large changes across the terrain and hilly areas of Cameron Highlands can be seen within a five-year period. Due to the continuous rise in both urban areas (9.40 km 2 ) and farmland (24.03 km 2 ), the decline of primary forests continues. Massive deforestation has occurred in the study area based on image interpretation and classification.  Illegal land clearing in small communities also contributes to the overall problem of deforestation. Although this may represent a small change, it is nonetheless an issue for the Land and Planning authorities of Cameron Highlands. This is because the areas of these activities are deep within the rainforest. Moreover, such uncontrolled activities might also be the result of landowners breaching the agreement of forest clearing, which is undertaken for plantations rather than the agreed purpose of urbanization. This may occur in other similar mountainous sites in addition to Cameron Highlands [49]. Predominantly in the areas of Kuala Terla, Tringkap, Brinchang, and Ringlet, in which numerous urban areas exist, it can be seen that farmlands are increasing in line with the decrease of urban housing. This could be due to the reduction in demand for housing areas and a rise in crop growth to meet the demand for western crops in the country. Towards the east, an area that was previously heavily forested has been cleared for land to build a road and a hydroelectric plant  Illegal land clearing in small communities also contributes to the overall problem of deforestation. Although this may represent a small change, it is nonetheless an issue for the Land and Planning authorities of Cameron Highlands. This is because the areas of these activities are deep within the rainforest. Moreover, such uncontrolled activities might also be the result of landowners breaching the agreement of forest clearing, which is undertaken for plantations rather than the agreed purpose of urbanization. This may occur in other similar mountainous sites in addition to Cameron Highlands [49].
Predominantly in the areas of Kuala Terla, Tringkap, Brinchang, and Ringlet, in which numerous urban areas exist, it can be seen that farmlands are increasing in line with the decrease of urban housing. This could be due to the reduction in demand for housing areas and a rise in crop growth to meet the demand for western crops in the country. Towards the east, an area that was previously heavily forested has been cleared for land to build a road and a hydroelectric plant and is now categorized as an urban area. Land cover change during the study period is shown in Table 4.
Rapid forest degradation continues to expand until 2019. Figure 6 shows that the North East region of Cameron Highlands experienced land cover change due to construction of a road network, in addition to agriculture plots adjacent to the highway. Primary forest (518.71 km 2 , accounting for 77.46% of total land area, a decline of 23.05 km 2 from 2014) was maintained as the highest class, followed by urban areas (61.10 km 2 , 9.12%) and farm and mixed vegetables (84.61 km 2 , 12.63%). The minor land cover classes of cleared lands and water bodies accounted for 1.10 and 4.17 km 2 , respectively. A noticeable change, of −6.89 km 2 , was found for the cleared land class, and a rise in water bodies of 0.76 km 2 ( Table 5). This change in water bodies is due to the operational status of the Ulu Jelai Hydroelectric dam. The trend of results obtained from this classification is consistent with the findings of Rendana et al. [26], Mohammadi et al. [50], and Razali et al. [51].   The accuracy of the OBIA classification for 2019 was validated using a ground-truthing assessment of GPS points collected from the study area. The results are applied in Figure 6. The producer's accuracy value was 94.60% and the user's accuracy value was 94.65%, resulting in an overall accuracy of 94.60%. Details of the confusion matrix table can be found in Table 6. Notably, only the classification  Tables 7 and 8, respectively), the confusion matrix was derived using Google Earth to validate the accuracy (see Figure 7).

Land Surface Temperature Change Analysis Based on Landsat Satellite Data Compared to MODIS
High temperature values were concentrated in the populous towns of Brinchang, Ringlet, and Kuala Terla (Figure 8, for 2009; Table 9), with typical values ranging between 23 and 26 • C in the months of the first quarter. The results of LST retrieval from Landsat shows a similar range in temperature during April. The monsoon season peaks in the period from the end of January to the start of March, hence there was variability in temperatures during this period. High and low temperatures exhibited differences resulting from the cloud patterns and wind dynamics in the highland areas. As a result, it was difficult to obtain a consistent average for climate measures in the study area. Based on the results of Landsat alone (Figures 8 and 9), it is difficult to ascertain the rise in temperatures relating to the deforestation pattern and the development of urban areas. Thus, the MODIS platform, in addition to air temperature data obtained from the Meteorological Department of Malaysia, was used as a guide to show the rising temperature trends.
Land 2020, 9, x FOR PEER REVIEW 16 of 26 sensors is significant, it is clear that the recorded temperature is correct due to the trend in rising temperatures.     The results are satisfactory with regards to a rising temperature, and this trend was observed after comparison with the multiple datasets. Although the rate of change is not as steep as previously thought, the change is clearly and steadily increasing. As shown in 2011, 2013, and 2019, due to cloud patterns and atmospheric conditions, it is difficult to comprehend the climate during this period; although other years show a better depiction of land surface temperature, this is a result of poor atmospheric correction effects [17,38,52]. By comparison, although MODIS shows lower variability of data between months, there are no absolute highs or lows, resulting in a more accurate depiction of the climate in Cameron Highlands. In 2014, the Landsat 8 sensor measured a maximum temperature of 40.9 • C, whereas MODIS recorded a maximum of 31.1 • C. A summary of the data can be seen in Table 9, which shows the relationship over time. Although the difference between the two sensors is significant, it is clear that the recorded temperature is correct due to the trend in rising temperatures.

Air Temperature against Landsat and MODIS.
Air temperature data were acquired from the Meteorological Department in Tanah Rata, Cameron Highlands. Maximum and minimum data for each month corresponding to the Landsat satellite data time frame are shown in Table 9. According to the data values from Met Malaysia, the highest recorded temperature was 26.4 • C, and the lowest temperature recorded was 11.8 • C. A negative difference denotes a higher estimation of LST compared to the recorded temperature, whereas a positive difference indicates a lower estimation compared to the recorded value. In the urban environment, the air temperature may vary in a broader range due to the urban heat index (UHI). For instance, the months of March and April display a higher temperature as the radiative energy emitted by the surface of pavements, vegetation, urban areas, soil, and water is at a maximum. This is due to an increase in the human influx to the region of Cameron Highlands, primarily due to an increase of visitors, and the passing of the monsoon season, which manifests as dry spells after the period of heavy rainfall [7,9]. The data are not relevant to the whole region but only to certain areas in which these data were taken from the weather retrieving nodes, because it is difficult to ascertain the "correct" temperature value for the region overall [53][54][55][56].
The year 2014 shows the highest deviation across all LST datasets with a minimum of −2 • C and maximum of −15.4 • C. The second highest deviation in temperature occurred in 2019, in which the minimum and maximum deviations were 13.7 and −5.7 • C, respectively. By comparison, a deviation between Met Malaysia and MODIS was found for 2014, with a minimum of −10.2 • C and a maximum of −5.6 • C (see Table 10). Clear skies are a major requirement for accurate LST retrieval. Thus, these unrealistic LST estimations are partly due to the erratic cloud coverage seen in 2011, 2013, 2018, and 2019 ( Figure 9). Because the base image was not trained with a specific emissivity value in specific areas, the resulting image is as shown in the figure [12,17,18,33,35,38].

Verification of LST
To ensure the correctness of the results, the LST must be verified. For this purpose, we obtained data derived from MODIS imagery and air temperature data from the Meteorological Department of Malaysia with links to the station in Tanah Rata, Cameron Highlands. Initially, MODIS was reclassified to the same resolution as the Landsat TIRS sensor, that is 100 m; it is a challenge to ascertain the coverage of the LST using a pixel size of spatial resolution of 100 meters. Other challenges of this method of verification relate to uncertainty in the data quality, namely the correct calibration of the image and apparent implications for the data due to the sensor of the satellite. Air temperature data obtained at ground level was used to address these challenges [42].
Due to the vegetation cover of the study area, much of the vegetation canopy temperature was approximately equivalent to the surface air temperature [24]. Hence, the Landsat and MODIS data for this particular area (primary forest) can be considered a good reference for the air temperature recorded by the weather station. Table 11 shows various datasets which were subjected to a comparison of user predictions for the calculation of the required values. The LST values obtained from Landsat for 2009 to 2019 were averaged, and the associated Root Mean Square Error (RMSE) and bias values were derived. The reference LST was taken from MODIS for air temperature because this was considered to be accurate. It can be seen that the comparison of MODIS LST with local meteorological data has the lowest values relative to the other datasets. With an RMSE of 4.53 and 5.56 • C for maximum and minimum values, respectively, the results are reasonable and consistent, given the nature of the dataset at a spatial resolution of 100 m.

The LULC Effect on LST
Because the study area is covered mostly by vegetation, the area that is most stable in terms of LST changes across the map is primary forest. However, in many other areas, particularly towns, urban zones blend with secondary forests (parks, plantations, and areas of human-made vegetative Land 2020, 9,372 18 of 23 recreation). These areas have a significant effect on the overall LST values because humans are carriers of heat and thus contribute to the anthropogenic heat distribution of a surrounding area. When this heat is massed, it tends to stick to nearby surfaces such as those of buildings. It is known that heat rises while cool air falls. Furthermore, the majority of buildings in the town areas of Cameron Highlands are living quarters with a minimum height of 10 stories and a maximum of 20 stories. The affected areas are Brinchang, Ringlet, and Tanah Rata, in which the population density is higher than that of other towns. These areas radiate an excess heat of 0.5 to 1.0 • C. Similarly, results from a study by Zhou et al. [22] showed that a similar urban environment in Austin, Texas, experienced a rise in temperature of 0.45 to 0.90 • C. The extent to which terrain factors affect the variation of LST depends on the intensity of external solar radiation on a given day. In his study, Zhao notes that the parks and forested regions of Austin have a cooling effect, in which the heat sinks of the area absorb heat and then release cold air.
In contrast, areas without parks display signs of overheating as a result of heat release, surface albedo, and an irregular complexity of landscape structures. Results obtained in the current study are similar to those of a study conducted in Oslo by Venter et al. [57] in which it was found that the temperature of the air was at least 0.6 • C higher in LST modelled using satellite imagery (in urban zones compared to the surrounding outer city). The authors used both Sentinel and Landsat as a means of comparison with meteorological data. The study clearly showed that the temperature of air increases through an impervious surface area paired with a decline of vegetation.
It has been found that the air temperature is at least 2-3 • C lower than the surface temperature depending on the altitude and surface type [48]. This is more evident in urban areas in which there are no pure vegetation pixels. Hence, the air temperature will be slightly lower (by 2-3 • C as previously mentioned) than the land surface temperature because ground radiation conduction is influenced by air particles [58,59]. In other regions of the world, studies have shown the usage of MODIS to estimate air temperature. In Northwest China, Lu et al. [60] conducted an experiment to determine the relationship between sensor temperature readings and on-site measurements. The results showed an RMSE between 2.39 and 3.05 • C; these results were dependent on seasonal and weather variability. Additionally, in the North Tibetan plateau of China, Zhu et al. [61] compared MODIS with air temperature to derive an RMSE of 7.45 • C and a bias of 6.21 • C. The cause of these higher values was attributed to the presence of cloud pixels and haze surrounding the study area. In Egypt, El Kenawy et al. [53] reported an average estimation of 5 • C according to seasonal change. These results are consistent with those of the present study, in which we derived RMSE values ranging between 4.53 and 8.13 • C for different sensor comparisons. Furthermore, Gomis-Cebolla et al. [62] studied the use of LST from MODIS to compare with locally extracted air temperature values in the Amazon rainforest. The study reported RMSE values between 1.58 and 1.85 • C. Because the study was carried out during the winter, the expected RMSE was below estimation. By comparison, if the study was conducted during the summer, the RMSE would be expected to be between 4 and 5 • C.
Various studies of Asian cities comparable in size and density have observed an excess of anthropogenic heating as high as 1.5 • C [21,48,[63][64][65]. A study by Wang et al. [66] conducted in Kuala Lumpur states that "air temperature is a good comparison against satellite-based data obtained temperatures, however, they do not reflect on the 'true' temperate climate of a surrounding area with a multitude of external implications such as wind speed, time, body induced heat radiation as well as heat created from vehicles". This statement is also supported by the study of Kong et al. [45], in which air temperature data also showed signs of outliers, mostly consisting of highs during the summer months of April to June, with estimated values of 0.3-0.5 • C. This could be a case of radiative heat being trapped amongst buildings or another case of climate-induced change as a result of increased influx into the region.
Cameron Highlands is known as a mountainous region with a climatology different from that of surrounding areas. Hence, multiple changes in temperature are frequent throughout the day. When paired with sudden changes in wind speed (which range on average from 2 to 3.4 kmph), this variation leads to uncertainty in the daily average temperature values. It should be noted that the undulating terrain of Cameron Highlands makes it difficult for heat to be radiated smoothly compared to the flatlands. This is more evident when various tall structures hinder the flow of air. Analysis of LST conducted in primary forests in urban areas reveal a higher pattern than LST of primary forests in a mixed zone of urban land cover. Portions of secondary forests are covered with glass roofing, whereas others have aluminum roofing. Furthermore, these roofs are either black in color or comprised of gray cloth to keep insects out. We classified these areas as secondary forests despite their signature showing as an urban zone due to surface reflectance. Notably, these areas have a similar pattern in their retention of heat to that of urban areas because plantations are situated below aluminum roofs [67]. Thus, heat will further radiate in these areas, thus contributing to an increase in the value of the overall LST, because these manmade objects do not regulate LST in the same way as natural objects [68,69]. A more significant surface area density of buildings has a proportionate relationship to high urbanization and, in turn, a weaker LST regulation ability. This results in an increase in the massing of local heat, in addition to increases in the urban heat Index (UHI) and LST. Therefore, it is vital to take note of these particular areas that might affect the overall distribution of LST.

Conclusions
This study explores the LULC and LST for the years 2009 and 2019 based on Landsat and MODIS data. Due to constant population growth in an environment of already cramped spaces for the expansion of urban housing, future trends for our tropical rainforests appear to be inevitable. Results obtained from this study show that the total urban area increased by almost 1.70% since 2009, and farmland and agriculture plots registered a significant 7.71% increase. Combined, these changes led to a decrease in primary forest cover of 8.88%. Particular areas of concern in the already crowded urban region are Tanah Rata and Brinchang due to the ongoing change in LULC. Due to the nature of the terrain and the historical record of landslides occurring in the area, it can be concluded that the area is highly susceptible to natural disasters. Irreversible actions and consequences may eventuate if further development and deforestation occur in the area. Changes in land use and land cover (LULC) over time have also affected the distribution of land cover areas, namely conversion of land cover types from either already deforested primary forests to urban areas, or towards further agriculture usage. Accuracy proved successful using on-site collected GPS points (94.60%). Further analysis is recommended with the use of high-resolution image sensors; Sentinel (10 m) and PlanetScope (3 m) could be considered. Using these image sensors, finer details can be extracted for an LULC classification, which would yield more classes and provide a more accurate description of the study site.
The land surface temperature showed a rapidly increasing trend, particularly in the urban areas of Cameron Highlands, in particular, Kuala Terla and Ringlet. During a ten-year period, the recorded temperatures showed a rise of about 2-3 • C across all channels, with an extreme of 5 • C on average. The use of Landsat 8 with the SWA method proved reasonable with the results obtained, although outliers arose partly due to the pre-processing stage. However, MODIS and weather data from Met Malaysia confirmed a rise in the region of LST, the primary focus of this study. It is recommended that the calibration of LST data be undertaken in future projects to avoid the uncertainty of erratic LST values as a result of cloud coverage. The results of this study suggest that LST is heavily dependent on seasonal radiation conditions mainly due to the monsoon season and its subsequent effects, which cause the overestimation of cooling values. Due to the undulating terrain of Cameron Highlands, which has a gradient of multiple forest altitude classes, future studies should calculate LST at different forest altitudes and take into account the forest cooling effect and the urban heat index. Population density also plays a role in the acquired LST values, in which concentrated portions of the study area displayed higher temperature values compared to smaller towns. Furthermore, thermal infrared remote sensing applications could be adapted for future studies to allow for better measurement of heat radiation in concentrated urban areas. Ground-based temperature data in towns can provide more accurate insights into urban heat indexes and their effect on air temperature. Correct usage of water vapor data for each zone should be taken into account in future studies, in addition to considering the use of Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) and Advanced Very-High-Resolution Radiometer (AVHRR) data. Research results provide a fuller understanding of the drivers of increased LST, which is mainly due to increased demand for urbanization and, in turn, is linked to the expansion of deforestation. However, it is a challenge to prevent urbanization without a change in the decision making of upper management. Thus, it is suggested that the state government authority considers enacting policies and strategies that support the control of deforestation in the highlands. This will address ecological and environmental concerns and help ensure that development takes place sustainably.