Highlights What are the main findings? Road TCE did not increase linearly with its grade, which was primarily determined by traffic flow (volume and speed). Urban main roads had higher TCE levels than other roads, while lower-grade roads (minor arterials) posed a higher thermal risk. The impact of TCE on RST was more pronounced on lower-grade roads, whereas the effect of the surrounding environment tended to be less noticeable. What are the implications of the main finding? An improved random forest (RF) model has been proposed to identify the nonlinear effects of segment-scale TCE and built environment factors on RST, which enriches our understanding of the urban thermal environment. Road traffic emissions contributed significantly to the local thermal environment, especially on lower-grade roads, providing a scientific basis for developing targeted strategies to reduce carbon emissions and mitigate urban thermal stress. Abstract Urbanization and rapid economic growth have exacerbated urban heat effects, increasing the frequency of heat-related disasters and intensifying human health risks. Urban traffic generates substantial carbon emissions and associated heat, which significantly alter roadside thermal environments and impact human activities. Numerous previous studies have investigated urban thermal environments and their influencing mechanisms. However, the relationships between road-level traffic carbon emission (TCE) and road surface temperature (RST) remain insufficiently explored. In this study, roadway segment-based TCE and RST were acquired by integrating hourly traffic flow information, localized vehicle carbon emission factors, high-resolution Landsat-8 remote sensing datasets, and the road network. Three commonly used linear regression models and an improved Random Forest (RF) model were utilized to assess the impact of TCE on RST for different grades of roads. The study showed that carbon emissions from road traffic exhibit a locally focused distribution pattern in space. Compared to other grades of roads, higher levels of TCE were observed in urban main roads. In summer, roads (e.g., minor arterials) with lower grades tended to have a higher thermal risk, with freeways having the lowest TCE and urban expressways experiencing the greatest TCE fluctuations. An improved RF model integrating the spatial weight matrix and Gaussian process could more efficiently identify the nonlinear effects of TCE on RST. The contributions of TCE to summer RST were 0.4, 0.37, 0.54, and 0.56 for freeways, urban expressways, main roads, and minor arterials, respectively. The relative impact of road TCE with lower grades on RST becomes more significant, while the impact of surrounding buildings and green areas tends to decrease. Our findings provide valuable insights for reducing urban carbon emissions and thermal risks. 1. Introduction These thermal anomalies not only exacerbate environmental risks but also pose severe threats to public health systems [ 9]. Epidemiological analyses have established a dose–response relationship between elevated urban temperatures and increased human morbidity/mortality rates [ 10, 11]. Empirical evidence from a Lancet study demonstrates that 2023 witnessed a 6% surge in heat-induced sleep deprivation relative to the 1986–2005 baseline, concomitant with a 167% increase in mortality among seniors over 65 compared to the 1990s levels [ 12]. Such findings underscore the critical necessity for developing comprehensive assessment frameworks to decipher thermal environment dynamics and their multiscale drivers—an imperative research frontier in urban climate and environment fields [ 13, 14]. LST, a critical indicator of urban heat, has been widely used in studies related to the urban thermal environment [ 15, 16]. In order to reduce urban thermal risks and hazards, numerous studies have investigated the factors influencing high surface temperatures and their driving mechanisms [ 17, 18]. These factors include urban spatial morphology [ 10], blue–green infrastructure [ 19], impervious surface or building coverage [ 20], and surface physical properties [ 21], etc. By reviewing the literature, most of the studies were executed by integrating multi-source remote sensing and other geospatial datasets. As the data or information related to energy consumption is difficult to obtain, such as electricity supply, oil consumption, equipment power parameters, etc., few studies have investigated the effects of urban carbon emissions on the thermal environment. In fact, energy consumption-induced carbon emissions have been identified as one of the major contributors to global urban warming and environmental degradation [ 22, 23, 24]. It has also been proven that reducing urban carbon emissions can effectively improve the thermal hazards and enhance public health [ 25]. There is a well-established academic consensus that high-concentration greenhouse gases (e.g., carbon dioxide and methane) trap more heat in the atmosphere, leading to a gradual rise in ambient temperature. Vehicular carbon emissions are inherently associated with anthropogenic sensible heat emissions, including engine exhaust and tire–pavement friction heat. Direct heat release and radiative trapping effects of carbon-based gases together contribute to localized road surface warming [ 19, 26, 27, 28]. Further, high-carbon-emitting areas in urban environments, such as industrial parks, commercial districts, and transportation operations (e.g., harbors), have been found to pose a higher thermal risk. According to the aforementioned investigations, previous studies have surveyed the surface temperature variations in different urban functional zones with high carbon emissions. However, the spatio-temporal variation in road surface temperature (RST) has lacked attention. Moreover, the investigation of the effect of roadway traffic carbon emission (TCE) on RST is remarkably scarce. Undoubtedly, roads are also an important part of urban functions [ 29], similar to residences, workplaces, and recreational spaces. Urban traffic is the largest contributor to anthropogenic carbon emissions in cities, accounting for more than 50% of neighborhood-scale CO 2 emissions [ 30, 31]. Notably, TCE not only directly represents greenhouse gas emission levels but also correlates significantly with traffic-generated waste heat [ 32, 33], making it a comprehensive indicator for studying the impacts of traffic activities on the thermal environment [ 34]. Higher TCE typically corresponds to denser traffic flow [ 35], longer vehicle operation times, greater engine intensity [ 36], larger vehicle tonnage, increased tire–road friction [ 37], and frequent start–stop cycles or low-speed idling in congested urban areas [ 38]. A large number of fuel vehicles emit waste heat and exhaust gases under driving conditions, which directly affect the thermal environment and air quality in the surrounding areas of roads [ 39, 40]. Unlike the other types of land in patches, relatively narrow roads of different grades act like a network, covering the dense urban built environment [ 41, 42]. Additionally, it should be emphasized that human activities (e.g., job-housing commuting, recreational activities, outdoor walking) are closely associated with the road network, which directly exposes urban inhabitants to thermal risks and health hazards [ 43, 44, 45]. Road traffic emissions induce pronounced thermal effects, substantially altering roadside thermal conditions, exacerbating urban heat stress, and exerting tangible impacts on the microclimate and human activities along road corridors. Therefore, it becomes urgent to investigate the spatio-temporal variations in TCE and RST and their relationships [ 46]. In this study, we first integrate hourly traffic flow information and localized vehicle carbon emission factors to identify the spatio-temporal characteristics of TCE for different grades of roads. High-resolution Landsat-8 datasets and road networks are then applied to capture RST variations based on roadway segments. Finally, the relationships between TCE and RST are systematically analyzed by linear regression models and an improved Random Forest (RF) model. This study will help to achieve a win–win situation: providing valuable insights for reducing urban traffic-related carbon emissions and improving the thermal environment of roads and nearby areas. 2. Study Area and Datasets 2.1. Study Area This study was conducted in Shenzhen, Guangdong Province, a megacity located in southern China with a population exceeding 17 million. It lies within the coordinates of approximately 113°46′ to 114°37′E and 22°27′ to 22°57′N. Shenzhen is characterized by its high-density urban layout, with a highly concentrated development pattern and a dense road network. Additionally, the city exhibits a high level of motorization, with a substantial proportion of the population relying on private vehicles [ 8, 47]. These attributes make Shenzhen an ideal case study for examining the interactions between TCE and RST, as well as their effects on human activities. The compact urban form, coupled with its advanced transportation infrastructure, provides a rich context for investigating the broader implications of urbanization, climate change, and human health. Figure 1 illustrates the study area and its road network. According to the 2023 Shenzhen Comprehensive Transportation Annual Evaluation Report compiled by the Shenzhen Municipal Transport Bureau in August 2024, in 2023, the cumulative road mileage in Shenzhen amounted to 8581.7 km. Specifically, freeways spanned 392.4 km, urban expressways measured 247.5 km, main roads extended for 1394.1 km, minor arterials covered 1089.0 km, and branch ways, along with other road types, totaled 5458.7 km. During peak traffic hours, the average speeds on freeways, urban expressways, and main arteries were 62.7 km/h, 40.1 km/h, and 27.1 km/h, respectively. All road mileage and peak-hour speed data presented in this paragraph are sourced from the aforementioned official report. 2.2. Data Collection and Preprocessing The key data used in this study include road network, traffic flow information, carbon emission factors, land surface temperature, building and green space, as shown in Table 1. (1) Road network data was originally obtained from the OpenStreetMap (OSM) platform. Since the traffic flow information in this study is linked to city roads of different grades (i.e., freeway, urban expressway, main arterial, minor arterial, and branch way), vehicular access roads in restricted or semi-restricted zones were not considered, such as gated communities, industrial parks, schools, and public parks. Each road segment was simplified to a vector polyline. The initial basic unit (i.e., road segment) in the network was defined according to road intersections. To address the issue of varying road segment lengths, especially the much longer stretches of freeways compared to branches, a segmentation approach was adopted. Specifically, the long segment was divided into short sections of approximately 200 m in length [ 48] to finely characterize TCE and RST. (2) Traffic flow data was sourced from the Shenzhen Road Traffic Emission Monitoring Platform established by the Shenzhen Urban Transport Planning Center Co., Ltd. (Shenzhen, China) (SUTPC). SUTPC is a professional institution directly under the Shenzhen Municipal Transport Bureau and serves the municipal transportation industry. The platform is an internal data middle platform of SUTPC, with the official website: https://transpaas.sutpc.com/#/tp/data-platform (accessed on 26 May 2025). It includes minute-level driving speeds and traffic volumes for each road segment. Among these, traffic volume refers to the number of different vehicle types, including fuel-powered or electricity-driven cars, buses, trucks, and taxis. It comes from the city-wide traffic sensing network deployed by the transportation and management authorities of Shenzhen. Details of the preprocessing process for traffic flow data were described in our prior research [ 49]. (3) Carbon emission factors were also developed by the SUTPC, which represents the carbon emissions generated per kilometer traveled by a vehicle. Initial laboratory calibration of these factors was completed in 2019, followed by operational application starting in 2020. The localized factors were further verified and updated in 2023 to reflect current vehicle emission characteristics. These time-calibrated carbon emission factors are derived from the internal data middle platform of SUTPC, as mentioned above. These localized factors consist of two components: static vehicle attributes and dynamic driving conditions. The static attributes are based on the “National Technical Guidelines for Motor Vehicle Air Pollution Emission Inventory in China”, which categorizes vehicles into 11 basic types. The dynamic driving conditions are the average speed, parking time ratio, and relative positive acceleration frequency of each type of vehicle. Each type of vehicle is divided into 25 service levels based on driving conditions in Shenzhen [ 49]. (4) Land surface temperature (LST) data were derived from Landsat-8 OLI/TIRS Collection 2 Tier 1 Level-2 imagery, which is publicly distributed by the United States Geological Survey (USGS). All processing was performed on the Google Earth Engine (GEE) platform, where LST values were extracted directly from the pre-calibrated ST_B10 surface temperature band. Cloud and shadow contamination were eliminated using the QA_PIXEL and QA_RADSAT quality bands, following the official specifications for Landsat-8 Collection 2 Level-2 data processing. (5) The primary building data for Shenzhen used in this study were obtained from the Smart City Sensing and Simulation Laboratory at Nanjing Normal University. This dataset, with a spatial resolution of 1 m, was developed using a deep learning-based semantic segmentation model trained on multi-source remote sensing data, and the relevant methods and products have been published in Scientific Data [ 50]. Building areas are typically characterized by concentrated human activities, reflecting the urban spatial layout, population density, and energy consumption. Additional datasets on building rooftops and the CNBH-10 m data [ 51] were utilized to ensure the accuracy of the building information [ 52]. The height data range is 0~600 m, and the resolution is 1 m. (6) Green space data was obtained from the Shenzhen Key Laboratory of Spatial Smart Sensing and Services. The data were derived using an RF classification method applied to high-resolution GF-2 satellite imagery, which also had a pixel size of 1 m [ 53]. The green spaces reflect the natural elements in the urban environment, which play crucial roles in regulating temperature, absorbing carbon dioxide, and mitigating urban heat effects. Details of the main datasets used in this study are summarized in Table 1. Data collection in this study covered the period from 2020 to 2023. Road network data and traffic volume data were continuously collected with good data integrity and availability. Considering that the carbon emission factors were revised and updated in 2023, potential deviations existed in the emission data of previous years. Therefore, 2023 was selected as the core research year, and summer was chosen as the typical research period. As illustrated in Figure 2, the 30 m resolution LST of the entire city on 14 August 2023 was adopted as the core temperature data, and winter temperature data were used for supplementary analysis. Landsat-8 OLI/TIRS imagery was selected for LST extraction due to its mature thermal infrared observation capability, 30 m spatial resolution suitable for urban road network analysis, and reliable atmospherically corrected Level-2 products. Multiple Landsat-8 images in August 2023 were compared, and a multiple cloud-free image acquired in August 2023 (local overpass time: ~10:30 a.m.) was finally adopted owing to its optimal data quality. Shenzhen is a coastal hilly city, and the solar exposure of road surfaces was relatively uniform at the satellite overpass time. To ensure spatiotemporal consistency, hourly carbon emissions corresponding to the satellite overpass hour were used for correlation analysis. In addition to the summer image, November and December were selected as representative winter months. Since this study focuses on the more severe thermal environment issues in summer, the seasonal comparison is designed to reveal the spatiotemporal distribution differences between the two seasons. Building data processing was conducted exclusively for areas along road segments, based on the 2020–2023 urban construction and development approval records provided by the Shenzhen Municipal Bureau of Planning and Natural Resources. Newly built or modified buildings alongside roads within this period were screened and checked, and approximately 160 inconsistent building records were removed to ensure data reliability for subsequent road-based analysis. The refined building data and green space data are presented in Figure 3a and Figure 3b, respectively. After the above screening, revision, and integration steps, the final number of valid road segment units in the study area was determined to be approximately 48,000. Road TCE was calculated on an hourly basis. The total hourly emissions were used to explore the overall spatio-temporal distribution characteristics of road TCE, while hourly emission data precisely matched with the satellite overpass hour of LST were applied for the correlation analysis between TCE and surface temperature. 3. Methodology The workflow of this study consists of the following three steps, as shown in Figure 4. The spatio-temporal characteristics of TCE for different grades of roads were first identified by integrating hourly traffic flow information and localized vehicle carbon emission factors. Then, Road segment-based RST variations were analyzed using high-resolution Landsat-8 datasets and the road network. Finally, the relationships between TCE and RST were explored by linear regression models and an improved RF model. 3.1. Traffic Carbon Emission (TCE) Calculation A bottom-up approach was employed to calculate the hourly carbon emissions for each road segment, considering traffic flow information from the road network and carbon emission factors [ 54]. The traffic flow data for each road segment included vehicle type information and driving condition details. Driving conditions are determined by both road grades and vehicle operating status, with parameters of vehicle operating status, including the average speed, the relative acceleration, and the ratio of stop time. Carbon emission factors were determined based on laboratory calibration for each vehicle type, representing the amount of carbon emissions produced per kilometer driven under different driving conditions [ 49]. The emissions were calculated according to the following formula: T E h , k = ∑ i T E h , k , i , (1) where T E h , k represents the carbon emissions for road segment k during time period h (unit: kg), and T E h , k , i represents the carbon emissions from vehicle type i traveling on road segment k during time period h (unit: kg). The emissions from individual vehicle types on each road segment were calculated as T E h , k , i = q h , k , i ( g ) ୍ଠ E F i ( g ) ୍ଠ L k , (2) where q h , k , i ( g ) represents the number of vehicles of type i traveling on road segment k during time period h with specific driving conditions (unit: vehicles); E F i ( g ) represents the carbon emission factor associated with vehicle type i under condition g (unit: kg/vehicle·km); and L k indicates the length of road segment k (unit: km). A total of over 6000 h of high-resolution, second-by-second exhaust emission and trajectory data were compiled for 11 distinct vehicle categories. Trajectory datasets encompassed vehicle identification numbers, timestamps, geographic coordinates (longitude and latitude), speed, and acceleration, with the ratio of stops estimated via indirect computational methods. These speed, acceleration, and the ratio of stop were incorporated as explanatory variables in a Random Forest model, trained against corresponding exhaust emission measurements. The model output was used to establish a queryable emission inventory, E F i ( g ) which accommodates variations in vehicle type and operational status. The model output is used to establish a queryable inventory with a total of over 1200 items, which takes into account variations in vehicle types and operating conditions. Taking several car models as examples, the output results of carbon emission factors are provided. Levels of Service (LOS) grades 1 through 5 represent traffic service levels for various road grades, ranging sequentially from free-flow conditions to severe congestion. For instance, Freeway/LOS 1 denotes a free—flow traffic state on freeways. The following Table 2 shows the carbon emission factors of some cars. Specific parameters for Freeway LOS 1 to LOS 5 are presented in Table 3, including average speed, relative acceleration, and stop-time ratio. These parameters were typical features comprehensively extracted based on the aforementioned trajectory information, and were calibrated with the data provided by the competent administrative authority. A vehicle tracking and monitoring method was employed to verify and calibrate the reliability of various types of vehicle exhaust emissions calculations. Over 60 h of such vehicle tracking tests were conducted, as shown in Figure 5, and the obtained data were utilized to update the emissions inventory. 3.2. Road Surface Temperature (RST) Retrieval The road network data and Landsat-8 remote sensing images are used to identify the RST for each road link. The detailed traffic flow information and localized carbon emission factors are integrated, considering both dynamic and static vehicle behaviors, to calculate carbon emissions for all links in the road network. A single-channel algorithm was used to retrieve the LST from the thermal infrared band of Landsat-8 imagery [ 55, 56]. The retrieval formula is expressed as follows: L S T = γ ⋅ ε − 1 ⋅ φ 1 ⋅ L s + φ 2 + φ 3 + δ , (3) where L s represents the radiance value of the sensor (unit: W/(m 2·μm·sr)); γ and δ are two parameters based on the Planck function; ψ 1 , ψ 2 , and ψ 3 are parameters related to the water vapor content in the atmosphere, and ε denotes the surface emissivity. To map the retrieved LST values to individual road segments and calculate their mean, the pixel values for each road segment were first extracted, and the average temperature of all pixels within each segment was then computed. Let L S T i j represent the LST of the j -th pixel on the i -th road segment, and let N i denote the total number of pixels on the i -th road segment. The mean RST of road segment i, denoted as RSTi, is expressed as R S T i = 1 N i ∑ j = i N i L S T i , j , (4) where RSTi is the mean RST of the i -th road segment; L S T i j is the LST of the j -th pixel on the i -th road segment; and N i is the total number of pixels on the i -th road segment. This method provides the representative surface temperature value for each road segment by averaging the LST of all pixels within the segment, thereby offering essential data for further investigation of RST variations across road segments. 3.3. Identifying Relationships Between TCE and RST The Getis–Ord Gi* was adopted to detect spatial clustering patterns, enabling the identification of “hot spots” (high-value clusters) and “cold spots” (low-value clusters) in spatial data. The Hot Spot Analysis (Getis–Ord Gi*) tool in ArcGIS v10.2 was used to conduct spatial clustering of vehicle carbon emissions and RST across Shenzhen’s road network. It identified high-value (hot spots), low-value (cold spots), and non-significantly clustered regions for both variables. For the interconnected road network, this method enables deeper spatial analysis. It delineates the continuity and connectivity of hot/cold spot clusters along linked road segments and intersections. It reveals spatial heterogeneity and coupling between emission hotspots and RST hotspots. By calculating the local spatial autocorrelation between each spatial unit and its neighboring units, it determines whether the unit exhibits statistically significant high or low clustering [ 57]. The formula is as follows: G i = ∑ w i , j x j − X ପ୍ତ ∑ w i , j S n ∑ w i , j 2 − ( ∑ w i , j ) n − 1 (5) where Gi represents the Getis–Ord Gi* statistic for the i-th spatial unit, wi,j denotes the spatial weight between units i and j, xj represents the observed carbon emissions or RST value for the j-th spatial unit. X ପ୍ତ is the mean of all observed values across spatial units, S is the standard deviation of all spatial unit values, and n is the total number of spatial units. Three commonly used linear regression models—ordinary least squares (OLS) model, spatial error model (SEM), and spatial lag model (SLM)—were applied to examine the simple relationships between TCE, building coverage, green space, and RST. Subsequently, the standard RF algorithm was improved by integrating spatial weight matrices and Gaussian processes to identify complex effects of TCE on RST. The enhanced model combines the traditional advantages of the RF algorithm with the integration of multi-source geospatial information. The spatial weight matrix was introduced, which took into account the geographical proximity of road connections and the potential impact of surrounding areas. The RF algorithm, as a powerful ensemble learning method, performs exceptionally well in capturing complex nonlinear relationships. However, the standard RF algorithm ignores the inherent spatial correlations present in spatially structured data. In geographical or spatially correlated data, features such as geographical coordinates (e.g., longitude and latitude) may become significant factors influencing the target variable. If only the standard RF algorithm is applied, the model might fail to capture the latent spatial effects, leading to a reduction in predictive accuracy [ 58, 59]. To effectively address the spatial correlations, a method that considers both spatial relationships and feature interactions must be designed. A combination of the spatial weight matrix and the Gaussian process was proposed as an improvement to the traditional RF algorithm. The core idea was to introduce spatial matrix adjustments to capture spatial relationships and to combine Gaussian processes to smooth the impacts of spatial features, thus enhancing the prediction accuracy of the model [ 60]. Specifically, the spatial weight matrix was used to adjust the data by representing the similarity between spatial locations. Gaussian process modeling was applied to further introduce spatial correlation. Finally, this improvement was integrated into the standard RF algorithm progressively. The initial parameters of the RF algorithm, such as the number of trees T, are set. The geographical coordinates of each observation point are extracted, and the distance matrix D is calculated for subsequent spatial weight calculation. (1) Adaptive-bandwidth Gaussian kernel weights (local). For each link Xi, the adaptive bandwidth hi is calculated as follows: h i = k 1 n i ∑ j ∈ N i X i j − X ପ୍ତ i 2 (6) where ni is the number of samples in the neighborhood of link Xi; Ni represents the set of its neighborhood links; Xij is the feature value of the j-th sample in the neighborhood. X ପ୍ତ i is the mean feature value of the neighborhood samples, and k is the adjustment coefficient. Based on this, the local spatial weight matrix is calculated as follows. W l o c a l X i = 1 d X i , X j + ε exp − d X i , X j 2 2 h i 2 (7) (2) Global spatial weights. The average distance between link Xi and all other links is calculated to obtain the global distance metric D g l o c a l X i . The global spatial weight matrix W g l o c a l X i is then constructed as follows. W g l o c a l X i = 1 D g l o c a l X i + ε (8) (3) Fused weight matrix. The fusion coefficient λ is determined through a cross-validation method. The local and global spatial weight matrices are fused to obtain the final spatial weight matrix W f i n a l X , as follows. W f i n a l X = λ W l o c a l X + 1 − λ W g l o c a l X (9) (4) Definition of the Gaussian process model. Let the Gaussian process be G P m x , k x , x ′ , where m x is the mean function and k x , x ′ is the kernel function (such as the commonly used squared exponential kernel function) to describe spatial correlations. It is assumed that the observed data are y = y 1 , y 2 , … , y n . Given the input X = X 1 , X 2 , … , X n , y follows a joint Gaussian distribution p y X = N m X , K X , X , where m X is the mean vector, K X , X is the covariance matrix, and K i j = k X i , X j . (5) Iterative learning process. The RF model F R F X is initialized, and the weighted features X ^ = W f i n a l X ⊙ X are calculated according to W f i n a l X . After obtaining the weighted features X ^ = W f i n a l X ⊙ X , iterative learning is performed using these features to update the RF model F R F ( X ^ ) . The final prediction is then given by: y ^ = F R F ( X ^ ) (10) Via this progressive process, the spatial weight matrix is first used to adjust the input features, and the Gaussian process is then applied to further enhance the modeling of the spatial correlation. This ultimately improves the predictive power of the RF model. The adjustment coefficient k in Equation (7) is a key parameter for controlling the smoothing level of the adaptive bandwidth. To determine its optimal value, a grid search (range: 0.5–2.0, step: 0.1) combined with spatial five-fold cross-validation was performed. Spatial stratified splitting was adopted to avoid spatial data leakage, and the optimal k = 1.0 was selected based on the highest R2 and lowest RMSE of the model. This parameter balances the local smoothing effect and the spatial heterogeneity of road network data, ensuring the stability of the spatial weight calculation. Considering the spatial autocorrelation of TCE and RST in road segments, and the introduction of a spatial weight matrix, spatial k-fold cross-validation was adopted for model training and validation instead of a completely random data split. Spatial k-fold cross-validation divides samples under spatial stratification rules. It ensures that training and test subsets are spatially non-overlapping and independent. Spatial stratification constraints were applied in the sample partitioning process. This approach effectively avoids spatial data leakage caused by adjacent similar road segments. It also guarantees the credibility and objectivity of model evaluation results [ 61]. The R 2 value, mean squared error (MSE), and root mean squared error (RMSE) are the metrics adopted to evaluate the performance of both the standard and enhanced models. R 2 measures the proportion of variance in the dependent variable that is explained by the independent variables in the model, with higher values indicating a better fit. The MSE quantifies the average squared differences between the predicted and actual values, indicating the accuracy of the model, with lower values representing better performance. The RMSE is the square root of the MSE, which offers a more interpretable measure of the average prediction error in the same units as the dependent variable. 4. Results 4.1. Spatial Distribution of TCE for Different Grades of Roads August 2023 was selected as a representative summer month, and hourly TCE for the entire month was calculated for each road link segment. The daily average TCE derived from this dataset was used to characterize the overall summer emission level. Subsequent correlation analysis with RST uses hourly TCE matched to the satellite overpass date and time, as detailed in later sections. Figure 6a depicts the spatial distribution of total TCEs by road segment in Shenzhen. According to statistics, the total daily TCE in Shenzhen was approximately 2.67 × 10 7 kg. The overall TCE level did not increase linearly with road grade, which was primarily influenced by traffic flow (volume and speed). Road mileage had a certain impact on the total TCE for roads of different grades. Figure 6b. Total length (left Y-axis, unit: km) and total TCE (right Y-axis, unit: kg) of different road grades in Shenzhen. Specifically, the total length of arterial roads is 1394 km, resulting in the highest TCE (1.37 × 10 7 kg), which is almost equivalent to the combined TCE of the other four grades of roads. This is because the main arteries, which have a large number of relatively long-distance commuters, generally exhibit low driving speeds and severe traffic congestion. It has been confirmed that frequent acceleration and deceleration of vehicles, as well as prolonged idling conditions, result in low energy utilization efficiency, thereby emitting a significant amount of exhaust pollutants and waste heat. Compared to main arteries, freeways and urban expressways with higher traffic volumes had relatively lower and the lowest TCE values, respectively. It is likely attributed to their shorter overall mileage and faster driving speeds. Previous studies have confirmed that vehicles in non-violent driving and idle conditions have lower carbon emissions, due to the high efficiency of internal combustion engines in utilizing fossil fuels [ 54]. Notably, the total mileage of secondary arteries slightly exceeds that of main arteries, but their TCE is only one-quarter of that of main arterial roads. This is likely due to the lower traffic volume and smoother traveling conditions on secondary arterial roads. Branches with the lowest segment-based average flow and slowest traveling speed have relatively low but not the lowest TCE. It is likely attributed to their longest mileage. 4.2. Spatial Variation in RST for Different Grades of Roads Figure 8 shows the spatial distribution of RST in Shenzhen on the selected representative summer day in August. The road facilities indicated by red lines were used to estimate the corresponding RST of the road infrastructure. Notably, areas such as Bao’an International Airport in the western part of the city, as well as coastal ports, exhibited significantly higher temperatures. Based on the carbon emission distribution results reported in the previous section, it is necessary to distinguish road grades for further research due to the basic attributes of roads of different grades, such as the number of lanes and the distance between intersections. These attributes play crucial roles in determining the traffic flow patterns and vehicle operation conditions on different roads, which, in turn, have significant impacts on carbon emissions. Complementing the summer RST pattern, Figure 9 presents the spatial distribution of summer–winter RST differences across Shenzhen. As observed in the summer RST results ( Figure 8), the city shows a distinct spatial gradient: northwestern inland areas recorded high summer temperatures, while southwestern coastal zones were relatively cooler. Central urban cores also exhibited elevated summer RST values, comparable to those in the northwest. However, the seasonal temperature difference reveals a different pattern: the northwestern zones showed only moderate summer–winter temperature variations, while central urban areas experienced significantly larger seasonal RST differences. These high-difference regions are concentrated in areas with intensive human activity and dense road networks, and exhibit clear spatial clustering characteristics. This clustered distribution provides a direct basis for the subsequent hot spot analysis of seasonal temperature variations. Figure 10 presents box plots of the RSTs for five grades of roads in Shenzhen during summer, while Table 4 reports the specific numerical values. The blue box plots represent the interquartile range, with the lower and upper bounds corresponding to the 25th and 75th percentiles, respectively. In summer, the maximum recorded RST reached 48.68 °C. The RST generally exhibited an increasing trend from highways to branch roads. The median RSTs for freeways, urban expressways, main roads, minor arterial roads, and branch roads in summer were 39.83, 40.81, 41.36, 41.46, and 41.75 °C, respectively. The difference in surface temperature patterns between highways and low-grade highways may stem from complex interactions between physical, environmental, and human factors. In summer, surface temperatures across five road grades—from expressways to branch roads—demonstrate an increasing trend. Specifically, both the mean and quantile values of expressways are lower than those of lower-grade urban roads. Potential causes include the fact that expressways feature more lanes and higher traffic volumes, alongside their extensive geographical exposure and intense sunlight. Conversely, major and minor roads are more susceptible to complex influences such as building clusters, greenery shading, and traffic exhaust emissions. Further investigation is needed to determine the contributions of these influencing factors. Figure 11 presents the Getis–Ord analysis results of RST, in which TCE hotspots, cold spots, and regions with no significant difference are highlighted in red, blue, and yellow, respectively. Getis–Ord analysis is a spatial statistical method used to identify clusters of high or low values, revealing significant hot spots and cold spots within geographic data. In comparison with the emissions shown in Figure 7b, the areas with high temperature aggregations are more dispersed. The Futian Central District, located in the southernmost center of Shenzhen, was found to have a cold spot in terms of RST. This may be attributed to the wide roads and large amount of surrounding greenery in this area. In Longhua, north of the Futian Central District, and in Nanshan District on the west side, the high RST areas were found to be fragmented. More importantly, the low-grade roads in some remote areas surprisingly exhibited the characteristic of high temperature aggregation. Significantly, it was found that certain Minor Arterials, which play a crucial role in alleviating the traffic burden on urban expressways, along with low-grade roads in peripheral regions, unexpectedly demonstrated a tendency of high temperature clustering. Narrower road widths combined with dense adjacent buildings create confined urban canyons, reducing convective cooling and trapping longwave radiation. While overall traffic volumes are lower, peak-hour congestion in these connector roads generates localized heat spikes from idling vehicles and tire–pavement friction. This phenomenon warrants further in-depth investigation, as it may have implications for understanding the complex relationships between road infrastructure and thermal conditions in urban areas. 4.3. Linear Modeling of the Impact of TCE on RST Figure 12 presents the spatial distribution of road network carbon emissions across different time periods (0:00–1:00, 7:00–8:00, 9:00–10:00, and 11:00–12:00). At night (0:00–1:00), overall emission levels were low, though some carbon emissions persisted on major roads. The 7:00–8:00 and 9:00–10:00 periods fall within the morning rush hour. Traffic volumes rose significantly during these times, leading to increased emissions. By 11:00–12:00, traffic volumes and corresponding emissions showed a gradual decline. Across all periods, the spatial distribution of emissions remained consistent with the overall daily pattern. High-grade roads consistently exhibited the highest traffic volumes and carbon emissions. TCE and RST data for the 9:00–10:00 period were adopted for subsequent analysis. The regression coefficients indicate that TCE and building area were positively correlated with RST, while the green space area was negatively correlated. The explanatory power of the OLS model, as indicated by the R 2 values, ranged from 0.04 to 0.27 across different road grades. Similarly, the SEM had R 2 values between 0.02 and 0.16. These relatively low R 2 values suggest that both the OLS model and SEM had limited capabilities in explaining the variations in RST, likely due to their inability to effectively capture spatial dependencies in the data. In contrast, the SLM demonstrated significantly improved performance, with R 2 values ranging from 0.63 to 0.76 across different road grades. This notable enhancement can be attributed to the ability of the SLM to account for spatial autocorrelation by incorporating a spatially lagged dependent variable, which captured the influence of neighboring observations on RST. This spatial consideration contributes to the better reflection of real-world dynamics where environmental factors are often spatially dependent. Furthermore, the performance metrics, including the Akaike information criterion (AIC) and log-likelihood score, revealed clear improvements in the SLM compared to the OLS model and SEM. The lower AIC values and higher log-likelihood scores of the SLM indicate a better model fit and greater explanatory power, reinforcing the importance of considering spatial effects when modeling RST. 4.4. Nonlinear Modeling of the Impact of TCE on RST Table 6 illustrates the core parameters involved in the process of tuning the RF model. Grid search combined with spatial five-fold cross-validation was utilized to systematically optimize the hyperparameters of the model. Spatial stratification constraints were imposed to ensure spatial independence between training and validation folds. The total number of parameter combinations was 1152. For each combination, the model was trained five times under spatial cross-validation, resulting in a total of 5760 trained models. The optimal parameter configuration was selected based on the R 2 score on the spatial validation set. Table 7 presents the model performance metrics, including R 2, MSE, and RMSE, for both the standard and enhanced RF models. The R 2 values for the standard RF model ranged from 0.41 to 0.66, the MSE values ranged from 2.50 to 3.69, and the RMSE values ranged from 1.58 to 1.92. In contrast, the enhanced RF model demonstrated significantly enhanced performance, with R 2 values ranging from 0.88 to 0.91, MSE values between 0.41 and 0.76, and RMSE values ranging from 0.64 to 0.88. These results collectively suggest that the model achieved adequate fit without evidence of overfitting. The substantial improvement in the performance of the enhanced RF model might be attributed to the incorporation of a spatial weight matrix based on longitude and latitude coordinates. This enhancement allowed the model to better capture spatial dependencies and heterogeneities that were not adequately addressed by the standard RF model. As described in Section 3.3, the spatial weight matrix was constructed using a distance-based approach, in which spatial proximity influenced the decision-making process of the model. By integrating this spatial information, the enhanced RF model was characterized by effectively reduced prediction errors and significantly increased explanatory power, as evidenced by the higher R 2 and lower MSE and RMSE values. Residual distribution analysis revealed that model residuals approximated normality (Shapiro–Wilk test, p = 0.07), with a mean residual of −0.02 ± 0.78 standard deviations and no significant heteroscedasticity. The use of spatial five-fold cross-validation further reduced the risk of overfitting caused by spatial data leakage. Figure 13 illustrates the contributions of different factors to RST. A clear upward trend could be observed in the contribution of road TCE to RST as road hierarchy decreases, with lower-grade roads showing a stronger influence compared to higher-grade roads. For freeways, the contributions of RST, surrounding building coverage, and green space area to LST were 0.40, 0.24, and 0.36, respectively; for urban expressways, the contributions were 0.37, 0.23, and 0.40, respectively; for main roads, the contributions were 0.54, 0.20, and 0.25, respectively; and for minor arterial roads, the contributions were 0.56, 0.22, and 0.22, respectively. These results highlight that as the road grade decreased, the relative impact of TCE on RST became more significant, while the influence of the surrounding building and green space areas tended to decrease. 5. Discussion and Limitations 5.1. TCE and RST Variations Based on Road Segments The research on the distribution of TCE based on road classification has more practical significance. From the perspective of total traffic emissions and spatial distribution, there are obvious differences and uneven distribution among different grades of roads. The total mileage of main roads is considerably longer than that of freeways and urban expressways. In terms of emission intensity per unit, main roads exhibit similar levels compared with high-grade roads. However, the width and number of lanes of the main roads are smaller than those of expressways and urban expressways. Under restricted road width, main roads are more prone to traffic congestion, frequent vehicle stop-and-go, idling, and low-speed cruising, which directly lead to substantial growth in carbon emissions. Although the mileage of trunk roads is only 1.28 times higher than that of secondary trunk roads, their total carbon emissions are several times higher. Therefore, the emission characteristics of different road grades are fundamentally different, and simple aggregation analysis will lose critical structural information. Compared to the extant research, which has typically focused on higher-grade roads [ 62], this study went further by considering lower-grade roads, an area often overlooked due to data limitations. From the spatial distribution and hot-spot/cold-spot analysis of RST, a series of contrasting and counterintuitive key characteristics have been identified. First, the spatial pattern of summer temperature increase shows a significant mismatch with the hot-spot and cold-spot distribution. From the temperature increase map, the central urban area (especially the Futian core district) presents the highest warming amplitude across the city. The large contiguous high-increase zone intuitively leads to a judgment that the central area must be a high-temperature aggregation region. However, the hot-spot analysis completely breaks this expectation. The central area does not form large-scale and significant high-temperature hot spots. Instead, the Futian area in the south-central part shows extensive low-temperature cold spots. This contradictory pattern may be attributed to the strong regulatory effect of the local urban environment on thermal conditions. As the earliest developed urban area, the Futian core district has high building density and intensive human and vehicle activities. Nevertheless, this district contains numerous parks, green spaces, and urban open areas. In addition, the building layout of the old town may form natural ventilation corridors. These factors may create a significant cooling buffer for the road surface thermal environment. As a result, road segment temperatures can hardly form contiguous high-temperature aggregations. This area finally presents a unique pattern: high overall warming amplitude, no significant high-temperature aggregation, and even local cold spots. Second, the hot-spot and cold-spot aggregation of RST is fundamentally different from that of TCE. Compared with carbon emission hot spots, RST hot and cold spots show completely distinct spatial features. Carbon emission hot spots are mostly large-scale, massive, and continuously extended along main traffic corridors. They may belong to activity-driven aggregation caused by dense traffic volume and frequent congestion. In contrast, RST hot and cold spots are fragmented, small-scale, and localized. They are mostly scattered in small clusters within the local road network and rarely form cross-regional contiguous aggregations. More importantly, the hot and cold spots of RST are almost completely spatially mismatched with those of carbon emissions. High-value hot spots of carbon emissions are highly concentrated on main roads, commuting and freight corridors with heavy traffic and frequent congestion. They are a direct reflection of traffic activity intensity. In contrast, RST hot spots mostly appear at the end of the road network in old urban areas. RST may be controlled by micro environmental factors such as street canyon geometry, building envelope structure, local ventilation conditions, and green space distribution. Therefore, RST exhibits stronger spatial heterogeneity and local dependence. This spatial mismatch clearly indicates that TCE and RST do not have a simple linear correspondence. The strong local spatial heterogeneity and non-stationarity of RST mean that traditional analysis methods ignoring spatial effects can hardly capture its real driving mechanism. 5.2. Relationships Between TCE and RST The introduction of a spatial weight matrix effectively address
