Effect of land cover change and elevation on decadal trend of land surface temperature: a linear model with sum contrast analysis

Land surface temperature (LST) is a significant factor in surface energy balance and global climatology studies. Land cover (LC) and elevation are two factors that affect the change of LST, and their effects depend on different geography. This study aims to demonstrate an alternative approach to examine the change of LST during 20 years (2001 to 2020) on Taiwan Island and to investigate the effect of LC change and elevation on a decadal trend of LST using a linear model that adjusting for each determinate factor. MODIS LST and LC data, as well as GMTED2010 elevation product, were downloaded available website. The natural cubic spline function was used to model annual seasonal patterns in LST. Linear regression model was used to estimate decadal change of long-term LST time series. Weighted sum contrasts linear regression was used to assess the effect of LC transformation and elevation on the decadal LST change by comparing adjusting mean of all factors. The adopted analysis method was an appropriate approach to assess categorical factors than those based on treatment contrasts, requiring specifying a control group to compare means and confidence intervals. Results showed that there was an increase in LST for most of the island. The average daytime and nighttime LST trends were 0.12 and 0.31 °C/decade, respectively. However, areas in the southern part of the north–south direction mountain range show a statistically significant increase in LST in both daytime and nighttime. The major landslides caused this noticeable change of surface temperature due to the catastrophic damage of typhoon Morakot in 2009. The results also revealed that the different pattern of LC change has a significant effect on daytime LST, but not on nighttime LST trends. The elevation above 600 m had affected both daytime and nighttime LSTs.


Introduction
Land surface temperature (LST) is a significant parameter for surface energy balance and urban heat island climatology studies. It is influenced by the land surface characteristics, for instance, land cover (LC) change of each type and elevation. The increasing urbanization and agriculture have caused significant changes in the LST (Khandelwal et al. 2018). Therefore, the warming temperature across the globe reveals an upward trend and shows the average temperature at the warming area more than the cooling area. According to NOAA's 2020, annual global climate of land and ocean temperature has increased at an average of 0.08 °C per decade since 1880. The rising of surface temperature caused by the transformation of LC and elevation (Guo et al. 2012). Dynamic of LC is another important factor causing changes in LST (Aboelnour and Engel 2018). Variation in LC can be detected from aerial images based on the difference in color in different areas and is often the result of the influence of human activities on natural ecosystems which have regional and global effects on the environment (Guo et al. 2012). Another significant factor influence LST is elevation. It is the height above or below a fixed geographic reference point, known as a geoid and is one of many factors that affects climate change. Some previous studies have examined the direct influence of elevation on LST from MODIS data and have found that elevation significantly affected LST across large areas (Phan et al. 2018).
Therefore, many studies have applied different statistical methods to evaluate the effect of LC indices on LST. According to Sahana et al. 2016, they used Landsat Satellite image to study the effect of land use/land cover (LULC) change on LST distributions in Sundarban Biosphere Reserve, India. Split window algorithm and spectral radiance model were used to retrieve LST and to assess the relationship between LULC change and LST distribution. They found an increase in non-evaporating surfaces and a decline in vegetation had increased LST. Moreover, Kayet et al. 2016 used three sets of Landsat image to investigate the effects of LULC changes on the LST distribution in Saranda Forest, Jharkhand. A remote sensing technique was used to detect the LU changes and their effect on LST and variation in mean LST from LC change hot spots. They found that the loss of vegetation area causes an increase in LST. LC classification system was used to classify LC variation from Landsat image and determine the impact of LU changes on LST with five LULC classes (Karakuş 2019). The results from this study indicated that the change in LST depends on the type of LC indices and that LST increased most in built-up urban areas, whereas it experienced a decreasing trend in rural areas. Furthermore, the study of the effect of LC heterogeneity on LST change in Phuket Island, Thailand (Wongsai et al. 2020a, b), using the Generalized Estimating Equations (GEE) model confirmed that LU had a significant effect on LST. The study results show that different LULC change patterns had a significant effect on long-term variation in LST.
Previous studies had evaluated the linear relationships between LC indices, elevation, NDVI and LST. Linear regression was used to investigate the correlation between LST and NDVI for each land use/land cover (LULC) in Guangzhou, China (Sun et al. 2012). They reported that the coefficients between LST and NDVI were great in cropland, grassland and forest. They also examined the correlation between LST and each LULC indices, including MNDWI, NDBI, NDVI, NDBAI and elevation. The result indicated that the correlations between NDBI, NDBAI and LST were significant positive (Sun et al. 2012). Moreover, Phan et al. also used a linear regression model to evaluate the relationship between LST and elevation in Northwest Vietnam. The aim of their study was to investigate the spatial and temporal correlation between the monthly average of daytime and nighttime. The results show that a stronger correlation was found at nighttime than daytime, and a stronger correlation was found in the hotter months during nighttime. The results also indicated a decrease of LST for each 1,000 m elevation increase (Phan et al. 2018). On the other hand, many studies had reported that LST increased with the increase of elevation (Shrestha and Aryal 2011;Dong et al. 2014;Phan et al. 2018).
The effects of LC, elevation and NDVI on LST with regression have been widely studied. Prasetya et al. used multiple linear regression to investigate the effect of elevation, LC and NDVI on decadal LST increase in Central Sumatra, Indonesia. They found that r-squared accounted for 37.6% of the LST change and it indicated 62.4% of the variation to other factors (Prasetya et al. 2020). However, the method in their study did not adjust for unequal sample size of those factors. Wongsai et al. used the analysis of variance (ANOVA) test to compare confidence interval (CI) of surface temperature change between 10 groups of LULC change patterns with insignificant number of sample in each group. They found that there was a significant difference in the mean of increasing temperature among the groups of LULC pattern. The comparative CI of the changing temperatures in the nighttime was around the average increasing temperature. The insignificant p-value from the ANOVA test indicated that there was weak evidence that the changing temperatures among the group of LULC patterns are different (Wongsai et al. 2020a, b). However, the ANOVA test of comparing mean in their study did not adjust for unequal sample size bias. It is considered robust to moderate departures from their assumption, but that is not true when the sample sizes are very different. According to Keppel (1991), there is no good rule of thumb for how unequal the sample sizes need to be for heterogeneity of variance to be a problem. Thus, if we have equal variances in our groups and unequal sample sizes, there is no significant problem. The only problem is if we have unequal variances and unequal sample sizes.
The main aim of this study is to introduce an alternative approach to investigate the effect of patterns of LC change and levels of elevation on decadal LST change during 20 years (2001 to 2020) where the number of sample size in each group of factors found in the study area is typically unequal. A linear model with weighted sum contrasts method adjusted for each determinate factor was adopted to demonstrate the comparative method of LC change and elevation effects on long-term LST change.

Study Area
Taiwan consists of a rugged mountain chain running north to south with a flat coastal plain to the west. The highest mountain in Taiwan rises to 3,952 m, and Taiwan is the fourth-highest island in the world. The total area of Taiwan island is 35,883 km 2 . The annual average temperature is approximately 22 °C, with an average of 21.7 °C in the northern area and 24.1 °C in the southern area. The coldest months, from January to March, with the lowest temperature around 10 °C. There may be occasional but unusual occurrences of frost or snow in high mountain areas in some years.
The hot season is between June and August, with the highest temperature of around 38 °C (Taiwan Weather 2019). Most lowlands were located along the western part of the island, where the urban and cultivation areas were established.

MODIS LST Products
LST data in comma-separated value (.CSV) file format were downloaded from the MODIS Land Product subset tool (ORNL DAAC 2018). The obtained MOD11A2 product, which provides a level 3 day MODIS/Terra LST and Emissivity at 8-day temporal and 1 km spatial resolution (Wan and Hulley 2015), contains daytime and nighttime LST (10:00 and 22:00 at a local over-pass time). Both daytime and nighttime LST data were used and obtained from Feb 18, 2000 to Feb 11, 2020, which is precisely 920 observations and cover a 20-year study period. We also generated polygons that match with the sinusoidal pixel size of the obtained 1-km MODIS LST data using projection/grid conventional formulas (Wongsai et al. 2020a, b;Abdulmana et al. 2021). Each polygon was linked to the corresponded LST time series obtained in.CSV file format. The created polygons were future used for aggregating land cover and extracting elevation data.

MODIS Land Cover Products
The MODIS land cover data for 2001 and 2019 were also downloaded from the United States Geological Survey (USGS) Data Pool website (USGS data 2020) . We used Terra MODIS MCD12Q1 version 6 product, which provides 500 m resolution global land cover types in 2001 and 2019. The acquired MODIS land cover data derived from six different classification schemes that offered different land cover classes (8 -51 classes). However, we use land cover type from the International Geosphere-Biosphere Programme (IGBP) classification, categorizing global land cover into 17 types. Two MCD12Q1 hierarchical data format files (.hdf) were required to cover the study area (tile h28v6 and h29 v6). Furthermore, we downscaled 500 m land cover data to 1 km resolution to match 1 km LST pixel by a majority (or maximum count) using Zonal Statistics function and the generated 1-km MODIS LST polygons as a polygon layer containing the zones in QGIS software. The 13 types of land cover found in the study area were re-classified into five major groups to simplify and reduce the distinction of land cover change pattern. The new groups were categorized according to their physical characteristics that might affect or relate to their thermal and biophysical properties, as detailed in Table 1.

Elevation Data
The Global Multi-resolution Terrain Elevation Data 2010 (GMTED2010) data set was used in this study (Danielson and Gesch 2011). Raster elevation data with 7.5-arc-second spatial resolutions (250 m) were retrieved from the United States Geological Survey (USGS) Earth Explorer website. The GMTED2010 product contains seven sets of elevation data for each of the 30-, 15-, and 7.5-arc-second spatial resolutions. However, the 30-arc-second resolution (1 km) product not accurately matched with 1 × 1 km 2 pixels of the obtained MODIS LST and land cover products. Thus we used 7.5-arc-second spatial resolutions (250 m) median elevation dataset and aggregating it into sinusoidal 1-km MODIS LST pixel areas by the median value. The process of elevation downscaling was also performed in QGIS software using Zonal Statistics function and the generated 1-km MODIS LST polygons. We grouped elevation into nine levels, as detailed in Table 2, for investigating the effect of elevation on change of LST in 20 years of study period.

Seasonally Adjusted and Time-independent Series
To estimate the change in 20-year LST time series, the data series need to be stationary and time-independent. First, a cubic spline function with annual periodic boundary condition (Wan and Hulley 2015) was used to extract seasonality from the LST time series. To do that, time series were re-ordered from 920 to 46 observations by Julian Days. The observed values within the same observed day that appear to be outliner from their distribution were omitted from this seasonal component extraction. We set eight knots at Julian Days 10, 50, 100, 150, 210, 260, 310 and 350 for an adopted cubic spline function to estimate its coefficients. Then, a linear regression model was used to fit the original LST data to the time point of data observation and the derived cubic spline coefficients (Wongsai et al. 2020a, b;Abdulmana et al. 2021). Finally, the seasonally adjusted time series were determined by subtracting the original time series from the obtained seasonal component and then adding the mean of the seasonal component.
As we used a regression model to assess the effect of LC change and elevation on long-term LST data, seasonally adjusted LST (LST s.adj ) time series need to account for autocorrelation. Time-dependent in time series data could affect the standard error of the estimated parameter in the regression-based analysis model, which could narrow its confident intervals, and the statistical assumption of independent and identically distributed errors term at two specified time points is violented. To eliminate autocorrelation in LST s.adj , a nonseasonal autoregressive integrated moving average (ARIMA) model with a maximum three lags was used to measure the maximum correlation of time series values between two different time lag (Wongsai et al. 2020a, b). Then the LST s.adj series were filtered to reduce autocorrelation by using the following formula: (1) where ϵ is the fitted residual from the selected ARIMA(p,0,0) model, σ is the term of population standard deviation, and μ is the term of population mean (Wongsai et al. 2020a, b).

The Weighted Sum Contrasts Linear Regression
In this study, the outcome was the estimated increase/ decreases LST in a decade resulting from the linear regression fitting of filtered LST s.adj data series, categories of elevation and groups of LC change pattern were the determinant of interest. When a determinant in a regression model is a categorical variable, it is called a factor, and it appears in the formula for the model as a set of k-1 parameters, where k is the number of different categories. The linear regression model for the study was now formulated as y = a + factor(x 1 ) + factor(x 2 ), where y is the LST trend in °C/decade, x 1 is the identity of the elevation group, and x 2 is the identity of the LC change group. Since factor(x 1 ) only has k-1 parameters, as well as factor(x 2 ), therefore the model was formulated as where b i and c i were the coefficients of x 1 and x 2 at identity i , respectively, and b 1 and c 1 = 0 in the case that no contrast option was selected when specifying the model. In general, if contrast option is not specified, the parameters are matched with each factor level in alphabetical order, but with the first parameter set to 0. A contrast in linear regression is a linear combination of determinants whose coefficients add up to zero, allowing comparison of different treatments (Casella 2008). The contrasts for this method are known as "treatment" contrasts because they were originally used in experiments designed to compare one or more treatment groups with a control group, where it makes sense to set the parameter associated with the control group to 0, so that the parameters associated with the treatments are the actual treatment effects. However, we have no control group in our study, and consequently, two categories in our determinant and covariate were omitted. To achieve the objective of our study in comparing the mean of all interested factors to assess their effect on the LST trend, we used weighted sum contrasts developed by Tongkumchum and McNeil. This method suggests constructing confidence intervals for comparing means that does not involve selecting a reference group and thus gives an informative confidence interval for comparing each means with the overall mean (Tongkumchum and McNeil 2009). When there is no control group, different contrasts called "sum" contrasts can be used. Sum contrasts constrain parameters associated with each factor level so that each gives a measure of its difference from the overall mean of the outcome. The formulation is similar to that for treatment contrasts, with additional terms ( b 1 and c 1 ), that is, For equal sample sizes, the sums of the sets of coefficients are 0, so But this is not the case in our study, where the sample sizes for the elevation and LC change group were unequal, as informed in Table 2. So that weighted sum contrasts method was adopted to measure mean in each factor level with unequal sample size.
During the process of fitting linear regression with weighted sum contrasts method, overall mean (average of all LST trends), overall r-squared (adjusted r-squared from regression model without contrast option), crude mean (mean of LST trends in each category of factors), adjusted r-square and p-value for each factor were calculated. All of the computed statistical model's parameters were investigated and compared. This method of weighted sum contrasts collaborating with linear and logistic regression had been used in many preceding studies. It used to compare blood lead levels among children in Pattani river, Thailand, 22 and investigated HIV mortality by age group and gender in Thailand between 2014 and 2015 (Tulu et al. 2020).
Moreover, the "democratic" confidence intervals for each category of factor were computed for assessing the variation in a mean. The confidence intervals from a two-sample t-test and analysis of variance (ANOVA) are valid for each sample individually (and use all the information in the data because they use pooled standard deviations), but they are not appropriate for comparing means. For comparing means, we use confidence intervals for the difference between the means, which is the same as the difference between each mean and the overall mean. They are consistent with the p-value and are "democratic" because they do not need to select a control group as a basis for comparison.

Land Cover Change and Elevation
Figures 1 shows map of Taiwan and 2 illustrates the distribution of land cover change in Taiwan between 2001 and 2019 and the elevation, which was classified into nine classes. Considering the conversion of land cove, the Unchanged DV area covered 58.64% of the island and found most in the mountain range around 1090.50 m above sea level by median (IQR: 527.75 -1839.25 m) of all unchanged DV areas. Unchanged CC, BB and WW lands contributed by 11.88%, 7.81% and 0.48%, respectively, of the island area. Most of them were located in the land whose elevation was lower than 50 m and covered approximately 21.05% of the island area. The overall elevation by the median for unchanged CC, BB and WW lands was around 27 m, 22 m c i x 2i and 2.5 m, respectively. Most of the unchanged SV land was found in the higher elevation (around 120.5 m by median) and contributed by 11.42%. Nearly two decades, land cover of Taiwan has changed only 9.81% of the entitled island, and the majority of the transformed areas were found in the land that its elevation was around 90.5 m by median (IQR: 24.5 -259.0 m). Areas with its land cover had transformed from SV to CC were the most land conversion as 34.51% of the LC change areas. LC change between SV to DV was the second conversion form which around 19.93%. Shrinkage of DV area that transformed to SV land and the increases of SV land converted from the cultivated land cover (CC) were contributed as the same amount as approximately 12.94% as described in Table 3.

Decadal LST Change
The change of LST in a decade, which estimated using linear regression model fitting for each filtered seasonally adjusted LST time series (each pixel), and its corresponding null hypothesis significance testing (p-value) are illustrated in Fig. 3. The daytime LST trends range from -1.24 to 1.62 °C/decade, and their average was 0.12 °C/decade. The daytime surface temperatures that increase more than 1.00 °C/decade with statistical significance (p-value < 0.05) were found in the southern part of the mountain range, as shown in Fig. 3(a). Some moderate increasing daytime temperatures (0.30 -1.00 °C/decade) with statistically significant also found in the southern mountain range and most in the cropland (39.64%) and urban area (42.64%) which located around the middle and southern part of western lowland (< 50 m above sea level). There were 29.35% of decreasing LST trend throughout Taiwan island. However, there was only 4.21% of all area that its surface temperature was decreasing with statistically significant, which found most in the valley and lowland along the earthen coastline of the island. The average of nighttime LST trends was 0.31 °C/ decade and ranged from -1.05 to 1.00 °C/decade. More than 72.14% of all area that its surface temperature was a statistically significant increase. The increasing surface temperatures were found throughout the island except for the south-earthen region, where the surface temperatures trends were declining (7.27% of the entire study area). However, only 21.06% of LST trends in this region was statistically significant decreased. Figure 4 illustrates the estimated mean of daytime and nighttime LSTs with comparative confident intervals after adjusting for elevation and land cover change with the overall mean for each factor. Graphs in the first column of Fig. 4 show the effect of elevation, whereas graphs in the second column show the effect of land cover transformation on decadal daytime and nighttime LST trends, graphs in the first and second row, respectively. The horizontal red lines indicate the overall mean by simple average LST trends from all pixels. The green dots indicate crude means which were the mean of LST trends in each group of elevation and land cover change. The label "r-sq:" shows adjusted r-squared from regression fitting with separated determinants. The label "Overall r-sq:" denotes adjusted r-squared from the regression model adjusting for elevation and land cover change factors. The label "n:" indicates the percentage of unequal sample size for each level of factors. The overall mean of daytime and nighttime LST trends was 0.12 and 0.31 °C/decade, which implies that the surface temperatures in the study area were escalating. Graphs in the first column clearly show an effect of high elevation (> 600 m) in both daytime and nighttime LSTs. In daytime, there was a sizable difference by comparing mean in each difference of LC conversion and the p-value, which for the difference between two or more population means, confirms that the hypothesis that the population means are equal can be rejected. The mean of LST trends in the unchanged BB area and mean of LST trends in areas where their land cover conversed from DV to SV land was above average as 0.24 °C/decade, which is two times higher than the overall mean. In contrast, mean of LST trends in conversion land from SV to DV land was lower than zero as -0.04 °C/decade and different from the overall mean by 0.16 °C/decade. Mean of LST trends in the other groups of LC change showed an insignificant difference from the overall mean. However, their confidence intervals were not overlapping each other or ride on the overall mean. The difference in LC change does not affect nighttime LST trends. The mean of LST trends in all different groups of LC change was very close to the overall mean. Some of their confident intervals overlapped each other and crossed the overall mean line, suggesting that there was insufficient evidence of the difference.

Confounding and Adjusted Comparative Mean
Effect of elevation on daytime and nighttime LST was similar. Means of daytime and nighttime LST trends at lowland (< 50 m) and land higher than 1000 m were above the overall mean, as displayed in Fig. 4. The sizable rise of an average of LST trends in the lowland was contributed by the built-up and cultivation land cover. Although, the urban and cropland area mainly not expanded or change as informed in Table 3. However, the increasing of built-up and impervious surface density in 1 × 1 km 2 during 20 years of study period may cause the rising trend of LST, especially in daytime. Surface temperature of land in high elevation generally lower than in the lowland due to the local climate, types of land cover and utility of land (Phan et al. 2018). If there was no change of land cover, LSTs trend to decline corresponding to higher elevation. From our result, 33.7%, 24.51% and 15.80% of unchanged DV, SV and CC land, respectively, (74.01% in total) were found within the areas that their elevations were between 50 and 350 m above the sea level (elevation level 2-4).
Moreover, the mean of LST trends in deforestation or degraded forest area (land converted from DV to SV), which   It implies that geographical relation between elevation and land cover type at local site affect the dynamic of LST. If the correlation between these two factors was disregard and their effect were analyzed separately, those results may not be accurate and comparable. In this study, elevation associated with LC change in the data and distorted the association between LC change and LST trends, which issued confounding. A confounding variable influences both a determinant and an outcome variable and can affect that association between them. The testing result from 10% rule of thumb informed that the parameter estimated for LC change in all levels has a percentage change of more than 21% from the unadjusted, or crude, estimate (from simple linear regression) to the adjusted estimate (from multiple linear regression). It confirmed that the association between LST trends and LC change is confounded by elevation. As mentioned in Sect. 2.6, a two-sample t-test (for the comparison between two groups of sample size) and ANOVA (for determining whether three or more populations are statistically different) are not appropriate for comparing factors that might affect each other. The regression model can be used to examine the association between multiple covariates and a numeric outcome. This model can be employed as multiple linear regression to see through confounding and isolate the relationship of interest (Pourhoseingholi et al. 2012). Moreover, there were unequal and unbalanced sample sizes of LC change in the study area. Since sample sizes in many levels of LC change (WW, SV to CC, SV to DV, DV to SV, CC to SV and others) were small and some levels were under-represented, which biased mean in the level of the factor of interest, as presented by the difference between crude mean and adjusted comparative mean with "democratic" confidence intervals. The weighted sum contrasts method corrects this bias. It gives the adjusted means of LST trends for each level of the factor of interest, which can be obtained simply by fitting the model using weighted sum contrasts and omitting the terms from all other factors. Graphs in Fig. 4 accurately show how LST trends vary with land cover change after accounting for elevation (as the confounding factor), which accurately reflects the population and reveals the effect of elevation and LC change on daytime nighttime LST trends.

Influence of LC Change and Elevation on LST Trends
Overall, the surface temperatures in the study area were escalating, especially at nighttime that up to 0.31 °C/ decade on average. Analysis of LST during the 20-year study period demonstrated that change of nighttime LST primarily increases all over Taiwan island compared to daytime LST. Over 73.67% of study areas show that its nighttime LST increases with statistically significant whereas daytime LST increases only in 26.72% of the study area. Considering Fig. 2(a) and Fig. 3(c), the distribution of increasing LST trends display roughly spatial relation with land cover change. However, LC type and the form of LC change had a negligible effect on variability of nighttime LST trends, which had insignificant differences, as shown in Fig. 4. Only mean LST in BB area (urban, built-up and barren) show a considerable difference from the overall mean. Using nighttime LST data may suitable for long-term regional and global climate change study because land coves and their change over time had a minor effect on LST increasing. This result supports the previous study in the smaller study area of Wongsai et al. 2020a, b where the ANOVA test in their study suggested that there was weak evidence that changing temperature amount groups of LC and their change patterns were different. Although, ANOVA analysis in their study was not adjusted for unequal sample sizes. On the other hand, daytime LST trends substantially varied in various LC type and their change pattern. Transformation of vegetation cover in both way, deforestation and afforestation, strongly impact daytime LST trend. There were around 58.65% of the study area covering with dense forest and more than 71.40% of that forestry area located above 600 m from sea level.
However, both daytime and nighttime LST in 20 years were increasing by overall. According to elevation graphs in Fig. 4, both daytime and nighttime LST trends in areas over 1,000 m above sea level, primarily unchanged DV land cover, appeared to escalate significantly. However, when the study area was divided into north and south regions by 23.60°N latitude, as shown in Fig. 1, areas over 600 m only in the south region advocate the overall daytime and nighttime LST trend. Figure 5 shows the effect of elevation and LC change in the north and south region on daytime LST trend. It reveals that not only the high elevation (> 600 m) area in the southern part of the island strongly contributed to the increase of LST trend, but also the class of land cover conversion. The pattern of increasing/ decreasing temperature in each level of elevation and LC change factors was similar to the daytime LST increase for the entire Taiwan island, as shown in Fig. 4.
Our preliminary investigation on the cause of increasing surface temperature in the southern highland found that there was apparently land cover change. Figure 6 reveals the soil erosion and massive landslide in many steep valleys around the southern mountainous region. The physical change of land cover in this area was noticeable since 2005, and this natural phenomenon distinctly visible from historical imagery in Google Earth planform after 2009. Typhoon Morakot had damaged Taiwan in 2009, especially in southern Taiwan. The deadliest typhoon to impact Taiwan in recorded history produced an amount of rainfall, mudslides and flooding. It damages vast agricultural area and causes many casualties (Chien and Kuo 2011;Wu 2013). The typhoon also was the cause of the transformation of the hillside and base of the valley covered by dense vegetation become the rocky and bare soil surface caused the increase of its surface temperature in this part of Taiwan. However, land cover around this region was unchanged from 2001 to 2019 based on MODIS land cover data. It might be that the total area of land cover change did not dominate the unchanged area within 1 × 1 km 2 area. However, a substantial quantity of change area, especially conversion from vegetation cover to built-up and bare land, and its spatial heterogeneity could significantly increase its surface temperature (Wongsai et al. 2020a, b).

Conclusion
In this study, we investigated the change of LST and LC during 20 years (2001 to 2020) and demonstrated the use of a linear model with weight some contrast for adjusting each determinate factor to assess the effect of LC change and elevation on a decadal trend of LST in Taiwan. Taiwan island is an example of diverse terrain with a high mountain range, gorges and sharp valleys and the flat to gently rolling plains, where 90% of the population lives. LC in Taiwan had changed only 9.81% of the whole island in the last two decades, and the majority of the transformed areas were found in the land that its elevation was around 90.5 m by median. This raises awareness when it comes to analyze and compare the average LST change in different groups of LC change and level of elevation as the factors of determinant that could have a significant difference in sample size and influence on each other. The presented method can compare several means with comparative confidence intervals based on weighted sum contrasts, which often more appropriate than those based on treatment contrasts requiring specifying a control group. The results reveal that the average daytime LST trends were 0.12 °C/decade while the average of nighttime LST trends was 0.31 °C/decade. The difference in LC change does not affect nighttime LST trends. However, the mean of daytime LST in urban, builtup and barren areas shows a significantly different overall mean. The effect of LC and elevation differed according to the difference in geography and sub-region. There was an effect of high elevation (> 600 m) in both daytime and nighttime LSTs. The area above 600 to 2000 m of elevation in Taiwan with their surface temperature trends to increase is worth investigating. Change of LC, elevation and other factors influencing the increase in LST should be investigated in other different geographical areas in further studies using the presented method. Fig. 6 a enlarged area of interest where high daytime LST increased with statistically significant and high-resolution airborne imagery from Goggle Earth platform in b 2001 and c 2019 that exposed the change of land cover in the southern mountainous region due to significant soil erosion in many steep valleys Author contribution Sahidan Abdulmana, Noppachai Wongsai and Sangdao Wongsai conceptualized this research idea, obtained the data, and performed statistical analysis, interpretation of the result and discussion. Final proofread was done by Noppachai Wongsai and Apiradee Saelim. All authors contributed to writing and editing the manuscript.

Data availability
The data and material used in this research will be available on request to the corresponding author.
Code availability The code used in this research will be available on request to the corresponding author.

Declarations
Ethics approval The authors confirm that this research is original and has not been published in any journal (in whole or part).

Consent to participate Not applicable.
Consent for publicatio All the authors have consented to publish this research.

Conflict of interest
The authors declare no competing interests.