A spatially-explicit approach to simulate urban heat islands in complex urban landscapes

Abstract. Mitigating urban heat islands has become an important objective for many cities experiencing heat waves. Despite notable progress, the spatial relationship between land use/land cover patterns and the distribution of air temperature remains poorly understood. This article presents a reusable computational workflow to simulate the spatial distribution of air temperature in urban areas from their land use/land cover data. The approach employs the InVEST urban cooling model, which estimates the cooling capacity of the urban fabric based on three biophysical mechanisms, i.e., tree shade, evapotranspiration and albedo. An automated procedure is proposed to calibrate the parameters of the model to best fit air temperature observations from monitoring stations. In a case study in Lausanne, Switzerland, spatial estimates of air temperature obtained with the calibrated model show that the urban cooling model outperforms spatial regressions based on satellite data. This represents two major advances in urban heat island modeling. First, unlike in black-box approaches, the calibrated parameters of the urban cooling model can be interpreted in terms of the physical mechanisms that they represent and can therefore help understanding how urban heat islands emerge in a particular context. Second, the urban cooling model requires only land use/land cover and reference temperature data and can therefore be used to evaluate synthetic scenarios such as master plans, urbanization prospects, and climate scenarios. The proposed approach provides valuable insights into the emergence of urban heat islands which can serve to inform urban planning and assist the design of heat mitigation policies.


and the spatial distribution of terrain features such as vegetation indices, without exploring how the observed patterns relate to the biophysical mechanisms that explain canopy-layer UHI. Such a limitation is important when models are used in simulations, for example to examine the effect of urban planning scenarios on air temperatures. As part of the Integrated Valuation of Ecosystem Services and Tradeoffs (InVEST) software, a suite of spatial models to quantify and value the goods and services from nature that sustain and fulfill human life (Sharp et al., 2020), an urban cooling model has been developed following recent 40 research on the effects of surface materials and vegetation cover on UHI (Phelan et al., 2015;Zardo et al., 2017). The aim of the urban cooling model is to simulate the spatial distribution of UHIs based on three key mechanisms, namely the shade provided by trees, the evapotranspiration of urban vegetation and the albedo of the urban surface. In a preliminary application of the model, Hamel et al. (2020) showed its capability to represent the spatial pattern of nighttime air temperature of the 2003 heatwaves in the Île-de-France region. 45 The main objective of this study is to extend such preliminary experiments by proposing a reusable computational workflow to apply the InVEST urban cooling model to predict the spatial distribution of air temperature in a given study area. The validity of the simulated results is optimized by calibrating some key parameters to best fit a set of air temperature measurements from monitoring stations. Additionally, the simulated spatial pattern of air temperature is compared with the one obtained with an alternative approach, namely a spatial regression over features extracted from satellite data.

Study area
Situated at the western end of the Swiss Plateau and on the shores of the Lake Léman, Lausanne is the fourth largest Swiss urban agglomeration with 420757 inhabitants as of January 2019 (Swiss Federal Statistical Office, 2018). As the second most important student and research center in Switzerland (after Zurich), the urban agglomeration of Lausanne has experienced substantial growth during recent decades, which has mostly occurred in the form of suburbanization .
A notable geographic feature of Lausanne is its elevation difference of about 500 m between the lake shore at 372 m.a.s.l. and the northeastern part of the agglomeration (see Figure 1 below). The area is characterized by a continental temperate climate with mean annual temperatures of 10.9 • C and mean annual precipitation of 1100 mm, and a dominating vegetation of mixed broadleaf forest. In line with urban economics and regional sciences, many works rely on administrative boundaries to define the spatial extent of the study. However, the way in which boundaries are constructed overlooks the characteristic scales in which landscape changes and environmental processes unfold, and might thus lead to equivocal results (Liu et al., 2014;Oliveira et al., 2014).
In consideration of such issues, the spatial extent for this study has been determined quantitatively by following the method 65 employed in the Atlas of Urban Expansion (Angel et al., 2012). The core idea is that a pixel is considered part of the spatial extent depending on the proportion of built-up pixels that surround it. In this study, a pixel is considered part of the spatial extent when more than 15 % of the pixels that lay within a 500 m radius are built-up. Additionally, in order to evaluate how temperatures change across the urban-rural gradient, the spatial extent has been extended by a 1000 m buffer. The above procedure has been applied to the rasterized LULC map by means of the Python library Urban footprinter (Bosch, 2020b). The 70 obtained spatial extent, displayed in Figure 1, has a surface of 112.46 km 2 .

Land use/land cover data
The land use/land cover (LULC) maps have been obtained by rasterizing the vector geometries of the official cadastral survey of August 2019 to the 10 m resolution. Such dataset is provided and maintained (i.e., weekly updated) by the cantonal ad-75 ministration of Vaud, and features the whole spatial extent of the canton of Vaud (Association pour le Système d'information du Territoire Vaudois, 2018). The classification distinguishes 25 LULC classes which are relevant to the urban, rural and wild landscapes encountered in Switzerland (Conference des Services Cantonaux du Cadastre, 2011). On the other hand, a 1 m binary tree canopy mask has been derived from the SWISSIMAGE orthomosaic (Federal Office of Topography, 2019), by means of the Python library DetecTree (Bosch, 2020a), which implements the methods proposed by Yang et al. (2009). The 80 tree canopy mask of the spatial extent of the study is shown in Figure 1.

Elevation data
The elevation map for the study area, which is displayed in Figure 1

Satellite data
The satellite dataset consists of the 8 Landsat 8 images in 2018 and 2019 which do not feature clouds over the study area and for days in which the maximum observed air temperature is over 25 • C (see the list of selected image tiles in Appendix A1).
Data from Landsat 7 has been excluded because of the scan line corrector malfunction.

Air temperature data 90
A dataset of consistent air temperature measurements in the study area has been assembled by combining data from 11 stations operated by various governmental and research sources, which are shown in Figure 1. The temporal resolution of the stations ranges from 10 minutes to 30 minutes. Given that the UHI effect in Switzerland reaches its maximal intensity around 9 p.m. (Burgstall, 2019), the remainder of this study evaluates it based on the air temperature observations of that hour.

Model description
The data inputs described above are used to compute the cooling capacity index, which is based on the physical mechanisms that contribute to cooling urban temperatures. More precisely, the cooling capacity index used in InVEST urban cooling model builds upon the indices proposed by Zardo et al. (2017), which are based on shading and evapotranspiration, and extends them 110 by adding a factor to account for the albedo. For each pixel i of the LULC raster map, the cooling capacity index is computed as in: where S i , AL i and ET I i respectively represent the tree shading, albedo and evapotranspiration values of pixel i as defined in the biophysical table, and w S , w AL and w ET represent the weights attributed to each component respectively. The values of 115 S i and AL i are retrieved from the biophysical table according to the LULC class of the pixel i (see Appendix A3). The tree shading coefficient corresponds to the average proportion of tree cover in all the 10 m pixels of each LULC class, which has been computed by coupling the 1 m binary tree canopy mask with the rasterized 10 m LULC map. The albedo coefficients are based on the local climate zone classification by Stewart and Oke (2012).
The evapotranspiration index ET I is computed as a normalized value of the potential evapotranspiration as in: where K c is the evapotranspiration coefficient, ET ref is the reference evapotranspiration raster for the period and area of interest and ET max is the maximum evapotranspiration value observed in the area of interest.
In line with the studies of Nistor et al. (Nistor and Porumb, 2015;Nistor et al., 2016;Nistor, 2016), the evapotranspiration coefficients are attributed to each LULC class by distinguishing four cases, namely the crop coefficient for single crops for 125 vegetation LULC classes, the water evaporation coefficient for surface water, the rock and soil evaporation coefficient for bare soils and rocks, and evaporation coefficients for artificial LULC classes (e.g., urban areas). The evapotranspiration coefficients attributed to the LULC classes of the Swiss cadastral survey are listed in subsection A3.
Following the recommendations of Allen et al. (1998), the daily evapotranspiration ET ref (in mm/day) has been estimated for each pixel using the Hargreaves equation (Hargreaves and Samani, 1985) as in: where T avg , T max and T min respectively correspond to the average, maximum and minimum T air (in • C) of each day and R a is the extraterrestrial radiation (in mm/day), which is in turn estimated for the latitude of Lausanne (i.e., 46.519833 • ) for each date following the methods of (Allen et al., 1998, Equation 21). The temperature values of each day have been extracted from the inventory of gridded datasets provided by the Federal Office of Meteorology and Climatology (MeteoSwiss), which 135 feature the minimum, average and maximum daily T air for the extent of the whole country at a resolution of 1 km. Such a dataset is obtained by interpolating 100 T air stations across Switzerland (including the MeteoSwiss Pully station of Figure 1) based on non-linear thermal profiles of major basins and non-Euclidean distance weighting that accounts for terrain effects (Frei, 2014).
In order to account for the cooling effect of large green spaces, the computed cooling capacity index of pixels that are part 140 of large green areas (> 2 ha) is adjusted as in: where g i is 1 when the pixel i is a green area and 0 otherwise (as defined in the biophyisical table), d(i, j) is the distance between pixels i and j, d cool is a parameter that defines the distance over which a green space has a cooling effect, and Ω i is the set of pixels whose distance to i is lower than d cool .

145
Then, a heat mitigation index is computed as: In order to simulate the spatial distribution of T air , the model requires two additional inputs. The first is the rural reference temperature T ref , where the UHI effect is not observed, e.g., in the rural surroundings of the city. The second is the magnitude of the urban heat island effect U HI max , namely the difference between the rural reference temperature and the maximum T air 150 observed in the city center. The two parameters are combined with HM i to compute the T i for each pixel i of the study area as in: Finally, the T air values of each pixel T no mix i are spatially averaged using a Gaussian function with a kernel radius r defined by the user.

Calibration and evaluation of the model
To compare the InVEST urban cooling model with the spatial regression based on satellite features, the urban cooling model is used to simulate the spatial distribution of T air for the same 8 dates used to train the spatial regression model, i.e., the dates of the selected Landsat images. It is implicitly assumed that no significative LULC changes have occurred throughout study period (i.e., from May 2018 to August 2019), and therefore all simulations depart from the same LULC raster, i.e., the 160 rasterized cadastral survey of the canton of Vaud as described above. Given the rugged terrain of the study area, the T ref has been set as the minimum average T air observed among the monitoring stations, while U HI max has been set as the difference between the maximum average T air observed among the monitoring stations and T ref . The values of T ref and U HI max for the 8 days considered in this study are displayed in Figure A1.
Although the documentation of the InVEST urban cooling model (Sharp et al., 2020) provides some suggested values 165 for several parameters of the model, their suitability depends strongly on the local geographic conditions of the study area.
Therefore, calibration of the parameters is required in order to better understand how the physical mechanisms beyond the emergence of UHIs take place in the context of Lausanne. Following the manual calibration approach drafted by Hamel et al. (2020), the target parameters are the weights attributed to the tree shading w S , albedo w A and evapotranspiration w ET , the distance over which green spaces have a cooling effect d cool and the T air mixing radius r. As an additional contribution, this 170 article implements an automated calibrated procedure based on simulated annealing optimization (Kirkpatrick et al., 1983) that aims at the minimization of the R 2 between the T air values observed in the monitoring stations and those predicted by the model 1 . The parameter values suggested in the documentation of the model are set as the initial state of the simulation annealing procedure, which corresponds to a T air mixing radius of r = 500 m, a green area cooling distance of d cool = 100 m, and weights attributed to tree shading, albedo and evapotranspiration of w S = 0.6, w A = 0.2 and w ET = 0.2 respectively. The 175 number of calibration iterations is set to 100. For both experiments, the three evaluation metrics are reported as their average over 10 runs.

Spatial regression of air temperature based on satellite data 185
The spatial regression to predict T air from features derived from satellite data is performed over a raster dataset on a per-pixel basis. A regression model is then trained to fit the observed T air measurements by minimizing the error at the pixels that correspond to the locations of the monitoring stations.
The regression operates in each pixel with the T air as the target variable, and the elevation, the LST and the normalized difference water index (NDWI) (Gao, 1996) as independent variables. Additionally, to account for the influence of temperature 190 and moisture surface conditions of each pixel, the LST and NDWI are spatially averaged over a series of circular neighborhoods with radii 200, 400, 600 and 800 m, thus reckoning 8 supplementary features. Based on previous research on the sensitivity of the landscape patterns-UHI relationships to the spatial resolution (Weng et al., 2004;Song et al., 2014), the target resolution has been set to 200 m.

195
The estimation of LST from Landsat 8 images follows the methods of Avdan and Jovanovska (2016). On the one hand, the data from the near-infrared (NIR) and red bands of Landsat 8 (i.e., bands 4 and 5 respectively) to compute the normalized difference vegetation index (NDVI), which is then used to estimate the ground emissivity ε λ . On the other hand, following the Landsat 8 data users handbook (Zanter, 2015), the data from the thermal band of Landsat 8 (i.e., band 10) is first converted to top of atmosphere spectral radiance L λ , from which brightness temperature BT is estimated (in • C) as in: where the K 1 and K 2 are band-specific thermal conversion constants embedded in the Landsat image metadata. Finally, the ground emissivity ε λ and the brightness temperature BT are used to compute the LST by inversion of Planck's Law as in: where ρ = 1.438 · 10 −2 m K is a constant computed as a product of Boltzmann constant and Planck's constants divided by the velocity of the light, and λ = 10.895 · 10 −9 is the average of the limiting wavelengths of the thermal band.
The NDWI is computed from the green and near infrared (NIR) bands of Landsat 8 (i.e., bands 3 and 5 respectively) as in:

Model selection and evaluation
Based on the work of Ho et al. (2014), three regression models have been considered, namely a multiple linear regression, 210 support-vector machine (SVM) and random forest. The accuracy of each regression model is assessed by means of a k-fold cross-validation procedure, where the regression samples are first shuffled and partitioned into 3 folds. Then, for each fold k, a regression model is trained using the other 2 folds and validated using the samples of such k fold. Finally, the model that shows the best validation score (i.e., the R 2 averaged over 10 repetitions of the k-fold procedure) is selected. Additionally, the MAE and RMSE are computed in order to evaluate the deviations between the observed T air and the predictions of each model.

215
On the other hand, the importance of each feature is evaluated by computing its permutation importance (Breiman, 2001), namely the average decrease of the regression accuracy when such feature is randomly shuffled. The training of the regression models, cross-validation and permutation feature importance described above have been implemented by means of the Scikitlearn library (Pedregosa et al., 2011).

Spatial regression of air temperature based on satellite data
When including all the samples, the R 2 for the linear regression, SVM and random forest are respectively 0.832, 0.014 and 0.960, with MAE of 1.198, 2.671 and 0.580 • C, and RMSE of 1.508, 3.652 and 0.738 • C respectively. The coefficients suggest that SVM is not well suited for such a regression in this study area, whereas the linear regression and random forest models obtain a very strong fit -the latter achieving the best performance. Nevertheless, the average cross-validation scores 225 suggest that the linear regression (average score R 2 = 0.733) is more robust to missing data and also less likely to over-fit the observations than the random forest regressor (average score R 2 = 0.658). F-test suggest that the significant variables for the linear regression are the NDWI when spatially averaged over a 800m, 600m and 400m radius (in decreasing order of significance). The following most significant variable is the NDWI spatially averaged over a 200m radius (p = 0.071) and without spatial averaging (p = 0.231), and the LST spatially averaged over a 400m radius (p = 0.277). With a p = 0.420, the does not appear to be significant in this particular regression. The low significance obtained for the LST features in this study might be attributable to the large time lag between the acquisition time of the Landsat images 235 (which ranges from 11:15 to 11:23 CET) and the time of the T air measurements (i.e. 21:00 CET).
The relationship between the predicted and the observed values is displayed in Figure B1. The MAE and RMSE of 1.198 and 1.508 • C respectively demonstrate a stronger fit than the 1.82 and 2.31 • C obtained in the study of Ho et al. (2014) in Vancouver.
The two plots of Figure 2 show the relationship between the elevation and T obs of each sample and the regression errors. While there is no discernable relationship regarding the elevation of the samples (i.e., the elevation of the monitoring stations), the 240 regression errors seem to be negatively correlated with T obs . This pattern, which was also noted by Ho et al. (2014), indicates that high temperature samples are systematically underestimated by the regression model whereas low temperature samples are consistently overestimated.
The series of predicted T air maps for the 8 dates as well as the prediction errors at the locations of the monitoring stations are displayed in Figure 3. While the range of temperatures exhibits important differences throughout the dates, the spatial 245 distribution of T air is seemingly consistent. The highest temperatures persistently occur in the most urbanized areas, whereas the lowest temperatures are take place in the higher elevations located east and north-east of the map. Finally, there seems to be no discernable pattern in space nor time regarding the prediction errors at the monitoring stations.

Simulation with the InVEST urban cooling model
The parameters of the model that result in the best fit of the station measurements are a T air mixing radius of r = 236.02 m, 250 a green area cooling distance of d cool = 89.21 m, and the weights attributed to tree shading, albedo and evapotranspiration of  low temperatures -the former being more prominent in this case, as noticeable from the asymmetry of the vertical axis in Figure 4.
The simulated T air maps for the 8 dates and the prediction errors at the monitoring stations are shown in Figure 5. As in the 265 spatial regression, the temperature ranges show important differences across dates yet the same spatial pattern of T air persists.
The simulated distribution of T air shows its highest values in the center of Lausanne and along the most urbanized (and hence less forested) zones along the main transportation axes, whereas the lowest temperatures are found in the forested areas located in the eastern and western extremes of the upper-half of the study area.

Model comparison 270
A comparison of the maps predicted by the spatial regression and the urban cooling model is displayed in Figure 6. In line with the temporal consistency of the spatial patterns predicted by the two approaches respectively, the comparison maps also show a spatial distribution of T air that persists throughout the dates. Such spatial pattern is strongly reminiscent of the elevation maps (see Figure 1 above) and reflects the fact that the elevation is explicitly considered in the spatial regression but not in the urban cooling model. The overall distribution of the T air pixel differences between the two approaches follows a normal 275 distribution that ranges from -9.620 • C to 11.929 • C (respectively reflecting lower and higher temperatures predicted in the spatial regression), which is considerably large range when compared to the small overall MAE and RMSE of both approaches.
Nonetheless, the way in which the histogram is centered around 0 • C suggests that the differences between the two approaches follow no particular correlation other than the spatial regression predicting more extreme T air values, which is not surprising considering that the range of T air is systematically bounded in the urban cooling model by the T ref and U HI max parameters.

Discussion
The results obtained of this study suggest that both the spatial regression based on satellite data and the InVEST urban cooling model are capable of predicting the spatial distribution of air temperature with a large degree of statistical determination.
Furthermore, the fact that a similar spatial pattern is predicted by both models suggests that the biophyisical mechanisms embedded in the urban cooling model are well represented. If that is the case, the urban cooling model presents two central 285 advantages with respect to the spatial regression.
Firstly, unlike regressions and black-box approaches, the fact that the biophysical mechanisms that drive the emergence of UHIs are represented explicitly allows for a physical interpretation of the parameters of the model. For example, in a comparative study of the relationship between the LST and the spatial configuration of trees in Baltimore and Sacramento, Zhou et al. (2017) suggest that the distinctive results observed in each city might be related to how the shading of trees and 290 evapotranspiration contribute differently to urban cooling in the climatic context of each city. More precisely, they suggest that in the dry climate of Sacramento, large patches of trees ameliorate the efficiency of the evapotranspiration, whereas with 13 https://doi.org/10.5194/gmd-2020-174 Preprint. Discussion started: 13 July 2020 c Author(s) 2020. CC BY 4.0 License. the humid climate of Baltimore, the gains from the tree shading are likely more important. The urban cooling model provides a suitable mean to quantitatively address such matters, i.e., by calibrating the model in the two cities, we can explore the weights obtained for each factor support such hypothesis. In the case study of Lausanne reported above, the weight attributed 295 to the tree shading w S = 0.59 is higher than the one attributed to the evapotranspiration w ET = 0.17. This is consistent with the local climatic conditions being more similar in Lausanne and Baltimore than in Sacramento, yet the weights obtained in this study might be partly determined by the initial solution provided. Nonetheless, to further understand this issue, validation and calibration of the InVEST urban cooling model in a broader variety of cities is required. Overall, the way in which the calibrated parameters differ from the recommendations in the documentation of the model are in consonance with the particular 300 characteristics of Lausanne. More precisely, the smaller mixing radius and cooling distances are consistent with the uneven relief of the study area.
The second major advantage of the urban cooling model is that once the model is calibrated for a given city, it can be used to evaluate synthetic scenarios such as those stemming from master plans, urbanization prospects or the like, and to spatially fact one of the central purposes of the InVEST suite of models (Tallis and Polasky, 2009). Statistical models like the spatial regression are not well suited to such a purpose since they rely on features such as the LST that are hard to obtain other than empirically.
The approach proposed in this article is nevertheless subject to some limitations that merit thoughtful consideration. On the one hand, as acknowledged in its user guide (Sharp et al., 2020), the design of the InVEST urban cooling model presents a 310 number of limitations, the most relevant to this study being the simplified and homogeneous way in which the air is mixed and the cooling effects of large green spaces. In complex terrains such as the Lausanne agglomeration, models with uniform weighting of space show considerable deviations from the observed distribution of air temperature (Frei, 2014). On the other hand, the relationship between the calibration parameters and the resulting R 2 is likely to define a complex optimization landscape with multiple local optima. As a metaheuristic that strongly depends on random decisions, the simulated annealing 315 procedure is susceptible to convergence to local optima, arbitrarily leading to different solutions in each run. A sensitivity analysis of the parameters of the urban cooling model as undertaken preliminarily by Hamel et al. (2020) for the Île-de-France region could serve as a basis to improve the simulated annealing procedure by careful design of appropriate neighborhood search and annealing schedule. Finally, the approach of the present study based on observations at the moment of maximal UHI intensity (i.e., 9 p.m. in Switzerland), however the factors that influence UHIs are likely to operate differently across the 320 diurnal UHI cycle. In fact, several studies point to distinct relationships between the spatial patterns of vegetation and daytime and nighttime UHIs (Anniballe et al., 2014;Sheng et al., 2017;Shiflett et al., 2017;Hamel et al., 2020). Considering the nature of the implications of UHI, e.g., energy consumption, work productivity or human health (Koppe et al., 2004;Santamouris et al., 2015;Zander et al., 2015), a sound understanding of the full diurnal UHI cycle becomes crucial towards the design of robust solutions.

325
Nevertheless, the limitations on how the urban cooling model represents the spatial air mixing and the cooling effects of green spaces seem hard to overcome with the current spatial sparsity of monitoring stations. Such major shortcoming, which contrasts with the growing availability of high-resolution LST datasets, is one of the main reasons why most of the UHI studies have focused on the latter (Jin and Dickinson, 2010;Zhou et al., 2019). As illustrated in this article, spatial regressions based on remote sensing features such as LST and NDWI do not necessarily replicate the air temperature measurements better than 330 biophysical models such as the InVEST urban cooling model. Therefore, improving the spatial density of the monitoring network becomes an imperative for further enlightening of the UHIs phenomena.

Conclusions
The present article presents a spatially-explicit approach to simulate UHIs with the InVEST urban cooling model, which is based on three biophysical mechanisms, namely tree shade, evapotranspiration and albedo. The proposed approach shows 335 how LULC and air temperature data can be combined to calibrate the parameters of the model to best fit measurements from monitoring stations by means of an automated procedure. The simulations performed for the urban agglomeration of Lausanne show that the InVEST urban cooling model can outperform spatial regressions based on satellite-derived features such as LST, NDWI and elevation. The way in which both approaches consistently predict the highest temperatures in the most urbanized parts of the agglomeration suggests that the enhancement of green infrastructure can be an effective heat mitigation strategy, 340 yet further exploration in other climatic contexts is required to fully understand this issue. To that end, the reusability of the computational workflow paves the way for further application of the urban cooling model to a broad variety of cities, which can serve to improve the understanding of the UHI phenomena and support the design of heat mitigation strategies.
Code availability. The code materials used in this article are available in a GitHub repository at https://github.com/martibosch/lausanne-heat-islands Appendix A: Data

355
The crop and water coefficients are based on Allen et al. (1998), while rock, soil and urban coefficients are derived from the results of Grimmond and Oke (1999) in the city of Chicago. Given that the evapotranspiration of the vegetation and crops is subject to seasonal changes in temperate zones such as Switzerland (Allen et al., 1998), the values that correspond to the midseason estimation (June to August) in Nistor (2016). The albedo values are based on the work of Stewart and Oke (2012). The shade column, which represents the proportion of tree cover of each LULC class, has been computed with a high resolution 360 tree canopy map of Lausanne and is therefore specific to the study area. Rows with a hyphen sign -in the shade column denote that the corresponding LULC class is not present in the study area. The source CSV file used in the computational workflow is available at https://github.com/martibosch/lausanne-heat-islands/blob/master/data/interim/biophysical-table-shade.csv.

B1 Spatial regression
The code of the spatial regression of air temperature from satellite data is available as a Jupyter Notebook (IPYNB) at https: //github.com/martibosch/lausanne-heat-islands/blob/master/notebooks/spatial-regression.ipynb.

B2 Simulation with the InVEST urban cooling model
The code of the spatial simulation of air temperature with the InVEST urban cooling model is available as a Jupyter Notebook (IPYNB) at https://github.com/martibosch/lausanne-heat-islands/blob/master/notebooks/invest-urban-cooling-model.ipynb

B3 Comparison
The code used for the comparison of the spatial regression and simulation of air temperature is available as a Jupyter Notebook 375 (IPYNB) at https://github.com/martibosch/lausanne-heat-islands/blob/master/notebooks/comparison.ipynb Author contributions. MB and ML designed the study and conducted the analysis with technical inputs of PH and RR and overall supervision of JC and SJ. MB developed the code and wrote the manuscript. All authors provided comments and contributed to the final version of the manuscript.
Competing interests. The authors declare no competing interests