Creating a Heat Risk Index for a Canadian City

This tutorial is a part of a larger collection of resources focused on climate change developed by  Esri Canada 's  Education and Research Group.  These resources use GIS to address climate change through mitigation and adaptation. For an introduction to our climate change resources, check out our  overview story map .

It is also a part of an initiative by Esri to develop and share resources demonstrating the diverse applications of GIS in climate action and resilience using examples from around the globe, including urban heat mapping workflows.

This tutorial walks you through the process of creating a heat risk index for a large Canadian city using the Calculate Composite Index tool in ArcGIS Pro. Municipal open data and satellite imagery from the ArcGIS Living Atlas of the World are used to derive values for different variables that contribute to heat risk. These variables are then combined into a single composite indicator and mapped to help identify areas of high risk across the city. 

 _{This tutorial was originally developed by Esri's Learn Team. You can find the official maintained version at this location: [ https://learn.arcgis.com/en/projects/customize-a-climate-resilience-index/ ]. You can find other tutorials in the tutorial gallery [ https://learn.arcgis.com/en/gallery/ ].}  

• Time required: 75 - 90 minutes.
• Materials required:
  • ArcGIS Pro 3.0+
  • ArcGIS Pro extensions: Spatial Analyst, or 3D Analyst + Image Analyst

• Last updated: October 2024.
• [Tested with: ArcGIS Pro 3.3.1]

By completing this tutorial, you will become comfortable with the following skills:

• Accessing image services from ArcGIS Online and filtering imagery for relevant scenes
• Calculating zonal statistics on raster data to summarize values within feature boundaries
• Summarizing feature attribute data within another layer’s boundaries
• Combining multiple variables into a single composite index 

•  City of Vancouver Open Data Portal 
•  City of Calgary’s Open Data Portal 
•  Metro Vancouver Open Data Portal 
•  ArcGIS Living Atlas of the World 

One of the many negative effects of climate change has been the increase in the frequency and intensity of extreme heat events. In cities, rising temperatures are exacerbated by the urban heat island effect – a phenomenon where urban areas experience significantly warmer temperatures than surrounding rural areas due to structures and materials in the dense urban environment that tend to trap heat, including waste heat generated through human activity. 

In response to a rise in heat-related illnesses and mortality, cities and organizations are developing local heat risk and vulnerability indices to identify areas and populations most at risk during extreme heat events. By quantifying and combining multiple factors or variables that contribute to heat risk into a single index value, cities can evaluate the varying level of risk within an urban area and determine where to prioritize adaptation and mitigation efforts.  

In this tutorial you will create a basic heat risk index for a Canadian city to map and identify high risk areas. The results can be used for local climate action planning, or may serve as a starting point for discussion with stakeholders that may lead to further analysis or refining the index. To build the index, you will incorporate four variables that you have found to be factors that contribute to heat risk – summer land surface temperature, percent tree canopy, population density, and percent impervious surface. 

For this tutorial, you can choose to create a heat risk index (HRI) for either the City of Calgary or the City of Vancouver.  

The two cities have unique and contrasting characteristics that make either of them an interesting choice for analysis. Calgary is a sprawling 825 km ²  inland city in the prairies with a population of over 1.3 million.

In contrast, Vancouver is a coastal city on the Pacific Ocean with a population over 660,000. While not as large in area or population as Calgary, Vancouver is more densely populated – in fact, the most densely populated municipality in Canada, with roughly 5,750 people per square kilometer.

After some thought, you decide that census dissemination areas will be suitable boundaries to use for your index – they are granular enough to show local variation, and they provide meaningful locational context to help target areas where future work can be done to address heat risk. Thus, your first step will be to add the census dissemination area (DA) boundaries for your city, along with census population data for those DAs. You will start by searching the  ArcGIS Living Atlas of the World , a living collection of geospatial data from around the world that is curated by Esri. 

1. Open ArcGIS Pro and sign in, if necessary. 
2. Under New Project, select Map, and name your project Heat_Risk_Index. Click OK. 
3. From the View tab in the Ribbon, select Reset Panes for Mapping from the Reset Panes menu. 
4. In the Catalog pane, select the Portal tab and click on the Living Atlas icon.  
5. In the search bar, type ‘2021 census dissemination’. 
6. Hover over the ‘Dissemination Area (DA) Boundaries – 2021’ dataset, and click the 'Path' URL in the pop-up that appears to view the item detail page. Then, do the same for the ‘Canadian Population and Dwelling Counts 2021’ dataset.
   1. Which of these is considered an authoritative dataset recommended by Esri Canada? (Hint: View the icons above the Description)
7. Right-click the authoritative dataset and click Add to Current Map.
8. Expand the group layer you just added in the Contents pane by clicking on the little arrow next to it.  
9. Zoom in towards your study area (i.e., Vancouver or Calgary). The layers make use of visible scale ranges, with progressively smaller census boundaries being displayed as you zoom in on the map. You only need a layer of the DAs in your city, so you will export a subset of the DisseminationAreas_21 layer. City attribute information is not available in the layer, so you will instead select the relevant DAs by location. 
10. Zoom in to the point where the census subdivision layer is visible (i.e., the layer checkbox is not greyed out). Alternatively, select 1:500,000 from the map scale list dropdown in the bottom left corner of the map view.  
11. Use the Select tool in the Ribbon Map tab to select the subdivision representing the city boundary. 
12. Zoom in further to activate the DAs layer. 
13. Select the DAs layer in the Contents pane, and click Select By Location in the Map tab. 
    1. Choose Have their center in as the Relationship type, and the census subdivisions layer as the Selecting Features. Ensure that Use the selected records is toggled on and click OK. 
    2. If you are working on the City of Calgary, this will have selected all the DAs except for one. To add it to the current selection manually, hold down the shift key and click on the DA while the Select tool is activated.

14. Right-click the DAs layer in the Contents and select Data > Export Features. a. Ensure that Use the selected records is toggled on. b. Name the output feature class ‘DAs_[CityName]’. Click OK.

Finally, you will remove the entire census group layer from the map, and open the attribute table of your exported DAs layer to explore the data.

15. Right-click the group layer in the Contents and select Remove. 16. Right-click the output DAs_[CityName] layer and select Attribute Table.

For a later phase in the analysis, you would like to have more significant figures than what is currently included in the 'Land area in square kilometres, 2021' field, so you will recalculate the area of each DA. You will also need to recalculate population density based on this new area measurement. 

Note that you may choose instead to use more recent population or population estimate data, but your output will differ somewhat from the image shown below.

1. Right-click any field header without an asterisk in the attribute table, and select Calculate Geometry to open the tool.
   1. Instead of overwriting the information in the currently selected existing field, you will create a new field by removing the current field name under Field (Existing or New) and typing LandArea2.
   2. For Property, choose to calculate the Area (geodesic).
   3. Set the units to Square Kilometers and click OK. 
2. Click Calculate at the top of the attribute table to open the Calculate Field tool. 
   1. Create a new field by typing PopDens2021 under Field Name. 
   2. Click anywhere in the dialog outside of Field Name, and set Field type to Float.  
   3. Complete the calculation Expression as below, then click OK:

PopDens2021 = !Population2021! / !LandArea2!

Hint: Double-click 'Population, 2021' in the Fields list. Select the division (/) operator, and double-click 'LandArea2'.

Finally, you will symbolize the DAs layer to visualize the variation in population density across DAs, and change the basemap.

3. Right-click the DAs layer in the Contents. a. Select Symbology to open the Symbology pane.   b. Choose Graduated Colors as the Primary Symbology symbolization method. c. Select PopDens2021 as the Field to symbolize. d. Change the color scheme to the Purples (Continuous) color ramp. Hint: from the dropdown, click Show names. 4. In the Map tab, select Light Gray Canvas from the Basemap dropdown menu.

This creates a choropleth map that distinguishes different classes or value ranges using different shades along a colour ramp. The classification Method and number of Classes can be changed to alter how the data is classified and displayed; you will accept the defaults. 

Which areas of the city appear to have the highest population densities? Which have the lowest? 

5. Click the Save button on the Quick Access toolbar in the top left to save the project. 

ArcGIS Living Atlas includes large collections of  multispectral satellite imagery  taken over multiple years, covering the globe. Preconfigured processing templates can be used in tandem with desired attribute filters to gather insights from multiple mosaiced images across time, including information like NDVI and classified land cover. In this next section, you will use imagery from Landsat 8 and 9 to calculate the median land surface temperature across DAs in recent summers.  

1. Add the Multispectral Landsat imagery layer from the Living Atlas to the map. 
2. Open the Layer Properties dialog for the imagery layer by right-clicking it in the Contents and selecting Properties, or by double-clicking the layer name. 
3. Click the Processing Templates tab, and then choose Band 10 Surface Temperature in Celsius as the template, and click OK. The map view updates, but the imagery displays as grey. To resolve this, you will update the raster symbology. 
4. Open the Symbology pane and select DRA from the Statistics dropdown menu. This dynamically applies a stretch to your imagery, based on the pixel values in the current extent.  
5. Change the Color scheme to a color ramp that helps to better visualize the difference in temperature, such as Yellow-Orange-Red (Continuous).
6. Toggle off the DAs layer in the Contents by clicking off the checkbox to better view the imagery. 

Now, you’ll apply spatial and attribute filters to select only the scenes you need from the millions of images available in the Multispectral Landsat image service.  

1. On the Definition Query tab in the Layer Properties, click the red x to remove the current definition query, then add a new one that matches the image to the right. This filters the dataset to only show imagery collected in the summer months of the last three years that has less than 5% cloud cover. 
   1. Click Apply, and then OK. 
2. Open the imagery layer attribute table. Thousands of records are still included; you will filter the table to show only those scenes that overlap the study area extent. 
3. If necessary, zoom back in to the extent of your study area. You can right-click the DAs layer in the Contents and click Zoom to Layer. 
4. Click the Filter by Extent button at the bottom right of the table.

Scrolling through the records, you note that the Acquisition time is always between 6PM and 7PM, approximately. While you may prefer to know the hottest mid-day temperatures, this still provides a good indicator of heat intensity, as there is less of a reprieve when elevated temperatures are sustained throughout the day and potentially into the night. 

Finally, you will update the mosaic operator, which currently mosaics the images together based on the first raster dataset listed. You will change the operator to instead use the average cell values in any areas of image overlap. 

5. In the Layer Properties Mosaic tab, specify Mean as the Mosaic operator. Click OK.  

You will run the Zonal Statistics as Table tool to calculate the median summer evening temperature within each DA. However, you are only able to run it on a small extent of the service imagery. Thus, you will first limit the processing extent to your local DAs layer. 

1. Click Tools in the Analysis tab to open the Geoprocessing pane.
   1. Search for and open the Zonal Statistics as Table geoprocessing tool. Alternatively, you can search for tools from the Command Search box at the top of the application window.
   2. Under the Environments tab, change the Processing Extent to the Extent of a Layer: choose your local DAs layer.  
   3. On the Parameters tab, point to the DA boundaries as the Input Feature Zones, and OBJECTID as the unique Zone Field. 
   4. For Input Value Raster, select the Multispectral Landsat layer. 
   5. Name the Output Table LST_[CityName].  
   6. Select the appropriate Statistics Type from the dropdown menu. Click Run. 
2. Open the table that is added to the Contents to view the results. You will perform a permanent join to append information in the output table to the DAs feature class so that you can symbolize and view the results. A common field with unique identifier values for each DA is required to match the records in each table – in this case, you will use the OBJECTID field. 
3. Remove the Multispectral Landsat layer from the map. 
4. Open the Join Field tool. 
   1. Choose the local DAs layer as the Input Table, and OBJECTID as the Input Field.  
   2. Point to the LST table as the Join Table, and choose OBJECTID_1 as the Join Field. 
   3. Select the field(s) you want to join to the DAs feature class from the Transfer Fields dropdown.  
   4. Click Validate Join. The Message window should say that the same number of records have been matched as exist in the input and join tables; if not, double check your input parameters.  
   5. Click Run. 
5. Toggle the DAs layer back on and use the MEDIAN field to symbolize it with graduated colors, using a colour ramp of your choice.  
6. Save the project. 

Note that your DAs layer and final overall HRI output may differ from what is shown here, as the Multispectral Landsat dataset is continually updated.

Finally, you will determine the percent tree canopy cover and percent impervious surface within each DA. While you could continue to work directly with satellite imagery for this, you decide that using more detailed, locally acquired or processed data may yield more accurate results. Thus for these two last variables, you will look for datasets provided by the City or larger Region. 

Continue scrolling if you have chosen the City of Calgary, or  click here  to jump to the workflow for the City of Vancouver.

Now that you have derived the values for individual variables, you will combine them to create a heat risk index for the city.  

Indices are calculated by weighting and aggregating relevant input variables following a mathematical method to obtain an overall index value. The variables are first preprocessed to standardize the inputs to a common scale, as they usually have incompatible units, direction, or value ranges and so are not directly comparable. If the data distributions are skewed, transformations may need to be applied for the data to follow a more normal distribution and avoid issues around data compression. 

You will begin by exploring the data distribution for each variable using  Data Engineering , and applying transformations where necessary. You will then use the Calculate Composite Index tool to standardize and weight the data, and calculate index values for each DA in the city. 

For more information on creating indices using this tool, see the  Best Practices guide . 

1. Right-click the DAs layer in the Contents and choose Data Engineering. 
2. Select the four variable fields of interest (hold down the Ctrl key to select multiple), right-click and choose Add to Statistics and Calculate. 

The statistics panel updates with a histogram preview of the data distribution, as well as key summary statistics describing the distribution. If you right-click on the preview, you will see you have options for standardizing, transforming, or reclassifying the data. You will standardize the fields when you create the index. For now, you will identify variables with skewed distributions that you will apply transformations to. Skewed distributions are not ideal for certain aggregation methods, as they can ‘compress’ the data to a small section of the distribution, which can in turn lower the impact of the variable on the overall index value due to low variance. 

For this analysis, you will transform data distributions that are at least moderately skewed, i.e. that have a Skewness less than -0.5 or greater than 0.5. Population density is heavily skewed, so you will start with that one. 

3. Right-click the Skewness value for PopDens2021. a. Select Transform Field. 

Various transformation functions are available to change the shape of the distribution; Box-Cox is the default. 

4. Leave the defaults and click OK. 

The transformed data is added as a new attribute field to the table.

5. Once again, choose to Add to Statistics and Calculate. 

You can see that the data now follows a more normal distribution, and the Skewness value has decreased. You can now use this transformed variable PopDens2021_BOX_COX when calculating the index, rather than the original.

6. Continue transforming all variables with a moderate or high level of skew. 

1. Open the Calculate Composite Index tool. 
   1. Select the DAs layer as the input table.
   2. Name the Output Features [City]_HRI. 
   3. Select your four final fields as the Input Variables. 
   4. All of your variables except for one have a direct relationship with heat risk; i.e., as the variable value increases, so does the risk. Check the Reverse Direction box next to the variable that has an inverse relationship to heat risk. Hint: Which one  mitigates heat ?

The next section provides options for scaling (standardizing) and combining the variables. Alternatively, a preset method can be selected that will automatically set these two parameters. The defaults are set to scale variables to an intuitive 0 to 1 scale using the range of values in the data, and then combine these values using the arithmetic mean. You will use the defaults. 

See the  documentation  for more information on the tool parameters. 

2. Expand the Variable Weights and Output Settings sections.

Here you can assign different weights to your variables, based on their relative importance. Because sources of heat are the primary drivers of risk for impacts on health, you will weight land surface temperature and population to be twice as important as percent canopy and impervious surface, which act to mitigate or exacerbate the heat from these sources. 

3. Assign a new weight of 2 to population density and LST.  a. Name the Output Index [City]_HRI.   b. Set the Minimum and Maximum Index Values to 0 and 10.   c. Under Additional Classified Outputs, check the boxes next to: Equal Interval, Quantile, and Standard Deviation. This creates additional classified layers based on HRI values. d. Leave the default number of output index classes.  e. Click Run. 

A group layer is added to the map, including charts on data relationships and distributions, and several layers visualizing HRI values across the city under different symbolization and classification methods. The layer that is initially toggled on displays the calculated HRI values using a continuous unclassed color scheme which reflects the continuous nature of the index values, with darker areas representing DAs with higher HRI values. 

1. Open the Distribution of [City]_HRI chart.
2. Change the basemap to Imagery and toggle off the DAs_[CityName]_UTM[]N_SummarizeWithin layer.
3. Use the cursor to drag a rectangle to select the two lowest index value bins in the histogram, and inspect the locations of the corresponding DAs in the map (you may want to toggle the group layer on and off). Then select the two highest index value bins, and view these locations in the map. 

Unsurprisingly, the lowest HRI values are associated with heavily treed areas with few to no residents, such as large parks and golf courses. The DAs with the highest HRI values are less predictable, but tend to occur in industrial or downtown areas, or other relatively built-up populated areas with little or no tree canopy. 

4. Open the Relationships of Scaled Variables and [City]_HRI chart. 

You can use the scatterplot matrix to investigate relationships in the data. The Pearson correlation coefficient (r) is shown in green boxes for each pair of variables, and between each variable and the index. A value that is closer to 1 indicates a high level of correlation, which can lead to unintentional weighting with correlated variables becoming more significant. In some cases, one or more correlated variables may need to be removed.  

R-squared values are also calculated for each variable-variable and variable-index pair and can be viewed by hovering over individual scatterplots. A low R-squared value in a variable-index pair indicates that a variable has low variance in values, and ultimately contributes less to the final index. A variable with very low variance may be removed, or an alternative preprocessing method may be applied. 

5. Examine the correlation coefficients between input variables.  

Why is the correlation between percent canopy and percent impervious surface so high? 

When there is high correlation between two factors, depending on the goals of your index it may be worth keeping both. You may choose to retain and sub-index the variables, or remove one of the two if the information is redundant. 

Are there any other input variables with relatively high levels of correlation with each other? If so, why do you think this might be the case? 

6. Examine the correlation coefficients and R-squared values between input variables and the HRI. 

Based on what you’ve seen, would you remove or sub-index any of the variables if you were to reperform the composite index calculation again? Why or why not? 

Finally, you will explore different ways to symbolize the HRI for your city. 

1. Toggle off [City]_HRI in the group layer, and toggle on the next sublayer symbolized by percentile. 

As with the previous one, this layer also uses an unclassed color scheme, but symbolizes DAs based on the percentile of index values, rather than the values themselves. Thus we see a more even distribution of colors in the map. 

2. Toggle the sublayers to view the equal interval classification. 

Here the data are divided into equally sized intervals along the range of values, allowing you to get a very general idea of the data distribution based on the prevalence of each colour.  

3. Toggle the sublayers to view the quantile classification. 

Here the data are divided such that each class contains an equal number of features, making this option ideal for showing intuitive divisions of equal percentages. For instance, with five classes, you can quickly see which DAs have index values within the top 20% of highest HRI values among the DAs. 

4. Toggle the sublayers to view the standard deviation classification. 

The final layer is symbolized using standard deviations of the HRI values distribution to visualize how much the HRI value deviates from the mean.  

The symbolization you ultimately choose depends on the purpose of the map and what you are trying to convey. In your case, you want to symbolize the data using five classes that represent intuitive, equally sized subranges from the range of HRI values. 

5. Toggle on the layer with the appropriate symbolization, based on the information above. a. Apply a graduated colors symbology. b. Use the Orange-Red (5 classes) color ramp. c. Select the correct field to be classified. 

Hint: the field name should reflect the name of the chosen symbolization layer. 6. On the Symbology pane Classes tab, select More > Format all symbols. a. From the Properties tab, change the Outline color to Arctic White.   7. Click the back arrow at the top of the pane. a. In the Classes section, change the numeric code class labels to qualitative ones: b. ‘Very Low’, ‘Low’, ‘Moderate’, ‘High’, ‘Very High’.   8. Select the layer in the Contents. a. In the Ribbon Feature Layer tab, increase Transparency to 35%. 

Can you use the two maps to compare the level of heat risk between DAs in the two cities? Why or why not? If not, how could you change your variable data or index parameters to make them comparable? 

Assuming the index values were comparable, which city do you think would have the higher HRI values in general? Why?  

This tutorial has demonstrated how to use satellite imagery and local data to derive values for individual variables contributing to heat risk. It has also covered how to preprocess and combine the variable data into a composite heat risk index that can be used to identify areas in a city that are at higher risk during extreme heat events. The index inputs and weights could be further refined or revised to incorporate feedback from stakeholders. 

In this tutorial you created a basic heat risk index. However, these indices can be much more complex, with numerous variables being aggregated through several tiers or dimensions of sub-indices. For example, a comprehensive heat vulnerability index could consider multiple variables that feed into larger sub-indices of heat exposure, population sensitivity, and adaptive capacity. 

While being able to identify areas of relative risk during heat events is valuable, it is only an important first step. Once at-risk areas have been identified, cities must determine actionable steps towards mitigation and adaptation – how can heat exposure be reduced? What cooling amenities should be established, and where? These questions may require further analysis, or revision of the original index. 

Interested in learning more about working with imagery? Check out this learning path on  Working with Imagery in ArcGIS . 

 For more tutorial resources on climate change, visit the  Higher Education Resource Finder .