⚠️This is the first operational version of the handbook, but it is still a work in progress and will be heavily updated during 2024!⚠️

Heat-Wave Risk assessment#

Based on the current climate scenario the occurrence of the heat-wave phenomenon should be more frequent in the future in Europe. The main problems connected with Heat-wave events are the overheating of the urban areas, which lowers the comfort of living or causes health issues[Integrated Assessment of Urban Overheating Impacts on Human Life], drought, and water scarcity. Nowadays, there are a lot of studies and methodologies on how we can mitigate the influence of these events. This toolbox wants to answer simple questions that are more frequently asked by crisis management local authorities, urban planners, or policymakers.

These questions are:

  1. What are the problematic areas? (most overheated areas)

  2. Who or What is exposed?

Import packages#

import zipfile        # working with the zip folders #
import os             # handling the current working directory
import rasterio       # netcdf and raster processing  #
from pathlib import Path # file system paths
import rioxarray as rxr # netcdf and raster processing
import rioxarray
import cdsapi         # API downloading
import numpy as np    # 2-3D array data handling
import pandas as pd   # data handling #
import xarray as xr   # 2-3D array data handling 
import cartopy.crs as ccrs # netcdf data projection 
import cartopy.feature as cf # netcdf data projection 
import cartopy.crs as ccrs # netcdf data projection 
import matplotlib as mpl  #  data plot #
import matplotlib.pyplot as plt #  data plot #
from rasterio.crs import CRS  #  raster dat handling 
from pylab import pcolormesh,show,colorbar,plot,title,legend,subplot,savefig
from xrspatial.classify import reclassify
from matplotlib import pyplot
from rasterio.plot import show
from glob import glob
import geopandas as gpd
import rasterio as rio
import plotly.graph_objects as go
from sklearn.linear_model import LinearRegression
from ipyleaflet import Map, DrawControl, Marker, LayersControl, LegendControl, GeoData, GeoJSON, WidgetControl
import ipywidgets as widgets
import leafmap.leafmap as leafmap
from localtileserver import get_leaflet_tile_layer, TileClient
from itertools import chain
import copy
from datetime import datetime, timedelta
from localtileserver import get_leaflet_tile_layer, TileClient
import json
import random
import requests
import rasterstats
import plotly.express as px

Create a directory structure#

# Define the directory for the flashflood workflow preprocess
workflow_folder = 'Heat_workflow'

# Define directories for data and results within the previously defined workflow directory
data_dir = os.path.join(workflow_folder,'data')
results_dir = os.path.join(workflow_folder,'results')

# Check if the workflow directory exists, if not, create it along with subdirectories for data and results
if not os.path.exists(workflow_folder):
    os.makedirs(workflow_folder)
    os.makedirs(os.path.join(data_dir))
    os.makedirs(os.path.join(results_dir))

Set your CDS API key#

URL = "https://cds.climate.copernicus.eu/api/v2"
KEY = "!!!!!!!!!!Put your KEY!!!!!!!!!!!!!!!!" ### put here your key

Observed 2m air temperature over Europe#

In this part of the notebook, we go to the smaller scale than NUTS2. These data give us information about the observed temperature over Europe. These data also help us find the periods with higher air temperatures for the selection of the Landsat8 images.

Observed temperature Europe 2011-2023#

# This takes time!!!
c = cdsapi.Client(url=URL, key=KEY)

zip_path = os.path.join(data_dir, 'eobs_airtemp_2011_2023.zip')

c.retrieve(
    'insitu-gridded-observations-europe',
    {
        'product_type': 'ensemble_mean',
        'variable': 'maximum_temperature',
        'grid_resolution': '0.1deg',
        'period': '2011_2023',
        'version': '28.0e',
        'format': 'zip',
    },
   f"{data_dir}/eobs_airtemp_EU_2011_2023.zip")
# Define zip file's absolute path
zip_path = os.path.join(data_dir, 'eobs_airtemp_EU_2011_2023.zip')
# Extract from zip file
with zipfile.ZipFile(zip_path, 'r') as zObject:
    zObject.extractall(path=data_dir)

2 Heat islands identification#

Heat islands

Heat islands are urbanized areas that experience higher temperatures than outlying areas. Structures such as buildings, roads, and other infrastructure absorb and re-emit the sun’s heat more than natural landscapes such as forests and water bodies. Urban areas, where these structures are highly concentrated and greenery is limited, become “islands” of higher temperatures relative to outlying areas [Heat islands] For the identification of the heat islands from the historical data we can use these data:

Data, needs to be manually downloaded from provided websites, it requires registration and manual choosing of your area and time period (except Euro-Cordex):

Historical data from satelite sensors:

  • Landsat8 land surface temperature (LST) for the 2016-2020 (15-30m spatial; 8-16 days tremporal). For the identification of the heat island we can use this product, which provide the data about calculated land surface temperature from the Landsat8 imagery. The LST product is available on the RSLAB website, or can be calculated from the L8 imagery bands. Download: ### [RSLAB LST data download] (only for 2015-2021) ## recalculation of the LST from Landsat8 imagery [Landsat8-9 raw images imagery]

  • Sentinel3 LST # [manual] # [Data Download]

Landsat 8, July 2016 land surface temperature 8 days composite#

We can use LST for the summer months (June, July, August) wherewe can expect the overheating of the urban areas. Based on the LST we can easily identified the heat islands (Dark Red areas).

  • for the best estimation of the overheated areas it is best to consult the selection of the days with local Meteorological authorities, they can provide the data about measured air temperature from past.

  • or you can look at the at the observed ait temperature in [CDS Observed temperature 12x12km]

for [Download of the precomputed LST] (for years 2013-2021) for [Calculation of the LST] for [Download of the precomputed LST] (for years 2013- till now)

Landsat8 Land surface temperature#

Before downloading the data create a new folder in the data_dir called LST, where you download data. Then you can proceed with this code.

  • Download all available [all available years for summer months June, July, August, in the south of the EU you also download a data form May and September] LST imagery and save it to the created LST folder. You do not need to unzip a data, because you can do it automatically with the following code (below)

  1. You can download the precomputed values from the RSlab web portal for years 2013-2021 # [LST data download]

    • How to use the RSLAB:

    • (1.) Draw a polygon on the map by clicking (not dragging) on the map to select the vertices of the polygon.

    • (2.) Select your preferred Landsat source. Please check the availability for each Landsat above. It is recommended to use MODIS and NDVI-based emissivity for Urban/Peri-Urban areas and ASTER emissivity for Natural/Isolated areas. Refer to the paper for more details-

    • (3.) Select your preferred emissivity source. Please check the important notes above about emissivity. For more details on this subject please refer to the paper

    • (4.) Click the “Calculate LST” button.

    • (5.) The results will be presented below the map. For each LST result you will find a “show” button that displays the LST as a layer on the map and calculates the mean/min/max LST value for the selected image.

    • (6.) For each LST result you will find a “download” button that downloads the LST as a .tiff image .to download all the images click the “download all” button at the end of the list. The downloaded .tiff is cut based on your polygon.

  2. Or you can calculate the LST directly from the Landsat imagery # [Manual for computation].

  • Save data to your Heat_workflow folder /data_dir/LST.

  • Then you can continue with the following code.

# This code unzips all downloaded LST data 
working_directory = os.path.join(data_dir, 'LST')
# Loop through all files in the directory
for file in os.listdir(working_directory):
    file_path = os.path.join(working_directory, file)
    # Check if the file is a zipfile
    if zipfile.is_zipfile(file_path):
        # Open the zipfile
        with zipfile.ZipFile(file_path, 'r') as zip_ref:
            # Extract all contents of the zipfile into the working directory
            zip_ref.extractall(working_directory)
from datetime import datetime
# This code extraxcts all dates from LST tif.files
# Create a list of file paths to your LST files using glob
L8list = glob(f"{data_dir}/LST/*LST*.tif")
# Initialize an empty list to store the datetime objects
lst_dates = []
# Loop through each filename
for file in L8list:
    # Extract the filename without the directory path
    filename = file.split('/')[-1]   
    # Extract the date and time part from the filename
    if "AllImages_LST" in filename:
        date_time_str = filename.split('.')[1]
        date_str = date_time_str.split('_')[0]
        time_str = date_time_str.split('_')[1]
    else:
        date_str = filename[4:12]  # Extract date part directly
        time_str = filename[13:19]  # Extract time part directly    
    # Combine date and time strings into a single string
    datetime_str = date_str + '_' + time_str    
    # Convert the combined datetime string to a datetime object
    datetime_obj = datetime.strptime(datetime_str, '%Y%m%d_%H%M%S')   
    # Append the datetime object to the lst_dates list
    lst_dates.append(datetime_obj)
# This code will create a raster stack from all downloaded data
# Load a data and crate a raster stack from all maps
L8list = glob(f"{data_dir}/LST/*LST*.tif")
#
with rasterio.open(L8list[0]) as src0:
    meta = src0.meta
#
meta.update(count = len(L8list))
# Save a data to working directory
with rasterio.open(f'{data_dir}/L8_raster_stack.tif', 'w', **meta) as dst:
    for id, layer in enumerate(L8list, start=1):
        with rasterio.open(layer) as src1:
            dst.write_band(id, src1.read(1))
# This code calculates mean values from the LST rasters which are not significantly influenced by clouds
# Open the raster stack
stack_path = f'{data_dir}/L8_raster_stack.tif'
raster_dataset = xr.open_dataset(stack_path)
# Count the number of pixels greater than 0 for each band
num_pixels_gt_zero = (raster_dataset > 0).sum(dim=('x', 'y'))
# Calculate the total number of pixels for each band
total_pixels = raster_dataset.sizes['x'] * raster_dataset.sizes['y']
# Calculate the proportion of pixels greater than 0
prop_gt_zero = num_pixels_gt_zero / total_pixels
# Filter bands where at least 70% of pixels are greater than 0
filtered_bands = raster_dataset.where(prop_gt_zero >= 0.5)
# Calculate the mean values for the filtered bands
mean_values_filtered = filtered_bands.mean(dim=('x', 'y'))
#mean_values_filtered = mean_values_filtered.fillna(0)
mean_values_filtered=mean_values_filtered['band_data']
mean_values_list = mean_values_filtered.values.tolist()
  • This code plots all LST rasters that we downloaded, it is important to check it and relise how many pictures are influenced with the clouds=dark blue.

# This code plots all rasters from LST raster stack 
from datetime import datetime
# Define the path to your raster stack
stack_path = f'{data_dir}/L8_raster_stack.tif'
# Open the raster stack
raster_dataset = xr.open_dataset(stack_path)['band_data']
# Get the number of bands (rasters)
num_bands = raster_dataset.band.size
# Calculate the number of rows needed based on the number of bands and desired number of columns (3 in this case)
num_rows = (num_bands + 2) // 3  # Ceiling division to get the nearest whole number
# Create a Matplotlib figure and axes
fig, axs = plt.subplots(num_rows, 3, figsize=(15, 5*num_rows))
# Flatten the axes array if num_rows = 1
if num_rows == 1:
    axs = axs.reshape(1, -1)
# Initialize variables to store global max and min values
global_max = float('-inf')
global_min = float('inf')
# Loop through each band and find global max and min
for i in range(num_bands):
    band_data = raster_dataset.isel(band=i)
    global_max = max(global_max, band_data.max().values)
    global_min = min(global_min, band_data.min().values)
# Loop through each band again to plot
for i, ax in enumerate(axs.flat):
    if i < num_bands:
        # Plot the raster for the current band with the same color scale
        raster_dataset.isel(band=i).plot(ax=ax, vmin=global_min, vmax=global_max)        
        # Set the date as the title
        ax.set_title(lst_dates[i].strftime("%Y-%m-%d"))
    else:
        # Hide the empty subplot
        ax.axis('off')

plt.tight_layout()
plt.show()
../../../_images/7e492b8f40488d5e642e1e35b37a75cfbe47f7677d7c5213dd8e692a965766a6.png

Reclass of the maximum LST values based on the LST temperature#

  • In this step, you reclassify the LST into 6 categories based on the temperature

  • You can change the treshold values for each category:

    1. No overheat - no data because of the water

    2. Very low < 20°C

    3. Low > 20°C

    4. Medium > 30°C

    5. High > 40°C

    6. Very High > 50°C

  • We need to realize that these values represent the Surface temperature, not the air temperature. The surface temperature reaches higher values than air temperature. Ground surface temperature has become an important component of weather. The temperature of the ground surface can be more than 10 C above air temperature on a sunny day and up to 10 C below air temperature on clear nights, when the surface loses heat by radiation to the cold infinity of outer space. Like most environmental sensors, infrared sensors have become more accurate and lower cost over the past 10 years. [Land surface temperature vs air temperature]

Vegetation respond to high temperatures#

  • If you are more interested in the effect on vegetation

  • For most species, actively growing tissue is damaged by brief exposure to temperatures above 45°C, while prolonged exposure can result in fatal injury. Temperatures between 30-40°C can be termed moderately high temperatures and result in reversible inhibition of metabolism (moderate heat stress). Temperatures above 40°C can be termed very high temperatures as they result in irreversible or prolonged inhibition of metabolism (severe heat stress) [source]. Based on this we suggest the values:

    1. No overheat - no data because of the water

    2. Very low < 20°C

    3. Low > 20°C

    4. Medium > 30°C

    5. High > 40°C

    6. Very High > 45°C

# This code calculates the maximum values from the raster stack of the LST and reclassifies the LST data into 5 groups based on the temperature
# Laod a data calculate maximum values from raster stack
L8 =f'{data_dir}/L8_raster_stack.tif'
L8 = xr.open_dataset(L8)
L8=L8.max(dim='band', skipna=True,  keep_attrs=True)
L8lst2016=L8['band_data']
#L8lst2016=L8lst2016.rio.clip_box(minx=18.67, miny=49.175, maxx=18.8, maxy=49.250)
# Reclasiffy data to the groups by the temperature 
lc_bins=[0,20, 30, 40, 50, 100] ### Effect on the human health
#lc_bins=[0,20, 30, 40, 45, 100] ### Effect on the vegetation
lc_values=[0, 1, 2, 3, 4, 5]
lc_class = reclassify(L8lst2016, bins=lc_bins, new_values=lc_values)
# Plot a data
fig, ax = plt.subplots()
oa = lc_class.plot(ax = ax, cmap='Reds', add_colorbar=False)
cbar = fig.colorbar(oa, ticks=[0, 1, 2, 3, 4, 5])
cbar.ax.set_yticklabels(['No overheat', 'Very Low', 'Low', 'Medium', 'High', 'Very High'], size=10)
ax.set_xlabel('lon [deg_east]')
ax.set_ylabel('lat [deg_north]')
plt.title('Overheated areas in Zilina city', size=16, pad = 10)
Text(0.5, 1.0, 'Overheated areas in Zilina city')
../../../_images/62a0781c526a144397cd29c886c893c789c6776e756a84d1f312d78ebf0a4559.png
lstbbox = L8.rio.bounds()
print("Bounding box:", lstbbox)
Bounding box: (18.67166284348108, 49.16416667916253, 18.846834323884384, 49.26253220277362)
# Save data on disk
lc_class.rio.to_raster(raster_path=f'{results_dir}/LSTZA.tif')

Calculation of the Landsat LST from Landsat 8 - 9 imagery#

Calculation of the Landsat LST from Landsat 8 - 9 imagery data are from european space agency. For the data download you need to choose the level of the data L1 or L2. L1 level contains all 11 bands L2 level contains bands 1-7 and band 10 We need to download the L8-9: [Manual], [Data], [Process]

You do not need to calculate you can download it from [RSLAB] for years 2013-2021#

  • but if you want to most recent data (2022-Now) you need to calculate it:

# Radiance
rad=0.00033420 * b10 + 0.1
# Brightest temperature
bt=1321.0789/np.log(774.8853/rad + 1)-272.15
# Normalized difference vegetation index
ndvi=(b5-b4)/(b5+b4)
ndvi_min=ndvi.min(skipna=True)
ndvi_max=ndvi.max(skipna=True)
# Proportion of the Vegetation 
pv=((ndvi +  ndvi_min)/(ndvi_max+ndvi_min))**2
# emisivity
emi=0.004*pv+0.986
# land surface temperature
LST=(bt+1)+10.8*(bt/14380)*np.log(emi)

Observed maximum 2m air temperature#

# This code crops the Observed temperature data to your selected BBox )lstbbox), and save a data to data_dir
# Read a data 
at = xr.open_mfdataset(f'{data_dir}/tx_*.nc', decode_coords='all') # to read a multiple datasets 
at.rio.write_crs("epsg:4326", inplace=True)
at=at['tx']
atEU=at
at_city=at.rio.clip_box(minx=lstbbox[0], miny=lstbbox[1], maxx=lstbbox[2], maxy=lstbbox[3])
# Delete the data from memory
# select a time window 
at_year=at_city.sel(time=slice("2013-04-01", "2023-06-30")) # Change the dates for specific year
at_year_m = at_year.mean(dim='time', skipna=True, keep_attrs=True)
at_year_m.rio.to_raster(raster_path=f'{data_dir}/at_m.tif')
at_year.rio.to_raster(raster_path=f'{data_dir}/at_y.tif')
# This code plots the map where we need to select the point for the plot of the graphs for computed indices
# First, create a tile server from local raster file
client1 = TileClient(f'{data_dir}/at_m.tif')
# Create ipyleaflet tile layer from that server
t1 = get_leaflet_tile_layer(client1, cmap='Reds', opacity=0.5, nodata=0, name='Heat-wave index')
#m = Map(center=client1.center(), zoom=client1.default_zoom)
m = Map(center=client1.center(), zoom=9)
m.add(t1)
control = LayersControl(position='topright')
m.add_control(control)
# Create lists to store point coordinates
point_lon_list = []
point_lat_list = []
# Create a DrawControl with point drawing enabled
draw_control = DrawControl(marker={'shapeOptions': {'color': '#FF0000'}})
# Add the DrawControl to the map
m.add_control(draw_control)
# Create a text widget to display coordinates
coord_output = widgets.Text(placeholder='Coordinates will appear here', disabled=True)
# Define a function to handle draw events
def handle_draw(self, action, geo_json):
    if action == 'created':
        if geo_json['geometry']['type'] == 'Point':
            # Extract coordinates of the point
            lon, lat = geo_json['geometry']['coordinates']
            # Update text widget with coordinates
            coord_output.value = f'Coordinates: ({lon:.2f}, {lat:.2f})'
            # Append coordinates to lists
            point_lon_list.append(lon)
            point_lat_list.append(lat)
            # Create and add a marker to the map
            marker = Marker(location=(lat, lon))
            m.add_layer(marker)
# Attach the event handler to the DrawControl
draw_control.on_draw(handle_draw)
# Arrange the map and text widget vertically
display(widgets.VBox([m, coord_output]))
at_year = xr.open_dataset(f'{data_dir}/at_y.tif') # to read a multiple datasets 
# This code extracts the data values for the selected pixel
lat = point_lat_list[0] 
lon = point_lon_list[0]
extracted_data = at_year.sel(y=lat, x=lon, method='nearest')
d=extracted_data['band_data']
air_temp = d.values.tolist()
# Sample data
start_date = datetime(2013, 4, 1)
end_date = datetime(2023, 6, 30)
date_list = [start_date + timedelta(days=i) for i in range((end_date - start_date).days + 1)]
values = air_temp  # Replace this with your actual data
# Create figure
fig = go.Figure()
# Add line trace with legend
fig.add_trace(go.Scatter(x=date_list, y=values, mode='lines', name='Observed 2 m air temperature', line=dict(color='blue'), showlegend=True))
fig.add_trace(go.Scatter(x=lst_dates, y=mean_values_list, mode='markers', name='LST mean', marker=dict(color='green', size=8), showlegend=True))
# Add title and axis labels
fig.update_layout(
     title={
        'text': "Observed Maximum air Temperature in the 2m for selected pixel",
        'x':0.5,
        'y':0.95,
        'xanchor': 'center',
        'yanchor': 'top'
    },
    xaxis_title="Date",
    yaxis_title="2 m Air Temperature",
    width=1100,  # Width of the graph
    height=600,
    legend=dict(
        orientation="h",  # Set legend orientation to horizontal
        yanchor="top",  # Anchor legend to the top
        y=1.1,  # Adjust y position
        xanchor="center",  # Anchor legend to the center horizontally
        x=0.5  # Adjust x position
    )
)
# Add slider
fig.update_layout(
    xaxis=dict(
        rangeslider=dict(
            visible=True
        ),
        type="date"
    )
)
# Add vertical lines for each date in date_list2
for date in lst_dates:
    fig.add_shape(
        type="line",
        xref="x",
        yref="paper",
        x0=date,
        y0=0,
        x1=date,
        y1=1,
        line=dict(color="red", width=1, dash="dash"),
        name='Highlighted Date'
    )
# Add custom legend item for vertical lines
fig.add_trace(go.Scatter(
    x=[None],
    y=[None],
    mode='lines',
    line=dict(color='red', width=1, dash='dash'),
    name='day with the LST from Landsat8'
))
# Show plot
fig.show()
  • When you take your Landsat8 LST data and look at the maximum observed 2m air temperature you will see if the LST temperature was reached during a heat wave or just during a hotter day.

  • Vertical red dash line represents the times for which we found Landsat8 LST pictures. You can zoom in / out in the picture with a slider.

  • The green dots represent the mean LST value for the LST image

  • If we cannot detect a green dot on the red line, that means that the LST image has values for less than 50% of the pixels (due to clouds)

  • The LST mean values for the LST images with cloud coverage lower than 50% will be plotted with the code below

import datetime
# Create a dictionary from your data
data = {'Dates': lst_dates, 'Values': mean_values_list}
# Create a DataFrame from the dictionary
import pandas as pd
df = pd.DataFrame(data)
# Plot using Plotly Express
fig = px.scatter(df, x='Dates', y='Values', labels={'Dates': 'Date', 'Values': 'LST Mean Values'}, color='Values')
fig.update_traces(marker=dict(color='green'))  # Set dot color to green
fig.update_traces(mode='markers')  # Display as dots
# Add title and center it
fig.update_layout(title_text='LST mean values (cloud coverage<50%)',
                  title_x=0.5,  # Center the title horizontally
                  title_y=0.9)  # Adjust vertical position if needed
fig.update_xaxes(tickformat='%d-%m-%Y')
fig.show()

3. Exposure and vulnerability of the population#

  • We can use the maps of the distribution of the population across the interested area and look if the most populated areas overlay with the most overheated parts.

  • For the Zilina pilot, we collected data from the Zilina municipality office about the buildings that are usually crowded with huge masses of people, e.g. hospitals, stadiums, main squares, big shopping centers, main roads, and bigger factories… You can apply your data in this part

  • If places like these are overheated, a huge number of people can be negatively influenced by the Heat. With that, the risk probability also rises, and thus, these areas are prioritized for Heat mitigation measures.

Critical infrastructure for Zilina city (overcrowded places)#

# This code loads the Critical infrastructure data, these data were created by the Zilina municipality office
ci=f'{data_dir}/ci_features_ZA.shp'
CI=gpd.read_file(ci)
CI_WGS=CI.to_crs(epsg=4326)

Vulnerable population#

  1. In the link [Source]:

    • Select the High Resolution Population Density Maps + Demographic Estimates

    • Select the Location

    • Select the Formats to GeoTIFF.

  2. Download the maps for the most vulnerable groups of the population elderly 60+ years and children under 5 years. Save the data in the data_dir where you create a folder population and save it there.

  3. When you download all these maps to the Heat-workflow data folder you can use this code for data handling:

    • in the first step we load all the maps of the critical population

    • then we calculate the sum of the vulnerable population from each of the maps

    • we classified the maps into 5 groups (equal intervals)

    • plot it next to a map of overheated areas.

# This code loads all population data and creates a raster stack from them 
poplist = glob( f'{data_dir}/population/*.tif')
#
with rasterio.open(poplist[0]) as src0:
    meta = src0.meta
#
meta.update(count = len(poplist))
#
with rasterio.open(f'{data_dir}/Population_raster_stack.tif', 'w', **meta) as dst:
    for id, layer in enumerate(poplist, start=1):
        with rasterio.open(layer) as src1:
            dst.write_band(id, src1.read(1))
pop=f'{data_dir}/Population_raster_stack.tif'
pop = xr.open_dataset(pop)
pop=pop.sum(dim='band', skipna=True,  keep_attrs=True)
pop=pop['band_data']
pop=pop.rio.clip_box(minx=lstbbox[0], miny=lstbbox[1], maxx=lstbbox[2], maxy=lstbbox[3])
# Calculate the number of bins (classes)
num_bins = 5
# Equal interval classification
min_value = np.nanmin(pop)  # Minimum population value
max_value = np.nanmax(pop)  # Maximum population value
bin_width = (max_value - min_value) / num_bins  # Width of each bin
pop_bins = [min_value + i * bin_width for i in range(num_bins + 1)]  # Define bin boundaries
pop_values=[0, 1, 2, 3, 4, 5]
pop_class = reclassify(pop, bins=pop_bins, new_values=pop_values)
# plot of the Landsat8 LST and population concetration over specified age 
fig, axes=plt.subplots(ncols=2, figsize=(18,6))
lc_class.plot(ax=axes[0], cmap='bwr')
axes[0].set_title('Landsat8 LST')
CI_WGS.plot(ax=axes[0], color='green', alpha=0.4) 
pop_class.plot(ax=axes[1], cmap='Reds')
axes[1].set_title('Population over 65 and under 5')
CI_WGS.plot(ax=axes[1], color='blue', alpha=0.4)
plt.draw()
../../../_images/ea6fc1126a12f75944707f9e7e4e3e9f68889d24fbe5c72d55c4adf61d50141a.png

Save data on the disk

# This code saves the data to results_dir
lc_class.rio.to_raster(raster_path=f'{results_dir}/risk_LST.tif')
pop_class.rio.to_raster(raster_path=f'{results_dir}/risk_pop.tif')

Load of the data

# This code creates a raster stack from risk_LST and risk_pop data, we need this step for the next processing of the data 
S2list = glob( f'{results_dir}/risk_*.tif')
#
with rasterio.open(S2list[0]) as src0:
    meta = src0.meta
#
meta.update(count = len(S2list))
#
with rasterio.open(f'{results_dir}/risk_raster_stack.tif', 'w', **meta) as dst:
    for id, layer in enumerate(S2list, start=1):
        with rasterio.open(layer) as src1:
            dst.write_band(id, src1.read(1))
# Plot a data
fig, ax = plt.subplots()
oa = lc_class.plot(ax = ax, cmap='Reds', add_colorbar=False)
cbar = fig.colorbar(oa, ticks=[0, 1, 2, 3, 4, 5])
cbar.ax.set_yticklabels(['No overheat', 'Very Low', 'Low', 'Medium', 'High', 'Very High'], size=10)
ax.set_xlabel('lon [deg_east]')
ax.set_ylabel('lat [deg_north]')
plt.title('Overheated areas in Zilina city', size=16, pad = 10)
# This code calculates a risk map by multiplying a risk_LST and risk_pop data 
risk=f'{results_dir}/risk_raster_stack.tif'
risk = xr.open_dataset(risk)
risk=risk['band_data']
risk=(risk[0])*(risk[1])
# Plot the risk map
fig, ax = plt.subplots(figsize=(12, 9))
im = risk.plot(ax=ax, cmap='RdYlGn_r', add_colorbar=False)
CI_WGS.plot(ax=ax, color='blue', alpha=0.3)
# Set the title
ax.set_title('Heat-wave risk identification in Selected area')
# Add a colorbar
cbar = fig.colorbar(im, ax=ax, shrink=0.95)
cbar.set_label('Risk Level (based on the risk matrix bellow)')  # Set the label for the colorbar
plt.show()
../../../_images/2dc7b26f73aeaca444dc8c4a5772509ea622fa1afba9a2cabc01fcfda98694f6.png

Save on disk

# This code saves the risk identification map in the results_dir
risk2.rio.to_raster(raster_path=f'{results_dir}/Heatwave_risk_identificationZA.tif')
  • Based on the risk interpretation map (above), we can identify the places most influenced by the Heat-waves (dark red), for better visualization we can load a map with leafmap (below).

  • In this notebook, we can also learn more about how to adapt to the heat-waves (part 6.), or how to interpret this risk map (risk matrix below)

  • The map of the infrastructure (blue) can give us better information about which areas we should prioritize, in case of adaptation or which buildings or squares, etc. are most exposed to heat

image.png [source]

Interactive mapping#

  • To see these maps on the interactive zoom in/out map with the Open Street Base map run the code bellow

# This code loads the Critical infrastructure map and covert it to EPSG:4326 crs, and saves it as Geojson
ci1 = gpd.read_file(f'{data_dir}/ci_features_ZA.shp', encoding='utf-8')
ci = ci1.to_crs(epsg=4326)
crs = ci.crs
#print("CRS of GeoJSON file:", crs)
ci.to_file(f'{data_dir}/ciZA4326.geojson', driver='GeoJSON')
# This code creates a tile client from risk maps
# First, create a tile server from local raster file
riskpop = TileClient(f'{results_dir}/risk_pop.tif')
riskLST = TileClient(f'{results_dir}/risk_LST.tif')
HWRI = TileClient(f'{results_dir}/Heatwave_risk_identificationZA.tif')
# This code creates ipyleaflet tile layer from that server
tpop = get_leaflet_tile_layer(riskpop, cmap='Reds',opacity=0.7, nodata=0, name='Risk population')
tLST = get_leaflet_tile_layer(riskLST, cmap='bwr',opacity=0.7, nodata=0, name='LST')
tHWRI = get_leaflet_tile_layer(HWRI, cmap='RdYlGn_r',opacity=0.7, nodata=0, name='Heat wave risk identification')
# This code loads the geojson Critical infrastructure map
with open(f'{data_dir}/ciZA4326.geojson', encoding='utf-8') as f1:
    geojson1 = json.load(f1)
geo_json1 = GeoJSON(data=geojson1, 
                    style={'opacity': 1, 'dashArray': '9', 'fillOpacity': 0.3, 'weight': 1},
                    name = 'Critical infrastructure for Zilina city')
# This code loads the geojson Critical infrastructure map
with open(f'{data_dir}/ciZA4326.geojson', encoding='utf-8') as f1:
    geojson1 = json.load(f1)
geo_json1 = GeoJSON(data=geojson1, 
                    style={'opacity': 1, 'fillOpacity': 0.3, 'weight': 1},
                    name = 'Critical infrastructure for Zilina city')
# This code plots all loades rasters and vectors on the ipyleaflet map
m = Map(center=riskLST.center(), zoom=riskLST.default_zoom)

control = LayersControl(position='topright')

m.add(tpop)
m.add(tLST)
m.add(tHWRI)

labels = ["Very low", "Low", "Medium", "Very high", "Extreme"]
colors = [(0, 104, 55), (35, 132, 67), (255, 255, 191), (255, 127, 0), (215, 25, 28)]

# Create legend HTML content with custom styles (smaller size)
legend_html = "<div style='position: absolute; bottom: 2px; left: 2px; padding: 10px; " \
              "background-color: #FFF; border-radius: 10px; box-shadow: 2px 2px 2px rgba(0,0,0,0.5); " \
              "width: 75px; height: 205px; '>"

# Append legend items (labels with colored markers)
for label, color in zip(labels, colors):
    color_hex = '#%02x%02x%02x' % color  # Convert RGB tuple to hex color
    legend_html += f"<p style='font-family: Arial, sans-serif; font-size: 14px; '>" \
                   f"<i style='background: {color_hex}; width: 10px; height: 10px; display: inline-block;'></i> {label}</p>"

legend_html += "</div>"


# Create a custom widget with the legend HTML
legend_widget = widgets.HTML(value=legend_html)

legend_control = WidgetControl(widget=legend_widget, position='bottomleft')
m.add_control(legend_control)

m.add(control)
m.add(geo_json1)

m

Conclusion#

  • You can add or remove a map by “click” on the “layer control” in the top right corner, or “Zoom in/out” by “click” on the [+]/[-] in the top left corner

  • We recommend first unclicking all the maps and then displaying them one by one, the transparency of the maps allows you to see which areas on the OpenStreetMap are most exposed to the heat, and in combination with the distribution of the vulnerable population, you can easily identify which areas should be prioritized for the application of the heat mitigation measures.

Vegetation characteristic NDVI#

  1. With NDVI you can detect a change in the plant’s health e.g. before and after heat-wave

  2. You can detect the areas with a lack of greenery

  3. You can see the effect of the greenery on heat-islands

  • In the simplest terms possible, the Normalized Difference Vegetation Index (NDVI) measures the greenness and the density of the vegetation captured in a satellite image. Healthy vegetation has a very characteristic spectral reflectance curve which we can benefit from by calculating the difference between two bands – visible red and near-infrared. NDVI is that difference expressed as a number – ranging from -1 to 1.

  • Normalized Difference Vegetation Index: Change Detection NDVI of a crop or a plant calculated regularly over periods of time can reveal a lot about the changes in their conditions. In other words, we can use NDVI to estimate plant health remotely.

  • A sudden drop in the NDVI values may be a symptom of crop health deterioration. The value drop can also correspond to normal changes, such as the time of harvesting, which is why NDVI should be counter-checked against other available data. Correct NDVI values interpretation can help agronomists raise healthier yields, save money on fertilizers, and take a better care of the environment. [Source]

image.png [source]

Download of the NDVI or Sentinel 2 imagery [download]#

  1. Register to the Copernicus browser website

  2. Open the Copernicus browser [copernicus browser]

  3. Zoom to your area

  4. Select: a Date -> cloud coverage to 10% -> select Sentinel-2 L2A

  5. From LAYERS: display NDVI

  6. Select a download from the panel on the right

  7. Select: Analytical -> Image format: Tiff 32bit -> Image resolution: Custom (10x10m) -> Coordinate system: WGS84 -> Layers: B04, B08 (do not download NDVI directly!!!) -> Download to your data_dir (and unzip it there)

  8. Download NDVI for a multiple years during the vegetation period in your country, to see a vegetation conditions

# This code loads all Sentinel2 bands and calculate the NDVI 
from glob import glob
B4aug = glob(f"{data_dir}/*B04*.TIFF")
B8aug = glob(f"{data_dir}/*B08*.TIFF")
## JULY
for i in range(len(B4aug)):
    B4p=B4aug[i]
    B8p=B8aug[i]
    B4=xr.open_dataset(B4p)
    B4=B4['band_data']
    B8=xr.open_dataset(B8p)
    B8=B8['band_data']
    NDVI=(B8-B4)/(B8+B4)
    #NDVI84=NDVI.rio.reproject_match(L8lst2016)
    target_crs = 'EPSG:4326'  # WGS84
    NDVI = NDVI.rio.reproject(target_crs)
    #ndvi=NDVI.rio.clip_box(minx=18.67, miny=49.175, maxx=18.8, maxy=49.250)
    p=Path(B4aug[i])
    # https://peps.python.org/pep-0428/
    x=str(p.stem)
    parts = x.split('_')
    c=parts[0]
    e=c.replace('-', '_')
    ndvi.rio.to_raster(raster_path=f'{data_dir}/ndvi_'+ e +'.tif') 
# Creates a raster stack from NDVI
# Corrected file path in glob to correctly format the year in the path
L8list = glob(f'{data_dir}/ndvi_*.tif')

if L8list:  # Check if the list is not empty
    with rasterio.open(L8list[0]) as src0:
        meta = src0.meta
        meta.update(count=len(L8list))  # Update metadata to reflect the number of layers

    # Updated file path and name formatting
    with rasterio.open(f'{results_dir}/ndvi_raster_stack.tif', 'w', **meta) as dst:
        for id, layer in enumerate(L8list, start=1):
             with rasterio.open(layer) as src1:
                 dst.write_band(id, src1.read(1))
# This code reclassifies the NDVI values to 4 classes based on the picture above
# Open the NDVI dataset
B4aug = glob(f"{data_dir}/*B04*.TIFF")
ndvi = xr.open_dataset(f'{results_dir}/ndvi_raster_stack.tif')
nd = ndvi['band_data']
# Define the bins and values for NDVI classification
ndvi_bins = [-1, 0, 0.3, 0.6, 1]  # These thresholds were based on the NDVI characteristics
ndvi_values = [0, 1, 2, 3, 4]
# Iterate over the years from 2016 to 2021
for i in range(len(B4aug)):
    p=Path(B4aug[i])
    x=str(p.stem)
    parts = x.split('_')
    c=parts[0]
    e=c.replace('-', '_')
    nd_year = nd[i]   
    # Reclassify the NDVI data
    ndvi_class = reclassify(nd_year, bins=ndvi_bins, new_values=ndvi_values)   
    # Save the reclassified NDVI data to a raster file
    ndvi_class.rio.to_raster(raster_path=f'{results_dir}/ndvi_reclass_plants_{e}.tif')
# This code plots the NDVI on the map
import os
import glob
import rasterio
import folium
import numpy as np
import matplotlib.pyplot as plt
from folium.raster_layers import ImageOverlay
from folium import LayerControl
from branca.colormap import linear
# Initialize the Folium map
m = folium.Map(location=[49, 18], zoom_start=5)
# Assuming results_dir is where your raster files are
# results_dir = '/path/to/your/data'  # Update with your data directory
# Find min and max for NDVI files
ndvi_files = glob.glob(os.path.join(results_dir, '*reclass_plants*'))
ndvi_min, ndvi_max = float('inf'), float('-inf')
for raster_file in ndvi_files:
    with rasterio.open(raster_file) as src:
        data = src.read(1, masked=True)
        ndvi_min = min(ndvi_min, np.nanmin(data))
        ndvi_max = max(ndvi_max, np.nanmax(data))
# Process and add each NDVI raster file to the map
for raster_file in ndvi_files:
    with rasterio.open(raster_file) as src:
        data = src.read(1, masked=True)
        normalized_data = ((data - ndvi_min) / (ndvi_max - ndvi_min) * 255).astype('uint8')
        colormap = plt.get_cmap('Greens')
        mapped_colors = colormap(normalized_data / 255)
        mapped_colors = (mapped_colors[:, :, :3] * 255).astype('uint8')
        bounds = [
            [src.bounds.bottom, src.bounds.left],
            [src.bounds.top, src.bounds.right]
        ]
        img_overlay = ImageOverlay(
            image=mapped_colors,
            bounds=bounds,
            opacity=0.7,
            interactive=True,
            cross_origin=False,
            zindex=1,
            name=os.path.basename(raster_file)
        )
        img_overlay.add_to(m)
# Add a color scale for NDVI
ndvi_color_scale = linear.Greens_09.scale(ndvi_min, ndvi_max)
ndvi_color_scale.caption = 'NDVI Scale'
ndvi_color_scale.add_to(m)
# Add layer control
LayerControl().add_to(m)
# Display the map
m

Heat-effect on plants#

Plants respond to high temperatures For most species, actively growing tissue is damaged by brief exposure to temperatures above 45°C, while prolonged exposure can result in fatal injury. Temperatures between 30-40°C can be termed moderately high temperatures and result in reversible inhibition of metabolism (moderate heat stress). Temperatures above 40°C can be termed very high temperatures as they result in irreversible or prolonged inhibition of metabolism (severe heat stress).[source])