URBAN HEAT ISLAND FOOTPRINT EFFECTS ON BIO-PRODUCTIVE RURAL LAND COVERS SURROUNDING A LOW DENSITY URBAN CENTER

The urban heat island (UHI) is a common effect caused by urbanization and has been studied to evaluate the thermal condition in cities worldwide. However, most previous UHI analyses are performed in major metropolitan cities. This study conducts a spatiotemporal analysis of UHI in a rapidly expanding low-density suburban centre and determines how bio-productive land covers react and the extent of the disturbance to each land cover based on time series land surface temperatures extracted from Landsat 7 ETM+ images. Two methods applied and compared are the single exponential decay method, which measures UHI footprint (UHIFP) on vegetation phenology, and the two dimensional Gaussian surface, which quantifies the influence based on distance from the local urban perimeter. Three spectral indices (Normalized Difference Vegetation Index (NDVI), Moisture Index (NDMI), and the Enhanced Vegetation Index (EVI)) were extracted and the residuals from the Gaussian model were compared based on these indices in order to better understand the thermal variations of each land cover within a UHI. The results show that the UHIFP of the studied low-density suburban centre is 1.4 times larger than the size of the urban centre, marginally smaller than previous analyses performed within high-density metropolises. All vegetated land covers experienced their maximum cooling effects before reaching the UHIFP perimeter while urban surfaces begin to diverge from the Gaussian model outside of the UHIFP. The residuals of sparse vegetation maintained strong correlations with each index throughout the growing season while NDMI retained the strongest relationships with every land cover. This study has helped us better understand the UHI effects of small communities with varied vegetation phonology based on the distribution of built-up pervious and impervious surfaces within the neighbourhood structure. The similar results from both methods indicate a strong urban cover influence overpowering the dominant distribution of agricultural surfaces throughout the growing season.


INTRODUCTION
The urban heat island (UHI) is the product of anthropogenic processes with urbanization which modifies atmospheric and surface properties and alters the energy balance and thermal environment (Yang, Huang, and Tang 2019). It has been studied as far back as early 19th Century (Howard 1818) and is a wellexplored phenomenon with urban centres of varying sizes (Oke 1973;Katsoulis and Theoharatos 1985;David R Streutker 2003;Krehbiel, Jackson, and Henebry 2016;Yao et al. 2019). Thermal infrared (TIR) data quantified from top of atmosphere radiances in satellite imagery is used to derive land surface temperatures (LST) (Tomlinson et al. 2011). UHI research, now commonly measured globally using LST, has accelerated since 2005 with the majority of regional and seasonal focus being in China and summer daylight hours (D. Zhou et al. 2019). Using the UHI intensity measured between urban temperatures and a referenced rural region, the footprint of the UHI effect is a new index quantifying the spatial extent of the rural area affected by the UHI .
Two commonly used approaches for determining the spatial influence of the UHI are the urban heat island footprint (UHIFP) and the footprint of surface urban heat island (SUHI) model developed by D. R. Streutker (2002;2003). The UHIFP, first introduced to measure the anthropogenic effects that urbanization has upon the warming of vegetation phenology (X. Zhang et al. 2004), has been adapted to determine the distance at which the UHI effects decay towards rural areas (D. Zhou et al. 2015). The footprint model of SUHI effect uses a Gaussian planar surface to fit the UHI effect and has been widely applied to study the UHI in many metropolitan areas (D. R. Streutker 2002;David R Streutker 2003;Martin, Baudouin, and Gachon 2015;Sobrino et al. 2012;Anniballe, Bonafoni, and Pichierri 2014;Anniballe and Bonafoni 2015;Yao et al. 2017;Qiao et al. 2019;D. Zhou et al. 2019). However, the footprints of both approaches have not been compared.
The majority of previous UHI studies using remote sensing are performed with mean annual data from daily products using coarse spatial resolution sensors such as the Moderate Resolution Imaging Spectroradiometer (MODIS) (Yao et al. 2017;Qiao et al. 2019). However, the smoothing effects with high-resolution thermal data are more appropriate for smaller scales despite having higher UHI magnitudes and spatial extents (Anniballe and Bonafoni 2015). Furthermore, the land cover products used in this analysis were also derived from Landsat which will help produce robust residual analyses based on each bio-productive surface.
The application of UHI studies in a cold and urban country such as Canada is increasingly necessary. Climate models indicate greater warming in colder regions, accelerated transformations of nitrogen oxides and volatile organic compounds into ozone, and the population's sensitivity to heat higher than 20˚C (Y. Wang, Berardi, and Akbari 2016). With most of the research focused on major metropolitan cities (Touchaei and Wang 2015;Y. Wang, Berardi, and Akbari 2016;Adamowski and Prokoph 2013), the UHI effect of low-density suburban centers, where the majority of Canadian urban growth occurs (Maoh and Kanaroglou 2007), are often not assessed individually.
This study examines the UHI effect of a low-density suburban centre near the Greater Toronto Area (GTA) from 2000 to 2019 using time series images of Landsat 7 Enhanced Thematic Mapper (ETM+). By comparing the footprint produced from both models, this analysis seeks to address three main questions: 1) what is the extent of thermal disturbance of the UHI effect to each bio-productive land cover in the surrounding rural areas?, 2) are the calculated footprints from UHIFP and SUHI consistent?, and 3) what are primary factors influencing the spatiotemporal pattern that a low-density suburban centre has upon differing land covers with varying distance from the urban perimeter?

Study Site
The location of this study is in Milton, Ontario, Canada. It is situated within the GTA, the most populous area throughout Canada, and is surrounded by a plethora of agricultural and urban land covers, as well as some forested and wetlands. Milton's population has risen 30.5% to 110,128 between 2011 and 2016 making it the most rapidly growing community in Eastern Canada (Hennessey 2017). Urban Milton (UM) and Rural Milton are defined based on population density due to their relatively homogenous neighbourhood structures and their stability for time-series analyses (Gordon and Janzen 2013). Figure 1 displays in boldened polygons the appropriate census tracts representing UM for which the urban surfaces within are representative of the urban centre. Gordon and Janzen (2013) defined urban census tracts as having a population density in excess of 150 people/km2 and exurban having a low gross population density and heavy reliance upon automobiles where over half of the labour force commutes to central cities for employment. This explains 13 of the census tracts; however, a fourteenth tract was selected due to its wide distribution of industrial urban surface components.

Land Cover Data
Version 3.0 of the Southern Ontario Land Resource Information System (SOLRIS) includes land cover maps in 2002, 2006, 2011, and 2016 in the Milton area with a pixel resolution as low as 15 m (Science and Research Branch of the Ministry of Natural Resources and Forestry 2019). Landsat imagery is associated with the nearest chronological land cover product (2004 imagery is analysed with the 2006 land cover product). The SOLRIS dataset employs methodology presented by Lee et al. (1998) to optimize the implementation of the classification techniques. It presents reliable land cover products with overall accuracy of 93% based on an assessment conducted within a polygon encompassing the study area (Sampson 2007

Remote Sensing Data Acquisition
In order to maintain data consistency, only LST data acquired solely from Landsat 7 ETM+ from 2000 to 2019 were used in this study. The ETM+ imagery has 7 main bands (along with an additional panchromatic band) with spatial resolutions ranging from 15 to 60 meters, and a 16-day temporal resolution (National Aeronautics and Space Administration 2020). The LST is extracted using the statistical mono-window algorithm (SMW) conceived by Ermida et al. (2020) from atmospherically corrected surface reflectance Landsat 7 ETM+ products in Google Earth Engine (Malakar et al. 2018). Its accuracy meets the threshold proposed by Guillevic et al. (2018) with overall precision (in RMSE) of 1.0 K determined by Ermida et al. (2020). A 20% local cloud coverage threshold was used to select appropriate LST images for UHI analysis.

Seasonality:
The annual timeframe used in this study was based around the growing season due to the complexity of environmental variables with the local climate. Soil moisture, snow cover, cloud cover, atmospheric moisture, precipitation, and anthropogenic gases are among the elements widely studied to separate the typical extremes of summer and winter measurements (Bernhardt and Carleton 2019;K. Wang et al. 2017;W. Zhou et al. 2014). Late autumn and winter data present far too many obstructions and interference caused by snow coverage. As a result, the annual calendar days included in the research consists of days 100 (April 9/10) to 290 (October 16/17). These days exhibit minimal frost and snow coverage along with increasing solar radiation and vegetated growth (Ministry of Agriculture, Food and Rural Affairs 2020). In order to differentiate the portions of the growing season for this analysis, April and May (referred to as 'Spring') are combined representing the end of frost coverage and the initial planting phase, June through August is the main summer period, and September and October (referred to as 'Autumn') represent the harvest and end of season.

UHIFP
UHIFP uses a single exponential decay model as the following to examine temperature trends from an urban center toward rural areas (D. Zhou et al. 2015). (1) The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII- B3-2021XXIV ISPRS Congress (2021 where ∆T = mean LST differences between the urban area and rural buffers A = maximum temperature difference using the midvalue amongst the three furthest rural buffers S = temperature decay rate D = distance away from the urban centre T0 = asymptotic value To determine the spatial properties of the UHI impact on rural surfaces, buffers are implemented as a factor of distance from the urban perimeter. In order to maintain a similar quantity of pixels within each buffer zone, twelve buffers surrounding the urban land cover area within UM census tracts are generated, each consisting of half the urban surface area (Zhou et al., 2015). Transportation surfaces are masked out to avoid the inclusion of roads in the rural parts of UM census tracts. Instead, a 15 m buffer around all other urban surfaces is used to regulate a centralized urban centre. As this study uses four land cover products addressing Milton's growth over 20 years, a different set of buffers are generated consistently maintaining a total rural area 6 times larger than the urban centre.
Previous research examined annual mean UHIFP based on daily MODIS products (X. Zhang et al. 2004;D. Zhou et al. 2015;C. Meng and Dou 2016). Due to the exchange of a 16 days revisit cycle of ETM+, a 5-years period is implemented with the central year determining the land cover product. With all 16 time periods (period 1 = 2000-2004), each set of results are inputted into the single exponential decay model.
All cloud cover pixels and elevations greater than 50 m from the highest point in UM were removed from the UHIFP model. The footprint determination is based upon the point at which each exponential model reaches 95% of their asymptotic values (X. Zhang et al. 2004). The final product represents the area surrounding UM affected by the local UHI within the proximate exurban buffers.

SUHI Gaussian Fit Model
The planar fit Gaussian surface was developed through the method described by D. R. Streutker (2002;2003) to quantify the UHI as a continuously varying surface. Its purpose is to present a quick and robust technique for analysing and parameterizing the spatial distribution of the SUHI for quantitative inter-city comparisons or single city time scale analyses (Anniballe and Bonafoni 2015). To access the Gaussian surface, the following equation is used: ( where x0, y0 = central location of the model's UHI a0 = magnitude ax, ay = spatial extent φ = orientation All cloud and water pixels must first be masked out from the LST image followed by the temporary removal of urban pixels within UM to produce datasets consisting entirely of rural pixels. After acquiring the mean rural temperature, it is then subtracted from the entire LST image to generate the SUHI Gaussian surface. Its results produce two elliptical thresholds representing the footprint of the SUHI in each image. One ellipse determines the distance from the central location at which the temperature decreases to a fraction of the magnitude, 61% (e -1/2 ) of the maximum value, and the second ellipse is a constant threshold where the SUHI planar fit surface temperature is greater than 1.0 K (Anniballe and Bonafoni 2015).

Residuals Analyses:
The SUHI planar fit model does not produce thresholds that align with the shapes of either UM or the rural buffers produced for the UHIFP model. Alternatively, the Gaussian surface produces residuals which can be analysed with the buffers and a land cover distribution map. Four analyses are made regarding the residuals' relationship with the planar surface; a ratio analysis, a buffer comparison, and a timeframe analysis.
The ratio analysis divided the quantity of positive residuals by the negative and the results with each land cover are displayed in a box-plot separated by months. Its purpose is to quantify the relationship that pixels within each land cover have with the planar fit model depending on whether the residuals are above or below the SUHI surface.
The Gaussian surface's residuals are also analysed using the buffer areas applied to the UHIFP model. The buffer analysis helps determine the SUHI model's spatiotemporal influence and provides a spatial analysis comparable with the UHIFP. Seasonal divisions are implemented based on vegetation growth and solar radiation separating the beginning, middle, and end of the growing season. The mean seasonal distribution of residuals for land covers within each buffer is compared.
The time frames used for each distribution of land covers are incorporated into a time frames analysis. Using Google Earth Engine, mean imagery for each year associated with a SOLRIS product was calculated and compared. General UHI magnitudes, spatial extents, and Gaussian surface residuals are analysed for each part of the growing season. Because of the nature of this analysis, a series of complications occur such as the inability to control local cloud coverage and a bias towards the few days with clear skies. These results are taken with caution and only conclusions regarding general UHI trends are made.

Vegetation Health and Moisture Correlation Analyses with Residuals:
Understanding the variables affecting the residuals of each bio-productive land cover in the SUHI model is achieved through correlation analyses between the residuals and vegetation and moisture indices. The multispectral Landsat 7 bands required from each date are blue, red, near infrared (NIR), and shortwave infrared (SWIR). Amongst the vegetation indices, NDVI is more commonly used due to its cancellation of topographic, cloud cover, atmospheric, and changing sun angle effects with vegetation monitoring (Alademomi et al. 2020). (3) Enhanced Vegetation Index (EVI) acts as an enhancement to it with its atmospheric correction, saturation in areas with high biomass, and soil reflectance reduction abilities (Alademomi et al. 2020). (4) The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B3-2021 XXIV ISPRS Congress (2021 edition) Normalized Difference Moisture Index (NDMI) behaves as an index to reflect the moisture content within the soil and vegetation canopy (Zhu et al. 2019). (5) The comparisons made with vegetation indices use Landsat 7 ETM+ Surface Reflectance Tier 1 products on the same dates as the LST acquisition. A cloud mask was installed, the residuals' pixels were converted to points, and all values from each vegetation/moisture index were extracted to the nearest point. The correlation coefficient is calculated with a range of 1.0 (positive correlation) to -1.0 (negative correlation).

UHIFP
The results from the UHIFP are determined based on the distance at which the exponential model reaches 95% of its asymptotic value. For Milton during the growing season, the mean UHIFP is achieved at 1.4 times the urban area ( Figure 2). As a result, the first three buffers are representative of the area within the local average UHIFP. The strongest outliers were the first two time frames (2000-2004 and 2001-2005) with marginally increased footprints (1.8 times the size of the urban centre). The size of the UHIFP within each set of mean sampled imagery was consistent after the first two time frames. The results experienced in Milton are much smaller than the conclusions drawn by X. Zhang et al. (2004), which were 2.4 times the area of an urban centre in eastern North America. Such results are related to the local conditions influencing the UHIFP analysis, such as population density, city size, proximity to a larger urban centre, and the surrounding vegetation distribution. Daily air temperature data from the nearest station to Milton, in Georgetown, is applied to separate the local UHI effects from regional climate change. Figure 3 presents the mean LST within the urban centre and the adjacent rural area 2 times the urban centre's size (first four buffers). The temperature difference between the urban centre and the first buffer exhibits minimal change while an increasing difference is evident between the first and second buffers after 2005. Following the second buffer, the mean temperature differences amongst buffer regions remain consistently small while most UHI effects have decayed. The data in Figure 3 presents mean LST increasing at a higher rate within the urban centre and the first buffer than the following buffers. As a result, the UHIFP of the very small low-density community, prior to its expansion, is larger relative to the urban centre's area.
The proportion of each land cover within the local UHIFP compared to its distribution throughout the area 6 times the size of UM can help determine the urban centre's impact on each type of surface. The portion of forested surfaces within the footprint was between 11% and 14% of its total. Wetlands and agricultural land covers contain slightly higher proportions within the UHIFP (19% -22% and 24% -27% respectfully). The majority of urban surfaces were contained within the footprint (between 65% and 73%); however, exclusive of the urban centre, only between 18% and 27% of the rural urban surfaces were located inside the footprint.

UHIFP Size Discussion:
The population arrangement and overall size of an urban centre are foundational in determining the footprint of a UHI. Population density can have a major effect on temperature differences within a city as it grows (Ramírez-Aguilar and Lucas Souza 2019). An increased distribution of impervious surfaces, coinciding with rising human activity frequency, contribute to the growing intensity of a UHI (F. Meng and Liu 2013). In comparison with previous studies applied to larger municipalities with denser population distributions, Milton's suburban stature is incapable of reproducing the increasing trends of sensible heat flux experienced in larger cities with less vegetated surfaces.
The proximity to much larger urban centres heavily affects the determination of the local UHIFP. All of the surrounding UHI footprints influence the referenced rural temperatures used for Milton's analysis. For example, the first two time periods (2000-2004 and 2001-2005) determined larger footprints due to Milton's smaller urban size and a reduced external influence from neighbouring UHI footprints. Milton's growth coincides with neighbouring communities' development and their influence on rural temperatures. As a result, an inevitable discrepancy is found between urban and rural temperatures in a suburban setting compared with a larger metropolis, being the primary influence on rural temperatures.

Vegetation Distribution within the UHIFP:
The vegetation distribution surrounding Milton is another major variable on its ecological impact upon the surrounding temperatures. While urban areas surrounding forested or highly vegetated areas produce a greater discrepancy in UHI, partly due to residential energy consumption for cooling during the summer (Imhoff et al. 2010), Milton's rural land covers consist largely of sparse vegetation. The Niagara Escarpment, a large nearby forested area, was omitted from the rural UHIFP due to the impact of elevation. Additionally, being an inland town with limited open water and wetlands surfaces in close proximity may result in lower evapotranspiration comparisons between urban and rural (D. Zhou et al. 2015;Li, Zhang, and Kainz 2012). As a result, urban centres lacking nearby dense vegetation may experience reduced UHIFP sizes.
The diminutive proportion of the study area's dense vegetation located within the UHIFP, as opposed to its distribution in the outer buffers, can have a major influence on the UHI effect on rural LST. Its cooling effects and spatial distribution have the ability to reduce the footprint encroachment on rural vegetation. Additionally, as the UHIFP contains a quarter of the local agricultural surfaces and the majority of urban surfaces, the spatial variability of vegetated land covers, and concentration of impervious surfaces, is a major contributor to a smaller footprint size.

UHIFP Limitations:
The limitations within this study for calculating the UHIFP with Landsat 7 ETM+ consists of an optimistic 16 days revisit cycle, only one time of day where imagery acquisition is possible, and the distribution of snow and ice surfaces in the winter. Previous research (D. Zhou et al. 2015) utilized MODIS products due to their near daily data availability, possibility for comparing daytime and night-time data, and much larger metropolitan sizes. In comparison, Milton's urban area is dwarfed by the far larger North American and Chinese cities. The resulting 5 years periods applied to determine mean pixel values for a UHIFP is unable to reproduce daily samples possible with MODIS analyses; however, this study's results were consistent with every period. Until a satellite is launched with similar spatial resolutions as Landsat and a major reduction in its revisit cycle, applying a 5 years average may presently be the optimal solution as uniform results were found within each 5 years period. With winter conditions, the model remains sufficient provided the urban area is a large enough size.

SUHI Model
The SUHI model was generated from a limited quantity of image acquisition dates due to a coarse temporal resolution and cloud interference. From the available results, the SUHI magnitude observes a rising trend (0.16 K/year). From Figure 4, it can be seen that the area with 1.0 K UHI threshold is increasing at a much more dramatic rate of 3626.7 pixels/year compared with the e -1/2 UHI threshold (a slope of 559.5 pixels/year). The surface thresholds indicated rising trends, and in the case of the 1.0 K temperature ellipse (Figure 4), greatly surpassing the area of UM. The rate at which the thresholds grow appears to align with the rate of expansion of the built-up environment. Based on the distribution and quantity of urban pixels within UM for each land cover product, the timeframe between 2006 and 2011 is when UM experienced its greatest degree of urban development. The declining vegetated surfaces along with increasing impervious surface, anthropogenic gases, and population densities are all highly correlated with higher surface temperature, especially during the urbanization process (Lu et al. 2015;Yu et al. 2019). The development of new built-up environments result in increased volume of ozone, faster pollutant production lofted to higher altitudes, and greater transport of the adverse effects on air quality (D. L. Zhang et al. 2011). Further evidence of these effects are displayed in the Table 2  Time series results incurred from the overall SUHI data, along with the time frames analysis with residuals data, must be taken with speculation due to the insufficient temporal resolution of the Landsat 7 ETM+ sensors. In order to accurately determine the SUHI effects, a spatial resolution similar to Landsat is required along with a 2 days revisit cycle and the overpass time immediately before sunrise (Sobrino et al. 2012). There are currently no spaceborne thermal sensors which satisfy these requirements. The available products from Landsat 7 meet the spatial resolution requirements; however they are heavily unqualified for the temporal obligation. Although trends are noticeable, conclusions based on magnitudes are incapable of being drawn due to insufficiencies with the products. Figures 5 and 6 list the ratio between quantities of positive and negative residuals in each month for each land cover type. The results from the ratio analysis indicate the distinctive impact of urban cool islands within the forested and wetlands surfaces surrounding UM. April and May data describe residual cooling behaviour less intense than in the summer and early autumn.  The most interesting results in this analysis are within the agricultural and urban surfaces. The agricultural residuals, being the dominant land cover throughout the study site, follows the Gaussian planar surface closely with only a marginally greater quantity of negative residuals. Spring and Autumn months exhibit results which better align with the SUHI model as opposed to summer months where a much more definitive distribution of negative residuals exist. The cooling effects climax in August with the average ratio being 0.72, as seen in Figure 5. The urban results contrast vegetated surfaces where positive residuals consistently outnumber negative residuals. As summer progresses, the distribution of positive residuals grows to an August peak.

Residuals Ratio Analysis:
It is well established that vegetated land covers behave as cooling zones within a UHI (Oke 1979;Sun, Wu, and Tan 2012;Li, Zhang, and Kainz 2012;Skelhorn, Levermore, and Lindley 2016;Q. Huang et al. 2019;Alademomi et al. 2020). The ratio analysis results in Figure 5 further exemplify that notion with only the urban land cover category obtaining a greater ratio of positive residuals compared to negative residuals. As summer progresses and thermal emissivity and capacity on impervious surfaces increases, the divergence urban land covers' residuals have compared to the vegetated covers increases (H. . The cooling effects from vegetated land due to evapotranspiration are largely removed when replaced with impervious materials (Oke 1979). The result is a surface with a heavily limited ability for evapotranspiration, and cooling properties, culminating in a warmer surface temperature compared to vegetation.
The crop surfaces in spring months (April and May) along with the autumn months of September and October (to a lesser degree) exhibit a reduction in its cooling abilities within the UHI. With the possibility of frost as late as May and cultivation occurring in September, the agricultural results in Figure 5 display considerable distributions of bare earth at the beginning and end of the growing season. The summer months of July and August, when sparse vegetation growth achieves its peak performance, provide its greatest ability of cooling below the SUHI surface.
Pixels within the other two land cover categories observed a much larger ratio of negative residuals throughout the growing season, especially in summer and autumn. Forests and wetlands, which usually have low LST due to its dense vegetation compared to crop lands which contain more sparse vegetation along with bare soil, are known to behave as cool areas within a UHI (Sun, Wu, and Tan 2012). In the spring, forested land covers deviated the least of any season (-0.27˚C compared to -0.54˚C in the summer) from the UHI temperature mean in Wuhan, China from 2005 to 2015 (Q. ). The densely vegetated land covers surrounding Milton behave in a similar manner with its most effective cooling abilities occurring after spring. Figure 7 shows the density of mean SUHI residual values for each land cover within the buffer regions. The buffer analysis indicates the spatial distribution of residuals from each SUHI image as distance increases from the urban centre perimeter. Applying the results from the UHIFP (rural area 1.4 times the size of the urban area), the first three buffers are within the local footprint of the UHI. Seasonality plays a major role in determining the cooling and warming abilities within a UHI with dense vegetation and urban lands. Summer is when impervious surfaces are warmest and deviate most from the Gaussian surface. Forests and wetlands, in contrast, provide their greatest cooling effect during summer and deviate most while moderately reducing its cooling abilities in autumn. Agriculture exhibits minimal seasonal effects within Milton's UHI.

Residuals Buffer Analysis:
The buffer analysis with both forested and wetlands pixels determined a cooling effect which is much more evident within the first three buffers and in the main summer months.
Throughout the year, the spatial influence which densely vegetated surfaces have on the UHI in Milton are most intense within the surrounding area 1.5 times the size of the urban centre (first 3 buffers), the rural area within the UHIFP. They remain definitive urban cool zones throughout the dataset with the highest mean residual value (-2.02 K) existing within the wetlands category at the fourth buffer zone during Spring. While green spaces within urban environments are known for the ability to improve local air quality, assist in the reduction of energy required for cooling, and optimize the natural ecological system (Ozyavuz 2012), the dense vegetation is most intense within the UHIFP. Agricultural land covers convey similar observations to the more densely vegetated pixels. Their respective cooling effects are limited to the first five buffers within every season and maximized within the first three (the area within the UHIFP). The mean residual values ranged between -0.67 K and -0.26 K within the first 5 (4 for autumn) buffers away from the urban centre before stabilizing with the Gaussian surface. Being the dominant land cover of the region, the distribution is more likely to contribute to the SUHI planar fit and produce residuals closest to zero, especially during the time before planting and after cultivation where its consistency is largely bare soil. The summer months contributed slightly more cooling as a result of its peak vegetated state (Figure 8), especially within the outer buffers.
The reaction within the first buffer zone for the urban land covers' residuals is due to this paper's methods of obtaining an urban perimeter. The urban centre was calculated excluding transportation within UM census tracts and joined with a 15 m buffer to replace the missing roadways and obtain a centralized urban centre. As a result, heavy traffic located on the major 401 highway along with the immediate outer roads are the main contributors to the first buffer's urban distribution.
Urban results provide the highest LST and, aside from Spring, reliably positive residuals. For the remaining buffers, the average urban residuals conform closely to the Gaussian planar fit until the final 6 buffers away from the urban centre where its positivity increases (especially during the summer and autumn months). The area located outside of Milton's UHIFP (buffers 4 -12) is representative of where the ecological footprint that Milton has on its surrounding vegetation phenology has decayed (X. Zhang et al. 2004). The positive residual values that urban pixels experience outside of the UHIFP indicates a localized warming divergence from the Gaussian surface, which has decayed.

Time Frames Analysis:
Mean LST for each time frame as each of the land cover datasets were used to determine general trends. Table 2 presents the SUHI magnitude, model thresholds, and the total quantity of available images within the region for each section of the growing season.
Results with SUHI magnitudes, thresholds, and measures of fit experience growth between the first and last time frame. An apparent rising SUHI magnitude is displayed within the first three time frames during the spring (3.93 K to 6.67 K) and summer (6.55 K to 8.09 K) of It should be noted that the quantity of possible observations is attributed to the number of products available within each tile, thus not every observation is incorporated into the area of the study area.

Vegetation and Moisture Index Correlations:
The results from the correlation analysis between the SUHI model's residuals and the moisture index ( Figure 9) provide the most consistently positive correlation results with all land cover categories. Wetlands and agriculture provided the strongest relationships throughout the growing season along with significantly positive results with urban pixels in the summer. In spring months, the strongest corrections occurs in the wetlands (0.52) while the weakest ones appears in urban land (0.31). Summer has agriculture being the strongest (0.74) with urban surfaces being second strongest (0.56). In autumn, agriculture and wetlands have stronger relationships (0.62 and 0.58 respectfully) and urban surfaces were the weakest again with moderate results (0.45).

Figure 9.
The correlation coefficient each land cover's residuals have with the NDMI product with the same Landsat 7 product. This comparison was determined to have the strongest correlation with the residuals of every land cover.
The results with EVI (Table 3) Table 3. The results from a correlation analysis with each of the indices. The summary statistics are similar to Figure 9 with only the inclusion of mean values from each set of months (M.).
Moisture was consistently the strongest variable affecting the SUHI residuals over every land cover. NDMI is considered as a cooling index useful for regulating thermal configuration with strong greening properties (Zhu et al. 2019). With a consistent moisture distribution throughout the study area, NDMI is capable of retaining a strong correlation with every bioproductive land cover throughout every season. Amongst all indices, its correlation with urban pixels was the strongest due to Milton's low-density residential distribution. Built-up pervious land in the form of lawns, parks, and recreational areas increase abilities for moisture retention.
The land cover which obtained the weakest correlation results with both vegetation indices was urban. This is due to the distribution of impervious surfaces, contrasting vegetation land covers and their abilities to retain moisture and reflect more infrared radiation, which are important in all three indices. The inclusion of built-up pervious land covers (in reference to urban recreation areas, such as golf courses or playing fields (Science and Research Branch of the Ministry of Natural Resources and Forestry 2019)) may be the primary contributor to the positive correlations. Additionally, the low-density structure of the residential neighbourhoods allow for a reduced distribution of impervious surfaces and increased abilities for infrared reflection, sparse vegetation distribution, and moisture retention. As a result, as noted in Table 3, urban surfaces correlate relatively well with NDVI in summer and autumn.
Agricultural surfaces dominated the correlations with NDVI and NDMI with mean correlation coefficients as high as 0.64 and 0.74 within the summer months respectfully. With EVI, the forested and wetlands surfaces obtained a marginally higher correlation, with April and May results deviating the most in Table 3. Considering that LST generally trends in contrast to NDVI and EVI (Alademomi et al. 2020), the stronger correlation with EVI is due to the denser vegetation and reduced soil reflectance while agriculture is distributed more dominantly and performs better with NDVI.
The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B3-2021 XXIV ISPRS Congress (2021 edition)

CONCLUSION
The two common methods for analysing the footprint of a UHI was adapted for a rapidly urbanizing low-density urban area in Canada using LST acquired from Landsat 7 ETM+. Both UHIFP and SUHI models were used for a comparison analysis and to determine the primary factors influencing the results. These techniques allowed for an enhanced comprehension of the effects that the local UHI has upon rural (or exurban) vegetation within differing bio-productive land covers with varying degrees of distance applied.
Milton, a secondary (or suburban) city, was determined to have an average rural UHIFP 1.4 times the size of the urban centre. It is smaller than previous research due to the study site's dwarfed size compared to previous studies, reduced population density, its proximity to larger urban centres at its periphery which influence the rural reference LST, and the rural vegetation distribution. Additionally, the use of 5 years mean imagery for Landsat 7 ETM+ provided an effective and consistent method for calculating the UHIFP with a 16-days temporal resolution.
Despite agricultural surfaces being the most dominantly distributed throughout the region, urban land covers dominated influencing the UHIFP and SUHI models. It mainly followed the SUHI Gaussian surface and exhibited increased warming effects outside of the UHIFP. Agricultural surfaces achieved their greatest cooling effects in the area 2.5 times the urban centre area. Alternatively, densely vegetated land covers maximized their cooling abilities within the rural area 1.5 times the size of the urban centre.
NDMI provided the strongest correlations with every land cover throughout every part of the growing season. As a result, moisture may be more important to the spatial pattern of the UHI and its footprint than vegetation health.
The results from this analysis describe the impact which small scale urban expansion poses for surrounding rural environments. The plethora of proximate agricultural surfaces and distribution of dense vegetation presents Milton's surroundings as an ideal representation for suburban communities across Canada and abroad. Conclusions made here are very useful for city planners, engineers, and geographers in their abilities to predict the impact of low-density urban expansion and its spatiotemporal influence. The strong correlations urban residuals found with NDMI and NDVI defines how the increased vegetation distributions in lowdensity neighbourhoods and urban centres influences the local UHI and rural vegetation.