Spatial Interconnections of Land Surface Temperatures with Land Cover/Use: A Case Study of Tokyo

: As one of the most populated metropolitan areas in the world, the Tokyo Metropolitan Area (TMA) has experienced severe climatic modiﬁcations and pressure due to densiﬁed human activities and urban expansion. The surface urban heat island (SUHI) phenomenon particularly constitutes a signiﬁcant threat to human comfort and geo-environmental health in TMA. This study aimed to proﬁle the spatial interconnections between land surface temperature (LST) and land cover/use in TMA from 2001 to 2015 using multi-source spatial data. To this end, the thermal gradients between the urban and non-urban fabric areas in TMA were examined by joint analysis of land cover/use and LST. The spatiotemporal aggregation patterns, variations, and movement trajectories of SUHI intensity in TMA were identiﬁed and delineated. The spatial relationship between SUHI and the potential driving forces in TMA was clariﬁed using geographically weighted regression (GWR) analysis. The results show that the thermal environment of TMA exhibited a polynucleated spatial structure with multiple thermal island cores. Overall, the magnitude and extent of SUHI in TMA increased and expanded from 2001 to 2015. During that time, SUHIs clustered in the compact residential quarters and redevelopment/renovation areas rather than downtown. The GWR models showed better performance than ordinary least squares (OLS) models, with Adj R 2 > 0.9, indicating that the magnitude of SUHI signiﬁcantly depended on its neighboring geographical setting, including land cover composition and conﬁguration, population size, and terrain. We suggest that UHI mitigation in Tokyo should be focused on alleviating the magnitude of persistent thermal cores and controlling unstable SUHI occurrence based on partitioned or location-speciﬁc landscape design. This study’s ﬁndings have immense implications for SUHI mitigation in metropolitan areas situated in bay regions.


Introduction
In the 21st century, tackling climate change is arguably considered the top priority to ensure human comfort and geo-environmental health [1]. The Intergovernmental Panel on Climate Change (IPCC) and the United Nations (UN) Framework Convention on Climate Change (FCCC) has emphasized the importance of addressing urban climate issues and enumerated the enormous consequences [2,3], such as global warming, weather extremes, etc. In particular, due to densified human activities and intensified urban sprawl, the climate crisis in the urban context is continuously amplified [4]. Urban areas have experienced overwhelming risks and pressures [5,6]. The underlying mechanisms and processes of the urban thermal-energy system are thereby altered, and the urban heat island (UHI) effect is gradually generated [7]. The UHI phenomenon is driven by surface temperature cliffs between urban and non-urban areas [8][9][10][11]. In light of this, extremely high temperatures in populated urban areas, in contrast with their surrounding rural spaces, have been a severe climatic threat [11,12]. The UN projected that more than twothirds of the world population would dwell in urban areas by 2050 [13,14]. These people would be exposed to high-temperature anomalies as well.
As one of the most vibrant, populated, and livable regions across the globe, the Tokyo Metropolitan Area (TMA) has been struggling to overcome the problems caused by urban expansion and extreme climate events since the adoption of the Kyoto Protocol (1997) [15][16][17][18][19]. As of 2015, the population of Tokyo reached 9.273 million, and the growth rate reached 14% compared with 2000 [20]. With a population density of 6169 persons per km 2 , Tokyo is home to 11% of Japan's total population but only covers less than 1% of the total area of the country [21]. The increase in population is mirrored by urgent demands for urban services and environmental supply. Thus, it is vital to foster a sustainable Tokyo, which contributes to a vigorous metropolitan area that embraces ecological harmony and economic prosperity. As Japan's central locus, Tokyo has always been struck by the UHI effect due to drastic metropolization, and it has acquired a series of problems, such as increased urban energy consumption, increased greenhouse gas emissions, and deteriorated urban air quality [22]. The centennial increment of Tokyo's annual average temperature was estimated at approximately 3°C by the Japan Meteorological Agency (JMA) [23][24][25][26]. Records of unseasonable high-temperature days/nights in Tokyo have been soaring [27,28]. Out-of-control UHI could wipe out hard-won gains for the urban environment and human health. Therefore, alleviating the UHI effect is highlighted as one of the most prominent challenges in Tokyo. Evaluating the evolutionary track of the urban thermal environment in Tokyo and UHI mitigation strategies yields experiences and lessons for other expanding metropolitan areas in bay regions.
Since the 2000s, the Tokyo Metropolitan Research Institute for Environmental Protection (TMRIEP) has set up massive meteorological observation frameworks in order to expand UHI investigations [25,29]. Plentiful and valuable findings on urban/suburban atmospheric heat conditions in Tokyo have been obtained from high-spatial-density monitoring stations [25,27,30]. This kind of canopy-layer UHI (CUHI) investigation based on isolated meteorological stations has contributed to UHI mitigation initiatives. Assessments of CUHI have identified the thermal impacts on specific locations over time using absolute atmospheric temperature. However, they ignored the spatial heterogeneity of UHI due to sparsely distributed in situ measurements [10]. In contrast, by estimating land surface temperature (LST), satellite-derived surface UHI (SUHI) monitoring at various spatiotemporal resolutions replenishes large-scale and long-term UHI observations. Satellite thermal infrared (TIR) remote sensing data, such as Landsat constellation and Moderate Resolution Imaging Spectroradiometer (MODIS) archives, have been extensively utilized for surface temperature measurement due to their short revisit intervals and broad spatial coverage [10,31].
Unlike in situ atmospheric temperature measurements, retrieved LST using Landsat data provides an in-depth evaluation of the thermal situation and surface parameter responses in the high spatial-sensitive context [8,10,31,32]. Together with various ecological, environmental, and socioeconomic parameters, the detailed information of radiative surface temperature can be interpreted to study underlying surface conditions and interactions between atmospheric circumstances and human activities [10]. Until recently, an incredible number of UHI-related scientific studies based on Landsat satellite data have explored the linkages between SUHI and its potential driving force from urbanization [10,31,[33][34][35][36][37][38][39][40]. Deilami et al. [41] identified the association between SUHI and land cover patterns in Brisbane based on spatial regression models at comparative cross-sectional and longitudinal perspectives. They found that the influences of population density and porosity vary significantly over space. Li et al. [42] depicted the spatiotemporal patterns of SUHI in Hangzhou and confirmed that increasing population and green space substantially affect changes in SUHI. Liu et al. [43] evaluated the impact of land cover/use on SUHI for 10 Chinese megacities and demonstrated that the land cover composition and the population had notable effects on altering the intensity of SUHI. Wang et al. [44] simulated land use/cover scenarios and impacts on land surface temperature in Sapporo and stated that urban expansion was the main driver affecting SUHI. Wang et al. [45] investigated the connections among urban landscape, population, and LST in three megacities along the Yangtze River (Chongqing, Nanjing, and Wuhan) and concluded that the effects of landscape composition and configuration on LST differed across cities.
However, few studies have depicted the spatiotemporal evolution of thermal gradients in Tokyo based on Landsat monitoring. The spatiotemporal interconnections between SUHI and its urbanization-induced influencing factors in Tokyo have also been scarcely interpreted due to the complicated urban setting. Tokyo's urban areas are contiguous with each other or border on urban fringes or suburban areas rather than rural areas. It is challenging to delineate and harmonize the contours of UHIs without knowing the distinct temperature differences between typical urban and rural sectors. Although several studies effectively identified the spatial patterns and characteristics of the thermal environment in Tokyo [46][47][48], the temporal dynamic component of land cover/use on LST in Tokyo was disregarded. Furthermore, the thermal environment mechanism is complicated [8,49,50] and impacted by mixed factors [10,31]. In light of this, interpreting the evolution of Tokyo's thermal environment requires a novel investigation into the dynamism of diverse underlying factors. In this context, a more comprehensive assessment of Tokyo's thermal environment is required in order to provide in-depth insights into the spatiotemporal dynamics of the influencing factors on urban thermal comfort. In addition, profiling how Tokyo has managed heat islands through urban planning can enrich the knowledge of urban sustainability and contribute to urban climate responses in the broader international context.
Hence, this study's chief objective was to conduct an integrative assessment of spatial interconnections between LST and land cover/use using multi-source long-term remote sensing data in Tokyo. The second objective was to capture and characterize the spatiotemporal variations and trends of the urban thermal environment for comfort and sustainability. We intend to provide a cross-sectional and longitudinal investigation that permits spatiotemporal evaluations of Tokyo's thermal environment based on land cover/use scenarios, coupled with remotely sensed data from multiple sources and spatial statistical analysis. This study's findings can support the validity of the analytical approaches undertaken and the universality of the research process for other metropolitan areas at the national and global scale, such as Osaka, Fukuoka, Shanghai, Hong Kong, Singapore, Los Angeles, etc.

Study Site
This study was conducted in the context of the Tokyo Metropolitan Area (TMA), Japan. We focused on the landscape and its influence on LST within a radius of 40 km from Tokyo Station (the city center; latitude: 35 • 40 51"N; longitude: 139 • 46 1"E), as shown in Figure 1. This area encompasses most of the urban-suburban-rural areas in Tokyo with many land details. More than 37 million people reside in the TMA [20]. Due to speedy urbanization, the TMA is a representative region of Japan and the world. The substantial settlement agglomerations and urban infrastructure constructions in TMA have given rise to urban warming and contributed to the climate crisis [25,26]. Under the climate scheme, TMA is situated in a humid subtropical zone with a complex aggregation of diverse land covers/uses. Consequently, TMA experiences hot and humid weather from May to early October. In recent decades, TMA has been subjected to frequent severe high temperatures in summer and autumn. The highest recorded temperature ranged between 34.2 • C and 39.5 • C from 2001 to 2015 [20]. This increased tropical weather resulted in a surge of heat-related health problems and destroyed the balance of the thermal environment.

Working Flowchart
This study used multi-source spatial data to examine the causal mechanism and associations between LST and land cover in TMA based on cross-sectional and longitudinal spatial regression models. The overall technical processing scheme is shown in Figure 2.

Land Cover/Use Specifications
In this study, land cover/use thematic maps of TMA were obtained from the global land cover (LC) products of the European Space Agency (ESA) Climate Change Initiative (CCI) [51]. The annual CCI-LC products at 300-m spatial resolution have been treated as persuasive data in a broad range of LC-related research [51,52]. The CCI-LC classification scheme is defined as 37 hierarchical categories, including 22 global and 15 regional classes based on the Land Cover Classification System (LCCS). The LCCS, developed by the UN (http://www.fao.org/docrep/003/x0596e/x0596e00.HTM), delivers great LC details and has flexibility in adjusting the thematic legend in order to serve diverse applications. The validation framework of CCI-LC products was established globally on international standards, using both the pixel-based uncertainty value and object-based expert verification of 2600 primary sampling units. The overall weighted-area accuracy of the CCI-LC 2015 product is assessed to be 71.1% [51].
The thematic details in TMA were extracted from global LC maps and rectified to northern Universal Transverse Mercator map projection zone 54 (UTM zone 54N). The LC thematic legend was rearranged into five categories considering TMA's urban setting and land use utilization together with SUHI profiling-agricultural land (AL), forest land (FL), mixed land (ML), urban fabric (UF), and water area (WA). The specifications of TMA's LC legend are shown in Table 1 reprocessed LC maps were then used to compute the landscape pattern metrics based on the landscape-level moving window in FRAGSTATS 4.2 developed by the University of Massachusetts, making it possible to quantify the relationship between the configuration of land cover/use and LST in the context of spatial heterogeneity [53].

LST Retrieval and SUHI Assessment
The multitemporal Landsat constellation images of TMA (Landsat collection 1 Level 2 (On-Demand); path: 107, row: 35) were utilized in this study. All Landsat scenes were geometric-corrected, radiometric-calibrated, and atmospheric-corrected by the United States Geological Survey (USGS; https://earthexplorer.usgs.gov/). The corresponding image products were also retrieved, including multispectral surface reflectance, spectral indices, and atmospheric brightness temperature in Kelvin (K). Considering the availability and comparability of data, time-series Landsat data in the same season (autumn in this study) or within three-month intervals were given priority. The selected images were cloudfree or had minimum cloud cover (<10%). In addition, the local hourly meteorological data archives from the JMA were used to confirm the existence of extreme weather factors in the 24 h before and after image acquisition. Based on these considerations, we utilized Landsat data from 24 September 2001, 16 October 2006, 17 September 2013, and 9 October 2015 in this study.
Initially, LST was the primary indicator to capture the spatiotemporal pattern of SUHI in TMA. Many studies have provided detailed descriptions of radiative transfer equation (RTE)-based LST retrieval using Landsat thermal bands [36,54,55]. The RTE can be expressed as Equation (1). Briefly, RTE-based LST retrieval involves three main procedures: (1) convert Landsat digital numbers (DNs) to top of atmosphere (ToA) radiance, (2) calculate at-sensor brightness temperature (BT), and (3) compute LST [56,57].
where L sensor,λ is Landsat spectral ToA radiance at the wavelength λ µm −1 , with W/m 2 ·sr·µm as the measurement unit; B λ (T) is Planck's law-defined equivalent blackbody radiance (W/m 2 ·sr·µm) at temperature T in Kelvin (K) and wavelength λ in µm −1 ; τ λ is the atmospheric transmissivity of Landsat sensor-land surface, unitless; ε λ is land surface emissivity at the wavelength λ µm −1 , unitless, and ε λ was estimated based on fractional vegetation proposed by Sobrino et al. [56] through the normalized difference vegetation index (NDVI) threshold method; L atm↑,λ is upwelling atmospheric radiance, with W/m 2 ·sr·µm as the measurement unit; and L atm↓,λ is downwelling atmospheric radiance, with W/m 2 ·sr·µm as the measurement unit. Three important parameters can be obtained from the module of atmospheric correction by NASA (https://atmcorr.gsfc.nasa.gov) [57][58][59]-atmospheric transmissivity (τ λ ), upwelling radiance (L atm↑,λ ), and downwelling radiance (L atm↓,λ ). The details of thermal image information and atmospheric parameters are shown in Table 2. LST-related studies are often highly climate-sensitive, and the investigated results are sensitive to the weather effect [42,43]. The responses of LST values are varied and complicated at different locations and dates. The evolution of the thermal environment in different phases should be investigated using the same statistical frame and harmonized legend. Although the selected Landsat images in this study were within the same season, it is meaningless and inappropriate to directly compare the spatiotemporal variations of the thermal environment for periods using absolute LST values due to enormous uncertainty and bias. Therefore, instead of absolute LST observations, spatiotemporal patterns of SUHI intensity were utilized in this study for analyzing and comparing the evolution of TMA's Remote Sens. 2021, 13, 610 7 of 26 thermal environment. The intensity of SUHI (SUHII) was measured by the estimated LST differences in urban/rural areas. The extrema of LST values were removed to eliminate uncertainty and errors in LST retrieval, making the temperature information comparable. The generated LST values were normalized to the range of 0-1 using the min-max rescaling method. For estimated LST values, the minimum and maximum values of LST were transformed into 0 and 1, respectively; the other LST values were rescaled into a decimal between 0 and 1 [43,54]. Herein, SUHII was computed as Equation (2), which subtracts the mean temperature value of non-urban area (non-UF category pixels) from the temperature value of each pixel in the study region [43,60,61] as follows: where LST urban,i and LST rural,i are the LST values of urban and rural pixels, respectively; LST i is the LST value of pixel i within the study region; and LST non−u f denotes the mean LST value of non-UF pixels. UF category pixels are substituted for urban pixels, while nonurban pixels contain the pixels of other land cover/use categories (AL, ML, FL, and WA).

Spatial Aggregation Pattern and Variations of SUHI
Since the adverse effect of SUHI is a consequence of urban sprawl and human activities over time [8,11,31,62], the relationship between the urban thermal environment and urbanization needs to be highlighted and interpreted to delineate the spatial clustering pattern of SUHII [42]. The variations in the spatiotemporal aggregations of SUHII in TMA mirror the evolution of the urban spatial layout from different angles. Herein, Moran's I (MI) [42,63] was computed to reflect the geographic pattern of TMA's thermal environment on a global scale (Equation (3)) as follows: where n is the total number of spatial analytical units, x i is the mean SUHII value for spatial unit i, w ij denotes the spatial weight between spatial units i and j, and x represents the mean SUHII value within the whole study area. The SUHI data were conceptualized and structured through spatial sparse weight matrix techniques. Using the inverse distance weighted strategy, the neighboring features were varied and influential, and spatial weights were computed to reflect the variations of the thermal environment. When the global thermal pattern of TMA is spatially autocorrelated, the occurrence of positive MI hints at spatial clusters of high or low SUHII values; conversely, negative MI implies a spatial dispersion of high and low SUHII values. In the case that MI has a value of 0, the global thermal pattern is spatially random. This study also adopted hot spot analysis to delineate the local spatial aggregations of SHUII in TMA based on Getis-Ord Gi* spatial statistics. The calculation of the Getis-Ord local statistic (G * i Z score ) is based on Equation (4) [63,64] as follows: where n is the total number of spatial analytical units, x j denotes the SUHII value for spatial unit j,X and S indicate the mean and variance of SUHII values within the whole study area, respectively, and w ij is the spatial weight between i and j.
The tool for hot spot analysis in ArcGIS 10.4, developed by Environmental Systems Research Institute (ESRI), was employed to estimate the Getis-Ord Gi* statistic for each spatial analytical unit of the SUHII data. The existence and locations of local clusters with statistically significant for a specified distance were identified spatially by the results of  (Table 3). This step was conducted by looking at each grid in the context of spatial neighbors. In terms of statistical significance, a larger and positive G * i Z score represents a significant hot spot that has more intense local clustering of high SUHII values, while a smaller and negative G * i Z score signifies a significant cold spot that has more intense local clustering of low SUHII values. When the G * i Z score approaches zero, it indicates the absence of apparent spatial clustering on the local scale. A confidence level of 90% or higher was utilized in further analysis. Additionally, the spatial centroid [65,66], together with standard deviational ellipse (SDE) [67,68] analysis, was implemented to explore the response of the thermal environment to the evolution of the urban landscape. The methodological details of centroid coordinates and SDE parameters can be found in Xu et al. [69] and Johnson et al. [70]. The spatial centroid movement of SUHI hot spots in 2001, 2006, 2013, and 2015 was calculated for a hint at the evolutionary trajectory of TMA's urban thermal environment. The UF centroids in the corresponding years were created to indicate the direction and track of urban expansion in TMA during the research periods. Based on the spatial pattern of SUHI hotspots, SUHII weighted SDEs across time series were also calculated to enhance our understanding of the spatiotemporal dynamic processes of TMA's urban thermal environment.

Geographically Weighted Regression Analysis
Linear regression analysis is intended to signify a dependent variable by a suite of independent variables with a linear weight combination [38,[71][72][73][74]. Referring to the issues of spatial non-stationarity and location independency, many studies, particularly those concerning geography and climatology [71,72,75], have shown that the geographically weighted regression (GWR) model is more precise and appropriate for quantifying the relations between dependent and independent variables than conventional regression models (e.g., ordinary least squares (OLS)) [41,72,76,77]. For instance, Zhao et al. [76] examined the relationship between SUHI and underlying biophysical factors in Austin and San Antonio, Texas, using (and comparing) OLS and GWR models, indicating the non-stationary spatial effect of relevant underlying biophysical factors on SUHI. Li et al. [72] assessed the relationships between diurnal/seasonal SUHI and underlying drivers in 419 major cities all over the world using the GWR model and compared the results to global regression models. Thus, recently, GWR is being utilized more to mirror spatial heterogeneity and geographical connections [10]. In GWR, the complicated spatially varying relationships of dependent and independent variables can be estimated and simulated using a locally linear form of a parameter matrix, which can be expressed mathematically as [41,43,72,74,76] where i indicates the location of a dependent variable (mean SUHII) for each spatial unit, labeled as the coordinate (u i , v i ); k represents the number of explanatory variables (driving forces) within a spatial unit, y i , x ik , and ε i stand for the values of the dependent variable, the kth explanatory variable, and the random disturbance for spatial unit i, respectively; and β i0 (u i , v i ) and β ik (u i , v i ) represent the estimates of the intercept parameter and the local regression coefficient for the kth explanatory variable for spatial unit i, respectively. The performance of the GWR model is sensitive to spatial kernel functions and the bandwidth statistical method. A set of evaluation parameters was computed based on ArcGIS 10.4 and RStudio for verifying the reliability and accuracy of the model, including the coefficient of determination (R 2 ), corrected Akaike information criterion (AICc), MI of the residuals, variance inflation factor (VIF), F-test, etc.

Selection of Spatial Determinants
The interconnections between SUHII and its influencing factors in TMA were clarified and highlighted using spatial regression analysis (OLS and GWR). Thus, the average SUHII of each spatial analytical unit was calculated as the dependent variable. Prior knowledge suggests that an array of potential explanatory variables can be considered [41,43,[76][77][78]. The candidates for explanatory variables are listed by type as follows: (1) Land cover/use composition variables-In terms of land cover/use composition, the proportion of urban fabric area (UFP), normalized difference built-up index (NDBI), forest proportion (FP), normalized difference vegetation index (NDVI), and water proportion (WP) were considered as candidate explanatory variables. UFP, FP, and WP came from the UF and FL categories of reprocessed CCI-LC layers. NDVI and NDBI were calculated using Landsat multispectral images; (2) Landscape metric pattern variables-Determining the size, morphology, and spatial arrangement of urban landscapes is vital to explain urban temperature anomalies. Here, four landscape metric parameters were chosen to quantify the characteristics of diversity, aggregation, and evenness in urban landscapes: Shannon's diversity index (SHDI), contagion index (CONTAG), patch density (PD), and patch richness (PR). Raster maps of landscape metrics were involved in further tests of explanatory regression; (3) Population variable-Population or population density indirectly affects the urban thermal field. Along with the vast production of anthropogenic heat, the concentration and overcrowding of the population impose serious pressure and demands on urban settlements and infrastructure constructions. This is why we incorporated the population as one of the essential determinants of SUHI formation. In this study, we used spatial demographic data at a high resolution of 100 m provided by the World-Pop Project, University of Southampton, UK (https://www.worldpop.org/) [79], to map the population distribution in TMA in 2001, 2006, 2013, and 2015. Numerous studies have gained valuable findings using this population data archive [80,81]; (4) Terrain variables-The urban terrain is an important influencing factor of the stark temperature difference between urban and rural zones. Fluctuations of the topography alter the intensity of solar radiation and the thermal properties of surface materials. Here, elevation, slope, and aspect were included as potential explanatory variables. The elevation data at 3 arc-second resolution were extracted from the digital elevation model (DEM) archives of USGS's Shuttle Radar Topography Mission (SRTM) [82]. Slope and aspect are calculated in ArcGIS based on the elevation.
Subsequently, for each analytical unit, the mean values of the dependent and explanatory variables were extracted and featured as spatial attributes. All data can be aggregated to perform spatial regression analysis.

Implementation of Spatial Regression Model
Regarding spatial regression, we generated a set of 1000 m × 1000 m grid cells as the appropriate analytical scale for observation. At an analytical scale of 1000 m, the spatial regression analysis can retain sufficient spatial information and weaken the scope of spatial dependence and autocorrelation, according to the findings of previous studies [43,83].
Initially, we carried out an exploratory regression analysis to determine the specific OLS model, examining all possible combinations of candidate explanatory variables that could best interpret the SUHI effect without a multicollinearity problem. The redundancy, model completeness, variable significance, bias, and OLS performance were all diagnosed and described in order to refine the initial OLS model. The final OLS model was determined and implemented with a range of independent variables, consisting of CONTAG, elevation, population density, UFP, and FP.
In this study, both cross-sectional comparison (OLS and GWR models) and longitudinal investigation (time sequence GWR models) were designed and performed. Ultimately, all spatial regression models were tested to determine the overall model fit, residuals, and improvement of GWR over OLS, with ANOVA. It is expected to verify the superiority of the GWR model in predicting the SUHI effect with underlying urban biophysical determinants over conventional regression models. On the other hand, time sequence GWR investigations provided a plausible explanation for accessing the relationship between land cover/use changes and variations in the magnitude of the SUHI effect. Intense human activities and land cover/use transformations alter the thermal properties of the urban land surface, which further causes the thermal imbalance of the urban system [11,50]. Figure 4 illustrates the spatiotemporal layout of the urban thermal environment across the TMA between 2001 and 2015. Figure 5 plots the curves for UFP, normalized LST (NLST), and SUHI variations based on urban-rural profiling. The spatial arrangement of UF in TMA between 2001 and 2015 exhibited typical urban-rural gradient distribution. The proportion of urbanized area gradually decreased from downtown to suburban areas. Spatially, radiating patterns of the thermal environment appeared in which the heat diffused from the highly urbanized inner-city area with high LST clustering outward to urban fringe areas and outskirts with low LST clustering. However, the peak value of LST was not located in the urban center despite the overall LST tendency of following the urban-rural temperature gradient pattern. We found that the values of NLST in 2001, 2006, 2013, and 2015 were similar. The curve of NLST arose from the origin (urban center) to the 10 km zone. In locations 10 km away from downtown, it climbed to the temperature peak and then gradually declined. The occurrence of NLST valleys and peaks was closely bound up with the accumulation of continuous urban areas. Accordingly, land cover/use exerted a considerable influence on the SUHI effect.

Interconnections of SUHII with Urban Development
This study examined SUHII hot/cold spots with statistical significance based on Getis-Ord local statistics, interpolated using the ordinary kriging method for visualization. Figure 7 shows the spatial locations of SUHII hot/cold spots and their variations. The hot spots were spatially distributed with a polycentric structure since 2001 (Figure 7a  From 2001 to 2015, SUHII hot spots I-III mainly emerged in the western part of TMA, between 10 km and 30 km from the urban core, and densely aggregated near low-rise residential buildings (Figure 7). The amount and extent of hot spot III in 2013 significantly increased compared to other years and those of hot spot II correspondingly decreased. The cold spot areas in TMA between 2001 and 2015 primarily occurred in areas with dense vegetation and water bodies, for example, Tokyo Bay, at the southeast edge of the study area. This finding suggests that abundant green space (forest, etc.) and large proportional water bodies would contribute to alleviating SUHI.
We combined the SUHII spatiotemporal aggregation patterns to yield the dynamics and variations of the thermal environment from 2001 to 2015, as shown in Figure 7e Table A1 of Appendix A. The SDE analysis suggests that the spatial distribution of SUHI hot spots had apparent orientational effects over the whole research period.

Spatial Relationships of SUHII and Land Cover/Use
The VIF values of auxiliary variables are estimated at a threshold of less than 7.5 (Table 4), proving no multicollinearity among auxiliary variables. Ultimately, the SUHI formation in TMA was evaluated and modeled based on the GWR model by the explanatory variables of the CONTAG, elevation, population, UFP, and FP models. The result and diagnostic parameters are outlined in Table 5. The validity of global OLS and local GWR models was judged by adjusted R 2 , AICc, and MI of the standardized residual value (StdResid). For each research period, the global OLS model assessed the statistically significant relationship between SUHII and all determinants (joint F and Wald statistic ρ The modified coefficient of determination (R 2 ) values for GWR models were all more than 0.9 (0.9082 in 2001, 0.9083 in 2006, 0.9322 in 2013, and 0.9132 in 2015). Compared to the OLS model, the MI of StdResid in the GWR model more closely approaches zero, indicating that the GWR model can provide better specifications of the spatial relationship between SUHII and its influencing variables than the conventional OLS model. The estimated F-values of GWR improvement by ANOVA also mirror a significant boost in simulating SUHII with the local GWR model. On account of spatially sensitive geographical variables, the GWR model enables the model parameters to spatially vary at a local scale, showing that corresponding spatial determinants have a non-stationary spatial effect on the SUHI phenomenon. Table 4. Summary of ordinary least squares (OLS) and geographically weighted regression (GWR) coefficients. SE, standard error; VIF, variance inflation factor; CONTAG, contagion index; UFP, urban fabric area; FP, forest proportion.

Variable
Signif. The localized spatial interconnection and differentiation between SUHII and underlying biophysical determinants in TMA from 2001 to 2015 are shown in Figure 9. The UFP has a significant positive effect on SUHII regardless of scale. We observed that the influence of UFP was more prevalent than that of other explanatory variables. The positive relationship constantly dominated the research region in terms of the UFP explanatory variable. However, the spatial relationship between SUHI and UFP was more intense in 2013 than in other years in which the area with higher coefficients (>2.5) significantly expanded. In the UF area from 2001 to 2015, the value of the mean UFP coefficient drastically increased from 2.33 to 4.11. However, the value peaked in 2013 and then dropped to 3.15 in 2015. Overall, the proportion of forest negatively correlated with SUHII between 2001 and 2015. The localized spatial relationship between SUHI and FP was noticeably diminished since 2001. The areas with high negative association had shrunk. In forest land, the values of mean FP coefficients varied in the range of −3.09 to −2.23. The peak of FP influence was in 2001, and the valley was in 2013. Interpreting the spatial relationship between SUHI magnitude in TMA and relevant exploratory variables using the GWR model, we can determine that there were many differences and linkages over time. The spatiotemporal variations of local R 2 ( Figure  10) mirror the spatial heterogeneity between SUHI and its exploratory variables over time. High values of local goodness-of-fit are measured in the region along the Tokyo Bay shoreline, while low values are concentrated in densely built-up areas or areas with IS coverage, such as downtown urban settlements. In addition, forces driving the formation of thermal differentiation in TMA were not static but varied from time to time as the urban environment and developments evolved. The spatially clustered distribution of StdResid derived from OLS and GWR models was identified ( Figure 10) and verified by the MI index. The spatial aggregation patterns of StdResid in the conventional OLS model are more pronounced and discernable than in the local GWR model. The performance of fit coverage, such as downtown urban settlements. In addition, forces driving the formation of thermal differentiation in TMA were not static but varied from time to time as the urban environment and developments evolved. The spatially clustered distribution of StdResid derived from OLS and GWR models was identified ( Figure 10) and verified by the MI index. The spatial aggregation patterns of StdResid in the conventional OLS model are more pronounced and discernable than in the local GWR model. The performance of fit between predicted and observed SUHII exceeds 86%. This indicates that SUHI formation in TMA during 2001-2015 can be simulated and investigated by these five auxiliary variables.

Land Use Policies and SUHI Magnitude
We clarified the spatial and temporal disparities and variations in the SUHI pattern of TMA between 2001 and 2015. The analysis showed that certain biophysical (proportion of urban fabric, forest ratio, etc.) and social (population) characteristics of TMA's urban environment favored the SUHI phenomenon. Spatial differentiation of SUHI formation and configuration by local non-stationarity GWR modeling has profound implications for urban comfort and sustainability in TMA.
Manley [84,85] proposed the term UHI and elaborated on the impacts of urban development processes on London's climate system under the hypothesis of well-defined utilization of urban land. Numerous studies exemplified the relationship between urban land utilization, urban forms, and thermal anomalies [10,31]. As embodied in the previous studies, urban land utilization should be a dominant driving force [8,31,86]. Thus, relevant policies and plans are significantly bound up with SUHI formation and configuration. In the previous analysis, urban thermal gradients between UF and non-UF categories in TMA were examined by a joint analysis of ESA-CCI land cover/use products and Landsatretrieved LST maps. Different from the results of classically defined SUHI investigations, the effects of SUHI in Tokyo were intensively clustered in residential and redevelopment areas, rather than in the central business district (CBD) ( Figure 5). Hirano et al. [87] and Tsunematsu et al. [46,47] assessed the impact of UHI on Tokyo from the perspective of urban energy. They demonstrated that the high-density residential areas presented a severe heat energy agglomeration phenomenon, in agreement with our results. As in TMA, in several Asian megacities (e.g., Hong Kong, Guangzhou, Shanghai, Singapore) [43,86,[88][89][90][91], the maximum SUHI magnitude was not detected in downtown or central commercial areas in recent years. The peak of diurnal thermal anomalies between urban and rural areas frequently emerged in high-density residential or industrial regions of these megacities. To cope with the urban thermal crisis, it is urgent to manage the thermal cores in the urban area of TMA.
The specific policy areas of development planning in TMA were designated for the upgrading and coordinating of regional development based on the National Capital Region Development Plan (https://www.toshiseibi.metro.tokyo.lg.jp/eng/) [92], as shown in Figure 11a. The official planning policies for the development of TMA resulted in a gradient urban-suburban-rural pattern in terms of urban infrastructure and settlement distribution. There are various types of urban land utilization in the built-up areas (BAs), specifically the TMA city proper, oriented for maintaining and improving the urban functions of TMA for socioeconomic activities, industrial manufacturing practices, and modern residential quarters. This phenomenon was difficult to explain using the spatial pattern of land cover/ and/or population shifts. In reality, the regional pattern and magnitude of SUHI in T probably varied due to the effects of land-sea breezes. That was confirmed by a coupl studies [48,[94][95][96]. TMA is situated at the border between land and ocean, in which thermal balance mechanism is highly climate-sensitive and unavoidably influenced by land-sea breeze system.  Green spaces and water bodies have a definitive role in delivering cooling effe Lemus-Canovas et al. [95] reported that the temperature of continuous urban fabric ar in the Barcelona Metropolitan Area was detected up to 2.5 °C higher than the green ur areas. In the context of Bangkok, Jakarta, and Manila, urban green spaces registe temperatures approximately 3 °C lower than impervious surface areas [33]. The nega daytime thermal intensity in suburban green zone preservation areas (SGZPAs) occur It is worth noting that the heat island magnitude of compact high-rise/middle-rise residential areas was less than that of low-rise areas in TMA (Figure 11b). Solar radiance was partially obstructed by tall buildings, even skyscrapers, which was related to the loss of surface temperature. This phenomenon implies that urban morphology and geometry affect the urban thermal field. In addition, the investment in sustainable infrastructure and greenlandscaped buildings has also contributed to relieving the high temperature of downtown TMA [27,93]. Analogously, Tsunematsu et al. [47] described the responses of different urban surfaces to thermal infrared (TIR) energy in downtown Tokyo at 2 m spatial resolution, confirming that the magnitude of heat islands in areas with commercial/office utilization was lower than in areas with compact residential buildings. Moreover, the intense SUHI effect did not occur frequently in the coastal regions of TMA. This phenomenon was difficult to explain using the spatial pattern of land cover/use and/or population shifts. In reality, the regional pattern and magnitude of SUHI in TMA probably varied due to the effects of land-sea breezes. That was confirmed by a couple of studies [48,[94][95][96]. TMA is situated at the border between land and ocean, in which the thermal balance mechanism is highly climate-sensitive and unavoidably influenced by the land-sea breeze system. Green spaces and water bodies have a definitive role in delivering cooling effects. Lemus-Canovas et al. [95] reported that the temperature of continuous urban fabric areas in the Barcelona Metropolitan Area was detected up to 2.5 • C higher than the green urban areas. In the context of Bangkok, Jakarta, and Manila, urban green spaces registered temperatures approximately 3 • C lower than impervious surface areas [33]. The negative daytime thermal intensity in suburban green zone preservation areas (SGZPAs) occurred between 2001 and 2015 due to the presence of urban greenery and open green landscapes. Large bodies of water (e.g., Tokyo Bay) have a notable influence on TMA's thermal environment as well, evidenced by their massive contribution to cold spots. Within the complicated urban system in TMA, such areas play a crucial role in balancing the urban surface energy system, fostering urban comfort with well-laid-out green spaces and water areas, in agreement with our results. Similar green areas and water spaces have also been developed in other megacities (e.g., Barcelona, Bangkok, Seoul, Hong Kong, etc.) [86,[95][96][97][98].
Based on the results of spatial regression modeling, we confirmed that compound effects from multiple driving factors influenced SUHI magnitude in TMA. Observing the spatiotemporal variations of model residuals, we noticed that overprediction of SUHI magnitude was mostly generated in the agricultural and mixed lands of Saitama, Ibaraki, and Chiba. In contrast, underprediction of thermal anomalies consistently occurred in TMA's downtown. This indicates that investigations of the SUHI effect in the downtown area should consider more potential driving factors due to its complexity. The ongoing vertical and horizontal urbanization of TMA [97,98] cannot be neglected in coordinating the urban environment. Thus, we plan to collect more parameters of TMA's urban environment (e.g., urban volume, sky view factor) to study and model the geographic processes of the complicated SUHI phenomenon with more advanced models and algorithms.
In the future, we expect to explore the interannual, seasonal, and diurnal thermal trends of TMA, fusing multisensor and multitemporal data. For example, satellite TIR images should be combined with synchronous data from microwave radiometer measure-ments and in situ observations to yield more reliable resultants. More detailed features of TMA's urban biophysical environment using local climate zones or unmanned aerial vehicles (UAVs) should be considered in further investigations.

Toward a Livable Urban Environment
In this study, we explored the geographic mechanism and process between land use/cover and the magnitude of the urban thermal effect in the urban aggregation of TMA based on spatial statistical analysis. The mechanism of SUHI formation in TMA is extremely complicated, influenced by more than the land cover/use composition and configuration. Nevertheless, we can partially study and interpret the confounding effects of SUHI using combinations of different variables. The detailed differentiation of whether the spatial mechanisms of heat islands in TMA were influenced by underlying physical factors was based on both global and local modeling. This can provide sufficient evidence for ameliorating the intensive thermal effect. Spatially, the geographic process of the SUHI effect in TMA manifested high heterogeneity according to the non-stationarity GWR modeling framework. Since the values of SUHII significantly vary over space, the thermal impact on TMA is context-sensitive and can be locally assessed based on spatial determinants from neighborhoods. That is mirrored in need for localized landscape planning and site-specific policy design to minimize the SUHI phenomenon. There are favorable evidence and references regarding the use of localized policy design to control the formation of SUHI across diverse cities around the world, such as Brisbane [41], Hangzhou [42], Las Vegas [77], Texas [76], and Ljutomer [74].
Furthermore, this study confirmed that SUHI barycenters and hotspot areas were detected in compact low-rise residential quarters through spatial aggregation analysis and spatial centroid movement. Tokyo is dominated by a large ratio of impervious surfaces and compact buildings. Due to scarce land assets and fierce competition for land resource development, it is hard to put into effect urban greening policies in severe thermal spots. In addition, it has been acknowledged that in Tokyo Metropolis there was a deteriorating trend in urban greenery; the area of emerging green coverage varied from 31 [17,20]. Impervious artificial landscapes and green spaces are dominant components in the urban setting. Thus, sustainable infrastructure and green retrofit, as the most appropriate and practical measures, should be increased to mitigate clustered SUHI.
We expect that local authorities in TMA could reinforce urban climate responses and management by adopting a group of mitigating measures, including expanding urban forests along the highways and in urban parks; cultivating and maintaining green corridors surrounding rivers, lakes, or ponds; and promoting small patches of green landscapes neighboring compact residences.
Japanese coastal megacities, such as Sendai and Fukuoka, have also built livable urban environments by adopting urban afforestation strategies [25,99]. In addition, the construction and enhancement of urban ventilation corridors have become important in UHI alleviation as well. The city temperature in dense IS areas can be decreased by enhancing the ventilation cooling effects and controlling the number of green spaces and buildings in SUHI hot spots. Asian coastal megacities such as Hong Kong, Shanghai, and Singapore have also implemented similar countermeasures for SUHI abatement [86,100].

Conclusions
In this study, we targeted the TMA as a crucial case study to evaluate the spatiotemporal response of land cover/use on surface temperature. This study's distinctiveness lies chiefly in its consideration of neighboring plans and localized landscape design on surface temperature in the metropolitan area of Tokyo over time and space based on a multitude of driving forces. In addition, the occurrence of SUHI depends on the compound effects of multiple drivers. Forces with different thermal impact magnitudes are not static but spatiotemporally varied, relying on the local geographical setting and urban development process. Hence, this study highlights the following: (1) The uniqueness and complexity of the SUHI effect in TMA: TMA's thermal environment profile showed an overall upward trend of LST and SUHII between 2001 and 2015. Based on the SUHI aggregation pattern, TMA was recognized as an area with multiple thermal cores. Dense thermal islands occurred in the compact residential quarters and redevelopment/renovation areas, rather than in downtown (CBD). The distribution of tall buildings shielding solar radiation and the contributions of water bodies and sea breezes might be two reasons for this phenomenon; (2) The identification of spatiotemporal interconnections between land cover/use and thermal environment based on global and local regression modeling: Land cover/use composition, landscape configuration, terrain, population density, and afforestation ratios have altered the balance of the urban thermal-energy system in TMA. Synergistic use of OLS (global) and GWR (local) models could effectively detect potential influencing factors and accordingly propose entire/partial SUHI mitigation strategies.
Overall, in order to create a sustainable and livable urban environment, local authorities in TMA should pay attention to offsetting the magnitude of SUHI cluster spots and carefully controlling the formation of SUHIs in non-significant areas. We suggest that urban afforestation be fine-tuned based on partitioned or location-specific landscape design by incorporating the green measures adopted by local communities and the private sector.

Data Availability Statement:
The data that support the findings of this study are available on reasonable request from the corresponding author.

Acknowledgments:
The authors highly appreciate the editors and anonymous reviewers for their valuable and constructive comments and suggestions to improve the manuscript.

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