Tidy and points-based hazard model

Signed-off-by: Joe Moorhouse <[email protected]>
os-climate · Oct 19, 2023 · 6de16a2 · 6de16a2
1 parent 3cc9e5b
commit 6de16a2
Show file tree

Hide file tree

Showing 4 changed files with 162 additions and 78 deletions.
diff --git a/src/test/data/hazard_model_store.py b/src/test/data/hazard_model_store.py
@@ -10,43 +10,13 @@
 
 
 class TestData:
-    longitudes = [
-        69.4787,
-        68.71,
-        20.1047,
-        19.8936,
-        19.6359,
-        0.5407,
-        6.9366,
-        6.935,
-        13.7319,
-        13.7319,
-        14.4809,
-        -68.3556,
-        -68.3556,
-        -68.9892,
-        -70.9157,
-    ]
-    latitudes = [
-        34.556,
-        35.9416,
-        39.9116,
-        41.6796,
-        42.0137,
-        35.7835,
-        36.8789,
-        36.88,
-        -12.4706,
-        -12.4706,
-        -9.7523,
-        -38.9368,
-        -38.9368,
-        -34.5792,
-        -39.2145,
-    ]
+    # fmt: off
+    longitudes = [ 69.4787, 68.71, 20.1047, 19.8936, 19.6359, 0.5407, 6.9366, 6.935, 13.7319, 13.7319, 14.4809, -68.3556, -68.3556, -68.9892, -70.9157 ]
+    latitudes = [ 34.556, 35.9416, 39.9116, 41.6796, 42.0137, 35.7835, 36.8789, 36.88, -12.4706, -12.4706, -9.7523, -38.9368, -38.9368, -34.5792, -39.2145 ]
 
     coastal_longitudes = [12.2, 50.5919, 90.3473, 90.4295, 90.4804, 90.3429, 90.5153, 90.6007]
     coastal_latitudes = [-5.55, 26.1981, 23.6473, 23.6783, 23.5699, 23.9904, 23.59, 23.6112]
+    # fmt: on
 
 
 def get_mock_hazard_model_store_single_curve():

diff --git a/src/test/kernel/test_hazard_model.py b/src/test/kernel/test_hazard_model.py
@@ -0,0 +1,78 @@
+from dataclasses import dataclass
+from typing import Dict, List, Mapping, NamedTuple, Sequence, Tuple
+
+import numpy as np
+from physrisk.kernel.assets import RealEstateAsset
+from physrisk.kernel.hazard_model import HazardDataRequest, HazardDataResponse, HazardEventDataResponse, HazardModel, HazardParameterDataResponse
+from physrisk.kernel.hazards import ChronicHeat, Wind
+from physrisk.kernel.impact import calculate_impacts
+from physrisk.vulnerability_models.real_estate_models import GenericTropicalCycloneModel
+
+import test.data.hazard_model_store as hms
+
+
+@dataclass
+class SinglePointData:
+    latitude: float
+    longitude: float
+    scenario: str
+    year: int
+    wind_return_periods: np.ndarray # years
+    wind_intensities: np.ndarray # m/s
+    chronic_heat_intensity: float # days over 35C
+    # etc
+
+class PointsKey(NamedTuple):
+    latitude: float
+    longitude: float
+    scenario: str
+    year: str
+
+class PointBasedHazardModel(HazardModel):
+    def __init__(self, points: Sequence[SinglePointData]):
+        """HazardModel suitable for storing relatively small number (<~ million say) of individual hazard
+        data points. 
+
+        Args:
+            points (Sequence[SinglePointData]): List of points.
+        """
+        self.points: Dict[Tuple[PointsKey, float, float], SinglePointData] = \
+            { self._get_key(p.latitude, p.longitude, p.scenario, p.year): p for p in points }
+
+    def _get_key(self, latitude: float, longitude: float, scenario: str, year: int) -> Tuple[float, float]:
+        return PointsKey(latitude=round(latitude, 3), longitude=round(longitude, 3), scenario=scenario, year=year)
+
+    def get_hazard_events(self, requests: List[HazardDataRequest]) -> Mapping[HazardDataRequest, HazardDataResponse]:
+        response = {}
+        for request in requests:
+            point = self.points[self._get_key(request.latitude, request.longitude, request.scenario, request.year)]
+            if request.hazard_type == Wind and request.indicator_id == "max_speed":
+                response[request] = HazardEventDataResponse(return_periods=point.wind_return_periods, intensities=point.wind_intensities)
+            elif request.hazard_type == ChronicHeat and request.indicator_id == "days/above/35c":
+                response[request] = HazardParameterDataResponse(point.chronic_heat_intensity)
+            # etc
+        return response
+
+
+def test_using_point_based_hazard_model():
+    # test that shows how data already present for a number of points can be used in a HazardModel
+    scenario = "rcp8p5"
+    year = 2080
+    assets = [
+        RealEstateAsset(lat, lon, location="Asia", type="Buildings/Industrial")
+        for lon, lat in zip(hms.TestData.longitudes[0:1], hms.TestData.latitudes[0:1])
+    ]
+    # fmt: off
+    wind_return_periods = np.array([10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 200.0, 300.0, 400.0, 500.0, 600.0, 700.0, 800.0, 900.0, 1000.0])
+    wind_intensities = np.array([37.279999, 44.756248, 48.712502, 51.685001, 53.520000, 55.230000, 56.302502, 57.336250, 58.452499, 59.283749, 63.312500, 65.482498, 66.352501, 67.220001, 67.767502, 68.117500, 68.372498, 69.127502, 70.897499 ])
+    # fmt: on
+    point = SinglePointData(hms.TestData.latitudes[0], hms.TestData.longitudes[0], scenario=scenario, year=year,
+                            wind_return_periods=wind_return_periods, wind_intensities=wind_intensities,
+                            chronic_heat_intensity=0)
+
+    hazard_model = PointBasedHazardModel([point])
+    vulnerability_models = {RealEstateAsset: [GenericTropicalCycloneModel()]}
+    results = calculate_impacts(assets, hazard_model, vulnerability_models, scenario=scenario, year=year)
+    impact_distrib = results[(assets[0], Wind)].impact
+    mean_impact = impact_distrib.mean_impact()
+    np.testing.assert_almost_equal(mean_impact, 0.009909858317497338)
diff --git a/src/test/kernel/test_hazard_models.py b/src/test/kernel/test_hazard_models.py
@@ -0,0 +1,78 @@
+from dataclasses import dataclass
+from typing import Dict, List, Mapping, NamedTuple, Sequence, Tuple
+
+import numpy as np
+from physrisk.kernel.assets import RealEstateAsset
+from physrisk.kernel.hazard_model import HazardDataRequest, HazardDataResponse, HazardEventDataResponse, HazardModel, HazardParameterDataResponse
+from physrisk.kernel.hazards import ChronicHeat, Wind
+from physrisk.kernel.impact import calculate_impacts
+from physrisk.vulnerability_models.real_estate_models import GenericTropicalCycloneModel
+
+import test.data.hazard_model_store as hms
+
+
+@dataclass
+class SinglePointData:
+    latitude: float
+    longitude: float
+    scenario: str
+    year: int
+    wind_return_periods: np.ndarray # years
+    wind_intensities: np.ndarray # m/s
+    chronic_heat_intensity: float # days over 35C
+    # etc
+
+class PointsKey(NamedTuple):
+    latitude: float
+    longitude: float
+    scenario: str
+    year: str
+
+class PointBasedHazardModel(HazardModel):
+    def __init__(self, points: Sequence[SinglePointData]):
+        """HazardModel suitable for storing relatively small number (<~ million say) of individual hazard
+        data points. 
+
+        Args:
+            points (Sequence[SinglePointData]): List of points.
+        """
+        self.points: Dict[Tuple[PointsKey, float, float], SinglePointData] = \
+            { self._get_key(p.latitude, p.longitude, p.scenario, p.year): p for p in points }
+
+    def _get_key(self, latitude: float, longitude: float, scenario: str, year: int) -> Tuple[float, float]:
+        return PointsKey(latitude=round(latitude, 3), longitude=round(longitude, 3), scenario=scenario, year=year)
+
+    def get_hazard_events(self, requests: List[HazardDataRequest]) -> Mapping[HazardDataRequest, HazardDataResponse]:
+        response = {}
+        for request in requests:
+            point = self.points[self._get_key(request.latitude, request.longitude, request.scenario, request.year)]
+            if request.hazard_type == Wind and request.indicator_id == "max_speed":
+                response[request] = HazardEventDataResponse(return_periods=point.wind_return_periods, intensities=point.wind_intensities)
+            elif request.hazard_type == ChronicHeat and request.indicator_id == "days/above/35c":
+                response[request] = HazardParameterDataResponse(point.chronic_heat_intensity)
+            # etc
+        return response
+
+
+def test_using_point_based_hazard_model():
+    # test that shows how data already present for a number of points can be used in a HazardModel
+    scenario = "rcp8p5"
+    year = 2080
+    assets = [
+        RealEstateAsset(lat, lon, location="Asia", type="Buildings/Industrial")
+        for lon, lat in zip(hms.TestData.longitudes[0:1], hms.TestData.latitudes[0:1])
+    ]
+    # fmt: off
+    wind_return_periods = np.array([10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 200.0, 300.0, 400.0, 500.0, 600.0, 700.0, 800.0, 900.0, 1000.0])
+    wind_intensities = np.array([37.279999, 44.756248, 48.712502, 51.685001, 53.520000, 55.230000, 56.302502, 57.336250, 58.452499, 59.283749, 63.312500, 65.482498, 66.352501, 67.220001, 67.767502, 68.117500, 68.372498, 69.127502, 70.897499 ])
+    # fmt: on
+    point = SinglePointData(hms.TestData.latitudes[0], hms.TestData.longitudes[0], scenario=scenario, year=year,
+                            wind_return_periods=wind_return_periods, wind_intensities=wind_intensities,
+                            chronic_heat_intensity=0)
+
+    hazard_model = PointBasedHazardModel([point])
+    vulnerability_models = {RealEstateAsset: [GenericTropicalCycloneModel()]}
+    results = calculate_impacts(assets, hazard_model, vulnerability_models, scenario=scenario, year=year)
+    impact_distrib = results[(assets[0], Wind)].impact
+    mean_impact = impact_distrib.mean_impact()
+    np.testing.assert_almost_equal(mean_impact, 0.009909858317497338)
diff --git a/src/test/models/test_wind_models.py b/src/test/models/test_wind_models.py
@@ -15,50 +15,8 @@ def test_wind_real_estate_model():
     year = 2080
     # mock some IRIS data for the calculation:
     store, root = hms.zarr_memory_store()
-    return_periods = [
-        10.0,
-        20.0,
-        30.0,
-        40.0,
-        50.0,
-        60.0,
-        70.0,
-        80.0,
-        90.0,
-        100.0,
-        200.0,
-        300.0,
-        400.0,
-        500.0,
-        600.0,
-        700.0,
-        800.0,
-        900.0,
-        1000.0,
-    ]
-    intensity = np.array(
-        [
-            37.279999,
-            44.756248,
-            48.712502,
-            51.685001,
-            53.520000,
-            55.230000,
-            56.302502,
-            57.336250,
-            58.452499,
-            59.283749,
-            63.312500,
-            65.482498,
-            66.352501,
-            67.220001,
-            67.767502,
-            68.117500,
-            68.372498,
-            69.127502,
-            70.897499,
-        ]
-    )
+    return_periods = [10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 200.0, 300.0, 400.0, 500.0, 600.0, 700.0, 800.0, 900.0, 1000.0]
+    intensity = np.array([37.279999, 44.756248, 48.712502, 51.685001, 53.520000, 55.230000, 56.302502, 57.336250, 58.452499, 59.283749, 63.312500, 65.482498, 66.352501, 67.220001, 67.767502, 68.117500, 68.372498, 69.127502, 70.897499 ])
     shape, transform = hms.shape_transform_21600_43200(return_periods=return_periods)
     path = f"wind/iris/v1/max_speed_{scenario}_{year}".format(scenario=scenario, year=year)
     hms.add_curves(