Relationship among land surface temperature and LUCC, NDVI in typical karst area

Land surface temperature (LST) can reflect the land surface water-heat exchange process comprehensively, which is considerably significant to the study of environmental change. However, research about LST in karst mountain areas with complex topography is scarce. Therefore, we retrieved the LST in a karst mountain area from Landsat 8 data and explored its relationships with LUCC and NDVI. The results showed that LST of the study area was noticeably affected by altitude and underlying surface type. In summer, abnormal high-temperature zones were observed in the study area, perhaps due to karst rocky desertification. LSTs among different land use types significantly differed with the highest in construction land and the lowest in woodland. The spatial distributions of NDVI and LST exhibited opposite patterns. Under the spatial combination of different land use types, the LST–NDVI feature space showed an obtuse-angled triangle shape and showed a negative linear correlation after removing water body data. In summary, the LST can be retrieved well by the atmospheric correction model from Landsat 8 data. Moreover, the LST of the karst mountain area is controlled by altitude, underlying surface type and aspect. This study provides a reference for land use planning, ecological environment restoration in karst areas.


LST inversion results.
Due to the lack of measured LST for verifying the inversion accuracy, the meteorological site temperature data were used for validation, which showed that the lowest, highest and average temperatures of Yinjiang County were 19 °C, 29 °C and 24 °C, respectively. The LST calculated by the atmospheric correction method was 14.75 °C to 40.59 °C, and the mean LST was 27.67 °C. Furthermore, the statistics showed that the LST of 99% of the study area was within 14.75 °C to 32.40 °C, thereby indicating that the LST retrieved by the atmospheric correction method was the same as the actual LST.
In Fig. 1a, the LST of the study area was high in the west and north and low in the east and south and showed an overall downward trend from northwest to southeast. Moreover, the LST formed a high-value area in Yinjiang County and other construction land-intensive areas, thus showing a regular urban heat island phenomenon (Fig. 1b). In the east of the study area, i.e. Fanjing Mountain, the LST formed a low-value area (Fig. 1c). The combination of the study area DEM elevation and land use data was used for superposition analysis, showing that the trend was mainly due to the combined effects of land use and topography. In addition, the combination of high-definition satellite maps showed high-temperature zones different from the traditional urban heat island in the county suburbs. Moreover, the highest LST was observed in the suburbs of the county.
Relationship between LST and elevation. The spatial distribution of the LST and altitude exhibited an opposite pattern. As shown in Supplementary Fig. S1, the LST colour changed from red to green from northwest to southeast, i.e. the LST gradually decreased with the increase in altitude. In the west-east and northwest-southeast directions, the terrain distribution alternated between valleys and hills. Correspondingly, the high-value peaks and low-value valleys of the LST were alternately distributed. Furthermore, in the area with relatively high elevation difference, the LST had a noticeably vertical variation rule with altitude change. Taking Fanjing Mountain and Langxi Valley as examples, Fig. 2a shows a low altitude in Langxi Valley area, where the LST was high. Figure 2b shows that the LST of Fanjing Mountain gradually decreased with the increase in altitude. However, the profile curves of the LST fluctuated more frequently than that of the corresponding elevation profiles. Meanwhile, in Fig. 2c, LST had a significant negative correlation with altitude, with a correlation coefficient of −0.782. The drop rate of the LST was 7.6 °C/km, which is 1.26 times that of the global average (6.5 °C/km). It is mostly probably because the underlying surface of the study area is complex, such as the hot valley, which is quite different from the area where the underlying surface is homogeneous.
Relationship between LST and land use. LST characteristics of different land use types. LSTs among different land use types were ranked as follows: construction land > unused land > cultivated land > water body > grassland > woodland ( Table 1). The LST of construction land and its standard deviation were the highest. Overall, the LST of construction land, unused land and cultivated land, which had high human activity intensity, was high; the LST of water body, forest land and grassland, which had low human activity intensity, was low. Therefore, land use types should be rationally planned and a cooling effect should be induced through green vegetation and water. Table 2 shows the post hoc test results of the mean differences in the LST among the different land use types in the study area 11 . The table has a 6 × 6 symmetrical matrix with 15 combinations. Except for that between the cultivated and unused lands, the differences among the pairs (14 pairs) were significant (sig. <0.01). The cultivated land was mainly dry land with low vegetation coverage, whereas the unused land was mainly bare rock in karst and sparsely vegetated grassland. Thus, the cultivated and unused lands exhibited only a slight difference in surface emissivity, thereby resulting in LSTs with no significant difference. Influenced by different heating capacity, heat conduction and thermal radiation properties, the LSTs of the woodland, grassland, water body and construction land exhibited significant differences, indicating that the different land use types contributed differently to the regional thermal environment effect.
LST features of different land use types based on slope direction and elevation. The solar radiation received by the earth is different because of the difference in slope directions of the earth's surface, which leads to the difference in LST 24 . In this study, the imaging time for Beijing was 11:15:48 a.m. and the mountain direction of the study area was southeast for the sunny slope and northwest for the shady slope. In the study, the mean LST of the sunny slope was 1.92 °C higher than that of the shady slope. In addition, in Supplementary Fig. S2, the average LST of each land use type per 100 m altitude range was obtained. Overall, in each altitude interval, namely at approximately  the same altitude, the LST of the different land use types exhibited differences, and their performances were ranked as follows: unused land > cultivated land > construction land > grassland > water body > woodland. Meanwhile, the range of LST change varied across different classes and altitudes. In addition, for all land use types, LST was negatively correlated with the elevation, and the linear goodness of fit (R 2 ) was more than 0.9, excluding that of unused land. As shown in Table 3, under sunny slope conditions at the same altitude, the LSTs of all land use types were ranked as follows: unused and construction lands > cultivated land > grassland > forestland > water body. By contrast, under shady slope conditions at the same altitude, the LSTs of all land use types were ranked as follows: unused land > construction land > cultivated land > grassland > water body > woodland. Moreover, for each land use class in the same altitude range, the LST under sunny slope conditions was higher than that under shady slope conditions. The maximum difference could reach 5.14 °C (Supplementary Table S1). Water body was excluded because the water bodies in lowlands, such as valleys, exhibited no noticeable difference between sunny and shady slope conditions. Figures 1   and 3 show that the LST and NDVI exhibited opposite spatial distribution patterns. In Fig. 4a and b, the LST showed the opposite trend with the corresponding NDVI in two different directions. At the macro level, the LST in Yinjiang County formed a high-value area, thus showing the city heat island phenomenon; moreover, the low-value zones of NDVI in this area exhibited a small effect, which are marked by red boxes in Fig. 4a and b. At the micro level, the high peaks of the LST value corresponded to a low 'valley' of the NDVI; meanwhile, the areas marked by blue boxes show the high peaks of the NDVI that exactly corresponded to the low 'valley' of the LST. Three points in the CD profile, i.e. A 1 , A 2 and A 3 , were characterised by relatively low NDVI and LST. The high-definition satellite map shows that these special points were in the Yinjiang River. The low LST and NDVI of the water body indicate its good cooling effect. The water temperature was relatively stable with vegetation coverage that was approximately zero, and the water LST was positively correlated with NDVI 25,26 . Meanwhile, the main part of the unused land was a non-vegetation area. Thus, the water bodies and unused land were not considered in the discussion about the relationship between NDVI and LST among the different land use types.

Relationship between LST and NDVI. Spatial distribution characteristics of LST and NDVI.
Quantitative relationship between LST and NDVI in single land use type. In the regression function shown in Table 4, Y and X denotes LST and NDVI, respectively, and the saliency coefficients are also listed. The table shows that the LSTs of the different land use types have negative linear relationship with the NDVI. In the regression functions, the absolute values of the first-term coefficients were ranked as follows: woodland > grassland > construction land > cultivated land, indicating that the LST of the woodland decreased fastest with the increase in NDVI, followed by grassland and construction land, and the LST of the cultivated land decreased relatively slowly. Therefore, the LSTs of the grassland and woodland were affected by vegetation cover more than the cultivated and construction lands, and their cooling effect was more noticeable than that of the other land use types. Thus, in urban lands, appropriately increasing the amount of green area in parks is conducive to improving the city thermal environment.
Relationship between spatial structure of land use types and LST. Because single land use is rarely observed, the relationship between LST and NDVI should be analysed in the spatial combination of different land use types. In Fig. 5a, the LST-NDVI feature space shows a unique obtuse-angled triangle shape ABC. Concretely, A denotes water body, which had low LST and approximately negative NDVI; B represents dry and bare soils (mainly construction land), which had high LST and low NDVI; and C denotes woodland, which had the lowest LST and highest corresponding NDVI. The area ratios of the six land use types constituting the feature space of the study area was ranked as follows: woodland > farmland > grassland > construction land > water body > unused land. In Fig. 5a, the different types intervened with one another in the LST-NDVI feature space; the construction land, mainly mixed with unused land, had high evaporation and low water content; the grassland was distributed in the transition zone between the cultivated land and the woodland, and its vegetation coverage was high, soil moisture was relatively abundant and evapotranspiration was small. The main distribution of the woodland was relatively concentrated and followed an overall northwest-southeast direction (Fig. 5b).

Discussion
Abnormal high-temperature zones. In summer, the urban heat island effect of urban construction and industrial lands was remarkable in the study area. However, in addition to those caused by the conventional urban heat island effect, many high-temperature zones with different shapes and sizes were observed in the suburbs of the county in the high-definition satellite map. No construction was observed, and the highest LST (40.59 °C) of the study area was found there. The abnormal high-temperature zones were distributed mainly in the urban suburbs without vegetation cover and with exposed rocks (Fig. 6). In terms of the special geological environment of the study area, the karst area had a vulnerable calcium environment 27 . Barren and shallow soil leads to low vegetation coverage 28 , and the infiltration of surface water, water loss 29,30 , human disturbance and destruction cause serious soil erosion 31 , thus resulting in the exposure of bare rocks and rocky desertification 32,33 . Therefore, the specific heat capacity of the surface substances and the thermal inertia were small. Eventually, the LST of the area increased rapidly in the daytime and was higher than that in the surrounding region, thereby forming an abnormal high-temperature zone. Thus, LST of the karst area were featured not only by the conventional urban heat island effect but also by the abnormal high-temperature caused by karst rocky desertification. In addition, LST can characterise rocky desertification in karst areas to a certain extent. In future research, we can use LST as an important factor and index for analysing and evaluating the extent of rocky desertification.
Main factors affecting the LST of the karst mountain area. Yinjiang County has prominent mountain features with mountain accounting for 76% of the county. In Supplementary Fig. S1, the terrain tilts from southeast to northwest, whereas the LST decreases from northwest to southeast and extends along the mountain range, thereby indicating a controlling role played by altitude in the LST of the study area. And the quantitative relationship shows a negative correlation between altitude and LST. This pattern is caused by two main factors. (1) The LST affected by air temperature also decreases with elevation. (2) In karst mountain areas, the surface vegetation is relatively good. Thus, the solar radiation received by the ground surface mostly spreads in latent heat form, and the LST is low. Meanwhile, the intensity and length of solar radiation differed among different slope directions. In this study, the LST also differed among different slope directions, which is consistent with the results of Wen L. J. 14 . Different land use types have different thermal capacity, thermal conductivity, roughness and surface albedo, thus leading to differences in the LST. In the study, at the same altitude in the same slope direction, the LSTs of the different land use types still differed. In summary, the LST in this study area is influenced by levels of elevation, slope direction and land use type, among which altitude plays a fundamental controlling role in the overall pattern of LST.

LST-NDVI feature space. The existing research results on the relationship between LST and NDVI vary
considerably under combinations of different underlying surface types (Table 5).    . The LST and NDVI data in Fig. 5a and b were extracted from LST and NDVI images randomly by using ArcGIS 9.3 (http:// www.esri.com).   In this study, the LST and NDVI scatter diagram exhibited an obtuse-angled triangle distribution under combinations of different land use types, which is consistent with the results of Zhou Y. but different from results of other researchers. This finding could be attributed to the differences in the research methods or that between qualitative and quantitative analyses or the different soil backgrounds and sensor angles in the study area. In addition, numerous studies show that LST is positively correlated with NDVI in water bodies. Therefore, to verify whether the water caused the difference, data from water body sampling sites were removed (Fig. 7). After the removal of data from water bodies, the LST and NDVI exhibited a negative linear correlation (Fig. 7a), which is consistent with the results of Cao L. and others. Meanwhile, the LST-NDVI feature space formed a new triangle BDC (Fig. 7b), which is consistent with the findings of Price. The point D represents saturated bare soil with low LST and NDVI, and the land use type was less in this area than in the others. Meanwhile, BD represents bare edge with low vegetation coverage, whereas BC represents dry edge, which is mainly reflected in urban bare soil and dry land types 34 .

Conclusion
In this study, the atmospheric correction algorithm was used to retrieve LST from Landsat 8 data and the LST spatial distribution patterns and influencing factors in a karst area were summarised by combining land use types with elevation data. Then, the differences in the LST among different land use types were determined, and the spatial distribution characteristics of the LST and NDVI in the study area and their quantitative relationship were discussed. The main conclusions were as follows. (1) The LST retrieved by atmospheric correction from Landsat 8 data was consistent with the actual temperature. Influenced by topography and underlying surface types, the LST of the study area showed an overall downward trend from northwest to southeast. In summer, the LST of the study area had abnormal high-temperature zones, perhaps caused by karst rocky desertification. (2) Multiple comparative analyses indicated that the LST difference among most of the land use types was significant, with construction land having the highest LST of 32.25 °C and forest land having the lowest LST of 25.04 °C. (3) For the single land use type, the water body had low LST and NDVI, while the LST and NDVI of forest land, grassland, cultivated land and construction land exhibited a negative linear correlation. (4) For the spatial structure of land use in the entire study area, the LST-NDVI feature space was an obtuse-angled triangle. After the removal of water body data, a significant negative linear relationship between LST and NDVI was observed in the quantitative analysis.

Materials and Methods
Study area. Yinjiang Autonomous County is in the northeastern part of Guizhou Province in China, west of Tongren area (Fig. 8), and its geographical coordinates are 108°17′52″-108°48′18″ E, 27°35′19″-28°20′32″ N. Yinjiang Autonomous County has a land area of 1,969 km 2 and a population of 445,000 people. Located in the northeastern area of Guizhou Plateau, Fanjing Mountain, the main peak of Wuling mountains, stands in the east of the county. In the middle are King Pier, Liangzi Slope, Eling Off and other mountains, back-like protruding terrain from south to north, such that Yinjiang County is high in the east and south, convex in the middle and low in the west and north, and its relative elevation is relatively large (Fig. 8c). The landform type is mainly karst landform, and the terrain is broken. Yinjiang County has a subtropical warm and humid monsoon climate. The average annual temperature is 15 °C to 16.8 °C, and the annual precipitation is 1,057 mm to 1,258 mm. The main vegetation types are subtropical evergreen broad-leaved forest, secondary coniferous and broad-leaved mixed forest and plantation. In addition, the soil erosion, karst, and rocky desertification areas account for 69.75%, 51.8% and 20.70%, respectively, of the study area.

Materials.
The data used in this study include Landsat 8 data, Yinjiang County Administrative Boundary, DEM data and the meteorological observation data, The Landsat 8 remote sensing images were acquired in summer, specifically on August 29, 2016, and were used for inverting the LST and calculating the NDVI. The imgaes are from the US Geological Survey and have multispectral bands with 30 m spatial resolution and a  WGS84 coordinate. Without clouds over the study area, the image quality is good and the data product Level 1 T, having been geometrically corrected on the basis of terrain, can be used directly. The concrete data sources are shown in Table 6 and the technical flow chart can be found in Supplementary Fig. S3.
Land surface temperature retrieval. At present, the common algorithms for LST inversion that are based on remote sensing images are radiation conduction equation (atmospheric correction method) 35 , single-channel algorithm 36 , single window algorithm 37 and split window algorithm 38 . The atmospheric correction method can not only be used to calculate LST with only one thermal infrared band but can also consider surface emissivity and atmospheric radiation effects. Thus, the atmospheric correction method is more comprehensive than the other methods. In addition, three parameters, i.e. transmittance, atmospheric upward radiance brightness and The LST and NDVI data in Fig. 7a and b were extracted from LST and NDVI images randomly by using ArcGIS 9.3 (http://www.esri.com). atmospheric downward radiation brightness in this study can be obtained by the atmospheric correction parameter calculator issued by NASA(http://atmcorr.gsfc.nasa.gov), which imports the imaging time, longitude, air pressure of the area and other related information of the input image. It makes the atmospheric correction method easy to use. Therefore, this study used the atmospheric correction method to invert the LST. The concrete method can be found in the "Method Details of LST Retrieval" section in Supplementary Information, in which surface emissivity was calculated by using the NDVI threshold method proposed by Sobrino 39 . NDVI. We used the significant differences in the red and near-infrared reflectance spectra of green plants to obtain the vegetation index NDVI (Fig. 3). The value of NDVI ranged from −1 to 1. Generally, NDVI > 0 in the growing season indicates vegetation cover. An increase in the NDVI value indicates an increase in green vegetation. NDVI > 0.5 shows good vegetation growth status and large coverage density. The formula for NDVI is expressed in Eq. (1) 26 .