Mandatory Exercises

Exercise 1 : First full impact calculation

We begin with a brief example of how to compute risk for winter storms in Switzerland. This example uses pre-computed data from the CLIMADA DATA API (and online database of data demonstrators).

As a first introductory exercise, please create a jupyter notebook impact_climada_YOUR_NAME.ipynb answering the different tasks below.

Tasks for this exercise:

1. Reproduce the code from this exercise and get familiar with it.

2. Reproduce all the plots, add (meaningful) titles to the plots.

3. For each plot, provide a brief explanation of what the plot shows (what are represented on the axis, the legend), and how it could be used. Explore the CLIMADA documentation if needed.

4. Save all the plots to .png files and send them along the notebook.

5. Answer all the Question n: by providing the answer in your impact_climada_YOUR_NAME.ipynb in written form prefixed by Answer n:.

Tip: the figure of a matplotlib.axes (typically named ax) is obtained as fig = ax.get_figure(). A figure can easily be saved to file.

Tip: files paths are most easily manipulated with the Pathlib library.

• Pathlib: documentation

• Pathlib: Article on how to use and why

Tip: Do not forget that you must have the CLIMADA environment climada_env activated when working with the following exercises. If working with jupyter lab select the climada_env kernel, see course material.

Importing the modules

# This code is to hide the warnings (such as deprecation notices).
# You may add this at the beginning of each script if you want to hide warnings.
import warnings

warnings.filterwarnings("ignore")

# This code is to set the resolution of printed graphs.
import matplotlib as mpl

mpl.rcParams["figure.dpi"] = 150

# Load the relevant climada modules
from climada.entity import ImpactFuncSet
from climada.entity.impact_funcs.storm_europe import ImpfStormEurope
from climada.util.api_client import Client

# API client module
client = Client()

The exposure layer

First, we load the exposures from the CLIMADA data API. The DATA API contains pre-computed or demonstrator data for different hazard and exposure layers.

Here we use the data produced by the LitPop module of CLIMADA, which combines nighlight satellite imagery with population census data to estimate the distribution of ‘physical assets’ (mainly buildings).

We retrieve this info for Switzerland.

Tip: For more details about the LitPop module see:

• Asset exposure data for global physical risk assessment

• LitPop documentation

assets_ch = client.get_litpop(country="CHE")

assets_ch is an Exposures object. The data it contains is saved in the form of a GeoDataFrame from GeoPandas (geopandas is a package to represented data as tables and is optimized for geospatial data).

This GeoDataFrame can be accessed via the .gdf attribute of the exposure. It contains information on the value of each data point, the numeric iso-3 code of the country in the region_id columns, and the geospatial coordinates in geometry.

In addition, the column impf_ sets the impact function label for each data point (see more about that below).

assets_ch.gdf

ax = assets_ch.plot_raster()
ax.set_title("Exposures")

Tip: do you obtain ‘deprecation notices? This is because Python packages evolve independently of one another. When someone wants to change something, say change the name exposures to exposure, developers do not simply change the name as this would break the code of everyone using the older naming. First, one makes a ‘deprecation’ warning to inform users that in future iterations, the name will change. Thus, users have time to adapt their code.

The hazard layer

Second, we need to load the hazard, in this case a probabilistic set of winter-storms, again from the CLIMADA DATA API.

Tip: For more details on the generation of winter-storms see:

• Comparing an insurer’s perspective on building damages with modelled damages from pan-European winter windstorm event sets: a case study from Zurich, Switzerland

storm_ch = client.get_hazard("storm_europe", properties={"country_iso3alpha": "CHE"})

Tip: checkout the iso3 country codes system.

• https://en.wikipedia.org/wiki/ISO_3166-1

• https://www.iso.org/iso-3166-country-codes.html

storm_ch is a Hazard object. It contains a probabilistic set of winterstorms. The data is stored in different arrays. Refer to the documentation to know more about its attributes.

# Number of storms
storm_ch.size

# Probabilities in terms of frequency per year
storm_ch.frequency

Question 1: What is the string value used to represent the hazard type of storm_ch, and how do you access it?

storm_ch.plot_intensity(event=-1)

Question 2: What does the ‘-1’ mean for .plot_intensity()? What other values are acceptable? What do they mean?

Tip: See the documentation of CLIMADA for the method plot_intensity. Or use help(storm_ch.plot_intensity)

The vulnerability, or impact function

Third, we need to relate the intensity of the hazard to the impact on the exposures, in others words we need to quantify the vulnerability. This is done with an ImpactFuncSet object.

For a real case study, this would have to be calibrated using existing impact data.

Here, we use the function that was calibrated for building damage in the canton of Zürich based on data from Welker et al. (2021) and that is directly available in CLIMADA.

This impact function is obtained with the method from_welker.

# Vulnerability
impf = ImpfStormEurope.from_welker()
impfset = ImpactFuncSet()
impfset.append(func=impf)
impfset.plot()

Tip: in CLIMADA, different exposure points can have different impact functions. Think for instance that a wind storm probably impacts differently industrial buildings and residential buildings. This is why the ImpactFunc is added to a set of impact functions, i.e. and ImpactFuncSet. In the present case, we considered a single impact function identical for all exposures points.

• For more info, see the documentation for ImpactFuncSet

Now, we have to assign to all exposure points the impact function defined above (with id=1).

For that we have to have a column impf_<Hazard_TYPE> in the GeoDataFrame of the exposure. The hazard type for Winter Storm is WS, and specify the id of the function to use for each exposure point.

# Associate the impact function 1 from WS to all assets in the exposures
assets_ch.gdf.rename(columns={"impf_": "impf_WS"}, inplace=True)
assets_ch.gdf["impf_WS"] = 1

Question 3: What do MDD, PAA and MDR mean?

Tip: See the CLIMADA paper or the CLIMADA documentation on impact functions.

Question 4: is the choice of the impact function appropriate for the chosen exposures?

Impact computation

We are now ready for the impact computation. The only thing to do is to combine the three elements - hazard, exposures, and impact function.

This is done by building an ImpactCalc object from the three elements and calling its .impact() method. The method returns an Impact object, which contains the results.

from climada.engine import ImpactCalc

impact = ImpactCalc(exposures=assets_ch, impfset=impfset, hazard=storm_ch).impact()

impact.plot_raster_eai_exposure();

Different metrics can now be extracted, such as:

# Average annual impact aggregated over all exposures points
impact.aai_agg

Question 5: What other metrics can be obtained from an Impact object?

Tip: See the documentation for the impact class

Tip: The full impact matrix is saved as an attribute and accessible via impact.imp_mat. The impact matrix contains the impact per exposure points (one column per exposure point) and event (one row per event).

Question 3: Can you calculate the average annual impact as a percentage of the total assets value at risk?

Tip: One way to get the total value of the assets at risk, is to use the method assets_ch.affected_total_value() together with the hazard object storm_ch. For a description, type:

help(assets_ch.affected_total_value)

Exercise 2 : Load NetCDF storm data into CLIMADA

In the previous exercise we carried out an end-to-end impact calculation using the hazard, exposure and vulnerability data available in the CLIMADA DATA API.

However, hazard, exposure and vulnerability data most often come from external sources. In this case, some pre-preprocessing is needed to load the data in CLIMADA.

Hazard and Exposure data often come as NetCDF files, since NetCDF is among the most commonly used data format for geographical data.

In this exercise, we will replicate what we have done in Exercise 1, but instead of using the hazard data from the Data API, we will load a NetCDF file that contains winter-storms in Europe.

The NetCDF file is named fp_eraint_1999122606_507_0.nc and can be found on the course website. Download it and save it into the folder /ip_python.

In this exercise, we will use the CLIMADA built-in methods to transform these data into a StormEurope CLIMADA object (a specific type of Hazard object), which will then allow us to carry out an impact calculation in exercise 3.

For this exercise, please save your script as a jupyter notebook with name wisc1_YOUR_NAME.ipynb.

Tasks for this exercise:

• Read the netcdf file for one of the storms using the xarray Python module. Print the xarray.

• Describe in words the data included in the file.

• When there is a question n, provide the answer in your wisc1_YOUR_NAME.ipynb in written form, prefixed by answer n:.

• Find the position (lon/lat) with the largest windspeed and print its value. Where is it (description of the place including the country, region, etc.)?

• For the impact calculation in the next exercise we will only use the windfields over Switzerland. Extract it from the whole dataset.

• Save the extracted data back to a netcdf file called storm_ch.nc.

Tip: NetCDF is a very useful format to store data with different coordinate systems, including information on how the data was made and what the data is (i.e. metadata).

• https://www.unidata.ucar.edu/software/netcdf/

• https://en.wikipedia.org/wiki/NetCDF

import xarray as xr

ds = xr.open_dataset("fp_eraint_1999122606_507_0.nc")

print(ds)

In this file, the maximum wind speeds at each latitude and longitude is saved in the variable max_wind_gust.

ds.max_wind_gust

The data can easily be plotted.

ds.max_wind_gust.plot();

Tip: To get the Get the axes object associated with the plot you can use ax = plt.gca()

Question 1: : What storm was this? When did it happen? Find information (in newspapers archives for example) about the impact of this storm in Europe? in Switzerland?

Tip: below is an example code that can be adapted to to select data for a specific area of interest.

min_lon, max_lon = 10, 30
min_lat, max_lat = 30, 45

mask_lon = (ds.longitude >= min_lon) & (ds.longitude <= max_lon)
mask_lat = (ds.latitude >= min_lat) & (ds.latitude <= max_lat)

ds_selection = ds.where(mask_lon & mask_lat, drop=True)

ds_selection.max_wind_gust.plot();

Tip: an xarray dataset can be saved to netcdf with ds.to_netcdf('filename').

Exercise 3 : Compute storm impacts to Switzerland

In this exercise we will compute impacts to Switzerland using the exposure and vulnerability data retrieved in Exercise 1 and the storm data generated in Exercise 2.

Save your script in a jupyter notebook name as wisc2_YOUR_NAME.ipynb.

Tasks for this exercise:

• Load the storm_ch.nc file directly using the CLIMADA StormEurope.from_footprints(filename) method. (For some data, CLIMADA comes with built-in functions that rely on xarray and allow reading the data into a CLIMADA object)

• Based on the code in Exercise 1, compute the impact from this storm.

• Plot the intensity of the event.

• Plot the impact.

• Extract two important metrics from the impact calculation:

  • The Aggregated Average Annual Impact aai_agg: this is the total impact that is expected each year on average.

  • The Total Value: this is the value of all affected assets.

• Find the name and date of this event. Provide a link to a newspaper article that describes the real damage from this storm. Compare the impact numbers to the ones obtained with CLIMADA.

Tip: the StormEurope class can be imported by doing:

from climada.hazard import StormEurope

Tip: The aai_agg and the tot_val are obtained directly form the impact object.

• The aai_agg is obtained from impact.aai_agg.

• The total value of assets affected from a hazard object can be obtained with the method affected_total_value() of an exposure object.

Exercise 4 : Compute impacts for different countries

In the following exercise, you will repeat the impact calculation done for Switzerland to four other European countries of your choice. This will require adapting the previous code to work for other countries too. For the hazard data, you can either use the storm data in the fp_eraint_1999122606_507_0.nc or just the DATA API.

Save your script to a jupyter notebook compare_countries_YOUR_NAME.ipynb.

Tasks for this exercise:

• Make a function storm_impact(country) that returns the impact to a given country.

• Compute the impact for four European countries of your choice.

• Save the aai_agg for each country. In addition, save the tot_val and compute the relative value of the aai_agg to the tot_val.

• Save these impact values in a pandas dataframe with columns (‘Country name’, ‘Country ISO’, ‘AAI AGG’, ‘Total Value’, ‘Relative AAI AGG’).

• Make a bar plot to compare the aai_agg and make another bar plot to compare the relative values.

• Save all plots to .png files.

• Provide a brief comment of your results.

Tip: a smart way to achieve the tasks above is to employ loops

Example code to construct the required dataframe and make the bar plot:

import numpy as np
import pandas as pd

cntry_list = ["Switzerland", "France"]
iso_list = ["CHE", "FRA"]
aai_agg = np.array([10e5, 16.3e5])
tot_val = np.array([53e7, 17.2e8])
rel_val = aai_agg / tot_val

df = pd.DataFrame(
    {
        "Country name": cntry_list,
        "Country ISO": iso_list,
        "AAI AGG": aai_agg,
        "Total Value": tot_val,
        "Relative AAI AGG": rel_val,
    }
)

df

df["Total Value"]

df.plot(x="Country name", y="AAI AGG", kind="bar");