From 6fedaf296fc52689a510cd314ace9b8facc77a7f Mon Sep 17 00:00:00 2001
From: Thomas Vogt <57705593+tovogt@users.noreply.github.com>
Date: Tue, 27 Jun 2023 11:43:11 +0200
Subject: [PATCH] TCTracks: improve hdf5 I/O (#735)

Use more efficient file format for cycling TCTracks objects.

This change is not backwards compatible. Users will have to
reload the original data and store the TCTracks objects again.

---------

Co-authored-by: emanuel-schmid <schmide@ethz.ch>
Co-authored-by: Lukas Riedel <34276446+peanutfun@users.noreply.github.com>
---
 CHANGELOG.md                                  |   1 +
 climada/hazard/tc_tracks.py                   | 189 ++++++++----------
 .../hazard/test/data/tctracks_hdf5_legacy.nc  | Bin 0 -> 89975 bytes
 climada/hazard/test/test_tc_tracks.py         |   6 +
 4 files changed, 93 insertions(+), 103 deletions(-)
 create mode 100644 climada/hazard/test/data/tctracks_hdf5_legacy.nc

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 9105b34c6..4f2c1ef9d 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -51,6 +51,7 @@ Removed:
 - Users can opt-out of the climada specific logging definitions and freely configure logging to their will, by setting the config value `logging.managed` to `false`. [#724](https://github.com/CLIMADA-project/climada_python/pull/724)
 - Add option to read additional variables from IBTrACS when using `TCTracks.from_ibtracs_netcdf` [#728](https://github.com/CLIMADA-project/climada_python/pull/728)
 - The `haz_type` attribute has been moved from `climada.hazard.tag.Tag` to `climada.hazard.Hazard` itself. [#736](https://github.com/CLIMADA-project/climada_python/pull/736)
+- New file format for `TCTracks` I/O with better performance. This change is not backwards compatible: If you stored `TCTracks` objects with `TCTracks.write_hdf5`, reload the original data and store them again. [#735](https://github.com/CLIMADA-project/climada_python/pull/735)
 
 ### Fixed
 
diff --git a/climada/hazard/tc_tracks.py b/climada/hazard/tc_tracks.py
index ccf959920..5e7b3916b 100644
--- a/climada/hazard/tc_tracks.py
+++ b/climada/hazard/tc_tracks.py
@@ -22,12 +22,10 @@
 __all__ = ['CAT_NAMES', 'SAFFIR_SIM_CAT', 'TCTracks', 'set_category']
 
 # standard libraries
-import contextlib
 import datetime as dt
 import itertools
 import logging
 from typing import Optional, List
-import pathlib
 import re
 import shutil
 import warnings
@@ -53,10 +51,6 @@
 from sklearn.metrics import DistanceMetric
 import statsmodels.api as sm
 import xarray as xr
-from xarray.backends import NetCDF4DataStore
-from xarray.backends.api import dump_to_store
-from xarray.backends.common import ArrayWriter
-from xarray.backends.store import StoreBackendEntrypoint
 
 # climada dependencies
 from climada.util import ureg
@@ -1366,21 +1360,27 @@ def write_hdf5(self, file_name, complevel=5):
             Specifies a compression level (0-9) for the zlib compression of the data.
             A value of 0 or None disables compression. Default: 5
         """
-        # change dtype from bool to int to be NetCDF4-compliant, this is undone later
+        data = []
         for track in self.data:
+            # convert "time" into a data variable and use a neutral name for the steps
+            track = track.rename(time="step")
+            track["time"] = ("step", track["step"].values)
+            track["step"] = np.arange(track.sizes["step"])
+            # change dtype from bool to int to be NetCDF4-compliant
             track.attrs['orig_event_flag'] = int(track.attrs['orig_event_flag'])
-        try:
-            encoding = {
-                f'track{i}': {var: dict(zlib=True, complevel=complevel) for var in track.data_vars}
-                for i, track in enumerate(self.data)
-            }
-            ds_dict = {f'track{i}': track for i, track in enumerate(self.data)}
-            LOGGER.info('Writing %d tracks to %s', self.size, file_name)
-            _xr_to_netcdf_multi(file_name, ds_dict, encoding=encoding)
-        finally:
-            # ensure to undo the temporal change of dtype from above
-            for track in self.data:
-                track.attrs['orig_event_flag'] = bool(track.attrs['orig_event_flag'])
+            data.append(track)
+
+        # concatenate all data sets along new dimension "storm"
+        ds_combined = xr.combine_nested(data, concat_dim=["storm"])
+        ds_combined["storm"] = np.arange(ds_combined.sizes["storm"])
+
+        # convert attributes to data variables of combined dataset
+        df_attrs = pd.DataFrame([t.attrs for t in data], index=ds_combined["storm"].to_series())
+        ds_combined = xr.merge([ds_combined, df_attrs.to_xarray()])
+
+        encoding = {v: dict(zlib=True, complevel=complevel) for v in ds_combined.data_vars}
+        LOGGER.info('Writing %d tracks to %s', self.size, file_name)
+        ds_combined.to_netcdf(file_name, encoding=encoding)
 
     @classmethod
     def from_hdf5(cls, file_name):
@@ -1396,16 +1396,32 @@ def from_hdf5(cls, file_name):
         tracks : TCTracks
             TCTracks with data from the given HDF5 file.
         """
-        ds_dict = _xr_open_dataset_multi(file_name, prefix="track")
-        track_no = sorted(int(key[5:]) for key in ds_dict.keys())
+        _raise_if_legacy_or_unknown_hdf5_format(file_name)
+        ds_combined = xr.open_dataset(file_name)
+
+        # when writing '<U*' and reading in again, xarray reads as dtype 'object'. undo this:
+        for varname in ds_combined.data_vars:
+            if ds_combined[varname].dtype == "object":
+                ds_combined[varname] = ds_combined[varname].astype(str)
         data = []
-        for i in track_no:
-            track = ds_dict[f'track{i}']
+        for i in range(ds_combined.sizes["storm"]):
+            # extract a single storm and restrict to valid time steps
+            track = (
+                ds_combined
+                .isel(storm=i)
+                .dropna(dim="step", how="any", subset=["time", "lat", "lon"])
+            )
+            # convert the "time" variable to a coordinate
+            track = track.drop_vars(["storm", "step"]).rename(step="time")
+            track = track.assign_coords(time=track["time"]).compute()
+            # convert 0-dimensional variables to attributes:
+            attr_vars = [v for v in track.data_vars if track[v].ndim == 0]
+            track = (
+                track
+                .assign_attrs({v: track[v].item() for v in attr_vars})
+                .drop_vars(attr_vars)
+            )
             track.attrs['orig_event_flag'] = bool(track.attrs['orig_event_flag'])
-            # when writing '<U*' and reading in again, xarray reads as dtype 'object'. undo this:
-            for varname in track.data_vars:
-                if track[varname].dtype == "object":
-                    track[varname] = track[varname].astype(str)
             data.append(track)
         return cls(data)
 
@@ -1541,88 +1557,55 @@ def _one_interp_data(track, time_step_h, land_geom=None):
             track_land_params(track_int, land_geom)
         return track_int
 
-def _xr_to_netcdf_multi(path, ds_dict, encoding=None):
-    """Write multiple xarray Datasets to separate groups in a single NetCDF4 file
-
-    Contrary to xarray's `to_netcdf` functionality, this only supports the "NETCDF4" format and the
-    "netcdf4" engine since the groups feature has been introduced by NetCDF version 4.
-
-    Parameters
-    ----------
-    path : str or Path
-        Path of the target NetCDF file.
-    ds_dict : dict whose keys are group names and values are xr.Dataset
-        Each xr.Dataset in the dict is stored in the group identified by its key in the dict.
-        Note that an empty string ("") is a valid group name and refers to the root group.
-    encoding : dict whose keys are group names and values are dict, optional
-        For each dataset/group, one dict that is compliant with the format of the `encoding`
-        keyword parameter in `xr.Dataset.to_netcdf`. Default: None
-    """
-    # pylint: disable=protected-access
-    path = str(pathlib.Path(path).expanduser().absolute())
-    with contextlib.closing(NetCDF4DataStore.open(path, "w", "NETCDF4", None)) as store:
-        writer = ArrayWriter()
-        for group, dataset in ds_dict.items():
-            store._group = group
-            unlimited_dims = dataset.encoding.get("unlimited_dims", None)
-            encoding = None if encoding is None or group not in encoding else encoding[group]
-            dump_to_store(dataset, store, writer, encoding=encoding, unlimited_dims=unlimited_dims)
-
-def _xr_open_dataset_multi(path, prefix=""):
-    """Read multiple xarray Datasets from groups contained in a single NetCDF4 file
-
-    The data is loaded into memory
+def _raise_if_legacy_or_unknown_hdf5_format(file_name):
+    """Raise an exception if the HDF5 format of the file is not supported
 
-    Contrary to xarray's `open_dataset` functionality, this only supports the "netcdf4" engine
-    since the groups feature has been introduced by NetCDF version 4.
+    Since the HDF5 format changed without preserving backwards compatibility, a verbose error
+    message is produced if users attempt to open a file in the legacy format. This function also
+    raises an error if the file is not supported for other (unknown) reasons.
 
     Parameters
     ----------
-    path : str or Path
-        Path of the NetCDF file to read.
-    prefix : str, optional
-        If given, only read groups whose name starts with this prefix. Default: ""
+    file_name : str or Path
+        Path to a NetCDF-compliant HDF5 file.
 
-    Returns
-    -------
-    ds_dict : dict whose keys are group names and values are xr.Dataset
-        Each xr.Dataset in the dict is taken from the group identified by its key in the dict.
-        Note that an empty string ("") is a valid group name and refers to the root group.
+    Raises
+    ------
+    ValueError
     """
-    # pylint: disable=protected-access
-    path = str(pathlib.Path(path).expanduser().absolute())
-    ds_dict = {}
-    with contextlib.closing(NetCDF4DataStore.open(path, "r", "NETCDF4", None)) as store:
-        groups = [g for g in _xr_nc4_groups_from_store(store) if g.startswith(prefix)]
-        store_entrypoint = StoreBackendEntrypoint()
-        LOGGER.info('Reading %d datasets from %s', len(groups), path)
-        for group in groups:
-            store._group = group
-            ds = store_entrypoint.open_dataset(store)
-            ds.load()
-            ds_dict[group] = ds
-    return ds_dict
-
-def _xr_nc4_groups_from_store(store):
-    """List all groups contained in the given NetCDF4 data store
-
-    Parameters
-    ----------
-    store : xarray.backend.NetCDF4DataStore
-
-    Returns
-    -------
-    list of str
-    """
-    # pylint: disable=protected-access
-    def iter_groups(ds, prefix=""):
-        groups = [""]
-        for group_name, group_ds in ds.groups.items():
-            groups.extend([f"{prefix}{group_name}{subgroup}"
-                           for subgroup in iter_groups(group_ds, prefix="/")])
-        return groups
-    with store._manager.acquire_context(False) as root:
-        return iter_groups(root)
+    test_ds = xr.open_dataset(file_name)
+    if len(test_ds.dims) > 0:
+        # The legacy format only has data in the subgroups, not in the root.
+        # => This cannot be the legacy file format!
+        return
+    # After this line, it is sure that the format is not supported (since there is no data in the
+    # root group). Before raising an exception, we double-check if it is the legacy format.
+    try:
+        # Check if the file has groups 'track{i}' by trying to access the first group.
+        with xr.open_dataset(file_name, group="track0") as ds_track:
+            # Check if the group actually contains track data:
+            is_legacy = (
+                "time" in ds_track.dims
+                and "max_sustained_wind" in ds_track.variables
+            )
+    except OSError as err:
+        if "group not found" in str(err):
+            # No 'track0' group. => This cannot be the legacy file format!
+            is_legacy = False
+        else:
+            # A different (unexpected) error occurred. => Re-raise the exception.
+            raise
+
+    raise ValueError(
+        (
+            f"The file you try to read ({file_name}) is in a format that is no longer"
+            " supported by CLIMADA. Please store the data again using"
+            " TCTracks.write_hdf5. If you struggle to convert the data, please open an"
+            " issue on GitHub."
+        ) if is_legacy else (
+            f"Unknown HDF5/NetCDF file format: {file_name}"
+        )
+    )
 
 def _read_one_gettelman(nc_data, i_track):
     """Read a single track from Andrew Gettelman's NetCDF dataset
diff --git a/climada/hazard/test/data/tctracks_hdf5_legacy.nc b/climada/hazard/test/data/tctracks_hdf5_legacy.nc
new file mode 100644
index 0000000000000000000000000000000000000000..a588e9c3e586d296f77c9240ccf6ac2023bab10c
GIT binary patch
literal 89975
zcmeHQ2YeJo7vB_O2%#8yO{k%hLJRfwHh};klu#82Aq4^<!4x_Y1O)^EQK}75s-UPS
zh=5_~T`4vYP^=Wij((Kyy`6cN&Ey6WxWjhy%jE6<?cB`l+c*8q?99mU4#j;c`M7&~
zE97o+HzwB2U<L6T11hW>-7PY_r@QY)_ZkMz!;9adpNQ1QT`}S)RSE<8R%gLG@Lwq&
z{Gx0X#S<<)+!Qea!YqOb_pu+u@ZWh=#0xi1mQO|T%1n=o9})acB)^~gOu#7D?_75y
z`^A0TDgVWIzzqrP;+0N}o(yC(Z~jYTCc$-NMDcBc7amFl#T)ZKGHz^aMpi~<Tykn+
zLhP92)P&fq)Z|R1Gz4<zenhHoBgirQs2uxId}1nOIwf{gdSXUKR(fLopNt&ht6~Hm
zN-4!xDWfQ9>B&Q56Gy`jW0O+ihAO4tG1-$^6n8gpm@9B1^=*uq^-xN%loR4I;~?3z
z(a8yk=|(0|8ko`KAt=X;*sP2=-&P7dnuk)9B^@7^nK(2pecWx+#Awj@1}I7q_<a$U
zW@_BXL?bCw91NqWQGGf^`__SX2cbM!sEp(UBT#W3C^#%EG`MB+=-|+h<{`fA73Dq%
z>jk0!VI(KSrlu+P2eh||plLCs=}piqd!kUzd3FnNn`DJwgJaj@-h1sErWcBv8^p+d
z2{|$xQ4e2~r&-iRN+8FhK8l_{ZZ=cX2e_U@@Hj}lOG?!n-!Q(MGEB^|(TSPy2}w=0
zw6yeu<WvxPhLJRjH^Ba*K?J@Hm~MG7IciaAOpGB(EPj~B@J?MLqGLLBi;j%}$wU5p
zfN64NBWoV-rbTw49D7T1tduwXEccS>-10Wp=V3HV9N!tjR+J}sFkgSnj7=w^)vggn
zCX2v*%w%RVD2bw!jG4QpLrhGBa$h9-q1!G+X@d!Lelx0eKtO0qrBBo8uNL6dmBANv
ze}KPH$tX&IQr=rBKN%t_;$ytR228FwlZ%x@?Cy@WmYWe8sT*i7lrSvDiif-V5nc(Q
zX5bG=R$fwBIq2;&X;$XaJc8Ll$91pA04*w4=5JU~P>RAur+C#slv`7lOpi-Q&dP|z
zIty#6No`qpcvy68zxH2<k@;GfLKS0$Ew)2)N=l!&l&nO%;<Q`hl$^=Fb=eE|Woj93
z(rOec$`c<VB`q}|B`y<{7^|ETP!>SGR2IUMM?$rL1(3b10D5qRw++0?egI-J9KJN>
z4*q$zXs&2#@qm5#T}+Qviv>)N>O@}3u#V3#6d3f{$ked3tn`e4jO5h#!~m?j8U?jz
z6x1{zIHXlj5d0}hWi9T-^M-sQFcY>aIW}7Bo}M-;F+DRGN<7L2UmMD1bYgl&a$0KJ
zkcO#7)zG$S<Ccws8xBiIYSuQmad1fEkPjNX_GYi15fP}b)M^r)2{G7p@M6r4Dmx}&
z^Bvc<O9CzlxFq0`fJ*`{3H*B`&{69aWpw$B0ebG35xQB1yN?giu@8by<kJXV<Wmnq
z<Wr|i<nvB#zW8}xn<Z{aU;KShTL^X`jB6Ae+M1oO8Nf9THO;Ps8NhW(0dsNLxcA$x
z*oB{fjlqM7XN)sVxUA_KFUtv79Bow)?i&*!CSE1eHC|R0a1{Yp6>v2HW7pQ~0?Yue
zs|y(Wt!7u#4B)z!fNKkwdZNVUeRhIMUE3AAs3+k10){?*9@j7rQWTh!$>U(EB#$#@
zdVoXaAfh|F4V`(ziTq#*KSaV0mGGNN_{}8z<`RAj3BRR;A12|$bdcf1=WQ+Fw~_GM
zO8D(0{Pq&QMZ#Ale5-`7N%-LszA<itd!igVNWyoN@FOLBJ`b5M{?3x{T_pS{3BRj^
zA1&c`lkj;l$d}I^lJGqx{1^$pmxSM&^Tl?ckAzQyeKCALN%;N}KFuzO@jNIAKS07C
zDB%wh_|J}Fj93YOu!J8c;SZ7U<0bqA&L@&63713xI2|hBVFFGT@NfZ-5O9isM+!Jq
zz-a>JVzhA|vR$!@bOC1wI8(q`0v;{kF#;Yd;Bf*TFW?CRo+#i+0)AM)lLeeD;3)!r
zM8H!8JWarl3ivSrPZ#hE0nZfh;{u)~;MoFxLcntbJXgT;1pK6c=L`5L0Y5F^1p;0u
z;6(ynEZ`*qUMk>c1pKUkpA+!&0)9cj%LM$QfR_vSC61|mi{=gL%M$(y34f)8ze>Vi
zE#a?`@YhQCuSodoB>eRf{;LxHYZCqj3IBBo{|yQMO$mRag#VU=|F(p`Ny6VO;lCr{
zZ;|ldmGHMp_}e7>?GpZb5`K=rr+)elflvCZOZYn_{9O|MZV7*nguhq9-zVYkm+%ir
z_y;BYL!3|Kp0A$Gtjm(F1H=>jJzhQQRh=b&gb`Z@{>az^OLLcq;<&nZpimzCf-N=q
zontWH?X;60XN1l9hb^4NbYxFMU0Y+dld%>=U0ZBCO%w{s17<-=;hKByDsoh(m|jXn
z2wVz$dGOxg?JWbk0pAD&1pCAl7x#PB+AA^~7UH3aDG8i}y4QH+>|<M%p8dg3a0gw-
z%hwCCNblH{-mqYn_+O-}-nCfrqOL|Q?%mCCI@ZrGv-0wE$HjL~7JGN_qGv~EpFL9K
z96k*zs=fHYg8G*ZGNqwf>lXGAhQ=#J-<fLbI@Wq{SW#BruD+F_iDPy6rLez9hkY&o
z7}J5+NG2mY+T^Z%)|a%YxMJnQRzLP~6<??z(7LH8H*APYNlZ<MW5!gaAx`W)gN0Q}
z;;78z_}HQ8Xgif0m#Qdl?|*xVQRJXv#RwczvG#`KCY6}DN32+_L)r7?R;*igtg~kO
zL;W&#(1W4LJrh^V`fYB_#$yKk(J^Gvx=A(4m;3PS<V%zL&Ae1~(aH7$zGx8G^?{#%
z+yB(5P2H;XYu)<j@h)F}_DtU<v0KO9xc14LtGYe^-jJUTJ+tw*FISxE_Wahjf|@l8
zIrc%Xg<r*w&$#&I=hr8U`Qh~TNf$pnTzS(E@9Db&SFC#}En)cZ&!;JKzE!FhQG$N&
zv_)pRkj=85@BceNd0ALTqr*Oh%ef@rl7LGBE(y3K;F7?<QvzcLy?g!J4_@i}+L^86
ze*XOUH!Zq-|Jjy@FCW`I>zgCbg-!bL!&kcbpV{Wtx!w58XLg*L@ae|+2_65485!NN
z+n-mTE4jHsnan@Z_ni!Azv=CmiF0bNKfZNa#a$~_eD&3W-p#TOUpjrd?Y2o1Cr-Ma
zopbi=(Y-s5M8(GR?=k3s2Tm`6`HPql;eBQOl~xjQC-Z(d-FEd?ND?mUl7LGBE(y3K
z;F7?<R{~=P#SdOqt>r`KJv`>lcw}Zr*qU{8W05X>)}L~BZ~fcI{>@$t*t0Wy%K4qs
zOBbtBwd!+MJeK)%RO$|?x?R@O2$u-!c&EHKTs@7v@^B%S1Y8nuNx&rmmjwQe5@1V=
z8w#+*n19ba-_zazm+_gqDr}0zJ9aijm+#DiV`l?{&;qH`H%0qGHwceYu(v}w|MS^1
zO<)<ndT(|qua_HJSihOk_b+*O1A4mN9iFyuyB+3NhA=tVStdRiGHN(`cf)U6?kd6o
z9cV6to~<K$D@BZh2Rs@e&v5+WJ^;o1Fc>?Oz;4;g;TgiRhrgq;S0vv)aM};<paTie
z)C0?2!T0qSYM(u`;OHE02t=_X^+DyDRCN3*dS6KIT8HN+i>v5WPOyyRuA=Yn_n^jj
z=Kqe*5pHVO@PeOVpyumPi-G<S^TWN{KGY!wOtF+*{4uz71J};nXY*9LD`S-aV6YG`
zc(Qwyco3qeHzsoXNyp0qH*CBf_GpYz+I-K52l-T3%kuL23`Qr`U^MSRJ|EU%+&PTg
zyMuhH**^^j<%qT-)}9f)*wHlXdj5%Mwvk~$C$9yDJ{G}n#UN-xT`S*<k*kl|7wy2d
zZtS2Z4!nE6Qr*|1$AXpTD!Kc&e`<5nMP8$eY#r?JM0WYq4qHB7m~^D|ynS9Rhn!J+
ze^h1L^%^~6>h=C~ahGN#&pzehcB$SwvN?ycTwo&U$eR%2)0_i#1zgT00ha_^5^zbt
zB>|TN{+$xQ)r|@V2D|rI_Cnor7T6*%&ZlW*Stp{YL<X~MU3?b6)rpV@%Efj`z$F2f
z1Y8nuN#I{30qjHsc-{yt-P6r|YoPZUWy<*FhlX7380l8wTHwoyZ}!@j-tLe6UZqNj
zX;ZdM^-$i(F7>6XkC6a!p=m6$h_3McV8o+7hHc`NZn&sR0xk)-B;b;OO9Czl{A(oO
z_ysAPw0LV{Ne_w{k1lu2c3wPFZ0_~mrP#vqEk90hyFKIe*-Bl`WVWVSdkTA@J@|X2
z6=P3luNbrBH`iq8M)Twuw+3ITv&{nGGMCxS3x8`^hn|t$`GN7KE*q_D@8ZGfivh->
z8uD;)&E$Mc*SER^so8$XdXoXrtDKfXzu7jf*8SFAQ5}>bK_^GLDf56|Q(63QR9&SG
z1TUvN65T{qP6M}!^8D43Rt1;O$}86gyM-%wz)wBJ_n8w9DGlMWhEly*$8pM4xU8TY
z-tgHIN+}}-|KHE+3d*Uf5?cPTv&sm#tgWma{p%9AaJj#M1dR9Z#k2#m!(OIMYuKyv
zpnS5|#|k6xn8FmHyG$p%2r|V%6A;5O-LS2wGF>W4cbU4_6r1T^;@r>=D1eJPM*>*y
zd;4MB9&Wwx{elbP>fn1xfzMcAdk}=@=Jbr{5YaOtT8oG^&dq5KT7s1reuyU{F*El7
zodoy*N-BQtLGi}>OkDzub9Lx_6V+|Hgz16{F#g;At*c?XG%iNf&Wn$$dGTDqpruPh
z-uo}d=f%@f^WuTOMa7^*Se1U6SuZGbTG}0bfe41^xA`_uNN|hj;GocuAZ9zku1BXX
zp+inhSsXiP;NJbkAL*`GRBMsZkNCWr^2p<xXMH-k(Y3K5KhGbtxTR0m<*{c+UHhiw
z!{^IayHfV-m8$Kpl-@LUZpo`NT2Iv5wmV-e`bwFjSE{UPxwF{TT@AMGsl0V}SmRMM
z8w@yKHC3-UqPK4nB}b*63N*blUmgHYJwj(F4^T3$s9<n&&uj*x5yMsIv=*m<oYv>G
z5vM_%hH~1RQ`&*rhU4~}S~;ax9m#PQPN|uu-a;?VM?)!e)L;)}P#MfAsY)WpB$#B5
zQ#kDq8E#eR3tY%sT5QDvH&3eRJn?%#bPl=v3pP{|+&roB_2J*RNsCXWV&UdVFU2qa
z*Se^zxOq~cDkXlMi;9n%XAz$Ha*3&<lhe~up~PmIKI28@%FWaE_VF=dGcprLQ5kdd
zq`JMxkU02S43#xEPtp_*_#6(EI5$tK+&$n!Ibylvx%o#D*1sKiTf@m)e8EZR-!_=a
z1VI{m&-%B%7G@D?&l)-Z=cnvjBg6NSqOjlTs9^*vmNlO1TiO_+w893Z@E>#|&APkK
zuPJNJbi>cDRl(@&KRyd_VzT7~(~I&Zy=cpGL9{R*Wt%9`38#LP%?&bK=#R1$@cU@f
z{*0X}v027s;oExC%6!0NIO|q2Eb!QodN~`jjOdkL7_r$KJJr_Tkfg;SrK>J0T?MyV
zF7VgLFq5q3bA37!4z-P@Llt)p{o36qiH><1htw4}=CJHyxlA5h{#b4oZrnz8hxg=^
z0(`15BqTICs6}wgX1qIQT*ud_#npzzXZoL%6s`WbzN>opP?UO9@1ll=c2WBSMO<b;
zR~7HC9Mv88F>2Pu-s;jN{neTS2EqLV^`*rl;r?j#*u;r|v(*n5PX+v#+BkHkswlJ7
zBa`Pq_$Sp8^%lbKpHbiH@&fQ*R!^*73;fsBhgNL@yjA_Z#CG@{^M{_;X7G}dcIt~h
z*{LTb?b2VYxLZGbXt&<A-yWcQ^lw`2)t4^Wt8f2xuO1q@Pk++7&!DkB`*cOwr_X4*
z55n!!hyJo3_y_fXpAJI!L;9exhatWrx^LH`z<*!AR_g-@e@rhn>lonU`jfp+!0%4#
zPmlfx@F#lm%})WJf-tA`(I1}%dPe`l+iE>`vz_7h{MOove%<AkhM)gsGpjv6?o^1?
z>|>8Lv4-z&Xl;;F@0K3JVSE_BIUP#ROetMT2Y*W*d325G+DjYnV>n8So)MpSr@xPO
zaCiyBAAieFLw`ysZ}=~|RMzbI3k&_VJpR$@+J(%T+NJo~xAYev9L9(7o7186%#_l#
zrGsgZ$MoqM(zTbiaUa4_TJ((gyoz!K^pLz_g{B028Snw%>w<4;$F~7q<Q;hY9ma?8
zV*Ca1m>0|w=8f`crpPzvnc_qqhT8{z1NddRK6_%gV&Pr^iauuo^VvY{c?<B(UcvC=
zA?$5DhRepV;LYhE7CnYTZ{c16N*;N1jp^D;8}DPd91LrE#%z2>aP#I`T=YL62gRGc
z=YJ5r-r|}Yl-Kqi`#f&IOnZ9~8+u!a=eb$j_&SA!#rFfBhWPu)<?i!`bsN=vSm$BA
z+8z8T@Co1tfX8|k!eM=oO?3mWuN169#zDOR)WYi%tUI`KMxTE5XcoDvPcMF97*wD4
z>5n|lWMrpLuZ4AN?~K>->=@T;vA-=h_qFt$A$=`xW<FkA+{QUF27iiQKBxJ#LEC4o
zH9j<6xd*#t`_On|-PY?GHOj`(5Dk8(WwO*3Gh&SH3$?4&I{0F|k=oj|EDBp{Pqch@
z+md$|taOGV%^S-SbCqv-`EDaG-^F&;(9K#<g^y+LSPL2~U%V}MkUjfR_I4KCiBd70
zjg6g3vtmXhPfy`-abj&ywj)!lrwT9`l>g)eelkRTnw3|StCMYIVjN54Tb7l3I3WaW
z;@|%82NVjM6v}yl$;nOw(x&fJn?4j*hT=*vI6p4XJNV*!OY;}!@A#bhNmXC$j%tc?
zm<6uoFw=|Y6ifAe&RGhV!atAZQn-=H35vW%yTLyTW=G-5A3gq^@n&KE@b14rK=X&V
zgc%vRzn$Y_E)Ov`+%!2&+y#G)x1zbt`Tq9LPJ+P7oMJZ@GX?B%+Zyi{XL8iXC3wT}
ziv$2v^EIBuZoSW!6=|W#sYh>M`EugX8@csWKY4xi_^)-=IQ}K$%D2zBLGEIo?eI@b
z%(*YQO!*YvTC+QfvJqao<8xFFqF-Ee3pC7QRwvGJ2Sln~2CmJxw1DpyFs|~O!xtBD
z2?3*ZfH^#BpBdu<v>D?%vl-*svKeE;ZN|99Y{otUE+$}Kjt!F)xZWXPUBEj9yi35l
z1-wVVdj-5t!21P!K)?q%ruT(vZGK+7He)QoX6(r^J?~2bep$dP1iVtfs|37Sz-t7I
z^%~}f(tkz3>jb=Bz^@7zZE(!-Z4mJ50)9il*d;Olexrcj67bss-X!470)9upTLg?>
z6f~#5RlwT>j9(5khksANIe<HA-J(!0NXEuBK;w5n=_`SRi4vdJ`uW16?w}LHW0^%K
z^2tU-<YO5|C-SKiF7nB4MC4=HMkn&IjH46zSk}>re5|w3iF~ZP(20Dqa1!~{q>6m3
zGtr5BvXBz_Sm&Y>`B*oj6ZzOZK}YhVuPb9Mk51%crvshH$94*x$gjxQ`S_J2{K^u(
zVIE;8{;HDjWFbcR`4$T&IDbZL`5}IQgkN34uOZ>rl<;eDK9SmzaCImEr-1^lE8uzp
zt}oyQ0&XbaMgndu;3fhN5^%79Lj)Wu;HCm@CgA1*ZXw{70wxPSVz;6TPFoAOjey$<
zxSfF83)muHRlrsOYXS}zaD;$62)LtwBL&=vV=C8uyE8^-3BQYkA0^>;mGGk_{B9CH
z&z7zJp?vm`gzqWg$4L0SB>dhIejf?HuY})E!tXEP)A*9!3w?F>K}q-l68=C5e~^S9
zE8!28@Z%)>ArgMPz^CU*knj^F{3Ho~sDwXE!cUg)hfDY)B>WW47u$=G5`L<LpC;jt
zlJFmr@Y5yy3<*C|!p{=;OHMJyXbFFegg;ioA1C3Dm+&WWK9Px%aFZwirw<EwvVgM%
zJVn5d2zaW1rwRB`0Y4_-=>nd?u~^<`O8AdU__HMZ*%JN}68;<sf3Ac-Pr`pv!k;hU
zKgIb}j!560mhcxy_zNZcMH2pE34e)%zf{71M#6tq;CKF$F`kp~pO^4oknopD_%BNM
z%Q>IOOOkLeQvgm^2zaG{R|$BvfY%6kt$<$<@Hzpn7x1eBeoep|1pK;y-w^Pd0^TU#
zw*>sQfHw(vvw+_b@D>5TE8wjH-X`Gf0)9`xIRf4xU|qmF1-whZy9K;Qz>c3UtOP0i
zxbU^+6!YEh2B!I|$`JDckA;GYQ#y{D!mkMFZ2szta>mrT?|Jcn-Gf}_^Z2`~%=_#<
z?ywREj;(ZE<C^Yn@L;;S3EuO7y{+E3oh|>~o7k+><V>YB1j052{fJcGMr?Jj9Q%>D
z7bpKuMh@|<Y5#!TqCC>*ZoOio)vgiv-Fa~A`hQcE;mL>296$~_4@OuJxW*r&^JT>N
zT0LI+aa@+u3Y_AbLRXCur8=jzI1S{qKBw@Qh6~~}l+)&%hH=`4)ApQNIn_9g;Isp$
zk(|;fD4Js$h0rL1Mgb4-`vW+o-uPgS<2g;@G?`QCO;T@=dUMnpqu$g69-ewPQ#hW+
zDfKd_w?M5uwbs<yQY%WWCbfFhYEi4RoX1P80M+hENp;l91M7gjcO5l$8lxK**9QCX
zT5A0IS_U0_fzhC<wZLw?rrIH~2H1;NS6`hIpn`3zdbE3hfs5{_X~gqZow^XN0sOY1
z`oD&a42st^!6Tji1fw<PHHO~@8UC4@usa;COIBg{Mrsp?tFgL!T|*--@89Yh^iRJA
z23@{`#eXu9#qnOTV8}zT`u%D~@%pR&p%6!Nb<v!b@SE0Z+86DBYT&EusbALwTk<;k
zZr}P~i{41j?iLI-=FRl#p<zJV>K#j}a6ep6O^VbNB}z}-)(yhN=$n@G1^S@AD=QYl
z#p`<xC4mk5aJ^*fR0x+2ybS%@Gg&~#K)7-GY`;lxovnxLo&w>f>Yp}$6!3JoH$y)X
zH4}KV^d(cD06bS8@!XR@pVD7lvH<8J{qYS;fIg$osqq|ydtQHR)(e1_>HS`N5%RWN
z5B%vRh--y@V*N_suY&Yf>))4N19Yu^{?J<By`rz2y$<MleO^dg>z9_c*88H{T9@{3
zYmG>3YyBg=t#wOwTkF?l!mRNlL#+Dv5bIYnLabQ}L#z+J9AdS+8Dfp!5n|oer?FLE
z7ieAbdZ0CHL!h<tnn3F>&j(spEef>upBHGI{#c;(%ZCH4#U}(>znT<iZ9Fs3T6#&K
z_2~y1TQ7VbVy*sbh&91G)H<y~sC8+bP;0Z_OKCrUQA&$_zm(>dQ%cKzt&~>n*;3k?
z8Ktzi{pB>D&-}G%C%_-{*Uswxn%5S8ZN>)htNgXCV`^w4-matdeY1{Ma$_BB`{p`Y
z^Y`j#XLi=nrtPbv?Kx6MTYkKbR`$a>TEvMuTA9Ojv|sf)TIZ}9+QBFMwPor4T0&oc
zZGpvKThq`Vs*6#ecQQW%?wRbd^IdO%`)IdaSdYKUAFq3*iQz$&Cx%B<nHU~cbz*pv
zY7@h|`%er{511I<dgp}jADT}HFWPcKcxc#!@E)xugpY1LA$&ud3E`i%oe;kD-1zV|
zy~c+x2md(uA7aLbR|G$~=lJlidW;V*2R^I&`0!fb2ZN6W{2bgX{`vUu8DLLY)jA>k
zO>IK>eGwDFgF0~gN=Msa51VT-s-eFry`SZ4na1eg(Ttv3#`C)e=+Y&hSc*0M#Grda
zPFW5gI&B#Y9_b&6XW;qHS=`&6GiWHM#gfh#_s@UG!hZzE!X_nsZRykT2Z-kn%gamL
zpq}zlcZC!Ky(p!obSw*XSw;2ti<P0itfH1FTGgNthpQO$DDaSvdqZ%a3GNf2JtSAx
z*R-B_Ilzjv)TTNHAMfLy3f!Mz?9HHkA8I}Aa!>>963|OM1C2cycpu*pGZrwPv(gsX
z;D~w@dav9s%UKM+Fk&-}a_UJdZjc!3hmFZ#=`Q2RGj7et7^IGct%t+dUpKxSU$Apl
z$b`{0(do|9;lw8|lp^R7+J1?dhZYuklQ(R`qZxi_zu(y>Dn?IXiDJJT9nmu<e>?3u
zxrNE(zxv;Ga?5}9&9mpmu`{a>%V#o9A1)C5J`!%gF9l=1GK!Lxo;)-*aWtIp5Sx?|
zH`Fk=V^e+>yPMcq6$S&w#$+{VjC%2$m)KG1X`_=964Q;bZZSv?E576*ndx!y8L?Rz
zab)UQlqDS>mzg*;Eq&Z=)1)J?$O6q<dXyI%rz?v-4mfI@yg{GL!DnB7y-}EO@&>Ju
z%TL}o$EWq}p1gr(S?F$g@<x}*6bM=Z!{s}9gY21ad-BHRsT7wn#0OB+cJjs`-tdVh
zZ=9V$zoD@{p1k3^hPc#@;mI2}#7Q|8=pXr?^?1xDZ%|K%Xu*3ttR!vN$^=1k45W@u
z@9+3p<fh>F`s59ova4#y;?%TCPh(X>Hlw~TG7*vSvhxp#q5{g&Wu}sJjg*QX@u;j3
zrXod{iV<NdK7^^@5EhHf={Fh~yD9GOERj;%7WTgNSt2~dwR0_IF-<q0B~p{;j%Z=d
z5^>wcek^LaQ_m8?hQe8AiPRarIb3p<2<w<Pu$~NLGBlqhLb4-T@Y;e3!|9iqBx%Ku
z>I<iyC1Oa@;*k2HF6$^b;Vco%B&$;&*ui8(^Rr;DPwtyZr5Df=z}=r1MD-R98Vj}K
z1M^xf_8RtyIz6`j#ry+lQV^nrSCnCV7WbYj3VR1@-eDp?!bq-Z9btR4OaDpUxt+|N
zH&~bs6x4%XcvZxU+zyf#b_Pd;9y-eITb^wlx$uPXmM#Ano!uzc**2|CcFvVbl@Hld
z@?NRpsIGUN+n*X38X6rO8q(Z$Y9L<E|DfI)W$E%Se{n`-_sd_U!w%{k-HXc3hnFxr
zmEESFc1%Xap|<1Omze6%+*+ZkEF7a%9Q4)NTuuCguTGQEMWLc&t)x$hrh~aU8|2uR
zi+edkk+N5^RPK3Y_zzN~1i2#Bj~(+0z4zkE^iP>E{ZJ*q(Foj&Ej#1FbcQin_lKMK
zo`<DnrDp`dfsOHr0l{Hmp^buCHVO_72nuRt{N3?p%dJ&vJ|J#yt8xspWZtd5YWDYD
z%1>8Y*I<fe-pxy`D^Y4?3%8rsms`bD;PSY4U9xa_^OpY+?4KTpvU2LD+;h8S+vN&2
z`Sn3Qs9=uySTl9YaGm|<t}D+vKIV}Svs%3BbmLrt10SFe;Z^KXls2dXYxWG&E`5}8
z(pd_(9uh+5QaF3Nd2<!60`CyoYZ@o}j4#-d4ZU64f(vW8*c*WVG}UEeNG?2CU1;28
z6-5Aa;i`Q^!S^<pe$`M^i@Z=5d|sxXUuNa{pRl;SAx3mAN^-zgNR`A0xHjYA0!E`9
z^ZlU$9wy*;0VfDJQNT$8MrAdp_n?5$*vNc;pnwMnI99-e1so^fAp-6s;LZX@V<~ff
zq6FMkz|jKkCgAP@#wjCn{5=I6Bj8>F?k!+kvM|So#$0CHPr&^Flk?ygUYye8w=c&f
zIsLmBx<&B(L~>eK!#30{bc2omih{b04t3UW*)6r}H(O{*b_>m2LE0Q}GoZ}?BW<eP
zycr5KR4aEQ1Zar%i+6|yhF99CM#0+QLqXcG(Lq{hXpr`^caWC!brY@Wn@u2G6Rlo+
z6W}+|S{-k!9gb?OwOZRqQ<O$p#jJ)9t^wR@0O9Iu&wN);OG>H-_rN*Y7=_|<LfvWB
zx%rTFGk#LQXwYZAkA{9`JV(HD1w2E*GX?y(fM*GqhE2eBxUpwaQXe(!_Xok?E?)hl
z>PS^l#;U)r&j#a`N7XH!vw=RTK0NwqpbNmDWv_np$9=k@9MF4RH~<E22lZzA4*@-_
z@7Z?*=utiE^Y`KYG2LwVf(NCN@r(HYS2Bcoa%OpJcxnmj=}zv}3!N{iRozdh^DpmJ
z*MGN7o%{7$>aH^{s(*g?gc>kvf_nJPacYNg<JFvY6V$TDryBh3pFE|eo_<EX@!1P%
zuX8KZmS4W6R{QoH<M$qycB=OLg&z*9yPm0P7~-H24(1E<h<QbW9hCLmI@SR-pXO1M
zM@f{8cnWfw2ep^AAO@OH*u9%Y7lg_Pjz&ICJ<S8ta^h(o7(JW9Oh1axl>4~`(=>AT
zlRxx!{HP)o5}R;+>mbL+oCYyB`}kar<6|BLF;{rtrR|Q7`5}loC9vK16tm-Y2BbsG
zT7``kqj>n96)pBBUV!`H)UO#`h}9#wycYXowlMB<jNChmeVA3j%?m}?Y)fj~$V6kZ
zr3mB$b+=Q`zMZ1EDZX)i`i)svF^>@y;(uoqg4tq^DgvT31(JWLQN#zQoIPw`j-X1Q
z)@qC$&JJh9RtY$D1v|Ht4&;Mj>{p^=Wq0Q$hOi|w9Cg(g3cKJo_0H_$O6V+xN{ra-
z<GO$P=Tet_92d%el(_8U{-ZCRKQ^{h)CH28lif=rH+mm;XE>6xta>aOj+Ce9@=!LJ
z)$pwndk#Kt*vHY1Hw;z$yB~`Q_HlHkirhY~;oB?$yY_LtH`7fIcA?n$GF8QQDG*G%
z87`lF+z@V(a$Ebjci*G9Y{Pz=eOyD{8i@9B>krXysNF&Pxb5!~*N0swXdhRNPt$SM
zKbo-qh1oui`V~YA-oIcCdE9v>5S#`=N9D`uy$fH9Bj#Zp^JOX;doioqd1t;dY*?rR
z_Um3lGPdOa=iA3o0cGhb3YAI1RQw22StCqEiZB%;!c=?+Q^6rD7MV6ZE$=hcGaIZG
z<H3m*#$fFYiP|~UV9f(E>@0(|fj3*EN(|O8Sw-2TO$}l)G#jju?1&b;_>y*R;KPH4
zMk=`l*r^6<h9oTxDZq7cSmgo6<%UUmd|Fz1LUJk?nX@JC`{1c?ECTg-09rsdU>}Ps
z0pX$LiH%Ooj890qh1^*sSL|CR96JkDTbqVWpb@a+4Az!HdWACAK!G;@=&TjU8{e&K
zlYbmZ3PQB-igNm|QWW+M9ykZj2JQ$W1v6Ny@#*zod^=>3hN^&$S)_$jT4o@HZc}H3
z-EKPRdv^QmJ7d?s>uJ1Y%RfeE6XiPFVoz^l)`ki8Fqs>%j^e#ESZfv%9UL6eoW84O
zT>Fe&lBxWX@#kn~RCc{4u{Y|f?8Gy+u>xWG@a`ohBRiFy-#e#&z+ymcNAdkE=xjy`
zS7l*Y4BlLK4o1jQUi%e`XVINR$6Cp!AAgKny|FmzE?AT43`NRzDbj6Q&9ax-JFufj
z8-J^;(R(khwETkwyeosX5HMJ4ewPMo$5W2Y5)IZ4UtyWcJ*vBN-P$!3%b5F~Vv1!p
zSfkdJD7CWrt-Wsb`H3aGSm9DPIvP_sze^S_Z{G4h!e3;tmRnP$$rWt=idRY)4R(|x
zzXpH@wgi9ic5TPU{4~Dj^y?mm!P*3_ayEmt21RH2Z`-@JfwL6uN0-7CXdJL#u5h1D
zN<309gEe})#<&DuJ9SF&8XV0%F(m&FomJZ}S=TK97J!caO24QWIrJw(Q7!U9oqrX{
zexc5pj0S5z4rdBzUa+Q`jwq|~3bSC1rzlZ}qSFDE`wjddsOEF(f;CnsPF%3gt*>hG
z@r0cQvFx_wx{j|^@L;Kd2mCS|AM-+pxkhTccO4%ycK#mve*0L*$GjL~uHF35DT>)~
z3)V{@=8K`-pErs}{srqs_p|c(lP@o)Wf;x7U|sifhy~npjQnjDtfM-GcT)L+wQ=qJ
z_PCmu$IHsqcUXmBUa+Q$fGBl=@>eyo1?wIMnFPc^np3Bgb4zK7i!AwQM);^G^9+Nt
zI9E7g#&6z>C+}>*`aM2Quf!uJ>epgX57&Y<N&CM^T?^L#&G$}EjV%>*f&Lx~)-*$U
zcMH~~-f#cV?G~&zwOV?;;%Q^Snhd^BEZ*~nY!w!)>5yvqg7s&onMCa_SbLqNo7e`R
zv$J4*fH%%=o+N>M3)ZE-WWT=cf^`bFTXXZwzhM0tZw<r+>p{Fxa`QAVSf~9!Nz$MJ
z7p#x|&ba@E|EzyuUa&64AB$+g`xk7I(a(o<zV7LQHENcx#S!B$$1YeuV+eEZ2(Oh`
z)zB#H4?YT``6uTuSW^Mz?^dBQNtlWsVJd5csYnr~VnmpV4`C`egvBEB^Xm5Jnd+Gr
ztb=%PqJ>$oj+dyNQx~kUu5#9b^#_g4Dv~K+Rz*HiS{un^XkM@;*%2*x@g?nCbCA7;
z9R4<_fpqHZs3A#<Lke(Rr!81xCRru7wFIl=?94RmExzk18Uf$^g7u*CETM2nx=;(&
znAhv~`Jt)zVUID7<R3?pf)FjdqD<gNaO6N5=<b=Kuy;`GCQJ4RBL!QqF8X}AnfP{S
zMVfxgR-}V?dv9B?Mx%0c>~^!1y6kOK9=RprGvjCp2Q@nPLTAgBl<Mq~b0cP;RJd|C
zW&!!@y_XBtEt=yAD^0`LbO^iFV%vS7^g8?g_B_Y)K3d&4iuY1{Z7q-I$Mt@HJxIH9
zpt<(bW9_t`tsS%r9U`?aI(ODSIvS-NIUcQT`7B0zVP{|Mq56Hbh$s7K{ZIDMMpo>r
zJ(k+n;J=d4SKBzax3+muFKzp&p4z@X-L+#MbTxi|CaQ~O&wsl&QoH2eh#eOUM+IZP
zF2FIu7a*^AtS~zQ`(ss6u7LNoDoTJIUl({y!AF800Ddg^+29v~Uj}|H_>JHXf&Uu(
zFW}uF?i)b=1b+oQ^3Q=k4gMJTec-o)$KSmQ9^<g*F%3)`^M(0Sz++x9&v1Nq4jg%%
z4IYldhK)MmO3rw|f>!|NjD!2&En~rt0Y4gi7WhoyXMj(K-=zT^31NrB^)T?s;Nd!F
z1n1+uRQL^s8wK<sxDUU>-{(lY0%>GJx|qI&=L_?QdBuFAqaBSdYWYj#;e?%}#&i6O
z<GopYGCrckW1q*3n>T-Ge|OL8_#Ld(S%#EaoCb1QpVLO1Qcok4<K~=(aoUDc{(xBd
z7%qa-NKU(ON)yO5(!`zb==w6EQ0HtQ$AdYI=QN4aWKL5!#m$1~u<wM9Ssg1%7K6%I
zPH|xk-NTG1Q#hT*>2yvX=XB|kdg`|EO<)&qQ`N6t7|=HAv(+rd=Qipkb$}hdoz$7@
zqoA(su4d~oK>Mk4=L`UToO=A1BnY3Pmgt!Qc$|86$HPFUsMhgQfj+AKZ~Y7iKU+;a
zGzajLYUhhj1755?du0jGXH})~3y{uB>dKp|Af9#Vp^h7XzfpCY{tocpgQN9z^?=82
zxIU`(c=IEmXVq&hzJ@Sx@O>h-0l$F%Yn?L&JFVT^jEBnMk9N=cu>YWrV@hczBPvg(
zeZW^Qj~s5E*_gz~j<mE?s38kZ;yk}$c}i|i^g+I$cSi<mp2jOp+w(1+E?{##qH2aa
zLy`XHQlwkVR<MD^9=ReFujI9l-g|K+FqCPWA0u~i!de&{vK9~=a+d~c!SxqT77f<6
z)@QlOHCTK0X*lYm_>Lo9OtH)cYt-furBQZ%Yp+`wti8)Wi<Zshaqns==rp^5^?xxq
z)xHea^LI#2N$C@pl9j03ACUhM?Ax|r(R1nvYgo{kf(;5}0<u$6t!~w5Cn}g@PFSO9
znYNI0@M~G)9UpT$h<W}CXJ!}%Ys0z9*$mdYe7`=_^2<lvouzQCJ((a~3TN+xwOobU
zELXT2`vc1q%wUb)uFc%c9p1gk(b|)VW%m-YW&0(2eF0E-w>ROHhEjg%oaKh1TI7W~
zzZHr7Ld^`;BEDh@Xf{}*nvN)|@d{(Gwj7?0CZh*&Qy!>;3`aGeQw`Se{W{TLEmwW~
z_<W+B`fB+%`=51ut%5_VAM0;i==hi|5cBdW)4wrVOD|N!lE6u*yR{)xzueYx#kKdV
z!!4yGweASOz2Cmo{iwfGx8iWnj9Ux1m4I6axTSz`m}!nLSim6y#$l>Cd{Y586EF^E
z&EabcxQ>8vm~0LYbHI6UJptoz+8n-tfU(bQzK_FhGj1&4CIZF+YYtyrz$FBX@7El@
zlz{ISu%CcS3mEzrdEzT8;Bo>kFW?FSt|;J20<J9JDgv%5;A#T)7jS@pQOV8nsUhH+
z0<HzPqt-3TC@&~XaIEaGYy6M#nH{ecZ7i%>%#PMCh(`tygmL?(*?F4*T$3>bvBjpS
z7~SC1SHP%4=qTS2-5Ew*LMQT3r_hOfE{=SB)G>5oc+@p?A|G`QoybSsLnrc42hoXq
z)J1e6A9WI)$Vc5oC-PB8(TRN2RdgaBbrzk-N8Lpy@==G;iG0*$bRr*h8lA`wFu3{o
zMn!HXd<{u>)OmDb{HXirL_XG`=tMr&rRYRHO+AWytXt8Ee5_;9iF~YU(TRMlbJ2->
zn(q+#H1#R+u`Wg@^07`vC-Q?BJ0Bn0Fmz&gY|GGze42U{`PjCh6ZzQ2p%eKubu03*
z4Mr#Ou`Na?@@eW<<YU{8PUK@B0G-IEsbi6k<1=(ruXX;DVH{nd6ZtgtEb=Wz==|YT
z3EwK=YZ89Agill7VtzUZ;lqq`x1s%z@FOMsP7;1+3BQYkAI14Zx=O-DQvgo8379tz
zHf|5w6}#vu;1~h-5^!$;_YrVk0rwMde*r%r;0FafK)?e9JV?N?0v;^jH~|k4aJ+yM
z1e_?~BmoZ<@Gt=<3wXGIM+i7Yz#|2mD&RDZ#d14J!hcA@PnYmBB>YSXKTE<NE#Z%m
z@W)E{<0Sm?68;1Uf1-pxNy2|v!k;YRXG{1~B>YDt{HYTDGztGv3I8z(f4YP}L&Bda
z;Xf|n&yw(GOZZPn_;V!uxf1?7&Zl~x>h~ul{P`08Qxg8u68-`Sf1!lGNWxz%;d8bU
z;6Ya|>XLv<0we*)&&y#C^^1N#{^t0Yv5#AH-}5af=DXkSNi%=gOJ1L`icK!(-|ktZ
zi<=>@&pR;sMYaKj^NByX$+*uk^0(RU8QraCSK97rTx;LXUx9h7d9uiQHl1$X?nxs*
zqO^dRf1+Nt-Lv(5OiVdk5HyAE)a{<0kQiH&dyj9<rLI5Y*1pbiLk+Aj(=>70xodZ}
z-E$PTa;?OlmgvW47RI}_dvc-tM~Q2@=YRCY^T)=PifX^V$97M8A9uOkb6w3Dx7+Ud
za*^|$+WQ*YJ?XdvO#G9O)a}A{Pnsi>Z}*(@KbCL1+dVh&nmzaPw05?8PVi=J(k-@o
z?((4^MYxp8C8mx}PESi6nV6axml8WFJuxFAD?PDlaSG_g0~WQNPdtaW2I6*4T6@g9
z-E&Yqipht^gxfv02Qw~b{i6x%UzoRhQon*|!TT4iA$JRB(xI6;Dql|TUHDoYVa4E>
z?Vi+ejWCpB!-5lTSk;gX!@EV8h{%}6`P)6IfRe6KnIufbk1&-r!c?RPQ!ye;#fLBz
z9KvFe@vr6Yy>!VoeagmdI_yNz1E1ZlCnde7CpF0dyhH!+&<=g+5?wzuXr~?;x=Wwi
zU^mb``il+r0^O(gi`ftN59nv_KLqr!zM|Ms`2G9(%ojg^Zyp`fU%2lCggdE!)9yo{
zAL-%MKZbZe(FdIV1j2u+FVju|J*A&`_!RI@>294rh2MXokC^%~@IM0HM|#T{AHwws
zeayyVM%*#lVc;Lo{i62*-L3EWt+nx8qH(7}jBgZ~zo&%XFG7#sQxZLW+X%mFgulh_
zGNH%#(4V{6&hUSDTdn9z`$ZVOVdY3G`VO_aSVMjex1ReU+<N>iKVv7+g@yjct|HuB
zWZV6OJC4x94yb&db|c}gDflh!q(YDJp+7jhgyDDk`)KIjJyFE)8)v#{==)9j(>he~
z(kd45g0a~?mF);Zt8nsYtdBc}%7Dk+L$r$tcN3vUjC{Or-dRMyza8Z}kc_l3{fTTx
zkG-8pxbp|SV#gyt8~Lys2!ETyJ;sOL!hJTz4?c%`WfI&6Z<!RHwQXW}6Yz$v|J!yE
zEuL#Jo{Q)*j)|@|@UmAI8T8x=M&T#k?05fR>JCz%b^xb!%^0VD&G-oc<0#mCAE%ei
z7^k_-7^k<*7^k((7^kz%7^nHoI9tF|1pJ7AaT?$J`)LBkX?yd1oUS)xT;?(3Q3A$g
z9rJx$&M{+L#xY}DzA<B5wlQN|t}$a=rZHn&o-t!wmN8>ojxl3ghB0GYelcTQb}?gG
zZUHtKF7gAua!qP{@g}&Ej3|DbmgTepr<FOa#%Xm<YjGOLX?;!`aY`M!P>!2(N=xpv
zWZs_ht(>wx4D3JUxC^H=e?uKTT56>w)CYLDft(KJlnCp3nvZ*+Cm~vx;~siSx}8ov
z?g0lP&N}X)%sU-ZBpWC&Sw)FC-}@qyq4~H6k{!{4ci*XZ`bSeYBYrH?=%@ql)Z-ou
zNm?Azfv<~2&ne$0z)Z41&MT2jMx^7W9Z(%+t3_!B+;PV};JRAFmr}Z?o5s0sa+T*s
zfqvUpu|>2=ZCU2vCx*K|;MEoR=aET4h!$Q^CN*Q}+;c@?@8C}hOa2Hwr$1fXf)429
zI~d;%i=s?>76q-sdsW4s+MUH%$%ArXABb#aoZ#c2LFc{D*&?0T!`ji=c^8juL8*Av
zW6#ej7)6mxu=sXlF?m&C!SheEOyX&I8I6`W*>Z!g*Lam@zxPH0&J(9;Xl48HIR{HH
zHvVt$#<A(nsO;RDvq!n-I{Xh(oCmu~6(@Pv|H8~O8)3Zd$l-6UP}@B{;YST0dg$BS
z^4Xa$=wa=Lu29i2C2lz)t||L5`|sCu(L`q`(xWa#x<zk=Rkv@+6=}<>b3Uf`UR?RO
zKhrorM($v{XK+xXppbyzmUp?`bHO(yYKz-FNA_U3%H6>7a<uxbQgm9qbf#G5?Vi=S
z?h&O{HveA9TW$9&+JYs?C&}+&mn>Y~CCmQ^_IHM%1f06vGq+n-GMdTBP7}A%J!}>#
zm}9ql&Vy&W*tg;`$H$C@T|Z6xh3pjpQL77=hMO6A{2ON}-2E+G3Mc<gS+2q*%M~u?
z({9-X+wNJJS187e0lr$CIWRd|dor=?UPAk6`z3pQ0nl4b>%*?$aF6jNUpExhA}`eW
ztw`(_YPQ{T-w~#O=Ix%<dESV!8m}PRJ!!Hm3SywaXvHv=`%?VLXfW#3?Vk94ow(gI
zSABNM)#s_EeSUO&t%7EuWv_-h@T`A&X=qsLSD|4S$LASMQ`6?}f5P!GqjBpCUN5ev
zm>p-ZhUT%sA3ZY&V<uI;$sIM~zK1}|V4|X=5`;%7^o-~b(K8}ii-?Wt6w^y-4hmli
zVo8Y0jLS&ORE%#8<3K3^#z?6d$!V#IH{Qoiv!VnjSWNFTt`&u-I1ykQ!}4&$t*mfy
zsMsbHI74i4XjDhxX&gtm$>jx&ScPCVSfh%7D0PAIS2Z$&wX{JjpVN3kR0((}BvX<4
zsu<J*<mzi2Q^7wwshK`jndf_KAUu`?){M;MmeN=Ga5b9W$vbyo>uVI{w;AzyzQ&$A
zGg$MU!eXez<0E?d*BP5#25WhQ{GS6ZgSG$ME6?-5=JupsYgSd++LP1pb$N%jCp0fd
z^j^0oCA(8)%!Dv+%7~SQ!5SHSBRk;t?C%AGHJU<{8?4cj+;Xb-mIiCocHY8ZjhfTj
zexHEaS$p3npf*=Fhnm-5jU?vgNrMJ7SWD}dSB`m9>K_JJ|H5prR-Hcq(ZcsHe*UZa
z7Y>>Ph5?h~R+@09)|KM5jZZF8*EL^Nb4gK3IJByv&J$58proD7?z&RJq4dNeV?Fm3
z+AO%7O9Czl{FM?&O6sFFQwG9$wDIcOhlas9wW;d-QJH`zs2OnnEuMY*aK)J_o_*VQ
z#2ldW)RJ$`gEMaDsb7OfdfoF$<4oN4bLYeTr_~wy0-%f4Qm-xr{&Q;CN0&i-FR77n
zD}b(2r(9eO>8w@%xB3<1tlXreGy0xC&wzF4S>318IbBgc*DG)O9IQmY(9e0C2mX1z
zedw29iTb6!XvkM^|7(5GD__I?Z}f64z6Jc9{>SR?z%uoMUSq`tz!&w&y)MG<F6t9c
kTm*bcUsdrE+`pvH8*&Nwm-G(~T>}1Py-Db0i088Ye<<WOP5=M^

literal 0
HcmV?d00001

diff --git a/climada/hazard/test/test_tc_tracks.py b/climada/hazard/test/test_tc_tracks.py
index b9eea67f8..580884fa9 100644
--- a/climada/hazard/test/test_tc_tracks.py
+++ b/climada/hazard/test/test_tc_tracks.py
@@ -44,6 +44,7 @@
 TEST_TRACK_CHAZ = DATA_DIR.joinpath('chaz_test_tracks.nc')
 TEST_TRACK_STORM = DATA_DIR.joinpath('storm_test_tracks.txt')
 TEST_TRACKS_ANTIMERIDIAN = DATA_DIR.joinpath('tracks-antimeridian')
+TEST_TRACKS_LEGACY_HDF5 = DATA_DIR.joinpath('tctracks_hdf5_legacy.nc')
 
 
 class TestIbtracs(unittest.TestCase):
@@ -339,6 +340,11 @@ def test_hdf5_io(self):
                 np.testing.assert_array_equal(tr[v].values, tr_read[v].values)
             self.assertEqual(tr.sid, tr_read.sid)
 
+        # attempting to read the legacy file format should fail gracefully
+        with self.assertRaises(ValueError) as cm:
+            tc.TCTracks.from_hdf5(TEST_TRACKS_LEGACY_HDF5)
+        self.assertIn("no longer supported by CLIMADA", str(cm.exception))
+
     def test_from_processed_ibtracs_csv(self):
         tc_track = tc.TCTracks.from_processed_ibtracs_csv(TEST_TRACK)