Forecast class¶
This class deals with weather forecasts and uses CLIMADA Impact.calc() to forecast impacts of weather events on society. It mainly does one thing: - it contains all plotting and other functionality that are specific for weather forecasts, impact forecasts and warnings
The class is different from the Impact class especially because features of the Impact class like Exceedence frequency curves, annual average impact etc, do not make sense if the hazard is e.g. a 5 day weather forecast. As the class is relatively new, there might be future changes to the datastructure, the methods, and the parameters used to call the methods.
Example: forecast of building damages due to wind in Switzerland¶
Before using the forecast class, hazard, exposure and vulnerability need to be created. The hazard looks at the weather forecast from today for an event with two days lead time (meaning the day after tomorrow). generate_WS_forecast_hazard is used to download a current weather forecast for wind gust from opendata.dwd.de. An Impact funtion for building damages due to storms is created. And with only a few lines of code, a LitPop exposure for Switzerland is generated, and the impact is calculated with a default impact function. With a further line of code, the mean damage per grid point for the day after tomorrow is plotted on a map.
[2]:
from datetime import datetime
from cartopy import crs as ccrs
from climada.util.config import CONFIG
from climada.engine.forecast import Forecast
from climada.hazard.storm_europe import StormEurope, generate_WS_forecast_hazard
from climada.entity.impact_funcs.storm_europe import ImpfStormEurope
from climada.entity import ImpactFuncSet
from climada.entity import LitPop
---------------------------------------------------------------------------
FileNotFoundError Traceback (most recent call last)
<ipython-input-2-8371ba553ad9> in <module>
2 from cartopy import crs as ccrs
3
----> 4 from climada.util.config import CONFIG
5 from climada.engine.forecast import Forecast
6 from climada.hazard.storm_europe import StormEurope, generate_WS_forecast_hazard
C:\shortpaths\GitHub\climada_python\climada\__init__.py in <module>
22 from pathlib import Path
23
---> 24 from .util.config import CONFIG, setup_logging
25 from .util.constants import *
26
C:\shortpaths\GitHub\climada_python\climada\util\__init__.py in <module>
20 """
21 import logging
---> 22 from pint import UnitRegistry
23
24 from .config import *
~\anaconda3\envs\climada_env_20210809\lib\site-packages\pint\__init__.py in <module>
12 """
13
---> 14 from .context import Context
15 from .errors import ( # noqa: F401
16 DefinitionSyntaxError,
~\anaconda3\envs\climada_env_20210809\lib\site-packages\pint\context.py in <module>
13 from collections import ChainMap, defaultdict
14
---> 15 from .definitions import Definition, UnitDefinition
16 from .errors import DefinitionSyntaxError
17 from .util import ParserHelper, SourceIterator, to_units_container
~\anaconda3\envs\climada_env_20210809\lib\site-packages\pint\definitions.py in <module>
11 from collections import namedtuple
12
---> 13 from .converters import LogarithmicConverter, OffsetConverter, ScaleConverter
14 from .errors import DefinitionSyntaxError
15 from .util import ParserHelper, UnitsContainer, _is_dim
~\anaconda3\envs\climada_env_20210809\lib\site-packages\pint\converters.py in <module>
10
11
---> 12 from .compat import HAS_NUMPY, exp, log # noqa: F401
13
14
~\anaconda3\envs\climada_env_20210809\lib\site-packages\pint\compat.py in <module>
161 # xarray (DataArray, Dataset, Variable)
162 try:
--> 163 from xarray import DataArray, Dataset, Variable
164
165 upcast_types += [DataArray, Dataset, Variable]
~\anaconda3\envs\climada_env_20210809\lib\site-packages\xarray\__init__.py in <module>
1 import pkg_resources
2
----> 3 from . import testing, tutorial, ufuncs
4 from .backends.api import (
5 load_dataarray,
~\anaconda3\envs\climada_env_20210809\lib\site-packages\xarray\testing.py in <module>
6 import numpy as np
7
----> 8 from xarray.core import duck_array_ops, formatting, utils
9 from xarray.core.dataarray import DataArray
10 from xarray.core.dataset import Dataset
~\anaconda3\envs\climada_env_20210809\lib\site-packages\xarray\core\duck_array_ops.py in <module>
13 import pandas as pd
14
---> 15 from . import dask_array_compat, dask_array_ops, dtypes, npcompat, nputils
16 from .nputils import nanfirst, nanlast
17 from .pycompat import (
~\anaconda3\envs\climada_env_20210809\lib\site-packages\xarray\core\dask_array_compat.py in <module>
3 import numpy as np
4
----> 5 from .pycompat import dask_version
6
7 try:
~\anaconda3\envs\climada_env_20210809\lib\site-packages\xarray\core\pycompat.py in <module>
53
54
---> 55 dsk = DuckArrayModule("dask")
56 dask_version = dsk.version
57 dask_array_type = dsk.type
~\anaconda3\envs\climada_env_20210809\lib\site-packages\xarray\core\pycompat.py in __init__(self, mod)
23
24 if mod == "dask":
---> 25 duck_array_type = (import_module("dask.array").Array,)
26 elif mod == "pint":
27 duck_array_type = (duck_array_module.Quantity,)
~\anaconda3\envs\climada_env_20210809\lib\importlib\__init__.py in import_module(name, package)
125 break
126 level += 1
--> 127 return _bootstrap._gcd_import(name[level:], package, level)
128
129
~\anaconda3\envs\climada_env_20210809\lib\site-packages\dask\array\__init__.py in <module>
1 try:
2 from ..base import compute
----> 3 from . import backends, fft, lib, linalg, ma, overlap, random
4 from .blockwise import atop, blockwise
5 from .chunk_types import register_chunk_type
~\anaconda3\envs\climada_env_20210809\lib\site-packages\dask\array\fft.py in <module>
5
6 try:
----> 7 import scipy
8 import scipy.fftpack
9 except ImportError:
~\anaconda3\envs\climada_env_20210809\lib\site-packages\scipy\__init__.py in <module>
134
135 # Allow distributors to run custom init code
--> 136 from . import _distributor_init
137
138 from scipy._lib import _pep440
~\anaconda3\envs\climada_env_20210809\lib\site-packages\scipy\_distributor_init.py in <module>
57 os.chdir(libs_path)
58 for filename in glob.glob(os.path.join(libs_path, '*dll')):
---> 59 WinDLL(os.path.abspath(filename))
60 finally:
61 os.chdir(owd)
~\anaconda3\envs\climada_env_20210809\lib\ctypes\__init__.py in __init__(self, name, mode, handle, use_errno, use_last_error, winmode)
371
372 if handle is None:
--> 373 self._handle = _dlopen(self._name, mode)
374 else:
375 self._handle = handle
FileNotFoundError: Could not find module 'C:\Users\ThomasRoosli\anaconda3\envs\climada_env_20210809\lib\site-packages\scipy\.libs\libbanded5x.UGR6EUQPIWHQH7SL62IWIXB5545VDNQZ.gfortran-win_amd64.dll' (or one of its dependencies). Try using the full path with constructor syntax.
[ ]:
#generate hazard
hazard, haz_model, run_datetime, event_date = generate_WS_forecast_hazard()
# #generate hazard with with forecasts from past dates (works only if the files have already been downloaded)
# hazard, haz_model, run_datetime, event_date = generate_WS_forecast_hazard(
# run_datetime=datetime(2021,3,7),
# event_date=datetime(2021,3,11))
[ ]:
#generate vulnerability
impact_function = ImpfStormEurope.from_welker()
impact_function_set = ImpactFuncSet()
impact_function_set.append(impact_function)
[ ]:
#generate exposure and save to file
filename_exp = CONFIG.local_data.save_dir.dir() / ('exp_litpop_Switzerland.hdf5')
if filename_exp.exists():
exposure = LitPop.from_hdf5(filename_exp)
else:
exposure = LitPop.from_countries('Switzerland', reference_year=2020)
exposure.write_hdf5(filename_exp)
[ ]:
#create and calculate Forecast
CH_WS_forecast = Forecast({run_datetime: hazard}, exposure, impact_function_set)
CH_WS_forecast.calc()
[ ]:
CH_WS_forecast.plot_imp_map(save_fig=False,close_fig=False,proj=ccrs.epsg(2056))
Here you see a different plot highlighting the spread of the impact forecast calculated from the different ensemble members of the weather forecast.
[ ]:
CH_WS_forecast.plot_hist(save_fig=False,close_fig=False)
It is possible to color the pixels depending on the probability that a certain threshold of impact is reach at a certain grid point
[ ]:
CH_WS_forecast.plot_exceedence_prob(threshold=5000, save_fig=False, close_fig=False,proj=ccrs.epsg(2056))
It is possible to color the cantons of Switzerland with warning colors, based on aggregated forecasted impacts in their area.
[ ]:
import fiona
from cartopy.io import shapereader
from climada.util.config import CONFIG
#create a file containing the polygons of Swiss cantons using natural earth
cantons_file = CONFIG.local_data.save_dir.dir() / 'cantons.shp'
adm1_shape_file = shapereader.natural_earth(resolution='10m',
category='cultural',
name='admin_1_states_provinces')
if not cantons_file.exists():
with fiona.open(adm1_shape_file, 'r') as source:
with fiona.open(
cantons_file, 'w',
**source.meta) as sink:
for f in source:
if f['properties']['adm0_a3'] == 'CHE':
sink.write(f)
CH_WS_forecast.plot_warn_map(str(cantons_file),
decision_level = 'polygon',
thresholds=[100000,500000,
1000000,5000000],
probability_aggregation='mean',
area_aggregation='sum',
title="Building damage warning",
explain_text="warn level based on aggregated damages",
save_fig=False,
close_fig=False,
proj=ccrs.epsg(2056))
Example 2: forecast of wind warnings in Switzerland¶
Instead of a fully fledged socio-economic impact of storms, one can also simplify the hazard, exposure, vulnerability model, by looking at a “neutral” exposure (=1 at every gridpoint) and using a step function as impact function to arrive at warn levels. It also shows how the attributes hazard, exposure or vulnerability can be set before calling calc(), and are then considered in the forecast instead of the defined defaults.
[ ]:
from pandas import DataFrame
import numpy as np
from climada.entity.exposures import Exposures
from climada.entity.impact_funcs import ImpactFunc, ImpactFuncSet
import climada.util.plot as u_plot
### generate exposure
# find out which hazard coord to consider
CHE_borders = u_plot._get_borders(np.stack([exposure.gdf.latitude.values,
exposure.gdf.longitude.values],
axis=1)
)
centroid_selection = np.logical_and(np.logical_and(hazard.centroids.lat >= CHE_borders[2],
hazard.centroids.lat <= CHE_borders[3]),
np.logical_and(hazard.centroids.lon >= CHE_borders[0],
hazard.centroids.lon <= CHE_borders[1])
)
# Fill DataFrame with values for a "neutral" exposure (value = 1)
exp_df = DataFrame()
exp_df['value'] = np.ones_like(hazard.centroids.lat[centroid_selection]) # provide value
exp_df['latitude'] = hazard.centroids.lat[centroid_selection]
exp_df['longitude'] = hazard.centroids.lon[centroid_selection]
exp_df['impf_WS'] = np.ones_like(hazard.centroids.lat[centroid_selection], int)
# Generate Exposures
exp = Exposures(exp_df)
exp.check()
exp.value_unit = 'warn_level'
### generate impact functions
## impact functions for hazard based warnings
imp_fun_low = ImpactFunc()
imp_fun_low.haz_type = 'WS'
imp_fun_low.id = 1
imp_fun_low.name = 'warn_level_low_elevation'
imp_fun_low.intensity_unit = 'm/s'
imp_fun_low.intensity = np.array([0.0, 19.439,
19.44, 24.999,
25.0, 30.549,
30.55, 38.879,
38.88, 100.0])
imp_fun_low.mdd = np.array([1.0, 1.0,
2.0, 2.0,
3.0, 3.0,
4.0, 4.0,
5.0, 5.0])
imp_fun_low.paa = np.ones_like(imp_fun_low.mdd)
imp_fun_low.check()
# fill ImpactFuncSet
impf_set = ImpactFuncSet()
impf_set.append(imp_fun_low)
[ ]:
#create and calculate Forecast
warn_forecast = Forecast({run_datetime: hazard}, exp, impf_set)
warn_forecast.calc()
The each grid point now has a warnlevel between 1-5 assigned for each event. Now the cantons can be colored based on a threshold on a grid point level. for each warning level it is assessed if 50% of grid points in the area of a canton has at least a 50% probability of reaching the specified threshold.
[ ]:
warn_forecast.plot_warn_map(cantons_file,
thresholds=[2,3,4,5],
decision_level = 'exposure_point',
probability_aggregation=0.5,
area_aggregation=0.5,
title="DWD ICON METEOROLOGICAL WARNING",
explain_text="warn level based on wind gust thresholds",
save_fig=False,
close_fig=False,
proj=ccrs.epsg(2056))
Example: Tropical Cylcone¶
It would be nice to add an example using the tropical cyclone forecasts from the class TCForecast. This has not yet been done.
[ ]:
[ ]:
[ ]: