diff --git a/src/test/data/hazard_model_store.py b/src/test/data/hazard_model_store.py index 0e253fc6..a377b66a 100644 --- 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 new file mode 100644 index 00000000..65538b06 --- /dev/null +++ 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 new file mode 100644 index 00000000..65538b06 --- /dev/null +++ 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 index 4909331a..870337ee 100644 --- 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(