Read legacy Exposures hdf5 metadata with affine 3

AUTHORS.md

@@ -39,3 +39,4 @@
 * Dahyann Araya
 * Giovanni Cozzolongo
 * Thomas Struys
+* rmz-oz

CHANGELOG.md

@@ -26,6 +26,7 @@ Code freeze date: YYYY-MM-DD
 
 ### Fixed
 
+- `Exposures.from_hdf5` reads files written by CLIMADA < 6.1 also with affine 3.x, whose `Affine` can no longer be restored from the pickled metadata by PyTables [#1311](https://github.com/CLIMADA-project/climada_python/issues/1311)
 - Preserve explicitly mapped dates in `HazardForecast.from_xarray_raster` [#1305](https://github.com/CLIMADA-project/climada_python/issues/1305).
 - Fixed asset count in impact logging message [#1195](https://github.com/CLIMADA-project/climada_python/pull/1195).
 - `Hazard.from_raster_xarray` now returns a sparse matrix instead of a sparse array [#1261](https://github.com/CLIMADA-project/climada_python/pull/1261).

climada/entity/exposures/base.py

@@ -23,7 +23,10 @@
 
 
 import copy
+import copyreg
+import io
 import logging
+import pickle
 import warnings
 from pathlib import Path
 
@@ -34,6 +37,7 @@
 import numpy as np
 import pandas as pd
 import rasterio
+from affine import Affine
 from deprecation import deprecated
 from geopandas import GeoDataFrame, GeoSeries, points_from_xy
 from mpl_toolkits.axes_grid1 import make_axes_locatable
@@ -83,6 +87,42 @@
 """MATLAB variable names"""
 
 
+def _reconstruct_legacy(cls, base, state):
+    """Replacement for copyreg._reconstructor when reading legacy metadata.
+
+    Up to affine 2.4, Affine was a subclass of tuple, and a protocol 0 pickle
+    of it (as written by PyTables) is restored with
+    ``copyreg._reconstructor(Affine, tuple, state)``. From affine 3.0 on,
+    Affine is no longer a tuple, so this raises a TypeError.
+    """
+    if base is tuple and issubclass(cls, Affine):
+        return cls(*state[:6])
+    # pylint: disable-next=protected-access
+    return copyreg._reconstructor(cls, base, state)
+
+
+class _LegacyMetadataUnpickler(pickle.Unpickler):
+    """Unpickler for Exposures metadata that PyTables could not restore."""
+
+    def find_class(self, module, name):
+        found = super().find_class(module, name)
+        if found is copyreg._reconstructor:  # pylint: disable=protected-access
+            return _reconstruct_legacy
+        return found
+
+
+def _unpickle_legacy_metadata(metadata):
+    """Return the metadata dict of an Exposures hdf5 file.
+
+    PyTables returns the raw pickle bytes when unpickling an attribute fails.
+    This happens for files written by CLIMADA < 6.1, whose metadata contains
+    an affine.Affine transform, when they are read with affine >= 3.0.
+    """
+    if isinstance(metadata, bytes):
+        return _LegacyMetadataUnpickler(io.BytesIO(metadata)).load()
+    return metadata
+
+
 class Exposures:
     """geopandas GeoDataFrame with metadata and columns (pd.Series) defined in
     Attributes.
@@ -1206,7 +1246,9 @@ def from_hdf5(cls, file_name):
         if not Path(file_name).is_file():
             raise FileNotFoundError(str(file_name))
         with pd.HDFStore(file_name, mode="r") as store:
-            metadata = store.get_storer("exposures").attrs.metadata
+            metadata = _unpickle_legacy_metadata(
+                store.get_storer("exposures").attrs.metadata
+            )
             # in previous versions of CLIMADA and/or geopandas, the CRS was stored in '_crs'/'crs'
             crs = metadata.get("crs", metadata.get("_crs"))
             if crs is None and metadata.get("meta"):

climada/entity/exposures/test/test_base.py

@@ -26,6 +26,7 @@
 import pandas as pd
 import rasterio
 import scipy as sp
+from affine import Affine
 from rasterio.windows import Window
 from shapely.geometry import MultiPolygon, Point, Polygon
 from sklearn.metrics import DistanceMetric
@@ -39,6 +40,7 @@
     INDICATOR_CENTR,
     INDICATOR_IMPF,
     Exposures,
+    _unpickle_legacy_metadata,
     add_sea,
 )
 from climada.hazard.base import Centroids, Hazard
@@ -504,6 +506,54 @@ def test_io_hdf5_pass(self):
             exp.data["geocol2"].geometry, exp_read.data["geocol2"].values
         )
 
+    # metadata as pickled by PyTables in files written by CLIMADA < 6.1, with
+    # affine 2.x (where Affine was a subclass of tuple)
+    LEGACY_METADATA = (
+        b"(dp0\nVdescription\np1\nVlegacy\np2\nsVref_year\np3\nI2018\n"
+        b"sVvalue_unit\np4\nVUSD\np5\nsVcrs\np6\nVEPSG:4326\np7\nsVmeta\np8\n"
+        b"(dp9\nVtransform\np10\nccopy_reg\n_reconstructor\np11\n"
+        b"(caffine\nAffine\np12\nc__builtin__\ntuple\np13\n"
+        b"(F0.5\nF0.0\nF19.75\nF0.0\nF-0.5\nF10.75\nF0.0\nF0.0\nF1.0\n"
+        b"tp14\ntp15\nRp16\nss."
+    )
+
+    def test_unpickle_legacy_metadata(self):
+        """metadata with an affine 2.x Affine can be read with any affine version"""
+        metadata = _unpickle_legacy_metadata(np.bytes_(self.LEGACY_METADATA))
+        self.assertEqual(metadata["description"], "legacy")
+        self.assertEqual(
+            metadata["meta"]["transform"], Affine(0.5, 0.0, 19.75, 0.0, -0.5, 10.75)
+        )
+        # metadata that PyTables could unpickle is passed through
+        self.assertIs(_unpickle_legacy_metadata(metadata), metadata)
+
+    def test_read_legacy_hdf5_pass(self):
+        """read an hdf5 file whose metadata holds an affine 2.x Affine"""
+        file_name = DATA_DIR.joinpath("test_hdf5_exp_legacy.h5")
+        with pd.HDFStore(file_name, mode="w") as store:
+            store.put(
+                "exposures",
+                pd.DataFrame(
+                    {
+                        "value": [1.0, 2.0],
+                        "latitude": [10.0, 10.5],
+                        "longitude": [20.0, 20.5],
+                    }
+                ),
+            )
+            store.get_storer("exposures").attrs.metadata = np.bytes_(
+                self.LEGACY_METADATA
+            )
+
+        exp_read = Exposures.from_hdf5(file_name)
+        file_name.unlink()
+
+        self.assertEqual(exp_read.description, "legacy")
+        self.assertEqual(exp_read.ref_year, 2018)
+        self.assertEqual(exp_read.value_unit, "USD")
+        self.assertTrue(u_coord.equal_crs(exp_read.crs, "EPSG:4326"))
+        np.testing.assert_array_equal(exp_read.value, [1.0, 2.0])
+
 
 class TestAddSea(unittest.TestCase):
     """Check constructor Exposures through DataFrames readers"""