More than half of the world’s population now lives in urban areas, and trends in rural-to-urban migration are expected to continue through the end of the century. Although cities create efficiencies that drive innovation and economic growth, they also alter the local surface energy balance, resulting in urban temperatures that can differ dramatically from surrounding areas. Here we introduce a global 1-km resolution data set of seasonal and diurnal anomalies in urban surface temperatures relative to their rural surroundings and use satellite-observable parameters in a simple model informed by the surface energy balance to understand the dominant drivers of present urban heating, the heat-related impacts of projected future urbanization, and the potential for policies to mitigate those damages. At present, urban populations live in areas with daytime surface summer temperatures that are 3.21°C (-3.97 - 9.24, 5th-95th percentiles) warmer than surrounding rural areas, such that 1.2 billion people are exposed to average surface summer temperatures in excess of 35°C that might put them at risk of heat-related illness. If design and infrastructure of cities remain unchanged, increased urban heat anomalies will add 0.19°C (-0.01, 0.47) to the daytime summer surface temperatures in urban areas in 2100 -- in addition to warming due to climate change. Such urban heating will increase the number of urban population living under extreme and potentially health-threatening temperatures by approximately 20% compared to current numbers. However we also find a significant potential for mitigation: 82% of all urban areas can optimize vegetation and/or surface albedo and reduce urban daytime summer surface temperatures for the affected population on average by -0.81°C (-2.55, -0.05).

Please see "Drivers and Projections of Global Surface Temperature Anomalies at the Local Scale" by Benz et al for details.

These are all data necessary to repeat the analysis shown in Benz et al. "Drivers and Projections of Global Surface Temperature Anomalies at the Local Scale"

The following describes the workflow in detail. To shorten the workload for researchers wanting to reuse the method, intermediate results are also made available. Specifics are noted in the description below.

Table of Content:
    1) Download all Data
    2) Determine global DT (results made available)
    3) Regression Model
    4) Mitigation Scenarios
    5) Future Scenarios

1) Download all necessary data
Most input data was exported from Google Earth Engine through the following codes:

LST, elevation and nighttime lights
https://code.earthengine.google.com/84c62f0f5640beacaa0915c42fc824dc
Aspect was determined from the dem.tif dataset in ArcGIS Version 10.5 using the Aspect (Spatial Analyst) tool.

Population density, NDVI and BSA
https://code.earthengine.google.com/1d8b929085c19a1f1938e5abcc288294

NDBI
https://code.earthengine.google.com/a4946431b797ccc4234a0a8501b611dd

Precipitation
https://code.earthengine.google.com/7c30c734ea363113b6fd56cc8c09fc06
(to download the different years variables year2 and year need to be adjusted, for different rcps the variable band needs to be adjusted)

Population Scenarios.
Downloaded from SEDAC at
https://sedac.ciesin.columbia.edu/data/set/popdynamics-1-km-downscaled-pop-base-year-projection-ssp-2000-2100-rev01
This data was later resampled to fit the resolution of everything else in ArcGIS Version 10.5, this was done manually, so there is no code to share.

Urban areas:
Global Human Settlement Layer - Urban Center Database was downloaded from https://jeodpp.jrc.ec.europa.eu/ftp/jrc-opendata/GHSL/GHS_STAT_UCDB2015MT_GLOBE_R2019A/
We are working with the V1-0 dataset, by now newer versions are available. The shapefile was converted to a Raster file in ArcGIS Version 10.5 in the same resolution and extend as all other data in order to identify urban pixels.

all data was saved in a 'data' subfolder and is called accordingly in the provided codes.

2) Determining global DT

Temperature Anomalies were determined globally in MATLAB.

First step was to identify all pixels for which the analysis was supposed to be run ( = all pixels with a LST value, and valid background LST). To do so we first identified all pixels with valid LST data, ran one analysis (e.g. winter daytime) for these and later run all other analysis only for pixels that gave us valid results in the first one. For your convenience row and col of these pixels are made available in allPixels.mat.

global Delta T are determined in three steps:
 a) for each season and day and night the code getDT.m is run. Due to memory constrains, this cuts the global dataset into several pieces.
 b) for each season and day and night all pieces are combined into a single image with combineDT.m. To generate this image MATLAB reads tidbits.mat which contains R1 the Geographic Cells Reference of the resulting image (=the Geographic Cells Reference of any input image).
 c) to account for the seasonal differences between Northern and Southern Hemisphere the code shiftSouthernHemisphere is run.
   
Similar analysis was run for NDVI, NDBI and BSA using Delta_Fluxes.m - The default code is running NDVI data, it needs to be adjusted for NDBI and BSA. This was only run for urban pixels based on the urban area’s dataset described above. For your convenience row and col of these pixels are saved in citypixels.mat.

For your convenience, all results are made available as images, a .csv file, and matlab .mat files:

The .tif files are global cover temperature anomalies in °C for all season, day and night (this includes rural pixels)

The .cvs file has one row per urban pixel showing again temperature anomalies (DT) in °C for all seasons day and night, all proxies of the surface energy balance (population, Delta NDBI, Delta BSA, and Delta NDVI), and the absolute BSA and NDVI values used to develop our mitigation scenarios.

DeltaT.mat file contains Delta T for all seasons as well as annual mean for day and night for all urban pixels; Fluxes.mat contains Population, Delta NDBI, Delta BSA and Delta NDVI; scenarios.mat contains absolute NDVI, absolute BSA, Precipitation (historic, 2045, 2065, 2085, and 2099); ssp2.mat, ssp4.mat and ssp5.mat contain populations of the respective SSP (2010, 2020, 2030 ... 2100)
These files are called in the codes described below where it is assumed that they are stored in a subfolder called 'results'.

3) Regression Models

The code regressionModels.m determines the result of the multivariate regression model (compare Table S3 in the manuscript). It reads DeltaT.mat and Fluxes.mat and creates betas.mat. Here all results are stored for day and night for summer (S), winter (W) and annual mean (no indicator). results are saved as an 1000x4 arrays displaying the coefficients of the 1000 samples for each of the 4 fluxes (log(Pop+1), D NDBI, D BSA, D NDVI)

4) Mitigation Scenarios

Mitigation Scenarios based on optimizing Vegetations are determined by Mitigation_NDVI.m, the ones based on optimizing albedo are determined by Mitigation_Albedo.m. Both of these read betas_all.mat and scenarios.mat. They append scenarios.mat with the coefficients (b) of the respective regression and the amount of Vegetation/Albedo that can be improved (Offset). They each export an image for summer day or rather summer night.
The Vegetation code reads a Precipitation (Pr) variable - this contains data of all nine possible timesteps. In reality only data for 2005 (historic), 2045, 2065, 2085, and 2099 (last possible) were prepared.

Both methods are combined in the code Mitigation_combined.m which exports images for summer day and summer night.

5) Future Scenarios

This analysis is performed in 2 steps.
    a) NDVINDI.m reads fluxes.mat to get the coefficient b between NDBI-NDVI and Population. The coefficient is added to betas_all.mat
    b) The code SSP_scenarios.m reads betas_all.mat, ssp_.mat, and (for the best-case scenario) scenarios.mat to determine Future changes in Temperature for different years, SSPs and scenarios. Results are median, 5th percentile, 95th percentile.