Hazard: winter windstorms / extratropical cyclones in Europe

Or: The StormEurope hazard subclass of CLIMADA

Auth: Jan Hartman & Thomas Röösli
Date: 2018-04-26 & 2020-03-03

This notebook will give a quick tour of the capabilities of the StormEurope hazard class. This includes functionalities to apply probabilistic alterations to historical storms.

%matplotlib inline
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = [15, 10]

Reading Data

StormEurope was written under the presumption that you’d start out with WISC storm footprint data in netCDF format. This notebook works with a demo dataset. If you would like to work with the real data: (1) Please follow the link and download the file C3S_WISC_FOOTPRINT_NETCDF_0100.tgz from the Copernicus Windstorm Information Service, (2) unzip it (3) uncomment the last two lines in the following codeblock and (4) adjust the variable “WISC_files”.

We first construct an instance and then point the reader at a directory containing compatible .nc files. Since there are other files in there, we must be explicit and use a globbing pattern; supplying incompatible files will make the reader fail.

The reader actually calls climada.util.files_handler.get_file_names, so it’s also possible to hand it an explicit list of filenames, or a dirname, or even a list of glob patterns or directories.

from climada.hazard import StormEurope
from climada.util.constants import WS_DEMO_NC

storm_instance = StormEurope.from_footprints(WS_DEMO_NC)

# WISC_files = '/path/to/folder/C3S_WISC_FOOTPRINT_NETCDF_0100/fp_era[!er5]*_0.nc'
# storm_instance = StormEurope.from_footprints(WISC_files)
$CLIMADA_SRC/climada/hazard/centroids/centr.py:822: UserWarning: Geometry is in a geographic CRS. Results from 'buffer' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.

  xy_pixels = self.geometry.buffer(res / 2).envelope


Let’s quickly see what attributes this class brings with it:

Type:           StormEurope
String form:    <climada.hazard.storm_europe.StormEurope object at 0x7f2a986b4c70>
File:           ~/code/climada_python/climada/hazard/storm_europe.py
A hazard set containing european winter storm events. Historic storm
events can be downloaded at http://wisc.climate.copernicus.eu/ and read
with `from_footprints`. Weather forecasts can be automatically downloaded from
https://opendata.dwd.de/ and read with from_icon_grib(). Weather forecast
from the COSMO-Consortium http://www.cosmo-model.org/ can be read with

ssi_wisc : np.array, float
    Storm Severity Index (SSI) as recorded in
    the footprint files; apparently not reproducible from the footprint
    values only.
ssi : np.array, float
    SSI as set by set_ssi; uses the Dawkins
    definition by default.
Init docstring: Calls the Hazard init dunder. Sets unit to 'm/s'.

You could also try listing all permissible methods with dir(storm_instance), but since that would include the methods from the Hazard base class, you wouldn’t know what’s special. The best way is to read the source: uncomment the following statement to read more.

# StormEurope??

Into the Storm Severity Index (SSI)

The SSI, according to Dawkins et al. 2016 or Lamb and Frydendahl, 1991, can be set using set_ssi. For demonstration purposes, I show the default arguments. (Check also the defaults using storm_instance.calc_ssi?, the method for which set_ssi is a wrapper.)

We won’t be using the plot_ssi functionality just yet, because we only have two events; the graph really isn’t informative. After this, we’ll generate some more storms to make that plot more aesthetically pleasing.

    method = 'wind_gust',
    intensity = storm_instance.intensity,
    # the above is just a more explicit way of passing the default
    on_land = True,
    threshold = 25,
    sel_cen = None
    # None is default. sel_cen could be used to subset centroids

Probabilistic Storms

This class allows generating probabilistic storms from historical ones according to a method outlined in Schwierz et al. 2010. This means that per historical event, we generate 29 new ones with altered intensities. Since it’s just a bunch of vector operations, this is pretty fast.

However, we should not return the entire probabilistic dataset in-memory: in trials, this used up 60 GB of RAM, thus requiring a great amount of swap space. Instead, we must select a country by setting the reg_id parameter to an ISO_N3 country code used in the Natural Earth dataset. It is also possible to supply a list of ISO codes. If your machine is up for the job of handling the whole dataset, set the reg_id parameter to None.

Since assigning each centroid a country ID is a rather inefficient affair, you may need to wait a minute or two for the entire WISC dataset to be processed. For the small demo dataset, it runs pretty quickly.

storm_prob = storm_instance.generate_prob_storms(reg_id=528)
2020-03-05 10:29:31,845 - climada.hazard.centroids.centr - INFO - Setting geometry points.
2020-03-05 10:29:32,248 - climada.hazard.centroids.centr - DEBUG - Setting region_id 9944 points.
2020-03-05 10:29:32,466 - climada.util.coordinates - DEBUG - Setting region_id 9944 points.
2020-03-05 10:29:33,506 - climada.hazard.storm_europe - INFO - Commencing probabilistic calculations
2020-03-05 10:29:33,620 - climada.hazard.storm_europe - INFO - Generating new StormEurope instance
2020-03-05 10:29:33,663 - climada.util.checker - DEBUG - Hazard.ssi not set.
2020-03-05 10:29:33,664 - climada.util.checker - DEBUG - Hazard.ssi_wisc not set.
2020-03-05 10:29:33,665 - climada.util.checker - DEBUG - Hazard.event_name not set. Default values set.
C:\shortpaths\GitHub\climada_python\climada\util\plot.py:311: UserWarning: Tight layout not applied. The left and right margins cannot be made large enough to accommodate all axes decorations.
Wall time: 2.24 s
<cartopy.mpl.geoaxes.GeoAxesSubplot at 0x1dafba69940>

We can get much more fancy in our calls to generate_prob_storms; the keyword arguments after ssi_args are passed on to _hist2prob, allowing us to tweak the probabilistic permutations.

ssi_args = {
    'on_land': True,
    'threshold': 25,

storm_prob_xtreme = storm_instance.generate_prob_storms(
    reg_id=[56, 528], # BEL and NLD

We can now check out the SSI plots of both these calculations. The comparison between the historic and probabilistic ssi values, only makes sense for the full dataset.

(<Figure size 1080x720 with 1 Axes>,
 <AxesSubplot:xlabel='Exceedance Frequency [1/a]', ylabel='Storm Severity Index'>)