Feature/change impact total value (#702)

chahank · Chahan Kropf · peanutfun · web-flow · commit 7a646c88bd17 · 2023-05-15T12:53:31.000+02:00
* Add first try affected value with threshold

* Remove test typos

* Update test affected value

* Fix affected value to correct results

* Ignore centroids -1 and exposure values &lt;=0

* Improve threshold affected value

Co-authored-by: Lukas Riedel &lt;34276446+peanutfun@users.noreply.github.com&gt;

* Deprecate Impact.tot_value attribute

* Replace attribute with getter/setter and declare private attribute.
* Adjust writers accordningly

* Code line cosmetics

* Add private function for legacy during deprecation period

* Update climada/entity/exposures/base.py

* Update climada/engine/impact.py

* Update climada/entity/exposures/base.py

* Format new code and improve docstrings

* Update CHANGELOG.md

* Update changelog for temporary method.

* Clarify docstring.

---------

Co-authored-by: Chahan Kropf &lt;chahan.kropf@posteo.com&gt;
Co-authored-by: Lukas Riedel &lt;34276446+peanutfun@users.noreply.github.com&gt;
Co-authored-by: emanuel-schmid &lt;schmide@ethz.ch&gt;
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -22,6 +22,8 @@ Removed:
 - `Impact.match_centroids` convenience method for matching (hazard) centroids to impact objects [#602](https://github.com/CLIMADA-project/climada_python/pull/602)
 - `climada.util.coordinates.match_centroids` method for matching (hazard) centroids to GeoDataFrames [#602](https://github.com/CLIMADA-project/climada_python/pull/602)
 - 'Extra' requirements `doc`, `test`, and `dev` for Python package [#712](https://github.com/CLIMADA-project/climada_python/pull/712)
+- Added method `Exposures.centroids_total_value` to replace the functionality of `Exposures.affected_total_value`. This method is temporary and deprecated. [#702](https://github.com/CLIMADA-project/climada_python/pull/702)
+
 
 ### Changed
 
@@ -36,6 +38,7 @@ Removed:
 - Tests with long runtime were moved to integration tests in `climada/test` [#709](https://github.com/CLIMADA-project/climada_python/pull/709)
 - Use `myst-nb` for parsing Jupyter Notebooks for the documentation instead of `nbsphinx` [#712](https://github.com/CLIMADA-project/climada_python/pull/712)
 - Installation guide now recommends installing CLIMADA directly via `conda install` [#714](https://github.com/CLIMADA-project/climada_python/pull/714)
+- `Exposures.affected_total_value` now takes a hazard intensity threshold as argument. Affected values are only those for which at least one event exceeds the threshold. (previously, all exposures points with an assigned centroid were considered affected) [#702](https://github.com/CLIMADA-project/climada_python/pull/702)
 
 ### Fixed
 
@@ -44,6 +47,7 @@ Removed:
 ### Deprecated
 
 - `Centroids.from_geodataframe` and `Centroids.from_pix_bounds` [#721](https://github.com/CLIMADA-project/climada_python/pull/721)
+- `Impact.tot_value`: Use `Exposures.affected_total_value` to compute the total value affected by a hazard intensity above a custom threshold [#702](https://github.com/CLIMADA-project/climada_python/pull/702)
 
 ### Removed
 
diff --git a/climada/engine/impact.py b/climada/engine/impact.py
@@ -80,8 +80,6 @@ class Impact():
         frequency of event
     frequency_unit : str
         frequency unit used (given by hazard), default is '1/year'
-    tot_value : float
-        total exposure value affected
     aai_agg : float
         average impact within a period of 1/frequency_unit (aggregated)
     unit : str
@@ -154,7 +152,7 @@ def __init__(self,
         self.at_event = np.array([], float) if at_event is None else at_event
         self.frequency = np.array([],float) if frequency is None else frequency
         self.frequency_unit = frequency_unit
-        self.tot_value = tot_value
+        self._tot_value = tot_value
         self.aai_agg = aai_agg
         self.unit = unit
 
@@ -246,7 +244,7 @@ def from_eih(cls, exposures, impfset, hazard,
                                  axis=1),
             crs = exposures.crs,
             unit = exposures.value_unit,
-            tot_value = exposures.affected_total_value(hazard),
+            tot_value = exposures.centroids_total_value(hazard),
             eai_exp = eai_exp,
             at_event = at_event,
             aai_agg = aai_agg,
@@ -257,6 +255,29 @@ def from_eih(cls, exposures, impfset, hazard,
                    }
             )
 
+    @property
+    def tot_value(self):
+        """Return the total exposure value close to a hazard
+
+        .. deprecated:: 3.3
+           Use :py:meth:`climada.entity.exposures.base.Exposures.affected_total_value`
+           instead.
+        """
+        LOGGER.warning("The Impact.tot_value attribute is deprecated."
+                       "Use Exposures.affected_total_value to calculate the affected "
+                       "total exposure value based on a specific hazard intensity "
+                       "threshold")
+        return self._tot_value
+
+    @tot_value.setter
+    def tot_value(self, value):
+        """Set the total exposure value close to a hazard"""
+        LOGGER.warning("The Impact.tot_value attribute is deprecated."
+                       "Use Exposures.affected_total_value to calculate the affected "
+                       "total exposure value based on a specific hazard intensity "
+                       "threshold")
+        self._tot_value = value
+
     def transfer_risk(self, attachment, cover):
         """Compute the risk transfer for the full portfolio. This is the risk
         of the full portfolio summed over all events. For each
@@ -860,7 +881,7 @@ def write_csv(self, file_name):
                          [self.tag['haz'].description]],
                         [[self.tag['exp'].file_name], [self.tag['exp'].description]],
                         [[self.tag['impf_set'].file_name], [self.tag['impf_set'].description]],
-                        [self.unit], [self.tot_value], [self.aai_agg],
+                        [self.unit], [self._tot_value], [self.aai_agg],
                         self.event_id, self.event_name, self.date,
                         self.frequency, [self.frequency_unit], self.at_event,
                         self.eai_exp, self.coord_exp[:, 0], self.coord_exp[:, 1],
@@ -901,7 +922,7 @@ def write_col(i_col, imp_ws, xls_data):
         data = [str(self.tag['impf_set'].file_name), str(self.tag['impf_set'].description)]
         write_col(2, imp_ws, data)
         write_col(3, imp_ws, [self.unit])
-        write_col(4, imp_ws, [self.tot_value])
+        write_col(4, imp_ws, [self._tot_value])
         write_col(5, imp_ws, [self.aai_agg])
         write_col(6, imp_ws, self.event_id)
         write_col(7, imp_ws, self.event_name)
@@ -1030,8 +1051,9 @@ def write_csr(group, name, value):
         with h5py.File(file_path, "w") as file:
 
             # Now write all attributes
+            # NOTE: Remove leading underscore to write '_tot_value' as regular attribute
             for name, value in self.__dict__.items():
-                write(file, name, value)
+                write(file, name.lstrip("_"), value)
 
     def write_sparse_csr(self, file_name):
         """Write imp_mat matrix in numpy's npz format."""
diff --git a/climada/entity/exposures/base.py b/climada/entity/exposures/base.py
@@ -36,6 +36,7 @@
 import contextily as ctx
 import cartopy.crs as ccrs
 
+from climada.hazard import Hazard
 from climada.entity.tag import Tag
 import climada.util.hdf5_handler as u_hdf5
 from climada.util.constants import ONE_LAT_KM, DEF_CRS, CMAP_RASTER
@@ -1019,11 +1020,15 @@ def concat(exposures_list):
 
         return exp
 
-    def affected_total_value(self, hazard):
-        """
-        Total value of the exposures that are close enough to be affected
-        by the hazard (sum of value of all exposures points for which
-        a centroids is assigned)
+    def centroids_total_value(self, hazard):
+        """Compute value of exposures close enough to be affected by hazard
+
+        .. deprecated:: 3.3
+           This method will be removed in a future version. Use
+           :py:meth:`affected_total_value` instead.
+
+        This method computes the sum of the value of all exposures points for which a
+        Hazard centroid is assigned.
 
         Parameters
         ----------
@@ -1043,6 +1048,39 @@ def affected_total_value(self, hazard):
         )
         return np.sum(self.gdf.value.values[nz_mask])
 
+    def affected_total_value(self, hazard: Hazard, threshold_affected: float = 0):
+        """
+        Total value of the exposures that are affected by at least
+        one hazard event (sum of value of all exposures points for which
+        at least one event has intensity larger than the threshold).
+
+        Parameters
+        ----------
+        hazard : Hazard
+           Hazard affecting Exposures
+        threshold_affected : int or float
+            Hazard intensity threshold above which an exposures is
+            considere affected.
+
+        Returns
+        -------
+        float
+            Sum of value of all exposures points for which
+            a centroids is assigned and that have at least one
+            event intensity above threshold.
+
+        """
+        assigned_centroids = self.gdf[hazard.centr_exp_col]
+        nz_mask = (self.gdf.value.values > 0) & (assigned_centroids.values >= 0)
+        cents = np.unique(assigned_centroids[nz_mask])
+        cent_with_inten_above_thres = (
+            hazard.intensity[:, cents].max(axis=0) > threshold_affected
+        ).nonzero()[1]
+        above_thres_mask = np.isin(
+            self.gdf[hazard.centr_exp_col].values, cents[cent_with_inten_above_thres]
+        )
+        return np.sum(self.gdf.value.values[above_thres_mask])
+
 
 def add_sea(exposures, sea_res, scheduler=None):
     """Add sea to geometry's surroundings with given resolution. region_id
diff --git a/climada/entity/exposures/test/test_base.py b/climada/entity/exposures/test/test_base.py
@@ -25,6 +25,7 @@
 from sklearn.metrics import DistanceMetric
 import rasterio
 from rasterio.windows import Window
+import scipy as sp
 
 from climada import CONFIG
 from climada.entity.exposures.base import Exposures, INDICATOR_IMPF, \
@@ -185,23 +186,29 @@ def test_assign_large_hazard_subset_pass(self):
         np.testing.assert_array_equal(assigned_centroids.lon, exp.gdf.longitude)
 
     def test_affected_total_value(self):
-        exp = Exposures.from_raster(HAZ_DEMO_FL, window=Window(25, 90, 10, 5))
-        haz = Hazard.from_raster([HAZ_DEMO_FL], haz_type='FL', window=Window(25, 90, 10, 5))
-        exp.assign_centroids(haz)
-        tot_val = exp.affected_total_value(haz)
-        self.assertEqual(tot_val, np.sum(exp.gdf.value))
-        new_centr = exp.gdf.centr_FL
-        new_centr[6] = -1
-        exp.gdf.centr_FL = new_centr
-        tot_val = exp.affected_total_value(haz)
-        self.assertAlmostEqual(tot_val, np.sum(exp.gdf.value) - exp.gdf.value[6], places=4)
-        new_vals = exp.gdf.value
-        new_vals[7] = 0
-        exp.gdf.value = new_vals
-        tot_val = exp.affected_total_value(haz)
-        self.assertAlmostEqual(tot_val, np.sum(exp.gdf.value) - exp.gdf.value[6], places=4)
-        exp.gdf.centr_FL = -1
-        tot_val = exp.affected_total_value(haz)
+        haz_type = "RF"
+        gdf = gpd.GeoDataFrame(
+            {
+                "value": [1, 2, 3, 4, 5, 6],
+                "latitude": [1, 2, 3, 4, 5, 6],
+                "longitude": [-1, -2, -3, -4, -5, -6],
+                "centr_" + haz_type: [0, 2, 2, 3, -1, 4],
+            }
+        )
+        exp = Exposures(gdf, crs=4326)
+        intensity = sp.sparse.csr_matrix(np.array([[0, 0, 1, 10, 2], [-1, 0, 0, 1, 2]]))
+        cent = Centroids(lat=np.array([1, 2, 3, 4]), lon=np.array([1, 2, 3, 4]))
+        haz = Hazard(
+            haz_type=haz_type, centroids=cent, intensity=intensity, event_id=[1, 2]
+        )
+
+        tot_val = exp.affected_total_value(haz, threshold_affected=0)
+        self.assertEqual(tot_val, np.sum(exp.gdf.value[[1, 2, 3, 5]]))
+        tot_val = exp.affected_total_value(haz, threshold_affected=3)
+        self.assertEqual(tot_val, np.sum(exp.gdf.value[[3]]))
+        tot_val = exp.affected_total_value(haz, threshold_affected=-2)
+        self.assertEqual(tot_val, np.sum(exp.gdf.value[[0, 1, 2, 3, 5]]))
+        tot_val = exp.affected_total_value(haz, threshold_affected=11)
         self.assertEqual(tot_val, 0)
 
 class TestChecker(unittest.TestCase):
