Exploring the Impacts and Temporal Variations of Different Building Roof Types on Surface Urban Heat Island

: This study examined the impact of different types of building roofs on urban heat islands. This was carried out using building roof data from remotely sensed Landsat 8 Operational Land Imager (OLI) and Thermal Infrared Sensor (TIRS) imagery. The roofs captured included white surface, blue steel, dark metal, other dark material, and residential roofs; these roofs were compared alongside three natural land covers (i.e., forest trees, grassland, and water). We also collected ancillary data including building height, building density, and distance to the city center. The impacts of various building roofs on land surface temperature (LST) were examined by analyzing their correlation and temporal variations. First, we examined the LST characteristics of ﬁve building roof types and three natural land covers using boxplots and variance analysis with post hoc tests. Then, multivariate regression analysis was used to explore the impact of building roofs on LST. There were three key ﬁndings in the results. First, the mean LSTs for ﬁve different building roofs statistically differed from each other; these differences were more signiﬁcant during the hot season than the cool season. Second, the impact of the ﬁve types of roofs on LSTs varied considerably from each other. Lastly, the contribution of the ﬁve roof types to LST variance was more substantial during the cool season. These ﬁndings unveil speciﬁc urban heat retention drivers, in which different types of building roofs are one such driver. The outcomes from this research may help policymakers develop more effective strategies to address the surface urban heat island phenomenon and its related health concerns.


Introduction
In recent decades, several developed and developing countries have undergone urbanization and industrialization. This rapid development in urban areas has led to an environmental issue known as the surface urban heat island (SUHI), which affects more than half of the urban population [1][2][3][4]. SUHI is a phenomenon in which urban areas have higher surface temperatures than their suburban/rural counterparts [5][6][7][8]. Land surface temperature (LST) is a significant parameter used to explore the exchange of surface matter, surface energy balance, and surface physico-chemical processes. Thus, it is one of the most essential sources used to characterize SUHIs [9]. LST represents the energy exchange between the land and atmosphere, an important geophysical parameter in ground-air systems [10].
Guangzhou and Foshan are cities in Guangdong Province, China (22 • 56 N to 23 • 93 N and 112 • 38 E to 114 • 05 E), both of which were selected as the study area ( Figure 1). They are located in southern China and the northern part of the Pearl River Delta, spanning a combined area of 11,231 km 2 . Guangzhou and Foshan experience a humid subtropical climate affected by the East Asian monsoon. Summers are hot and humid, while winters are warm and less humid. Guangzhou and Foshan are core regions of the Guangdong-Hong Kong-Macao Greater Bay Area (GBA), with the highest population density in Guangdong Province and being the national center of South China. By the end of 2018, the population of these two cities exceeded 22 million.

Data Sources
Imagery from the Landsat 8 Operational Land Imager (OLI) and Thermal Infrared Sensor (TIRS) with a spatial resolution of 30 m were used as the data sources in this study ( Figure 2). The data were recorded during the hot and cool seasons; the specific dates and corresponding temperature ranges are listed in Table 1 [20]. No cloud-free images were available during the hot seasons in 2014 and 2020. Thus, images recorded during the hot season from previous years (i.e., 2013H and 2019H) were used as substitutes for 2014H and 2020H (Table 1). High spatial resolution imagery acquired from Google Earth historical images (for 2014, 2016, 2018, and 2020) and GaoFen-2 imagery (for 2016 and 2019) were used to assist in sample selection. Building outlines containing building height information were used to obtain data on the building density and determine the average building height. The building outline data were supplied from a geographic information system (GIS) company (ZhongKe Tuxin Co., Ltd. Location: Suzhou Science and Technology City Chinese Academy of Sciences Geographic Information and Cultural Science and Technology Industry Base; http://www.tuxingis.com//index.html, accessed on 12 November 2020). As these building outlines were only available for urban areas, there were no samples taken from suburban or rural areas. All remotely sensed data were pre-processed with calibration and atmospheric correction, and all data were re-projected to the Universal Transverse Mercator (UTM) projection in zone 49N.

Method
This section is made up of three parts. Section 3.1 discusses the data preparation, including the definition of the finer ISA class, LST retrieval, and sample selection. Section 3.2 focuses on data variance analysis using charts and analysis of variance (ANOVA). Finally, Section 3.3 discusses the process of establishing the multivariate linear regression including variables such as building density, average building height, distance to city center, dummy variable construction, and the standardization of variables. A flowchart of the research design is shown in Figure 3.

Definitions of Types of Building Roofs and Natural Land Covers
Following careful investigations using high spatial resolution imagery and field trip surveys, we subdivided the land surfaces of the study area into five ISA classes on the basis of different types of building roofs and three natural land cover classes ( Figure 4). First, five ISA classes, namely, WS, BS, DM, DR, and RR, were divided on the basis of their color and material properties. The WS consisted of white building roofs made of polished white galvanized steel, while the BS included a steel panel painted in blue. DM refers to a rusty steel panel that is brown or dark in color, while the DR contains all dark-colored roofs that usually expand over large areas. RR represents residential roofs that consist of loosely distributed residential buildings within a neighborhood. Due to the limited spatial resolution of data sourced from Landsat imagery, the RR reflectance may be affected by shadows and vegetation. As such, the RR was assigned as an independent class for two reasons: (1) its land-use type is residential, where each patch has a small area, and the building heights are generally low, making sample selection difficult, and (2) it is a common land-use type scattered across urban areas with a significant impact on the LST distribution. Types of vegetation cover were classified as tree (TR) or grass (GR) categories; agricultural land was assigned to the latter class. The water (W) class, which contains surface waters such as lakes and rivers, was also used, while soil was excluded because its heat property is affected by moisture.
To elucidate the difference between classes, we collected data on the physical characteristics from previous studies, including the thermal conductivity coefficient and albedo value ( Table 2) [21][22][23][24][25][26][27]. The thermal conductivity coefficient illustrates that metal-like materials, such as WS, BS, and DM, have the highest values. This indicates that these materials absorb heat more rapidly and are thus quick to heat up. In other finer ISA classes, such as DR and RR, the thermal conductivity coefficients were found to be much lower than the WS, BS, and DM. The thermal conductivity coefficient of W was similar to that of RR (0.606). Albedo was an additional factor recognized to affect LST and was related to color; a white color has the highest albedo, while a black color has the lowest. A higher albedo value implies that more energy is reflected, and less energy is absorbed by the materials. As a result, the temperature of a high-albedo material is generally lower than a low-albedo material. Table 2 shows that the albedo of WS was the highest, whereas the albedo values of DM, DR, and RR were close to 0.1. The albedo of BS albedo was between 0.1 and 0.75, in terms of the color gradation. Generally, the albedo of GR was higher than TR, as its lowest albedo was much larger than TR; however, TR had the highest albedo according to the different types of trees. The albedo of W was almost zero, indicating negligible reflectance in this class.

Land Surface Temperature Retrieval
The original LST was retrieved using a method of atmospheric correction, the radiative transfer equation, in the single-channel algorithm with Band 11 in Landsat 8 TIRS [28][29][30]. This study adopted the calculation method developed by Jimenez-Munoz [30], which describes the detailed calculation process.

Sample Selection and LST Extraction
Samples of different building roofs and natural land covers (one pixel: 30 m × 30 m) were collected from the LST imagery, with the assistance of high spatial resolution images. Then, we used three windows to present the high spatial resolution, LST, and Landsat imagery using ENVI, the commercial software (https://www.ittvis.com/envi/, accessed on 1 November 2018). These three windows were geographically linked such that when moving the focus on the high spatial resolution image window, the corresponding location will change accordingly in the other two windows. We then investigated the entire study area in the high spatial resolution image window and labeled the building roofs that covered more than 9 pixels in the LST imagery. The pixel of the LST image, which was in the center of the building roof, was selected as the sample (see Figure 5, rectangle outlined by green dashed line). Finally, the value of the selected pixel in the LST image was extracted and stored in a table for statistical analysis. The selection of the natural land cover type used the same method as the building roof sample selection. The hot-and cool-season readings for different years were taken from the same samples in the LST imagery; the number of samples for each class is listed in Table 3.

Between-Class Comparison Analysis
Bar charts were used to illustrate the mean LST for each class; a comparison between classes may be easily interpreted by comparing the height of each class in the bar chart. In addition, one-way ANOVA was applied to compare the LST means to test whether there were any statistically significant differences between the different land cover types. In this study, ANOVA combined with a post hoc test was applied to LST data for each season of every year to determine whether there was a statistically significant difference between different land cover types during the hot and cool seasons.

Seasonal Variation Analysis
Seasonal variation was analyzed by comparing the variance of the mean LST in hot and cool seasons. The bar chart describes different values using the height of the bar and enabled an easy comparison of the differences between the two values. Thus, the mean LST of each class was assembled in a bar chart to illustrate the variation between seasons and years. The outline of a building along with building height information were utilized to calculate the building density (B_D) and average building height (B_H) in each sample. This information was downloaded from the open street map (OSM). Building density was computed using the total area of the building within a sample (pixel) divided by the total area of the sample (1 pixel, 30 m × 30 m = 900 m 2 ). The determination of B_D and B_H were based on Equations (1) and (2), respectively:

Regression Analysis
where B i,a and A s are the area of building, i, within a sample and the total area of the sample, respectively; B i,h is the height of building i; and N is the number of buildings in a sample. We also calculated the distance between the sample and city center. As this study was based on two cities, we calculated the distances between samples and city centers using the Euclidean distance; the shortest distance to the city center was used for the analysis. The distance was calculated using Equation (3): where DIS is the distance between the sample and city center, x 1 and y 1 are the sample coordinates, and x 2 and y 2 are the city center coordinates.

Dummy Variables Construction
Dummy variables were applied for attributes with more than one distinct category in regression analysis [31]. These variables operate as switches, such that they filter nonrelated cases in a regression model. Dummy variables were limited to two specific values: 1 or 0. Typically, 1 indicates the presence of a qualitative attribute, while 0 indicates otherwise. In detail, an dummy variable was created to represent an attribute with two or more distinct classes [31]. For each class, all missing values were recoded to "0" to minimize the impact of missing values on model results. The dummy variables were created to indicate whether "0" is a missing value or is the average (good) performance [32]. For example, the regression equation obtained with quantitative variables using the least square method is often difficult to fit; therefore, it is challenging to accurately reflect the internal relationship between elements in the city. Thus, dummy variables were introduced to filter irrelevant factors in the model, which may enhance model robustness [33][34][35]. The building roof types (WS, BS, DM, DR, and RR) and natural land covers (TR, GR, and W) were recoded as dummy variables (i.e., WS_D, BS_D, DM_D, DR_D, RR_D, TR_D, GR_D, and W_D, respectively). For example, if a sample belonged to WS, a value of 1 was assigned to WS_D, while the remaining values of WS_D were assigned a 0.

Variable Standardization
As the values of the collected variables were at different magnitudes, applying them directly to the regression model may generate biases. Therefore, it was necessary to standardize variables prior to their application in the regression model. The z-score was used to standardize the variables (Equation (4)): and V z is the standardized variable, x i is the value of sample i, M is the number of samples, x is the mean of the sample, and s is the standard deviation of the sample. LST, DIS, B_D, and B_H were standardized as LST_Z, DIS_Z, B_D_Z, and B_H_Z, respectively.

Multivariate Linear Regression
A multivariate linear regression model was used to explore the correlation and impact scales of each independent variable. Among the statistical models, linear regression analysis showed promise because of its reasonable accuracy and relatively simple implementation when compared to other methods [36]. Multi-linear regression (MLR) analysis is a means to evaluate the relationship between independent and dependent factors [37]. As such, this model has been widely used to analyze the relationship between urban environmental physico-chemical parameters and LST. MLR is a generalization of the simple linear regression model; this model allows for more than one predictor variable.
The standardized variables and land cover dummy variables were set as independent variables, and the standardized LST was selected as the dependent variable in the MLR model. The least squares model was applied to estimate the regression parameters. The regression equation may be expressed as per Equation (5): where LST z is the dependent variable standardized LST; xz i is the standardized independent variable, i (DIS_Z, B_D_Z, and B_H_Z); m is the number of standardized independent variables; xd j is the dummy variable, j, which may be WS_D, BS_D, DM_D, DR_D, TR_D, GR_D, and W_D; n is the number of dummy variables; and b 0 is the intercept of the regression model. The reference dummy variable was assigned as RR_D.

LST Distribution
Eight scenes of LST images encompassing hot and cool seasons for 2014, 2016, 2018, and 2020 were retrieved using the atmospheric correction method with Landsat TIR imagery. Figure 6 depicts all LST results using the same legend; these results illustrate that the LSTs during the hot seasons were much higher than those during the cool seasons. Most LSTs in the study area were higher than 30 • C in the hot seasons and were lower than 25 • C in the cool season. Urban area LSTs were higher than those in mountainous areas, which were covered by forests. Areas with high LST during the hot seasons were larger than those in the cool seasons. Extremely high LST areas (>35 • C in the hot season [38] and >25 • C in the cool season) were concentrated in eastern Foshan and eastern and northwestern Guangzhou. LSTs in the water areas were lower than those in other regions, where the shapes of rivers were clearly identifiable in the LST imagery. LSTs over areas covered by clouds were much lower than those in cloud-free regions.

Comparisons of Mean LSTs by Building Roof Types and Natural Land Classes in 2018
There were apparent differences in LST between various types of building roofs and other natural land surfaces. In general, seasonal changes in LSTs between various land classes were much more apparent when compared with their annual changes. The seasonal changes in the LSTs of each class in 2018 are shown in Figure 7; in descending order, the mean LST was DR > BS > DM > WS > RR > GR > TR > W during the hot season. However, this pattern slightly differed during the cool season, in which the descending order was DM > BS > DR > WS > GR > RR > TR > W. The LSTs of BS, DM, and DR were similar, with an average of approximately 45 and 25 • C during the hot and cool seasons, respectively. The LSTs of WS during the hot and cool seasons were lower than those of BS, DM, and DR, averaging at 42 • C during the hot season and 21 • C during the cool season. The mean LSTs of TR and W were much lower than the finer ISA classes; this was approximately 10-15 and 5 • C cooler in the hot and cool seasons, respectively. The mean LSTs for GR were lower than the finer ISA classes, albeit higher than TR and W; the means were also higher than those of RR during the cool season. To test whether there was a statistically significant difference in the mean LSTs between building roof types and natural land cover classes, we carried out post hoc tests using the Games-Howell model. Table 4 presents the results of these tests, whereby the p-values are displayed in the lower left panel, and the mean differences of LSTs are shown in the upper-right region. The entries were diagonally symmetrical, showing mean LST differences and corresponding p-values. The mean LST difference was calculated using the mean LST of the class in the second row subtracted by the mean LST of the class in the corresponding column.
WS was statistically different from the other classes, with the exception of RR in 2018 ( Table 4). The mean LST of WS was significantly lower than the finer ISA classes of BS, DM, and DR; however, it was significantly higher than RR and other pervious surface classes during the hot season. In terms of the hot and cool seasons, the differences between WS and other classes were larger during the hot season than the cool season.
The mean LST of BS was slightly higher than DM, whilst being lower than DR during the hot season of 2018; these differences were +0.5 and −0.08, respectively. Conversely, this pattern changed as the mean LST of BS was lower than DM during the cool season. The mean LST of DM was slightly larger than DR and significantly larger than the previous surface classes. The difference in the mean LST between DM and RR was approximately 4.2-5.2 • C during the hot season, and approximately 3.6-6.3 • C during the cool season; this indicates higher variance during the atter. The mean LST of DR was significantly different from that of RR, with this difference increasing during the hot season; it reached approximately 3.9-5.1 • C during the hot season and ≈2.8-4.8 • C during the cool season. The LST of RR was significantly higher than that of previous surface classes during the hot season; however, this difference was lower than the mean LST of other ISA classes. During the cool season, the LST of RR was lower than that of GR, implying a cool spot compared to GR during the cool season. The mean LSTs of TR, GR, and W were significantly different.

Seasonal Variability of LSTs by Years
By comparing the four years (see Figure 8 and Table A1), the pattern for most years was similar to that of 2018. This meant the mean LSTs were significantly lower than those of BS, DM, and DR, while being higher than those of RR, TR, GR, and W. Generally, differences were larger during the hot season compared to the cool season. The mean LSTs of BS were slightly lower than those of DM, although they were slightly higher than the mean LSTs of DR for most years, ranging from 2014 to 2020. The differences were not statistically significant, except during the hot season of 2013; the differences between the mean LSTs of DM and DR were <1 • C. This was particularly the case in 2016, when the mean differences in LSTs were near 0.1 • C. The differences in the mean LSTs between DM and the pervious surface classes were larger compared to BS in 2016, 2018C (2018 cool season), and 2020; however, they were lower in 2014 and 2018H (2018 hot season). The mean LST differences between the DR and pervious surface classes were smaller than the differences between DM and BS in 2014, 2018C, and 2020. The mean LST of RR was significantly lower than the other ISA classes, although it was higher than the previous surface classes. The differences during the cool season were much lower than the differences during the hot seasons.

Model Summary
Multivariate linear regressions were carried out using LST, standardized variables (B_D_Z, B_H_Z, and DIS_Z), and dummy variables (WS_D, BS_D, DM_D, DR_D, GR_D, TR_D, and W_D) from 2014 to 2020. The results indicate that these variables fitted well with the LSTs in the regression models; this was because the R 2 values of all four models were >0.74, reaching 0.87 for the hot season of 2014. This means that >74% of the variance in LSTs may be explained by the independent variables in the models (see Figure 9 and Table A2). Given that all p-values in the regression models were <0.01, the models were statistically significant (Table A3). Moreover, the R 2 values during the hot seasons were higher than their cool season counterparts, indicating that the hot season regression models were slightly better performing (Figure 9).

Regression Coefficients
RR was the reference group for all regression models; thus, the coefficients of the other dummy variables may be interpreted in comparison with the reference dummy variable, RR. The dummy variables of the finer ISA classes (WS, BS, DM, and DR) had positive coefficients from 2014 to 2020 in the cool and hot seasons (Table 5). This indicates that these finer ISA classes generated higher positive LSTs compared to RR. The scale patterns of these coefficients were consistent during the cool season; the descending order of these coefficients was DM > BS > DR > WS, in 2014, 2018, and 2020 (Table 5). In 2016, the coefficient of DR was higher than those of DM and BS. The WS coefficients were lower than BS, DM, and DR during the hot and cool seasons for four years, implying that WS had a lower contribution to LST. Generally, the coefficients of WS, BS, DM, and DR during the cool seasons were higher than those in the hot seasons. This demonstrates a stronger heating effect on LST during the cool season than the hot season.

Discussion
ISAs are human-made features in which water cannot infiltrate the soil; this term focuses on the physical characteristics of water infiltration. Different compositions of the ISA alter the urban surface energy budget [19]. The assumption that all ISA classes are the same material does not accurately reflect the ISA heating properties; as such, specific drivers of urban heat retention within individual urban areas cannot be quantified [19].
The results of this study demonstrate that various finer ISA classes showed different heat characteristics and thus differed in their contribution to LST. The LSTs of BS and DM were much higher than those of WS, RR, and the pervious surface classes. BS and DM had similar thermal conductivity coefficients around 10-60 (W/m K) at 25 • C ( Table 2) [27]. This means they may heat up rapidly when they absorb energy from the sun. The LSTs of DR and DM were much higher than those of WS and RR; DM had the highest impact on high LST values, whereas WS showed the least heating effect as per the corresponding regression coefficients. These results may be due to their low albedo characteristics (close to 0.1); this pattern aligns with the color of absorbing heat, where black > white. Black materials are able to absorb more heat than white materials under the same environmental conditions, as white-colored materials reflect more solar energy than materials with other colors [22,25,39].
The influence of different surface materials on LST varied during the hot and cool seasons. The impact magnitudes of finer ISA classes were more substantial during the cool season, although they were weaker than those of vegetation and water during the hot season ( Figure 10). Some studies have concluded that the radiation intensity of an ISA increases in the summer compared to other seasons, resulting in an increased albedo and a corresponding reduction in LST [39][40][41]. By contrast, rural areas experienced an increase in temperature that was up to 0.27 • C during the summer [42]. An additional factor that affected the LST contribution of finer classes was the moisture content of air between seasons. Air was less humid during the cool season than the hot season, reducing the surface heat capacity and increasing the surface heat radiation during the former [43,44]. Therefore, the contribution of finer ISAs to LST was weaker during the hot season than the cool season.
Parameters such as building density, average building height, and distance to the city center were used in the regression analysis. Most research to date has demonstrated that the horizontal and vertical differences in three-dimensional cities have an impact on LST. Many researchers have often linked building height to the urban shadow; the higher the building, the lower the temperature in this area due to the effect of shadows [45,46]. However, there was also a significant positive correlation between building density and LST during the hot season, and a negative correlation during the cool season. This is related to the monsoon climate, building ventilation, and human activities during different seasons [47,48]. There were four key limitations in this current study. First, only five types of building roofs were discussed because of the limited sample collection within the study area. Finer ISA classes, such as brick, stone, and plastic, should be collected for analysis in future studies. Second, the analysis of the impact of finer ISA classes on LST was limited to daytime imagery, particularly 11:00 a.m. local time. These impacts may vary during the afternoon or nighttime. Moreover, this research focused on the relationship between finer ISA classes and LST; however, the factors that affect urban LST are complex, including urban three-dimensional space, urban functional area distribution, and heat emissions from human activities. At the same time, LST was closely related to surface energy balance (SEB); thus, it is recommended that future research analyze the changes and driving factors of urban LST within a multi-dimensional context. Third, the analysis method used only accommodates multiple regression; as such, a new method of comparative analysis is required in future research to fully reflect the relationship between each finer roof class and the LST. Finally, shadows may have an important effect on the LST. During the sampling process, we selected samples with no or limited shadowed areas. While shadow areas were not removed during the data pre-processing stage of this study, this step should be considered in future studies.

Conclusions
This study discussed the impacts of finer ISA classes on LST and their temporal variations. Bar charts and ANOVA were used to demonstrate the different LST characteristics of each class. Multivariate regression models were applied to selected variables (where RR was set as the reference) during the hot and cool seasons for 2014, 2016, 2018, and 2020. The key findings from the analysis were: (1) The mean LST of different types of building roofs were statistically different from each other; these differences were stronger during the hot season than the cool season. The LST of WS was significantly different from that of BS, DM, DR, and RR. The mean LSTs of BS, DM, and DR were significantly higher than those of WS and RR. (2) The impacts of building roof types on the LST differed. WS, BS, DR, and DM were positively correlated with LST, while RR was negatively correlated with LST during the cool season. The impact of WS on LST was significantly different from that of BS, DR, DM, and RR. The impacts of the finer ISA classes of BS, DR, and DM on LSTs were similar; however, they demonstrated a consistent descending order pattern (i.e., DM > BS > DR >WS), during the cool season. (3) The contributions of finer ISA classes to LST variance were stronger during the cool season than the hot season, as their standardized coefficients were larger in the former than the latter.
This study provides useful information to mitigate the UHI effect by unveiling specific urban heat retention drivers within urban areas. This includes different types of building roofs and pervious surface classes. This information may assist policymakers in developing effective measurements to deal with the UHI phenomenon and related health concerns.