Determining the Mechanisms that Influence the Surface Temperature of Urban Forest Canopies by Combining Remote Sensing Methods , Ground Observations , and Spatial Statistical Models

The spatiotemporal distribution pattern of the surface temperatures of urban forest canopies (STUFC) is influenced by many environmental factors, and the identification of interactions between these factors can improve simulations and predictions of spatial patterns of urban cool islands. This quantitative research uses an integrated method that combines remote sensing, ground surveys, and spatial statistical models to elucidate the mechanisms that influence the STUFC and considers the interaction of multiple environmental factors. This case study uses Jinjiang, China as a representative of a city experiencing rapid urbanization. We build up a multisource database (forest inventory, digital elevation models, population, and remote sensing imagery) on a uniform coordinate system to support research into the interactions that influence the STUFC. Landsat-5/8 Thermal Mapper images and meteorological data were used to retrieve the temporal and spatial distributions of land surface temperature. Ground observations, which included the forest management planning inventory and population density data, provided the factors that determine the STUFC spatial distribution on an urban scale. The use of a spatial statistical model (GeogDetector model) reveals the interaction mechanisms of STUFC. Although different environmental factors exert different influences on STUFC, in two periods with different hot spots and cold spots, the patch area and dominant tree species proved to be the main factors contributing to STUFC. The interaction between multiple environmental factors increased the STUFC, both linearly and nonlinearly. Strong interactions tended to occur between elevation and dominant species and were prevalent in either hot or cold spots in different years. In conclusion, the combining of multidisciplinary methods (e.g., remote sensing images, ground observations, and spatial statistical models) helps reveal the mechanism of STUFC on an urban scale. Remote Sens. 2018, 10, 1814; doi:10.3390/rs10111814 www.mdpi.com/journal/remotesensing Remote Sens. 2018, 10, 1814 2 of 18


Introduction
With the acceleration of urbanization, more and more urban forest (definition: all the trees within the administrative boundaries of the city) distributed in highly populated environments has become an important part of urban landscapes [1].Urban forest mainly includes urban forest land, trees in parks, trees along streets, and trees dispersed in residential land.In addition to their use for landscape decoration, residents appreciate the services provided by the branches and leaves of urban forest [2,3].For example, the dense shade that limits heat re-emitted from pavement, walls, and adjacent objects, and the direct radiant heat from the sun.Water transpiration through leaves further reduces the radiant heat, and this cooling effect mitigates the urban heat island effect and improves the quality of the urban ecological environment [4].A growing number of countries have started to place large-area monitoring of urban forest on their political agenda.Moreover, the target of forest-resource surveys in China, India, Sweden, and the UK has also expanded from traditional forests to urban forest [5,6].
The large-area monitoring of urban forest currently relies mainly on remote sensing images, combined with ground-survey methods to evaluate the various functional types and different spatial configurations [5,7,8].This method is promising for large-area monitoring of forest resources [9,10].However, the method uses remote sensing images as essential data input and ground surveys as auxiliary verification data instead of integrating multisource data [1].Furthermore, the monitoring studies of urban forest ecological functions are limited to national and global scales, with quantitative studies on the urban scale being rare (the urban scale often means the area within city-level administrative boundaries).The transpiration cooling function of urban forest on the urban scale is strongly heterogeneous in both the spatial and temporal domains, and the quantitative study of spatiotemporal heterogeneity poses a new challenge for ecologists [2,11].As various environmental factors interact (indirect effect), the diversity, complexity, nonlinearity, and heterogeneity of these factors influence the entire process of urban forest growth [12], and interactions between these environmental factors vary in strength [13].The indirect effect of interactions compared with the direct effect from a single factor may have a greater impact on the urban forest cool island effect (the "cool-island effect refers to the air temperature of the surrounding area being cooler than that of the central area" [14]).However, interactions are difficult to quantify by using empirical equations or routine experimental methods.Environmental factors affecting the surface temperatures of urban forest canopies (STUFC) include biotic factors (e.g., dominant tree species, forest age, and canopy density) and abiotic factors (e.g., area, altitude, slope gradient, and site quality) [4,15].Knowledge about interactions between these factors can improve simulations and predictions of the spatial patterns of urban cool islands.
Clarifying the mechanisms that influence the STUFC requires applying a combination of research methods to quantify the spatial heterogeneity of STUFC [16,17].Remote sensing at various resolutions, ground surveys, and spatial statistical models have all been used by researchers to study how STUFC relates to environmental factors such as forest area, canopy density, vegetation richness, and vegetation index [2,[18][19][20].Overall, these studies indicate that the transpiration cool island effect of urban forest is complex and is influenced by many environmental factors.However, only a few studies have comprehensively investigated the quantitative relationship between multiple environmental factors and STUFC [21].The influence of these multiple environmental factors on STUFC at the urban scale and the strength of the interactions between these factors remain unclear [22].Therefore, a quantitative comprehensive study of STUFC spatial heterogeneity and its determining factors is in order.
The present study aims to clarify how multiple environmental factors interact and influence the STUFC.We focus on the interactions that occur in urban forest land and trees in park aggregations because many scholars have confirmed that urban forest land and trees in parks are the most important components that affect the urban cool island effect [23][24][25].Combining Landsat-5/8 Thematic Mapper (TM) images with forest-management planning inventory (FMPI) data, we obtain the spatial distribution of STUFC on an urban scale.To accomplish these aims, we fuse multisource data (including remote sensing imagery and ground observations) and use spatial statistical analysis methods to determine the mechanisms.As a case study, we tested this method by using the urban forest of Jinjiang as a Chinese city representative of one that has undergone rapid urbanization in recent years.We address two main questions: (1) Which direct environmental factors most strongly influence the STUFC?(2) How do interactions between multiple environmental factors influence the STUFC?

Overview
To reveal the factors that influence the STUFC in cluster regions (hot and cold spots, the area with statistical significance that was calculated by using a spatial statistical method), we developed a hierarchical spatial statistical analysis method that could be widely applied elsewhere.The method consists of three main steps (Figure 1): First, we built up a raster and vector database for the integration and fusion of multiple sources of spatial data (including forest-management planning inventory (FMPI), digital elevation models, population, and remote sensing images).Second, we calculated the spatial distribution of STUFC from thermal remote sensing images and then applied hot-spot analysis with an optimal spatial analysis threshold to locate the hot and cold spots of the STUFC.Third, focusing on clustering regions, we qualitatively analyzed the interaction mechanisms related to multiple environmental factors as they affect STUFC.For an overview of the data sources, and other research methods, please see Appendix A in the Supplementary Materials.

Study Area
Jinjiang is located in southeastern China; its geographical coordinates are longitude 118 • 24 -118 • 43 , northern latitude 24 • 30 -24 • 54 (Figure 2).The geographical location of Jinjiang is in the subtropical zone.This city is subject to the oceanic monsoon climate of the southern subtropics, which is characterized by high air temperature, rich light and heat, and abundant precipitation.The mean annual temperature is 18.3-21.3• C; relative humidity is 76%, and mean annual precipitation is approximately 1130-1820 mm.The terrain is generally smooth and consists mainly of plains and hills.The elevation varies from sea level to 518 m.

Creation of a Multisource Spatial Database
Three types of data were used in this study, which covers the years 2004 and 2014.The first type was the FMPI data with 90% sampling accuracy, which was acquired from field surveys and observations and included the graphical data and the corresponding attribute database in 2004 and 2014.FMPI data, which are compiled by the national forestry bureau and its affiliations in China, provide forest ground survey data on a regional scale.The attribute database contains data on urban forest stands, soil, topography, and other environmental factors.The FMPI data used herein are extracted from the provincial FMPI data based on the administrative boundaries of city and mainly include urban forest land and trees in parks, but not trees along streets and water courses or trees in home gardens.The FMPI data for Jinjiang County were provided by the Forestry Department of Fujian Province, China.
By using large-scale sampling methods, this forest resource inventory collects detailed information about the characteristics and conditions of each type of forest.In each forest patch, only trees with a minimum diameter at breast height (DBH) of 8 cm were sampled.The forest inventory investigation recorded the forest data and graphical data from each forest patch to establish forest files for resource management, reasonable management, and continuous use.FMPI data contain (1) stand data (age, patch area, plantation density, DBH, tree height, and volume), (2) soil data (soil depth, humus depth, and site index), and (3) topographical data (elevation, slope degree, slope direction, and slope position).All factors were measured for all units within each forest patch, with the average being used as the factor value for each patch.
The accuracy of forest patch attributes was tested based on differences in volumes by combining systematic sampling and stratified sampling.According to the data, forest patches varied in size from 0.05 to 1400 ha in 2004 and from 0.01 to 1029.18 ha in 2014 and the average area went from 9.64 ± 63.08 ha to 3.90 ± 30.65 ha over the same period (mean ± SD).In 2004, 410.821 ha of a total 10,896.937ha in patches were in small patches (<1.44 ha); in 2014, 887.352 ha of a total of 9565.680ha in patches were in small patches.The total amount of Jinjiang urban forest were 2777 × 10 4 trees in 2004 and 2331 × 10 4 trees in 2014.The stem density was 2471 trees•ha −1 in 2004 and decreased to 2082 trees•ha −1 in 2014.The average DBH grew from 5.73 cm in 2004 to 10.30 cm in 2014.The main species of urban forest include Pinus massoniana, Acacia confusa, Eucalyptus grandis, Casuarina equisetifolia and Dimocarpus longan.From 2004 to 2014, the main species were Acacia confusa and Eucalyptus grandis (the other species decreased in number).On average, in 2004 and 2014, 97.73% and 82.45%, respectively, of all urban forest in Jinjiang was composed of small trees (with DBH < 15 cm) (Figure 3).
The second data type that characterizes anthropogenic activity contains the population density distribution raster data and the percent of impervious area represented by the buffer circles around the forest patches [26], which were calculated by using statistic year books and land-use vector data, respectively.The kernel density tool in the Arctool box of Arcmap 10.1 was used to interpolate the population density.In terms of the patch location, a 100 m buffer was appropriate, and the impervious surface area was calculated by using the spatial calculation tool of Arcmap.
The third data type was inverted STUFC maps from Landsat-5/8 TM remote sensing images in summer.The summer season (July to September) is the optimal time to study the urban forest cool island effect.The dense canopy provides the largest temperature gap between STUFC and non-forest land surface temperature.Moreover, people focus on the cooling effect of forest patches during summer more than during the other seasons.Because the heat island affects human health relatively little in winter, spring, and autumn, the few tree leaves that remain in these seasons provide only limited shade.Thus, images acquired during the summer season are ideal to discern STUFC variations and determine the contributions of various environmental factors over different years.For more detailed information, see Appendix A in the Supplementary Materials.

Inversion of Surface Temperature of Urban Forest Canopies
The STUFCs were inverted from the thermal infrared band (10.40-12.50µm) of Landsat-5 TM images (Row 119/Path 43) acquired on 26 July 2004, 11 August 2004, and 12 September 2004 and of Landsat-8 images acquired on 6 July 2014, 22 July 2014 and 7 August 2014.The six images were acquired at approximately 10:00 a.m.local time, and the study area was always under clear atmospheric conditions.The average summer STUFC patches were calculated by using a single channel algorithm in terms of the three images of the same year [27].
First, TM and thermal infrared (TIR) data were applied to correct radiometric and geometrical distortions to ensure accurate calculations.The root mean square errors of the comparison with the ground control points were less than 0.50 pixels for each of the six images.Second, we used atmospheric correction methods for the data from the TM and TIR to eliminate how view angle and atmosphere affect surface reflectance.Third, we used the cubic convolution interpolation method to resample the thermal infrared band of Landsat 8 resolution from 100 to 30 m and of Landsat 5 resolution from 120 to 30 m so as to satisfy the resolution of the visible spectrum.Fourth, the STUFC was estimated by using the mono-window retrieval algorithm.Finally, the STUFC was normalized to evaluate the interactions among the multiple factors influencing urban forest temperature.Nondimensionalization can eliminate the influence of viewing angle and of the different image-acquisition times.A relative STUFC index T R , proposed by Zhao et al. [28], was used to normalize STUFC as follows: where STUFC is the surface temperatures of urban forest canopies of each pixel and T b is the background surface temperature, which is defined as the average temperature of the study-area images.For details of the processes, see Appendix A in the Supplementary Materials.

Accuracy Assessment
The STUFC is linearly related to the near-surface air temperature [29], which provides a basis for testing remote sensing retrieval of STUFC.Thus, we used the daily average land surface temperature data measured at eight meteorological stations in Jinjiang to evaluate the accuracy of remote sensing retrieval of STUFC.We used the Pearson correlation coefficient to evaluate the linear relationship between observed daily-average land surface temperature data and STUFC estimated from remote sensing images and FMPI.The results show that STUFC correlate significantly with land surface temperature from meteorological stations (Pearson correlation coefficient of 0.903; p < 0.01; see Figure 4).

Determination of Optimal Threshold Distance
To determine the optimal distance for distinguishing the surrounding patches, we calculated the local Getis-Ord G i * statistics with threshold distances varying from 250 to 3500 m in 250 m intervals.The hot and cold spots and the number of neighboring FMPI polygons under these different threshold distances were used to determine the optimal threshold distance for spatial autocorrelation analysis [30,31].The total area of hot and cold spots was relatively stable and the minimum number of neighbors for each patch always exceeded eight, which satisfies the normal distribution assumption of autocorrelation.Based on the G i * normality assumption and the variation at every local neighborhood threshold, we determined the optimal threshold to be 2250 m for both 2004 and 2014 (see Table 1 and Figure S1 in Appendix C in the Supplementary Materials).The global Moran's I and the local Getis-Ord G i * indexes were used to determine the statistical significance and to identify the spatial positioning of hot and cold spot cluster regions of STUFC on an optimized threshold distance [32,33].The geographic analysis and extension module of Arcmap 10.1 was used to calculate both of these indexes.
The most commonly used index in spatial autocorrelation analysis is global Moran's I [34], which uses the following formula: where n is the number of observed space units (polygonal region or point), X is the statistical mean, X i and X j are the observed values in regions i and j, respectively, and W ij is the weight of region i relative to region j.This study uses the spatial statistical method of local Getis-Ord G i * to estimate the ratio of local values with different weights to the overall value based only on distance measurement standards, without considering the condition of boundary adjacency: where d is the region size, W ij is the weight matrix of the sample points and region j, and x quantifies importance.
To compose a graph of STUFC clusters, the local Getis-Ord G i * index was used with an optimal threshold distance of 2250 m.Regions with an absolute value Z greater than 1.65 (p < 0.1) calculated by Arcmap 10.1 [35] were considered significant.The significant areas were divided into hot spots and cold spots according to their Z score (the Z score is the standard deviation), and the remaining regions were considered nonsignificant [32].Hot spots (Z ≥ 1.65) are areas of statistically significant clustering of high STUFC (90% confidence interval).Cold spots (Z ≤ −1.65) are locations of statistically significant clustering of low STUFC (90% confidence interval).Nonsignificant areas are defined as locations showing no significant spatial correlation (1.65 > Z > −1.65).See Appendix A in the Supplementary Materials for more information.

GeogDetector Modeling
The objects of an eco-environmental study are always characterized by spatial heterogeneity or autocorrelation.Heterogeneity may be detected by comparing the sum of the spatial variance of sub-areas with that of the total study area.The spatial correlation is based on the similarity of the spatial distributions between two factors.The GeogDetector model (http://www.geodetector.cn/)uses heterogeneity and autocorrelation to evaluate the contribution of the factors based on their q values [36].It quantifies the interaction between factors X1 and X2 by comparing the sum of their q values with that of a new factor X3 = X1 ∩ X2, which it constructs [37].Unlike basic statistical and geographic modeling tools that constrain data input, GeogDetector processes both categorical and numeric variables in parallel.To summarize, we identified the most dominant environmental factors affecting the STUFC and their impacts and also quantified the interactions between them.
To begin the modeling process, we assumed that multiple environmental factors jointly influence the STUFC and that the spatiotemporal distribution reflects these environmental factors.First, the average values of the attributes of the factors were calculated in terms of the forest patches, and then the values were classified according to the classification criteria of the FMPI documentation (e.g., SDe) [38].The factors and STUFC values that were not classified in the manual were classified by using the Natural Breaks (Jenks) method (e.g., patch area) [39].Finally, we regarded the normalized STUFC value as the explained variable and the classified environmental factor values as the explanatory variables and put them into the GeogDetector model to do the analysis.For additional details, please see Appendix A in the Supplementary Materials.

The Surface Temperatures of Urban Forest Canopies, Stand Structure, and Anthropogenic Activity
The spatial distribution of STUFC reveals characteristics of spatial autocorrelation, which departs from the basic assumption that the traditional statistical methods are independent of each other.Therefore, we use herein spatial statistical methods that integrate spatial information into classical statistical analysis to reveal spatial distributions, patterns, processes, and relationships involving STUFC.This approach also clarified the spatial correlation between location-related STUFC and multiple influencing factors.The aggregated STUFC areas (absolute value of Z score ≥ 1.65 represents statistical significance) identified through the spatial statistical method were more suitable than the non-aggregated areas (1.65 > Z Score > −1.65) to determine the relationship and mechanisms.
Landsat 5 and 8 images from 2004 and 2014 were integrated with the corresponding FMPI data (there were 3300 stand patches in 2004 and 5963 stand patches in 2014).For each of these two years, a multisource dataset was generated from the same boundaries and coordinate system, in which the STUFC of patches corresponded to certain environmental factors.In 2004 and 2014, the relative average STUFC was −0.020 and −0.061, respectively (Figure 5).From 2004 to 2014, as the urban forest structure changed, the standardized STUFC changed accordingly.Specifically, upon using the global Moran's I, our analysis of STUFC was spatially autocorrelated, and its spatial clustering decreased (Moran's I = 0.476 in 2004 and 0.178 in 2014; see Table 2).The total area of urban forest stands decreased by 1331.26 ha from 2004 to 2014 (10,896.94ha in 2004; 9565.68 ha in 2014), while the average forest age increased by 5.09 years during the same time period (15.03 ± 7.91 years in 2004; 20.12 ± 10.86 years in 2014), and the average canopy density increased by 0.046 (0.31 in 2004; 0.39 in 2014).Urban forest stands became more fragmented between 2004 and 2014; the shape index of urban forest patches decreased from 301.5 to 182.4.The areas of dominant tree species such as Dimocarpuslongana, Casuarina equisetifolia, Acacia confusa, Pinusmassoniana, Eucalyptus grandis, and other broadleaved trees also changed.For example, the total stand areas of D. longana, C. equisetifolia, and P. massoniana decreased by 14.89%, 0.94%, and 0.96%, respectively, and the total urban forest areas of A. confusa and E. grandis increased by 4.69% and 18.13%, respectively.See Figures S2 and S3 in Appendix C of the Supplementary Materials for the variation of the urban forest characteristics.To quantify the effect of anthropogenic activity on STUFC population density and proportion of impervious area were used to represent anthropogenic activity.The average citywide population density of Jinjiang in 2004 was 1099.0 ± 511.4 people ha −1 versus 1187.3 ± 602.0 people ha −1 in 2014 (Figure 6).The aggregation degree (Moran's I) of the STUFC decreased from 0.731 in 2004 to 0.642 in 2014.The impervious surface area surged.The proportion of impervious surface area within the 100 m buffer circle increased from 28.13 ± 24.85% to 42.45 ± 15.95%.These changes are very significant in our statistical analysis (p < 0.001, t-test).

Mechanisms that Influence Surface Temperatures of Urban Forest Canopies
The mechanism influencing STUFC contains two parts: direct and indirect effects.Overall, although the dominant environmental factors affecting the direct effects on STUFC in different periods varied in the two periods, forest-stand groups remained the dominant environmental factors.Specially, forest stand and topography groups were the dominant environmental factors that influenced STUFC in 2004.In 2014, the dominant environmental factors were transformed to forest stand and anthropogenic activity groups.Specifically, the patch area and dominant tree species, which, compared to other environmental factors, were the main factors contributing to STUFC clustering (Figure 7, Table 3).The GeogDetector model revealed the mechanisms by which environmental factors impact the STUFC.
To further clarify the mechanisms in the hot and cold spots of STUFC, we used local Getis-Ord G i * statistics.The patch area and dominant tree species were the main factors contributing to STUFC clustering (Figure 7, Table 3).The dominant species not only exerted a strong direct effect, but also interacted intensely with other environmental factors in cold and hot spots in the two periods.Among the different cluster regions and during the different time periods, both linear and nonlinear interactions were observed between the 12 selected factors we studied: q value of factor 1 ∩ q value of factor 2 > the sum of the q values of factors 1 and 2 (effect of factor 1 and of factor 2); q value of factor 1 ∩ q value of factor 2 > q value of factor 1 or q value of factor 2 (effect of factor 1 or of factor 2). See Figure 8 and Figures S4-S6 in Appendix C of the Supplementary Materials.
Meanwhile, environmental factors had varying impacts on STUFC in different clustering regions.For example, in 2014, the main influential factors in hot spots were dominant tree species (q value represents power of the determinant indicator of Geodetector model, q = 0.019), patch area (q = 0.013), and ISP100 ("ISP100" refers to the percent of impervious surface area within a 100 m radius of the patch edge.)(q = 0.012).The main influential factors in cold spots were ISP100 (q = 0.045), dominant tree species (q = 0.032), stand age (q = 0.016), and slope position (q = 0.016).Furthermore, the environmental factors with interactions varied during different time periods.For instance, for the cold-spot area, the major interaction factors changed from DS in 2004 to ISP100 in 2014 (see Figure 7 and Figures S4-S6 in Appendix C of the Supplementary Materials).

Mechanisms that Impact the Surface Temperature of Urban Forest Canopies
The quantitative analysis of the GeogDetector spatial statistical model shows that the dominant species factor explained more STUFC heterogeneity than the other environmental factors, which indicates that differences in transpiration and shading between different tree species are the dominant reasons for STUFC on the city scale.The above result is consistent with the conclusion of green roof-temperature-reduction services [40].The differences between transpiration and shading is mainly related to stomatal conductance and biological traits of leaves [41].Studies mainly focus on the variational characteristics of stomatal conductance of different plant leaves and the influence of environmental factors on stomatal conductance [42][43][44][45][46].However, the physiological mechanisms of leaf stomatal behavior are still not fully understood [46].The comprehensive effects of multiple environmental factors on stomatal conductance have not been quantified, resulting in the use of fairly simple model assumptions and climate parameters in ecological-process models to calculate the transpiration rate of vegetation types [47].
Future research should thus consider the interaction between environmental factors and their cumulative effect on stomatal conductance [43,48].We quantified the effects of tree species and multiple environmental factors on the heterogeneity of STUFC and inferred that the stomatal conductance of individual plants differing from the synergistic effects of multiple environmental factors, might result in differences in the transpiration of different tree species.These differences are cumulatively reflected on the regional scale and have larger impacts on the heterogeneity of STUFC than other environmental factors.
Furthermore, environmental factors not only directly affect the stomatal conductance of plant leaves, but also indirectly affect the shading by changing the biological characteristics of the plant's morphology, structure, physiology, and phenology [49].Studies have confirmed that the effect of shading by trees depends on the canopy structure, such as canopy size, tree shape, canopy density, branching pattern, and leaf density (e.g., leaf area and average leaf dip angle), the optical properties of the leaves, such as size, length, width, thickness, texture, and brightness, and the topography (slope and aspect relative to direction of sunlight) [50][51][52][53].These biological traits on the regional scale correlate significantly with environmental factors [54].However, the research methods are influenced by three main factors, including the size of the study area, the comparability of trees, and the availability of spatial attribute data.Most regional-scale surveys are based on the same environmental gradient consider the different species as a general object to study the relationship between biophysical traits, environmental factors, and shading effects, ignoring differences in the biological traits of various tree species [55].Tsukaya emphasized that the biological characteristics of different tree species led to different response magnitudes and mechanisms in terms of changes in environmental factors [56].For example, the blade size is the result of a combined effect of multiple environmental factors.Smaller blades reduce the boundary-layer resistance and heat loss, which is beneficial to maintain proper blade temperature in arid environments.The analysis of factors through the GeogDetector model, which breaks the basic assumptions of sample independence used in traditional statistical methods and insert spatial factors into the research model.We analyzed how biological differences affect shading on a regional scale.Such an analysis reveals that, compared with other environmental factors, the biological heterogeneity of tree species makes a major contribution to the STUFC heterogeneity.
Meanwhile, we should clearly understand that, at the regional scale, the formation of heterogeneity of STUFC is driven by the major factors and the combined effects of environmental conditions and various biological processes.The GeogDetector model indicates that the spatial interaction between various environmental factors and their linearly or nonlinearity enhanced STUFC in hot and cold spots.Thus, the synergetic interaction between two factors can have a greater influence than the direct effect of a single environmental factor, and the combined effect may even be greater than the sum of the direct effects of the two factors.For instance, in cold spots, the influence of DS and ISP100 impacts the STUFC in the two years.The strongest interaction between the two factors was DS ∩ ISP100 in 2004 but, in spite of the strong interaction ISP100 ∩ DS, this changed in 2014 to ISP100 ∩ SA.As seen in Figure 6, the impervious surface area in Jinjiang city expanded significantly from 2004 to 2014, which might explain the enhanced effect of ISP100 not only in cold spots but also in the hot spots, nonsignificant areas, and in the entire study area.As urban ecological processes continue, the evaluation of such interactions will become even more important [57].Most often, previous studies of urban land surface temperature and its mechanisms have used classical multivariate statistical analyses to quantify the correlations between anthropogenic activity and land surface temperature [4,20].Attributes such as landscape composition (e.g., vegetation coverage, impermeability, water proportion), spatial configuration (e.g., shape index, Shannon diversity index), and land surface temperature have been investigated, but these studies have ignored the complex effects and interactions of other environmental factors (e.g., soil and topography) [58].Moreover, spatial analysis is often neglected in these statistical analyses, whereas nonspatial characteristics and attributes are included.The GeogDetector model uses spatial variance to evaluate geospatial data, in contrast with general statistical analysis methods.Therefore, this model is able to clarify the effects of geographical spatiotemporal changes caused by, and multiple environmental factors combined with, anthropogenic activity on STUFC.
In summary, our results indicate that, in 2004 and 2014, for the STUFC cluster regions (hot spots and cold spots), the interaction between dominant factors that have direct effects and other factors (e.g., dominant species or patch area ∩ other factors) affecting STUFC produces a stronger effect than do the other environmental factors (e.g., canopy density ∩ stand age).The interaction intensity changes with the progress of urbanization, which shows that, after considering the direct and indirect interactive effects, we should focus on the dominant factors that cause direct effect in the future planning of urban green space and the structural design ecological process models.These are the most important factors affecting the STUFC, whereas the interactions between them would have strong spatial and temporal heterogeneity.This is consistent with the hypothesis of the variability of multi-factor interaction strength [54].Specifically, the type of interaction may vary by category (such as linear, nonlinear, and independent), intensity (such as attenuation and enhancement), and symmetry (such as symmetric and asymmetric) based on the characteristics of species or individuals.However, previous studies suffer from knowledge deficits caused by their neglect of interaction between influencing factors that cause spatial-temporal heterogeneity.Recent quantitative studies of interactions have used control experiments and model simulations at the plot and regional scales [59][60][61][62].The key challenge to reveal the mechanisms of macro-ecosystem ecological function is to quantify the interaction process of influential factors across temporal and spatial gradients, and to clarify the temporal and spatial heterogeneity based on quantitative results.

Recommendations for Further Study
Considering our research results, we advance two recommendations: First, STUFC research should combine the respective advantages of remote sensing, ground surveys, and spatial statistical analysis to solve regional-scale urban environmental problems and to elucidate how interactions between multiple environmental factors affect STUFC information pertaining to macro-spatial patterns and micro-forest stands must be integrated in new models, and a multi-source database at a unified spatiotemporal scale must be built.Second, the influential factors in cluster regions of STUFC must be identified in land surface temperature studies.As the dominant factor exerting direct and indirect effects on STUFC, tree species should be considered in future urban forest planning.Future quantitative studies on heterogeneity should be carried out on different spatiotemporal scales, in different cities and seasons, and with remote sensing imagery at various resolutions.This study reveals several environmental factors that are correlated with STUFC in urban areas (e.g., dominant tree species, stand area).However, the limited environmental factors selected in this study cannot fully express the driving mechanisms behind STUFC.Besides, the study is aimed at urban forest land and trees in parks and does not cover trees along streets and dispersed trees in residential land.In fact, although the number of trees along streets and dispersed trees in residential land is much lower than that in urban forest land and trees in parks, it still plays a role in the urban scale cold island effects.However, the resolution in the thermal band for Landsat images is 100-120 m, which is much greater than a street width or a residential green space.The mixed-pixel problem causes the contribution of the PA factor to the STUFC to be underestimated.The smaller the patch, the higher will be the proportion of these mixed pixels relative to pure pixels.Because of this resolution effect, the smaller patches might appear warmer than the large patches.Future research may consider investigating how these components of urban forest influence the cool island effect and use higher-resolution images (e.g., thermal images from unmanned aerial systems [63]).

Conclusions
Understanding the interaction between landscape patterns and processes should be reconsidered by using methods of modeling and analysis that use multiple sources of observation data and multi-disciplinary approaches.
This study advances traditional methods of analysis in this field by developing a novel modeling method that combines remote sensing imagery, FMPI data, spatial statistical modelling, and multiple environmental factors.It then integrates them into a spatial dataset with a unified scale to model urban forest regions with strong spatial heterogeneity.The results demonstrate that the degree of spatial clustering of STUFC decreases during rapid urbanization.Although different environmental factors exert different influences on STUFC, in two periods with different hot spots and cold spots, the patch area and dominant tree species are the main factors contributing to land surface temperature clustering in urban forest.The interaction between multiple environmental factors increases STUFC, both linearly and nonlinearly.
In conclusion, the combination of multidisciplinary methods (e.g., remote sensing images, ground observations, and spatial statistical models) helps reveal the mechanism of STUFC on an urban scale.This is done by quantitatively expressing nonlinear interactive systems on the urban scale.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2072-4292/10/11/1814/s1.Appendix A: The detailed information about the study, Appendix B: The references used in Appendix A, and Appendix C: Figures S1-S6.
Author Contributions: S.Z.interpreted the data and wrote the manuscript; S.D. visualized the data.X.S. provided the data resources.C.X. and Y.L. conceived of the methodology of Geodetector.W.C. took charge of the project administration.Q.C. was the tutor for land surface temperature inversion.Y.L. and J.T. did the investigation of the observation data.W.M. helped to visualize the data.Y.R. conceived the methodology and supervised the data analysis.

Figure 1 .
Figure 1.Flow chart describing the proposed methodology.The data were processed by combining Landsat-5/8 TM remote sensing images, FMPI, and the GeogDetector model.

Figure 2 .
Figure 2. Map of study area.The inset shows the location in China.

Figure 3 .
Figure 3. Changes in FMPI data in Jinjiang City between 2004 and 2014.

Figure 4 .
Figure 4. Comparison of STUFC between on-site observations and Landsat data.

Figure 5 .
Figure 5. Spatial distribution of temperature in 2004 and 2014 obtained using a 2250 m threshold and local Getis-Ord Gi* statistics.Red (GiZscore ≥ 1.65) indicates a hot spot, blue (GiZscore ≤ −1.65) indicates a cold spot, yellow (1.65 > GiZscore > −1.65) indicates a nonsignificant area, and gray indicates altitude.The numbers indicate the main regions of mutual conversion of cold spots, hot spots, and nonsignificant areas.

Figure 6 .
Figure 6.Population density in different cluster regions classified with local Getis-Ord G i * statistics, which indicate population changes between 2004 and 2014.

Table 1 .
Number of neighboring patches found around each hot and cold spot.
2.5.2.Statistical Analysis Based on Global Moran's I and Local Getis-Ord G i *
Note: p-values of all areas are <0.001 in 2004 and 2014; ** indicates high significance.

Table 3 .
Degree of influence of urban forest attributes, soil, topography, and population on STUFC in hot spots (HS), cold spots (CS), and nonsignificant regions (NS) in Jinjiang, China.Influence (q value) of interactions in cold spot areas in 2004 and 2014.