Assessing the Impact of Land Cover Changes on Surface Urban Heat Islands with High-Spatial-Resolution Imagery on a Local Scale: Workﬂow and Case Study

: Low-altitude remote sensing platform has been increasingly applied to observing local thermal environments due to its obvious advantage in spatial resolution and apparent ﬂexibility in data acquisition. However, there is a general lack of systematic analysis for land cover (LC) classiﬁcation, surface urban heat island (SUHI), and their spatial and temporal change patterns. In this study, a workﬂow is presented to assess the LC’s impact on SUHI, based on the visible and thermal infrared images with high spatial resolution captured by an unmanned airship in the central area of the Sino-Singapore Guangzhou Knowledge City in 2012 and 2015. Then, the accuracy assessment of LC classiﬁcation and land surface temperature (LST) retrieval are performed. Finally, the commonly-used indexes in the ﬁeld of satellites are applied to analyzing the spatial and temporal changes in the SUHI pattern on a local scale. The results show that the supervised maximum likelihood algorithm can deliver satisfactory overall accuracy and Kappa coe ﬃ cient for LC classiﬁcation; the root mean square error of the retrieved LST can reach 1.87 ◦ C. Moreover, the LST demonstrates greater consistency with land cover type (LCT) and more ﬂuctuation within an LCT on a local scale than on an urban scale. The normalized LST classiﬁed by the mean and standard deviation (STD) is suitable for the high-spatial situation; however, the thermal ﬁeld level and the corresponded STD multiple need to be judiciously selected. This study exhibits an e ﬀ ective pathway to assess SUHI pattern and its changes using high-spatial-resolution images on a local scale. It is also indicated that proper landscape composition, spatial conﬁguration and materials on a local scale exert greater impacts on SUHI.


Introduction
Against the background of extensive urbanization, cities and conurbations experience rapid expansions, especially in developing countries. Vast areas of rural land are changed into impervious surfaces such as buildings and roads, which generally absorb and re-radiate solar radiations effectively [1][2][3][4][5]. These easily-stored solar radiations and artificial heat generated from productions, transportations, and civilian life lead to severe urban heat islands (UHI) [6][7][8][9]. This phenomenon, as a chief corollary of urbanization, critically affects not only the urban environment and ecology but also human health [10][11][12]. Consequently, how to mitigate these impact of the UHI effect has become one of the hot research issues since the inevitability of urbanization [13].

Platform and Data Acquisition
In this study, an unmanned airship was used as a platform (Figure 2), which could lift a payload of up to 10 kg for a flight duration of up to 3 hours. It can automatically maintain the altitude and cruise according to a predefined route. The TIR imaging camera (TH9260) was responsible for taking TIR images that have one band with a spectral range of 8 to 13 µm. The digital SLR camera (Canon 450D with 18-55 mm F/3.5-5.6 IS) was used to obtain visible images. They were both mounted on a double three-axis platform. The favorably large payload and long duration enabled the simultaneous capture of visible and TIR images in one single flight. See Reference [35] for the integration of the whole system.  Outdoor data acquisition included the airborne observation and synchronous field measurement. The airship observing campaign at a height of 600 m was conducted between 12:00 and 12:30 on 23rd August 2012 and between 13:00 and 13:30 on 21st June 2015. The height is corresponding to the spatial resolution of visible and TIR images within 0.15 m and 0.8 m, respectively, as shown in Table 1. The airship observing method and principles were adopted from Reference [35]. To verify the accuracy of the LST retrieved, the typical surface temperatures were measured in synchrony with the airship observing campaign. In this field campaign, the contact temperature loggers (Ibutton and BES-1) with a circular planar contact sensor with an accuracy of 0.5 • C were used to measure the surface temperatures of trees, lawns, paddy water, water bodies, concrete pavements, asphalt pavements, building roofs, and soils at one-minute intervals.

Methodology
This study aims to determine the effective method to assess SUHI patterns based on LCTs and their changes on a local scale. The workflow is as follows (Figure 3). First, the visible and TIR images obtained by observations are mosaiced and orthorectified to visible and TIR orthophotos. Subsequently, the orthophotos are spatially co-registered to a high-resolution non-offset optical Google map. This is followed by the classification of the LCT from the visible map. Then, an LST map is retrieved from the TIR map based on the algorithm of TH9260. Finally, statistical analyses on LC, LST, and SUHI are performed. The following sections detail these phases and present the relevant data.

Pre-Processing
The hundreds of visible and TIR images need to be pre-processed, with regard to aspects including the mosaicking, geometric correction, and spatial co-registration. The mosaicking for both kinds of images was performed by the Agisoft PhotoScan software, which has been widely used in UAV image mosaicking [48]. For the visible orthomosaic, the processing with PhotoScan was twofold. Firstly, the software was deployed to identify and match the feature points between images with the Scale Invariant Feature Transform (SIFT) algorithm, and to reconstruct the scene geometry with Structure from Motion (SfM) techniques [49]. Then, it applied a dense multi-view reconstruction to the aligned images to establish the scene geometry, and finally exported an orthophoto [50]. The whole processing was automatic and involved no human interactions. The TIR images obtained from the TH9260 camera can be identified only by its companion software "Image Processor Pro II"; accordingly, the ASCII format file with the temperature information of each pixel had to be exported and converted into the Tiff format image through the IDL program in the ENVI software. For the TIR images, the subsequent mosaicking was performed in a manner similar that for the visible images [51].
The other pre-processing step is to perform spatial co-registration for the visible and TIR orthophotos, which is necessary for subsequent analyses on LCT and LST. We chose the non-offset optical Google map with geographic coordinates in the study area as the reference. The visible and TIR orthophotos in 2012 and 2015 were loaded into ArcMap and co-registered to the Google map with twenty as possible as evenly-distributed, clear, and discernible feature points e.g. intersections of roads and corners of houses in all maps. The second-order polynomial transformation was selected for this process in ArcMap.

Classification and Accuracy Assessment
Given the central role of the algorithm in ensuring successful classification, five supervised classification algorithms were deployed in the ENVI software [52]: Maximum likelihood, parallelepiped, minimum distance, Mahal distance, and neural network classification algorithms [53]. Based on land-use characteristics in the study area and spectral characteristics of the visible image, LC was identified and classified into six types: Dry farmland (including crops and lawns), paddy field, woodland, bare land, built-up land, and water body. After the classification, accuracy assessments were performed by comparing the consistency of validation data (over 500 checkpoints for each LCT) and the classification results for the different algorithms. Overall accuracy and kappa coefficient, as well as user's accuracy and producer's accuracy, were used as quantitively assessing indexes [54]. Finally, each LCT was extracted from the classified maps for discerning the changes in the areas during the period of interest in ArcMap.

LST Retrieval and Accuracy Assessment
Due to functional limitations of the TH9260, traditional satellite remote sensing methods cannot be directly applied. The LST retrieval has to be performed based on the corresponding function in the "Image Processor Pro II" software, which provides a retrieval algorithm based on the parametric correction of atmospheric transmittance, mean atmospheric temperature, and emissivity. This process is fourfold. First, we retrieved LST from the brightness temperature (BT) in the images (including a representative LCT) by correcting these aforesaid parameters. Next, we extracted the BT and LST of some corresponding sample pixels to regress their linear relationships. Then, we used the mask to extract the TIR maps for each LC according to the previous LC classification. Finally, we engaged spectral calculation to obtain the LST maps of each LCT based on the linear relationship between BT and LST, and then mosaicked them to the whole LST maps for 2012 and 2015.
The atmospheric correction was performed with the MODTRAN 4 model. In this paper, the settings of input parameters were in place [55]: the model atmosphere was selected as tropical; the atmospheric path was the slant path between two altitudes; the aerosol model was set to Rural Extinction; the default visibility was 5 km; the seasonal aerosol profile was set to Spring-Summer; the weather was fine, and no cloud or rain. The range of wavelengths was from 8 µm to 13 µm, and the wavelength space was set to 0.1 µm. The atmospheric temperature at the first boundary took the observed mean value during the airship's campaign at flight altitudes (29 • C in 2012 and 28.5 • C in 2015). Assuming the spectral response functions on one broad spectral band of the TH9260 according to a normal distribution, the mean atmospheric transmittances were 0.64 in 2012 and 0.68 in 2015.
Given the lack of in situ measurement of the emissivity of LCTs and non-mixed pixels in our study, reference values released by the MODIS UCSB Emissivity Library were used to retrieve LST [56,57]. The mean values of emissivity (Table 2) are as follows: dry farmland (0.95), paddy field (0.97), woodland (0.96), barren land (0.96), built-up land (0.95), and water body (0.99). After LST retrieval, the surface temperatures recorded synchronously with the in-situ measurements were used to evaluate the accuracy of the retrieval results.

Thermal Indexes Based on LC
In this section, the thermal effect contribution index (H i ), thermal pixels proportion index (D 1 ), and regional thermal pixels proportion index (D 2 ) were used to discuss the impact of various LCTs on the urban thermal environment.
H i indicates the impact on the regional mean LST caused by LCTs. For calculating H i , we first aggregated the differences between the LSTs of the pixels on each LC higher than the regional mean LST as the initial thermal effect contribution index (H i '), and then normalized H i ' to compare them among various LCTs. H i ' and H i are respectively calculated by Equation (1) and Equation (2), (1) where T ij represents the LST of pixel j on the land cover i more than the regional mean LST, T a0 the regional mean LST of the study area, n i the pixel numbers on the LCT i more than the regional mean LST. A greater H i implies a more significant thermal effect of the LCT. On the other hand, D 1 and D 2 are respectively represented by: where D 1 represents the ratio of the area with LST above the mean in an LCT, D 2 the ratio of the area with LST above the mean in the total study area, and N i the pixel numbers on the LCT i, and N the total pixel numbers in the study area.

Normalized LST Level Map
Given differences in climate background, it will be challenging to directly compare LST maps across different years [58]. Accordingly, the normalized LST (NLST) could maintain the spatial distribution of LST and eliminate the influence of varied climate background, providing an objective comparison [59]. The LST maps need to be normalized by Equation (5) [60]: where LST i ' is the normalized LST of i pixel, LST i the LST of i pixel, LST max and LST min respectively the maximum and minimum LST values in a map. Such normalization uniformly rescales the distribution of LST to the 0-and-1 range, as well as reduces the climatic differences. Then, we classified the NLST into different thermal field levels (from 3 to 6) through the equalization method and the mean value-standard deviation (STD) method (Table A3), and chose the preferred method and numbers of corresponding thermal field levels based on the results (Figures A1-A4) [61]. Subsequently, the area of each level and SUHI pattern were calculated and investigated during the period of interest. The primary purpose of the classification of the NLST was to understand the changes in the statistical distribution of temperature data and mitigate SUHI.

SUHI Ratio Index Computation
In order to further study SUHI variations between 2012 and 2015, the surface urban heat island ratio index (SUHIRI) is computed, as follows: where m represents the level number in the previous normalized step, n the number of the levels higher than the mean temperature, Wi the weight of the levels higher than the mean temperature (the order value of the level in m, for example , the weight value of high and extra-high levels are 5 and 6, respectively), and P i the area proportion of the SUHI level of i in the study area [59,61].

Results of Mosiacking and of Spatial Co-Registration
With the methods in Section 2.1, we obtain the co-registered visible and TIR mosaic maps in the study area for 2012 ( Figure 4) and 2015 ( Figure 5). For each mosaic map, objects such as various constructions and even their material characteristics are apparent and sharply-defined in the visible map, while the edges of the corresponding objects and their TIR distributions are obvious in the TIR map.

250
To assess co-registration accuracy, we checked the errors of some feature points between the 251 Google map and the mosaic maps. The co-registration error and the root mean square error (RMSE)

252
of each map are summarized in Table 3. The values of the RMSE for the visible maps were 0.74 pixel 253 in 2012 and 0.80 pixel in 2015, while those for the TIR maps were correspondingly 2.33 pixels and 254 1.86 pixels. It is reasonable that the co-registration accuracy of the visible map with higher spatial 255 resolution and more obvious object features was higher than that of the TIR map. The co-registration 256 accuracy of both kinds of maps is regarded as satisfactory for LC classification and further analyses 257 of the thermal environment.

250
To assess co-registration accuracy, we checked the errors of some feature points between the 251 Google map and the mosaic maps. The co-registration error and the root mean square error (RMSE)

252
of each map are summarized in Table 3. The values of the RMSE for the visible maps were 0.74 pixel 253 in 2012 and 0.80 pixel in 2015, while those for the TIR maps were correspondingly 2.33 pixels and 254 1.86 pixels. It is reasonable that the co-registration accuracy of the visible map with higher spatial 255 resolution and more obvious object features was higher than that of the TIR map. The co-registration 256 accuracy of both kinds of maps is regarded as satisfactory for LC classification and further analyses 257 of the thermal environment.  To assess co-registration accuracy, we checked the errors of some feature points between the Google map and the mosaic maps. The co-registration error and the root mean square error (RMSE) of each map are summarized in Table 3. The values of the RMSE for the visible maps were 0.74 pixel in 2012 and 0.80 pixel in 2015, while those for the TIR maps were correspondingly 2.33 pixels and 1.86 pixels. It is reasonable that the co-registration accuracy of the visible map with higher spatial resolution and more obvious object features was higher than that of the TIR map. The co-registration accuracy of both kinds of maps is regarded as satisfactory for LC classification and further analyses of the thermal environment.  We used the visible maps in 2012 and 2015 to assess the classification accuracy from different classification algorithms. From Table 4, it is evident that the maximum likelihood method has yielded the best result on the LC classification maps. For 2012 and 2015, the values for the overall accuracy were respectively 92.88% and 90.32%, while those for the Kappa coefficient were respectively 0.91 and 0.88 (Appendix A: Tables A1 and A2 detail the classification results in 2012 and 2015 from the maximum likelihood method).

Patterns of LC Classification in 2012 and 2015
The results of the LC classification for 2012 and 2015 from the maximum likelihood algorithm are shown in Figure 6. It is readily discernible that areas of the built-up land and bare land have expanded rapidly from 2012 to 2015, at the expense of other LCT, especially woodland. Other observations are that the changes were mainly concentrated in the southwest and northeast areas, and that the changed areas were generally less than the unchanged ones.
The area proportion of each LCT was calculated, as shown in Table 5. From Figure 6 and Table 5, the order of the areas of the different LCTs did not change though the areas did vary drastically. The most dominant LCT was woodland, which also exhibited the most significant decline from 55.66% of the study area in 2012 to 35.49% in 2015. Built-up land and bare land were the two fastest-growing types, respectively expanding from 9.36% to 18.19% and from 17.87% to 29.60%, reflecting the robust and rapid development during these three years of analysis. Furthermore, the substantial growth of bare land earmarked for future constructions implies that the pace and extent of development would only continue and intensify in the future. For dry farmland and paddy field, their areas have witnessed only slight changes, thereby underlining the excellent protection bestowed upon agricultural land. Note that the increased proportion of dry farmland (3.58%) does not translate into the reclamation of new agricultural land. Two reasons underlie this, the first of which concerns the growth of lawns, which had hitherto not been classified separately due to their small total area and similar spectral properties in visible imageries. For example, the artificially-grown lawns near the Phoenix Lake (a large artificial lake) in 2015 were identified as dry farmland. The other reason concerns the variety of crop types and harvest. A little field without any crops but reflective water on the northwest side of the study area (Figure 6a) was classified as built-up land in August 2012; however, it was later re-classified as dry farmland with crops in June 2015. Additionally, the Phoenix Lake which appeared in 2015 did not lead to an increase in the area of the water body (from 4.84% to 3.05%) due to the disappearance of small-spread water bodies elsewhere.  To further characterize the temporal and spatial evolutions of LC in the study region, change 296 detection analysis was used to determine the quantities of conversion from a particular LCT to 297 another. A cross matrix (Table 6)      To further characterize the temporal and spatial evolutions of LC in the study region, change detection analysis was used to determine the quantities of conversion from a particular LCT to another. A cross matrix (Table 6) provides information regarding the transition between LCTs. For the woodland, around 28.48% in 2012 was converted into bare land, 11.84% into built-up land, and 10.12% into dry farmland, illustrating that the major source of the developed and developing areas was founded on converted woodland. For the built-up land, around 21.20% was converted into bare land, meaning that existing constructions were demolished and that planning was underway for earmarking the land for new uses. Additionally, 14.61% and 13.61% of the built-up land were turned into dry farmland (most of them were lawns) and woodland respectively, thereby highlighting the importance of ecological protection during the development. For the bare land (representing the area that had been planned for construction in our case), it was converted into built-up land (24.3%), woodland (14.15%), and dry farmland (12.92%), suggesting a balance between urban development and ecological preservation.

Assessment of LST Retrieval
The measured real surface temperature as reference points and the corresponding LST values for comparison are shown in Table 7. The temperature deviations between the retrieved and measured temperatures varied from -0.45 • C to 3.08 • C, and the mean of absolute deviation and RMSE were 1.63 • C and 1.87 • C, respectively. In general, water body, woodland, dry farmland, and paddy field have higher temperature accuracy than built-up land (building roof, asphalt pavements, and concrete pavement), given the relatively stable emissivity of the former four LCT and their temperature ranges that were similar to the ambient temperatures. The temperature deviations of built-up land ranged from 1.65 • C to 3.08 • C, mainly because the emissivity of the man-made impervious surfaces was set to be the same in this study. The broad emissivity of some pavement made by concrete may in fact vary from 0.85 to 0.96 in the database published by the MODIS UCSB Emissivity Library [56,57]. The temperature deviation of water in the Phoenix Lake was -0.45 • C, while the value of pool (small water body) was 1.63 • C, which is because of the much difference in their depth and volume. LST maps and their associated descriptive statistical results in 2012 and 2015 are shown in Figure 7 and Table 8, respectively. The results ( Figure 7) reveal a generally increasing trend in the LST distributions from 2012 to 2015, as is consistent with the expansion of the built-up land in the study area (Figures 4 and 5). The LST of the impervious areas such as the built-up land, especially roads, is noticeably higher than that of the pervious areas, such as the water body, woodland and so on, implying the impacts of construction activities of human being on the thermal environment.    anthropogenic activity [11,62]. Dry farmland with the second-highest LSTmean-which was even 342 slightly higher than that of bare land and much higher than that of woodland-implied that dry  ). Secondly, the LST of water body depends on its size, depth, and fragmentation 351 [63]. In 2012, the water bodies were small, shallow, and scattered. (Table 6)  From Table 8, the maximum LST mean in 2012 was registered by the built-up land, followed by dry farmland, bare land, paddy field, water body, and finally by woodland, while the maximum LST mean in 2015 was registered by the built-up land, followed by dry farmland, bare land, woodland, paddy field, and finally by water body. We noted that the LST mean orders were the same for all LCTs except woodland, which exhibited the lowest LST mean (26.99 • C) in 2012 but the third lowest one (36.69 • C) in 2015. The values of LST mean on built-up land in 2012 and 2015 were both much higher than those of others because of its impervious characteristics such as thermal properties and anthropogenic activity [11,62]. Dry farmland with the second-highest LST mean -which was even slightly higher than that of bare land and much higher than that of woodland-implied that dry farmland with low coverage did not contribute to cooling despite its vegetation. The values of LST mean of woodland in 2012 and of water body in 2015 were the lowest due to active transpiration and to evaporation, respectively. The difference in the order was attributed to two reasons. Firstly, the LST of woodland depends critically on transpiration and evaporation. For these two natural processes, the appreciable variations in the study area could be attributed to: Differences in the transpiration intensity (caused by the variation on the leaf area index (different years and months)); and the diverse water contents of soil in the woodland (caused by the atmosphere (rainy weather two days before observation in 2012)). Secondly, the LST of water body depends on its size, depth, and fragmentation [63]. In 2012, the water bodies were small, shallow, and scattered. (Table 6). Conversely, in 2015, most of them had been converted into other LCTs. And the vaster and deeper Phoenix Lake appeared, leading to the lowest LST. Moreover, the difference in the values of the LST mean of all LCTs up to 10.29 • C between 2012 and 2015 implied the necessity of normalization for the LST change pattern analysis.
From LCTs, the built-up land, bare land, and woodland had more significant STDs due to more complex compositions in their respective fields, whereas the water body and paddy field were contrary. The dry farmland showed much different STDs (1.87 • C in 2012 and 2.67 • C in 2015) due to the massive increase in the areas of lawns in 2015, whose properties somewhat differed from crop plants. From 2012 to 2015, STD of each class varied with the area, besides for water body and built-up land. The STD of water body increased from 1.23 • C to 1.42 • C, while the area decreased from 4.84% to 3.05% (Table 5). This was contributed to the appearance of the Phoenix Lake, which differed prominently from the previously small and shallow water bodies. The STD of built-up land decreased with the increase of its area, resulting from the rise in the complexity of built-up land to a certain extent by the heavy rain two days before observation in 2012.

Thermal Contribution from Various LCTs
The results of H i , D 1 , and D 2 were calculated ( Table 9). The highest H i was registered by built-up land, followed by bare land, dry farmland, and woodland. The values of H i of built-up land up to 48.11% and to 60.35% were much higher than those of the above three types, meaning that its contribution to the thermal environment was much more significant than the LST pattern shown. Conversely, the lowest values of H i were registered by paddy field and water body (2.49% and 2.48% in 2012, and 0.18% and 0.33% in 2015, respectively), illustrating that almost all the LSTs of their pixels were below the regional mean LST, which was positive to reduce SUHI. This result also differed somewhat from the finding from the LST pattern: the LST of woodland (26.99 • C) was significantly lower than that of paddy field (29.55 • C) and slightly lower than that of water body (27.71 • C) in 2012. It implied that the use of H i minimizes the impact of weather-related factors and is more suitable for assessing the effects of LCT on the thermal environment than the LST pattern. Furthermore, in contrast to the increasing of H i with the area growth of the other LCTs, the decreasing H i of bare land from 24.66% to 16.13% with its growing area from 17.87% to 29.60%, illustrating that the LCT contributed more than the area. The values of D 1 of built-up land were highest in all LCTs, followed by bare land and dry farmland, showing dominant effect exerted by built-up land on the thermal environment. The values of D 1 of built-up land was nearly similar in 2012 and 2015, whereas those of bare land and of paddy field drastically dropped from 94.38% to 68.47%, and from 57.08% to 7.56%, respectively. The main reason was that the sharp growth of built-up land by 8.83% of the study area directly increased the regional mean LST, which was higher than the LST in some of the bare land and most of the paddy field. Additionally, for paddy field, another reason was noteworthy: the LST of woodland with dominant areas (55.66% in 2012 and 35.49% in 2015) was lower than that of paddy field in 2012, but the reverse was observed in 2015.
The value of D 2 of bare land was the highest, instead of built-up land. All values of D 2 had positive relationships with their LC areas, illustrating that D 2 is more sensitive than H i for the area.

Applications of NLST Level on a Local Scale
We used the LST map in 2015 to assess the applicability of NLST classification on a local scale. Figure 8 depicts the proportion of the area and the intensity of SUHI with the different methods (Appendix A: Table A3 defines the relationship between level numbers (LNs) and SUHI; Figures A1-A4 detail the NLST level maps with the different methods). The SUHI area (from 2.46% to 7.72% along with 3-6 levels) classified by the equalization method was much smaller than that by the mean value-standard deviation method, which obviously did not match the reality by checking the NLST level maps (Figures A1-A4), indicating that the equalization method is inappropriate for the high-spatial LST map on the local scale. For the mean value-standard deviation method, the area and intensity of SUHI as LNs 4 and 6 were the same whenever the STD multiple was x = 1, 1.5, or 2, while the results of LNs 3 and 5 were fluctuated and were sensitive. From Figures A1-A4, with increasing LNs, the details of the NLST level map were more abundant, and could more comprehensively reveal the distribution characteristics of SUHI. Accordingly, LN 6 was more suitable for our case. For LN 6, when the STD multiple was x = 1, the extra-high temperature areas (accounting for 30.08% of the study area) consisted of all built-up land and some soil land, while the low temperature areas (31.46% of the study area) consisted of all water body, paddy field and most of the woodland. Apparently, it is unreasonable that most areas were concentrated on the lowest and highest two levels. And it was inconsistent with the LST distribution range of the various LCTs. When the STD multiple was x = 2, the demarcation between the low and sub-medium temperature areas was within the LST range of the Phoenix Lake; this finding likewise was unreasonable. Lastly, when the STD multiple was x = 1.5, the low-temperature areas comprised water body and dense woodland, whereas high-temperature areas comprised mainly built-up land. Additionally, the profile of the object was clearer than when x = 1. Part of the bare land belonged to the sub-medium temperature area, and the SUHI area included most of the built-up land, part of the bare land and dry farmland, as was consistent with the actual situation. Therefore, we chose the NLST level map classified by the mean value-standard deviation method with six levels and x = 1.5 to analyze the SUHI of the study area ( Figure 9).
(Appendix A: Table A3 defines the relationship between level numbers (LNs) and SUHI; Figures A1-

397
A4 detail the NLST level maps with the different methods). The SUHI area (from 2.46% to 7.72% 398 along with 3-6 levels) classified by the equalization method was much smaller than that by the mean 399 value-standard deviation method, which obviously did not match the reality by checking the NLST 400 level maps (Figures A1-A4), indicating that the equalization method is inappropriate for the high-

427
To quantitatively analyze the spatial structure of the SUHI, the area Proportion of each level and 428 the SUHI intensity were calculated for both years (Table 10). According to

Patterns of NLST Level and Changes in SUHI
The NLST level maps with LN 6 for 2012 and 2015 are shown in Figure 9. By inspection, two findings are evident: the SUHI region-comprising extra-high, high, and sub-high temperature areas-were located on built-up land, some bare land, and a little dry farmland in both years; and the extra-high temperature area witnessed a noticeable expansion from 2012 to 2015.
To quantitatively analyze the spatial structure of the SUHI, the area Proportion of each level and the SUHI intensity were calculated for both years (Table 10). According to Table 10, the area proportions of the low, sub-medium and medium temperature levels had changed little during these three years, while the proportion of the SUHI area increased slightly from 43.81% to 44.94% of the study area. However, the SUHI intensity increased by 3.32 • C (from 7.22 • C to 10.54 • C), which was attributed to the growth of the extra-high temperature area from 14.76% to 24.62%, with the expansion of built-up land.   Based on Table 10 and the cross matrix of NLST levels (Table 11)   Based on Table 10 and the cross matrix of NLST levels (Table 11), the inference can be drawn that the increase of area of the extra-high temperature level in 2015 arose mainly from the conversion from areas of the sub-medium, medium, sub-high and high temperature levels (each with a conversion ratio exceeding 20%) in 2012. Such a finding appears to imply that the development of built-up land would inevitably increase the extra-high temperature area and lead to environmental deterioration. However, mutual transformation is noted between SUHI areas and non-SUHI ones, as reflected by the conversions of extra-high temperature areas (24.05% into medium temperature areas and 18.07% into sub-high temperature) and of high temperature areas (29.58% into medium temperature areas and 18.99% into sub-high temperature). In other words, we can reasonably mitigate the SUHI intensity or control its direction and range through judicious urban planning, including to control the proportion among various LCTs and their spatial distributions. For instance, the Phoenix Lake as an excellent cold source was built in 2015 from the bare land in the central region instead of the disappeared sporadic and small water bodies [63].
For Equation (6), the levels of SUHIRI were as same as the NLST levels: Low, sub-medium, medium, sub-high, high, and extra-high. Among these, sub-high temperature, high temperature, and extra-high temperature areas were considered as the SUHI regions, meaning that n = 3 and W i = 4, 5, 6 respectively. The weight of each level equaled the row number of the level. P i indicated the area proportion of level i. The SUHIRI was used to analyze changes of SUHI intensity during the period of interest, which were 0.35 for 2012 and 0.39 for 2015, revealing the increasing trend of the SUHI.

Discussion
The primary purpose of this study is to explore the LCT's impact on SUHI on a local scale by a visible camera and a TIR camera. We first used the visible imagery (RGB channel) to classify six LCTs and obtain the sufficient classification accuracy in 2012 and 2015, according to the recommended 85% to 90% overall accuracy [53]. Limited by the property of the RGB channel, it is hard to classify more LCTs without other wavebands and geographic data [64]. For instance, woodland (with its relatively dark color, texture, and shadow) and dry farmland (with its light and uniform color) can be easily distinguished. Conversely, dry cropland and lawns are hardly distinguishable due to their similar texture and the same spectral feature on the RGB channel. They were therefore considered in one LCT in our study.
From Tables 7 and 8, the difference of LST mean between heat and cooling sources is more significant than the mean of absolute deviation and RMSE. It is very reasonable that the retrieved LST is are acceptable for assessing LST patterns. Moreover, the problems of LST retrieval in our study are different from in satellite remote sensing field. Firstly, LST on the map with the higher spatial resolution usually exhibits more variabilities, especially for the built-up land [35,39]. Therefore, the error control of LST retrieval in our study tends to be relatively challenging. Secondly, more LCTs with accurate emissivity are able to obtain more accurate LSTs in a high-spatial-resolution map. However, it is inconvenient to measure the real value of emissivity for every surface in practice and hard to classify more LCTs only using the RGB channel. Fortunately, our results are acceptable based on the data from the MODIS UCSB Emissivity Library [57]. Finally, the atmospheric transmittance calculated from MODTRAN is easily influenced by atmospheric parameters such as the absorption coefficients of atmospheric molecules and of atmospheric aerosols in satellite-based LST retrieval [65,66], while the impact of atmospheric transmittance is much smaller at the shooting height of around 600 m.
We introduced several indexes widely used in satellite remote sensing to assess the thermal environment and SUHI, and their application on a local scale. A decrease in the LST was registered by built-up land, dry farmland or bare land, paddy field, and water body or woodland (the descending order was in agreement with most results from satellites [67][68][69]). While the LSTs between different LCTs exhibited more significant differences than the result from the satellites due to the presence of mixed pixels [70,71]. Similarly, the higher spatial resolution shows more details in surface materials and their properties, revealing why the RMSE of LST in each LCT derived from satellites are generally smaller in summer daylight [39,72]. For an LST map with high spatial resolution, the LST is determined mainly by the thermal properties of materials such as thermal conductivity and heat capacity. The range of its fluctuations is closely related to the heterogeneity of the materials. The built-up land usually consists of numerous surface materials (such as concrete, asphalt, and metal) and has greater thermal capacity and conductivity. Therefore, the highest mean LSTs and the larger fluctuations of LSTs (most in 2012 and second most in 2015) are found on built-up land. The wider gap of LSTs between different LCTs shows the potential that the local thermal environment can be improved by utilizing the heat transfer between the different LCTs [73]. It is not clearly illustrated by satellites. Accordingly, compared with the application of satellite remote sensing on an urban or city scale, low-altitude remote sensing is more suitable for assessing the impacts of landscape composition, spatial configuration, and even materials on a local or micro scale. For example, paying attention to the importance of small cold sources can improve the micro thermal environments [74,75]. More specifically, the proper compositions between built-up land and water body or woodland need to be ensured and the cool pavements or roofs need to be used [76][77][78]. In the meantime, the difference of LSTs and their fluctuated ranges caused by different resolutions demonstrates that it is improper to directly compare their SUHI intensities, even with the same time and place. This has been presented by Sobrino (STD of LST is from 4.4 • C to 0.8 • C with spatial resolution from 4 m to 1000 m in the same district) [34].
The LST and LC classification maps depict the patterns of LST and LC and their relationship but do not accurately appear the thermal contribution of the various LCTs. Then, LCT-based indexes i.e. H i , D 1 , and D 2 are demonstrated to be applicable to our local map. H i highlights the LCT's impact on the SUHI: the built-up areas of 9.36% and 18.19% contributed to SUHIs of 48.11% and 60.35%, respectively. To compare the impacts of LC changes on SUHI in the situation with high spatial resolution, we attempted to classify the LSTs with the equalization method and the mean value-standard deviation method, as well as different STD multiple (x = 1,1.5, and 2). We finally determined the most suitable combination (the mean value-standard method with x = 1.5 and LN 6) in our study area. Similar to the LST map, the NLST level map is closely related to the LCs and their shapes, which indicates the importance of LCs together with their compositions and materials. In addition, the role of the woodland and water body as the high-quality cold source is especially highlighted. Woodland and water body were both located on the non-SUHI areas, whereas some of the dry farmland as a kind of vegetation was situated on the high temperature area, even the lawns beside the Phoenix Lake in the extra-high temperature area in 2015. It indicates the differences between woodland and lawns in mitigating SUHI. It also allows us to elucidate the principle for regulating the microclimate on the local scale: more trees with intense evaporation such as arbors, rather than large areas of lawns. The change of landscape composition, spatial configuration, and materials tends to be crucial for the improvement of the civilians' thermal comfort and reduce their heat stress.

Conclusions
We studied the surface thermal characteristics and its spatial and temporal changes in the Sino-Singapore Guangzhou Knowledge City in 2012 and 2015. Meanwhile, we examined the workflow and assessed the accuracy of LC classification, LST retrieval, and applicability of some indexes for evaluating SUHI on a local scale.
Our results reveal that the common digital camera and TIR camera mounted on an unmanned airship can conveniently and affordably capture fine and high-spatial-resolution images, thereby enabling detailed characterization of visible and TIR images for complex landscape composition and spatial configuration.
Several detailed conclusions are drawn: (1) The supervised maximum likelihood method yielded satisfactory results in classifying the high-spatial-resolution imagery; the overall accuracy was 92.88% and 90.32% in 2012 and 2015, as well as the Kappa coefficient was 0.91 and 0.88. Whereas, cropland and lawns were indistinguishable in the RGB channel. (2) The deviation of the retrieved LST and RMSE indicate that the retrieval results are acceptable for assessing LST distributions and further SUHI patterns.
(3) The LSTs demonstrated greater consistency with LCTs and more fluctuation within an LCT on a local scale than on an urban scale. Built-up land contributed to the highest SUHI intensity, while woodland and water body had the most cooling contribution. Additionally, more significant discrepancies in LST were found between built-up area and woodland on a local scale, indicating that proper landscape composition is an excellent way to mitigate SUHI and improve thermal comfort in the local microclimate. (4) Indexes from the satellite field such as the H i and SUHIRI are applicable for a local scale, but the application of the NLST level map necessitates the judicious selection of parameters. (5) Low-altitude remote sensing data can provide a reference for local regulator planning, in comparison to satellite remote sensing data, which provides a reference for urban planning. More specifically, for the local scale, landscape composition, spatial configuration, and even materials are of greater concern.
Furthermore, the following works are envisioned: 1) To continue observing this study area to check the SUHI change whether the aim of urban planning is achieved upon completion of the development; 2) and to further investigate the relationship between the parameters (e.g., landscape composition, spatial configuration, and materials) and SUHI. Acknowledgments: The authors would like to thank Jianer Jiang for his help in the control of our unmanned airships, Yufeng Zhang for his help in writing advice, Xu Yuan for his help in the data acquisition, and Qianlong Qi for his help.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
Tables A1 and A2 show the accuracy assessment of LC classification by maximum likelihood algorithm in 2012 and 2015.   Sub-high -T mean ≤ T< T mean + STD T mean + 0.5 STD≤ T< T mean + STD T mean ≤ T< T mean + STD High T≥ T mean + 1.5 STD T≥ T mean + STD T≥ T mean + STD T mean + STD≤ T< T mean + 1.5 STD Extra-high ---T≥ T mean + 1.5 STD the mean value-standard deviation method: x=1.5 Low T< T mean -1.5 STD T< T mean -1.5 STD T< T mean -1.5 STD T< T mean -1.5 STD Sub-medium --T mean -1.5 STD≤ T< T mean -STD T mean -1.5 STD≤ T< T mean -STD Medium T mean -1.5 STD≤ T< T mean + 1.5 STD T mean -1.5 STD≤ T< STD T mean -STD≤ T< T mean + STD T mean -STD≤ T< T mean Sub-high -T mean ≤ T< T mean + 1.5 STD T mean + STD≤ T< T mean + 1.5 STD T mean ≤ T< T mean + STD High T≥ T mean + 1.5 STD T≥ T mean + 1.5 STD T≥ T mean + 1.5 STD T mean + STD≤ T< T mean + 1.5 STD Extra-high ---T≥ T mean + 1.5 STD the mean value-standard deviation method: x = 2 Low T< T mean -2 STD T< T mean -2 STD T< T mean -2 STD T< T mean -2 STD Sub-medium --T mean -2 STD≤ T< T mean -STD T mean -2 STD≤ T< T mean -STD Medium T mean -2 STD≤ T< T mean + 2 STD T mean -2 STD≤ T< STD T mean -STD≤ T< T mean + STD T mean -STD≤ T< T mean Sub-high -T mean ≤ T< T mean + 2 STD T mean + STD≤ T< T mean + 2 STD T mean ≤ T< T mean + STD High T≥ T mean + 2 STD T≥ T mean + 2 STD T≥ T mean + 2 STD T mean + STD≤ T< T mean + 2 STD Extra-high ---T≥ T mean + 2 STD