diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index f71df736..84d0af62 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -3,7 +3,6 @@ repos: - repo: https://github.com/pre-commit/pre-commit-hooks rev: v4.6.0 hooks: - - id: trailing-whitespace - id: check-merge-conflict - id: end-of-file-fixer - id: name-tests-test @@ -15,28 +14,32 @@ repos: - id: check-yaml - id: check-symlinks - id: detect-private-key - - id: check-ast - - id: debug-statements - repo: https://github.com/Lucas-C/pre-commit-hooks rev: v1.5.5 hooks: - id: remove-tabs - - repo: https://github.com/psf/black - rev: '24.4.2' + - repo: https://github.com/astral-sh/ruff-pre-commit + rev: v0.5.7 hooks: - - id: black - - id: black-jupyter + - id: ruff + args: [--fix, --exit-non-zero-on-fix, --config=pyproject.toml] + - id: ruff-format - - repo: https://github.com/s-weigand/flake8-nb - rev: v0.5.3 + - repo: local hooks: - - id: flake8-nb - additional_dependencies: - - pep8-naming - # Ignore all format-related checks as Black takes care of those. - args: - - --ignore=E2, W5, F401, E401, E704 - - --select=E, W, F, N - - --max-line-length=120 + - id: mypy-cache + name: "create mypy cache" + language: system + pass_filenames: false + entry: bash -c 'if [ ! -d .mypy_cache ]; + then /bin/mkdir .mypy_cache; fi; exit 0' + + - repo: https://github.com/pre-commit/mirrors-mypy + rev: "v1.11.1" + hooks: + - id: mypy + verbose: true + args: ["--show-error-codes", "--install-types", "--non-interactive"] + additional_dependencies: ["pytest", "types-requests"] diff --git a/docs/conf.py b/docs/conf.py index e19e907f..f3c65a48 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -60,7 +60,11 @@ autosummary_generate = True # Docstrings of private methods -autodoc_default_options = {"members": True, "undoc-members": True, "private-members": False} +autodoc_default_options = { + "members": True, + "undoc-members": True, + "private-members": False, +} # -- Options for HTML output ------------------------------------------------- # https://www.sphinx-doc.org/en/master/usage/configuration.html#options-for-html-output diff --git a/pyproject.toml b/pyproject.toml index 28ce31f9..e78068da 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -78,21 +78,11 @@ test = [ "pytest-cov", "sphinx-pyproject" ] -lint = [ - "isort", - "black", - "pyproject-flake8", - "flake8", - "flake8-nb", - "mypy", - "pre-commit", - "tox" -] -"black[jupyter]" = [] -pandas = [] dev = [ - "pandas>=2.0.3", - "geopandas>=0.13.2" + "mypy", + "pre-commit", + "ruff", + "tox", ] [tool.pdm.scripts] @@ -113,40 +103,83 @@ addopts = "-v" warn_unreachable = true ignore_missing_imports = true -[tool.isort] -profile = "black" -multi_line_output = 3 -include_trailing_comma = true -force_grid_wrap = 0 -use_parentheses = true -ensure_newline_before_comments = true -line_length = 120 - [tool.coverage.run] source = "src" omit = "tests/*" relative_files = true -[tool.yapf] -blank_line_before_nested_class_or_def = true -column_limit = 88 - -[tool.black] -line-length = 120 -exclude = ''' -/( - \.git - | \.tox - | \venv - | \.venv - | \*.env - | \build - | \dist -)/ -''' - -[tool.flake8] -max-line-length = "120" -extend-ignore = [ - "E501", +[tool.ruff] +# Exclude a variety of commonly ignored directories. +exclude = [ + ".bzr", + ".direnv", + ".eggs", + ".git", + ".git-rewrite", + ".hg", + ".ipynb_checkpoints", + ".mypy_cache", + ".nox", + ".pants.d", + ".pyenv", + ".pytest_cache", + ".pytype", + ".ruff_cache", + ".svn", + ".tox", + ".venv", + ".vscode", + "__pypackages__", + "_build", + "buck-out", + "build", + "dist", + "node_modules", + "site-packages", + "venv", +] + +[tool.ruff.lint] +extend-fixable = [ + # Instead of trailing-whitespace + "W291", "W293" + ] + +extend-select = [ + # Instead of pydocstyle + "D", + #Instead of flake8 + "E", "F","B", + # Instead of pep8-naming + "N", + # Instead of flake8-debugger or debug-statements + "T10", ] + +ignore = [ + "E203", + "E501", + # Avoid incompatible rules + "D203", + "D213", + + # Ignore this rules so that precommit passes. Uncomment to start fixing them + "B006", "B008", "B904", "B012", "B024", + "D", + "F401", +] + +[tool.ruff.lint.extend-per-file-ignores] +# Ignore `D` rules everywhere except for the `src/` directory. +"!src/**.py" = ["D"] + +[tool.ruff.lint.pycodestyle] +max-line-length = 120 + +[tool.ruff.format] +quote-style = "double" +indent-style = "space" +skip-magic-trailing-comma = false +line-ending = "auto" +docstring-code-format = true +docstring-code-line-length = "dynamic" diff --git a/src/physrisk/api/v1/common.py b/src/physrisk/api/v1/common.py index 5915cd7e..c79387b8 100644 --- a/src/physrisk/api/v1/common.py +++ b/src/physrisk/api/v1/common.py @@ -1,7 +1,7 @@ from typing import Dict, List, Optional, Union import numpy as np -from pydantic import BaseModel, Field +from pydantic import BaseModel, ConfigDict, Field class TypedArray(np.ndarray): @@ -23,25 +23,33 @@ class Array(np.ndarray, metaclass=ArrayMeta): pass -class Asset(BaseModel, extra="allow"): - """Defines an asset. An asset is identified first by its asset_class and then by its type within the class. +class Asset(BaseModel): + """Defines an asset. + + An asset is identified first by its asset_class and then by its type within the class. An asset's value may be impacted through damage or through disruption disruption being reduction of an asset's ability to generate cashflows (or equivalent value, e.g. by reducing expenses or increasing sales). """ + model_config = ConfigDict(extra="allow") + asset_class: str = Field( description="name of asset class; corresponds to physrisk class names, e.g. PowerGeneratingAsset" ) latitude: float = Field(description="Latitude in degrees") longitude: float = Field(description="Longitude in degrees") - type: Optional[str] = Field(None, description="Type of the asset //") + type: Optional[str] = Field( + None, description="Type of the asset //" + ) location: Optional[str] = Field( - None, description="Location (e.g. Africa, Asia, Europe, Global, Oceania, North America, South America)" + None, + description="Location (e.g. Africa, Asia, Europe, Global, Oceania, North America, South America)", ) capacity: Optional[float] = Field(None, description="Power generation capacity") attributes: Optional[Dict[str, str]] = Field( - None, description="Bespoke attributes (e.g. number of storeys, structure type, occupancy type)" + None, + description="Bespoke attributes (e.g. number of storeys, structure type, occupancy type)", ) @@ -81,7 +89,8 @@ class IntensityCurve(BaseModel): intensities: List[float] = Field([], description="Hazard indicator intensities.") return_periods: Optional[List[float]] = Field( - [], description="[Deprecated] Return period in years in the case of an acute hazard." + [], + description="[Deprecated] Return period in years in the case of an acute hazard.", ) index_values: Optional[Union[List[float], List[str]]] = Field( [], @@ -99,37 +108,40 @@ class IntensityCurve(BaseModel): class ExceedanceCurve(BaseModel): """General exceedance curve (e.g. hazazrd, impact).""" + model_config = ConfigDict(arbitrary_types_allowed=True) values: np.ndarray = Field(default_factory=lambda: np.zeros(10), description="") - exceed_probabilities: np.ndarray = Field(default_factory=lambda: np.zeros(10), description="") - - class Config: - arbitrary_types_allowed = True + exceed_probabilities: np.ndarray = Field( + default_factory=lambda: np.zeros(10), description="" + ) class Distribution(BaseModel): """General exceedance curve (e.g. hazazrd, impact).""" + model_config = ConfigDict(arbitrary_types_allowed=True) bin_edges: np.ndarray = Field(default_factory=lambda: np.zeros(11), description="") - probabilities: np.ndarray = Field(default_factory=lambda: np.zeros(10), description="") - - class Config: - arbitrary_types_allowed = True + probabilities: np.ndarray = Field( + default_factory=lambda: np.zeros(10), description="" + ) class HazardEventDistrib(BaseModel): """Intensity curve of an acute hazard.""" - intensity_bin_edges: np.ndarray = Field(default_factory=lambda: np.zeros(10), description="") - probabilities: np.ndarray = Field(default_factory=lambda: np.zeros(10), description="") + model_config = ConfigDict(arbitrary_types_allowed=True) + intensity_bin_edges: np.ndarray = Field( + default_factory=lambda: np.zeros(10), description="" + ) + probabilities: np.ndarray = Field( + default_factory=lambda: np.zeros(10), description="" + ) path: List[str] = Field([], description="Path to the hazard indicator data source.") - class Config: - arbitrary_types_allowed = True - class VulnerabilityCurve(BaseModel): """Defines a damage or disruption curve.""" + model_config = ConfigDict(arbitrary_types_allowed=True) asset_type: str = Field(...) location: str = Field(...) event_type: str = Field(description="hazard event type, e.g. RiverineInundation") @@ -139,10 +151,9 @@ class VulnerabilityCurve(BaseModel): intensity: List[float] = Field(...) intensity_units: str = Field(description="units of the intensity") impact_mean: List[float] = Field(description="mean impact (damage or disruption)") - impact_std: List[float] = Field(description="standard deviation of impact (damage or disruption)") - - class Config: - arbitrary_types_allowed = True + impact_std: List[float] = Field( + description="standard deviation of impact (damage or disruption)" + ) class VulnerabilityCurves(BaseModel): @@ -154,9 +165,13 @@ class VulnerabilityCurves(BaseModel): class VulnerabilityDistrib(BaseModel): """Defines a vulnerability matrix.""" - intensity_bin_edges: np.ndarray = Field(default_factory=lambda: np.zeros(10), description="") - impact_bin_edges: np.ndarray = Field(default_factory=lambda: np.zeros(10), description="") - prob_matrix: np.ndarray = Field(default_factory=lambda: np.zeros(10), description="") - - class Config: - arbitrary_types_allowed = True + model_config = ConfigDict(arbitrary_types_allowed=True) + intensity_bin_edges: np.ndarray = Field( + default_factory=lambda: np.zeros(10), description="" + ) + impact_bin_edges: np.ndarray = Field( + default_factory=lambda: np.zeros(10), description="" + ) + prob_matrix: np.ndarray = Field( + default_factory=lambda: np.zeros(10), description="" + ) diff --git a/src/physrisk/api/v1/exposure_req_resp.py b/src/physrisk/api/v1/exposure_req_resp.py index c4507718..25567173 100644 --- a/src/physrisk/api/v1/exposure_req_resp.py +++ b/src/physrisk/api/v1/exposure_req_resp.py @@ -11,7 +11,8 @@ class AssetExposureRequest(BaseModel): assets: Assets calc_settings: CalcSettings = Field( - default_factory=CalcSettings, description="Interpolation method." # type:ignore + default_factory=CalcSettings, # type:ignore + description="Interpolation method.", ) scenario: str = Field("rcp8p5", description="Name of scenario ('rcp8p5')") year: int = Field( @@ -29,7 +30,9 @@ class AssetExposureRequest(BaseModel): class Exposure(BaseModel): category: str value: Optional[float] - path: str = Field("unknown", description="Path to the hazard indicator data source.") + path: str = Field( + "unknown", description="Path to the hazard indicator data source." + ) class AssetExposure(BaseModel): @@ -40,7 +43,9 @@ class AssetExposure(BaseModel): description="""Asset identifier; will appear if provided in the request otherwise order of assets in response is identical to order of assets in request.""", ) - exposures: Dict[str, Exposure] = Field({}, description="Category (value) for each hazard type (key).") + exposures: Dict[str, Exposure] = Field( + {}, description="Category (value) for each hazard type (key)." + ) class AssetExposureResponse(BaseModel): diff --git a/src/physrisk/api/v1/hazard_data.py b/src/physrisk/api/v1/hazard_data.py index 92b54367..9c1b696b 100644 --- a/src/physrisk/api/v1/hazard_data.py +++ b/src/physrisk/api/v1/hazard_data.py @@ -10,13 +10,15 @@ class Colormap(BaseModel): """Provides details of colormap.""" min_index: Optional[int] = Field( - 1, description="Value of colormap minimum. Constant min for a group of maps can facilitate comparison." + 1, + description="Value of colormap minimum. Constant min for a group of maps can facilitate comparison.", ) min_value: float = Field( description="Value of colormap minimum. Constant min for a group of maps can facilitate comparison." ) max_index: Optional[int] = Field( - 255, description="Value of colormap maximum. Constant max for a group of maps can facilitate comparison." + 255, + description="Value of colormap maximum. Constant max for a group of maps can facilitate comparison.", ) max_value: float = Field( description="Value of colormap maximum. Constant max for a group of maps can facilitate comparison." @@ -51,7 +53,9 @@ class Period(BaseModel): """Provides information about a period, which currently corresponds to a year, belonging to a scenario.""" year: int - map_id: str = Field(description="If present, identifier to be used for looking up map tiles from server.") + map_id: str = Field( + description="If present, identifier to be used for looking up map tiles from server." + ) class Scenario(BaseModel): @@ -70,7 +74,10 @@ class HazardResource(BaseModel): """Provides information about a set of hazard indicators, including available scenarios and years.""" hazard_type: str = Field(description="Type of hazard.") - group_id: Optional[str] = Field("public", description="Identifier of the resource group (used for authentication).") + group_id: Optional[str] = Field( + "public", + description="Identifier of the resource group (used for authentication).", + ) path: str = Field(description="Full path to the indicator array.") indicator_id: str = Field( description="Identifier of the hazard indicator (i.e. the modelled quantity), e.g. 'flood_depth'." @@ -83,14 +90,23 @@ class HazardResource(BaseModel): indicator_model_gcm: str = Field( description="Identifier of general circulation model(s) used in the derivation of the indicator." ) - params: Dict[str, List[str]] = Field({}, description="Parameters used to expand wild-carded fields.") + params: Dict[str, List[str]] = Field( + {}, description="Parameters used to expand wild-carded fields." + ) display_name: str = Field(description="Text used to display indicator.") - display_groups: List[str] = Field([], description="Text used to group the (expanded) indicators for display.") + display_groups: List[str] = Field( + [], description="Text used to group the (expanded) indicators for display." + ) description: str = Field( description="Brief description in mark down of the indicator and model that generated the indicator." ) - map: Optional[MapInfo] = Field(None, description="Optional information used for display of the indicator in a map.") - scenarios: List[Scenario] = Field(description="Climate change scenarios for which the indicator is available.") + map: Optional[MapInfo] = Field( + None, + description="Optional information used for display of the indicator in a map.", + ) + scenarios: List[Scenario] = Field( + description="Climate change scenarios for which the indicator is available." + ) units: str = Field(description="Units of the hazard indicator.") def expand(self): @@ -122,7 +138,9 @@ def expand_resource( deep=True, update={ "indicator_id": expand(item.indicator_id, key, param), - "indicator_model_gcm": expand(item.indicator_model_gcm, key, param), + "indicator_model_gcm": expand( + item.indicator_model_gcm, key, param + ), "display_name": expand(item.display_name, key, param), "path": expand(item.path, key, param), "map": ( @@ -132,7 +150,13 @@ def expand_resource( item.map.model_copy( deep=True, update={ - "path": expand(item.map.path if item.map.path is not None else "", key, param) + "path": expand( + item.map.path + if item.map.path is not None + else "", + key, + param, + ) }, ) ) @@ -154,7 +178,8 @@ class InventorySource(Flag): class HazardAvailabilityRequest(BaseModel): types: Optional[List[str]] = [] # e.g. ["RiverineInundation"] sources: Optional[List[str]] = Field( - None, description="Sources of inventory, can be 'embedded', 'hazard' or 'hazard_test'." + None, + description="Sources of inventory, can be 'embedded', 'hazard' or 'hazard_test'.", ) @@ -168,7 +193,9 @@ class HazardDescriptionRequest(BaseModel): class HazardDescriptionResponse(BaseModel): - descriptions: Dict[str, str] = Field(description="For each path (key), the description markdown (value).") + descriptions: Dict[str, str] = Field( + description="For each path (key), the description markdown (value)." + ) class HazardDataRequestItem(BaseModel): @@ -176,7 +203,9 @@ class HazardDataRequestItem(BaseModel): latitudes: List[float] request_item_id: str hazard_type: Optional[str] = None # e.g. RiverineInundation - event_type: Optional[str] = None # e.g. RiverineInundation; deprecated: use hazard_type + event_type: Optional[str] = ( + None # e.g. RiverineInundation; deprecated: use hazard_type + ) indicator_id: str indicator_model_gcm: Optional[str] = "" path: Optional[str] = None diff --git a/src/physrisk/api/v1/hazard_image.py b/src/physrisk/api/v1/hazard_image.py index 217b3d51..afc9c687 100644 --- a/src/physrisk/api/v1/hazard_image.py +++ b/src/physrisk/api/v1/hazard_image.py @@ -17,7 +17,9 @@ class Tile(NamedTuple): class HazardImageRequest(BaseHazardRequest): - resource: str = Field(description="Full path to the array; formed by '{path}/{id}'.") + resource: str = Field( + description="Full path to the array; formed by '{path}/{id}'." + ) scenario_id: str year: int colormap: Optional[str] = Field("heating") diff --git a/src/physrisk/api/v1/impact_req_resp.py b/src/physrisk/api/v1/impact_req_resp.py index bb5e753c..700f92f9 100644 --- a/src/physrisk/api/v1/impact_req_resp.py +++ b/src/physrisk/api/v1/impact_req_resp.py @@ -1,14 +1,22 @@ from enum import Enum from typing import Dict, List, NamedTuple, Optional, Sequence -from pydantic import BaseModel, Field, computed_field - -from physrisk.api.v1.common import Assets, Distribution, ExceedanceCurve, VulnerabilityDistrib +from pydantic import BaseModel, ConfigDict, Field, computed_field + +from physrisk.api.v1.common import ( + Assets, + Distribution, + ExceedanceCurve, + VulnerabilityDistrib, +) from physrisk.api.v1.hazard_data import Scenario class CalcSettings(BaseModel): - hazard_interp: str = Field("floor", description="Method used for interpolation of hazards: 'floor' or 'bilinear'.") + hazard_interp: str = Field( + "floor", + description="Method used for interpolation of hazards: 'floor' or 'bilinear'.", + ) class AssetImpactRequest(BaseModel): @@ -16,19 +24,31 @@ class AssetImpactRequest(BaseModel): assets: Assets calc_settings: CalcSettings = Field( - default_factory=CalcSettings, description="Interpolation method." # type:ignore + default_factory=CalcSettings, # type:ignore + description="Interpolation method.", + ) + include_asset_level: bool = Field( + True, description="If true, include asset-level impacts." + ) + include_measures: bool = Field( + False, description="If true, include calculation of risk measures." + ) + include_calc_details: bool = Field( + True, description="If true, include impact calculation details." + ) + use_case_id: str = Field( + "", + description="Identifier for 'use case' used in the risk measures calculation.", ) - include_asset_level: bool = Field(True, description="If true, include asset-level impacts.") - include_measures: bool = Field(False, description="If true, include calculation of risk measures.") - include_calc_details: bool = Field(True, description="If true, include impact calculation details.") - use_case_id: str = Field("", description="Identifier for 'use case' used in the risk measures calculation.") provider_max_requests: Dict[str, int] = Field( {}, description="The maximum permitted number of requests \ to external providers. This setting is intended in particular for paid-for data. The key is the provider \ ID and the value is the maximum permitted requests.", ) - scenarios: Optional[Sequence[str]] = Field([], description="Name of scenarios ('rcp8p5')") + scenarios: Optional[Sequence[str]] = Field( + [], description="Name of scenarios ('rcp8p5')" + ) years: Optional[Sequence[int]] = Field( [], description="""Projection year (2030, 2050, 2080). Any year before 2030, @@ -61,9 +81,12 @@ class RiskMeasureDefinition(BaseModel): class RiskScoreValue(BaseModel): - value: Category = Field("", description="Value of the score: red, amber, green, nodata.") + value: Category = Field( + "", description="Value of the score: red, amber, green, nodata." + ) label: str = Field( - "", description="Short description of value, e.g. material increase in loss for 1-in-100 year event." + "", + description="Short description of value, e.g. material increase in loss for 1-in-100 year event.", ) description: str = Field( "", @@ -72,18 +95,23 @@ class RiskScoreValue(BaseModel): ) -class ScoreBasedRiskMeasureDefinition(BaseModel, frozen=True): - hazard_types: List[str] = Field([], description="Defines the hazards that the measure is used for.") - values: List[RiskScoreValue] = Field([], description="Defines the set of values that the score can take.") +class ScoreBasedRiskMeasureDefinition(BaseModel): + hazard_types: List[str] = Field( + [], description="Defines the hazards that the measure is used for." + ) + values: List[RiskScoreValue] = Field( + [], description="Defines the set of values that the score can take." + ) underlying_measures: List[RiskMeasureDefinition] = Field( - [], description="Defines the underlying risk measures from which the scores are inferred." + [], + description="Defines the underlying risk measures from which the scores are inferred.", ) # for now underlying measures defined directly rather than by referencing an ID via: # underlying_measure_ids: List[str] = Field( # [], description="The identifiers of the underlying risk measures from which the scores are inferred." # ) - # should be sufficient to pass frozen=True, but does not seem to work (pydantic docs says feature in beta) + # It is not enough to pass frozen=True, since not all attributes are hashable def __hash__(self): return id(self) @@ -124,7 +152,9 @@ class AcuteHazardCalculationDetails(BaseModel): hazard_exceedance: ExceedanceCurve hazard_distribution: Distribution vulnerability_distribution: VulnerabilityDistrib - hazard_path: List[str] = Field("unknown", description="Path to the hazard indicator data source.") + hazard_path: List[str] = Field( + "unknown", description="Path to the hazard indicator data source." + ) class ImpactKey(BaseModel): @@ -136,6 +166,7 @@ class ImpactKey(BaseModel): class AssetSingleImpact(BaseModel): """Impact at level of single asset and single type of hazard.""" + model_config = ConfigDict(arbitrary_types_allowed=True) key: ImpactKey @computed_field # deprecated: use key instead @@ -161,9 +192,6 @@ def year(self) -> str: description="""Details of impact calculation for acute hazard calculations.""", ) - class Config: - arbitrary_types_allowed = True - class AssetLevelImpact(BaseModel): """Impact at asset level. Each asset can have impacts for multiple hazard types.""" @@ -173,7 +201,9 @@ class AssetLevelImpact(BaseModel): description="""Asset identifier; will appear if provided in the request otherwise order of assets in response is identical to order of assets in request.""", ) - impacts: List[AssetSingleImpact] = Field([], description="Impacts for each hazard type.") + impacts: List[AssetSingleImpact] = Field( + [], description="Impacts for each hazard type." + ) class AssetImpactResponse(BaseModel): @@ -196,12 +226,18 @@ def __init__(self, risk_measures: RiskMeasures): def _key(self, key: RiskMeasureKey): return self.Key( - hazard_type=key.hazard_type, scenario_id=key.scenario_id, year=key.year, measure_id=key.measure_id + hazard_type=key.hazard_type, + scenario_id=key.scenario_id, + year=key.year, + measure_id=key.measure_id, ) def get_measure(self, hazard_type: str, scenario: str, year: int): measure_key = self.Key( - hazard_type=hazard_type, scenario_id=scenario, year=str(year), measure_id=self.measure_set_id + hazard_type=hazard_type, + scenario_id=scenario, + year=str(year), + measure_id=self.measure_set_id, ) measure = self.measures[measure_key] asset_scores, asset_measures = ( @@ -212,11 +248,14 @@ def get_measure(self, hazard_type: str, scenario: str, year: int): measure_ids = self.measure_definition.asset_measure_ids_for_hazard[hazard_type] # measure definitions for each asset measure_definitions = [ - self.measure_definition.score_definitions[mid] if mid != "na" else None for mid in measure_ids + self.measure_definition.score_definitions[mid] if mid != "na" else None + for mid in measure_ids ] return asset_scores, asset_measures, measure_definitions - def get_score_details(self, score: int, definition: ScoreBasedRiskMeasureDefinition): + def get_score_details( + self, score: int, definition: ScoreBasedRiskMeasureDefinition + ): rs_value = next(v for v in definition.values if v.value == score) return rs_value.label, rs_value.description diff --git a/src/physrisk/container.py b/src/physrisk/container.py index 36fcc257..49c92311 100644 --- a/src/physrisk/container.py +++ b/src/physrisk/container.py @@ -28,11 +28,16 @@ def __init__( self.store = store self.reader = reader - def hazard_model(self, interpolation: str = "floor", provider_max_requests: Dict[str, int] = {}): + def hazard_model( + self, interpolation: str = "floor", provider_max_requests: Dict[str, int] = {} + ): # this is done to allow interpolation to be set dynamically, e.g. different requests can have different # parameters. return ZarrHazardModel( - source_paths=self.source_paths, store=self.store, reader=self.reader, interpolation=interpolation + source_paths=self.source_paths, + store=self.store, + reader=self.reader, + interpolation=interpolation, ) @@ -48,7 +53,9 @@ class Container(containers.DeclarativeContainer): inventory_reader = providers.Singleton(InventoryReader) - inventory = providers.Singleton(_create_inventory, reader=inventory_reader, sources=config.zarr_sources) + inventory = providers.Singleton( + _create_inventory, reader=inventory_reader, sources=config.zarr_sources + ) source_paths = providers.Factory(create_source_paths, inventory=inventory) @@ -56,11 +63,15 @@ class Container(containers.DeclarativeContainer): zarr_reader = providers.Singleton(ZarrReader, store=zarr_store) - hazard_model_factory = providers.Factory(ZarrHazardModelFactory, reader=zarr_reader, source_paths=source_paths) + hazard_model_factory = providers.Factory( + ZarrHazardModelFactory, reader=zarr_reader, source_paths=source_paths + ) measures_factory = providers.Factory(calc.DefaultMeasuresFactory) - vulnerability_models_factory = providers.Factory(DictBasedVulnerabilityModelsFactory) + vulnerability_models_factory = providers.Factory( + DictBasedVulnerabilityModelsFactory + ) requester = providers.Singleton( Requester, diff --git a/src/physrisk/data/colormap_provider.py b/src/physrisk/data/colormap_provider.py index 8208bcd6..0f6e983c 100644 --- a/src/physrisk/data/colormap_provider.py +++ b/src/physrisk/data/colormap_provider.py @@ -781,7 +781,12 @@ # for test cases map_test = dict((str(i), [i, 0, 0, 0]) for i in range(256)) -colormaps = {"flare": map_flare, "heating": map_heating, "heating_2": map_heating_2, "test": map_test} +colormaps = { + "flare": map_flare, + "heating": map_heating, + "heating_2": map_heating_2, + "test": map_test, +} def colormap(id: str): diff --git a/src/physrisk/data/geotiff_reader.py b/src/physrisk/data/geotiff_reader.py index 51ad2913..840560d3 100644 --- a/src/physrisk/data/geotiff_reader.py +++ b/src/physrisk/data/geotiff_reader.py @@ -37,7 +37,13 @@ def dataset_read_bounded(dataset, longitudes, latitudes, window_half_width=0.01) offsets = [[0, 0], [-hw, -hw], [-hw, hw], [hw, hw], [hw, -hw]] points = [] for offset in offsets: - points = chain(points, [[lon + offset[0], lat + offset[1]] for (lon, lat) in zip(longitudes, latitudes)]) + points = chain( + points, + [ + [lon + offset[0], lat + offset[1]] + for (lon, lat) in zip(longitudes, latitudes) + ], + ) samples = np.array(list(rasterio.sample.sample_gen(dataset, points))) samples.resize([len(offsets), len(longitudes)]) diff --git a/src/physrisk/data/hazard_data_provider.py b/src/physrisk/data/hazard_data_provider.py index 77f029a2..4b3e457d 100644 --- a/src/physrisk/data/hazard_data_provider.py +++ b/src/physrisk/data/hazard_data_provider.py @@ -31,7 +31,12 @@ class SourcePath(Protocol): """ def __call__( - self, *, indicator_id: str, scenario: str, year: int, hint: Optional[HazardDataHint] = None + self, + *, + indicator_id: str, + scenario: str, + year: int, + hint: Optional[HazardDataHint] = None, ) -> str: ... @@ -56,7 +61,9 @@ def __init__( ValueError: If interpolation not in permitted list. """ self._get_source_path = get_source_path - self._reader = zarr_reader if zarr_reader is not None else ZarrReader(store=store) + self._reader = ( + zarr_reader if zarr_reader is not None else ZarrReader(store=store) + ) if interpolation not in ["floor", "linear", "max", "min"]: raise ValueError("interpolation must be 'floor', 'linear', 'max' or 'min'") self._interpolation = interpolation @@ -94,14 +101,18 @@ def get_data( path: Path to the hazard indicator data source. """ - path = self._get_source_path(indicator_id=indicator_id, scenario=scenario, year=year, hint=hint) + path = self._get_source_path( + indicator_id=indicator_id, scenario=scenario, year=year, hint=hint + ) if buffer is None: values, indices, units = self._reader.get_curves( path, longitudes, latitudes, self._interpolation ) # type: ignore else: if buffer < 0 or 1000 < buffer: - raise Exception("The buffer must be an integer between 0 and 1000 metres.") + raise Exception( + "The buffer must be an integer between 0 and 1000 metres." + ) values, indices, units = self._reader.get_max_curves( path, [ @@ -109,7 +120,9 @@ def get_data( Point(longitude, latitude) if buffer == 0 else Point(longitude, latitude).buffer( - ZarrReader._get_equivalent_buffer_in_arc_degrees(latitude, buffer) + ZarrReader._get_equivalent_buffer_in_arc_degrees( + latitude, buffer + ) ) ) for longitude, latitude in zip(longitudes, latitudes) diff --git a/src/physrisk/data/image_creator.py b/src/physrisk/data/image_creator.py index 8cbba013..c31db74f 100644 --- a/src/physrisk/data/image_creator.py +++ b/src/physrisk/data/image_creator.py @@ -50,7 +50,9 @@ def convert( bytes: Image data. """ try: - image = self._to_image(path, colormap, tile=tile, min_value=min_value, max_value=max_value) + image = self._to_image( + path, colormap, tile=tile, min_value=min_value, max_value=max_value + ) except Exception as e: logger.exception(e) image = Image.fromarray(np.array([[0]]), mode="RGBA") @@ -95,14 +97,18 @@ def _to_image( tile_size = 512 # data = self.reader.all_data(tile_path) if len(data.shape) == 3: - index = len(self.reader.get_index_values(data)) - 1 if index is None else index + index = ( + len(self.reader.get_index_values(data)) - 1 if index is None else index + ) if tile is None: # return whole array data = data[index, :, :] # .squeeze(axis=0) else: # (from zarr 2.16.0 we can also use block indexing) data = data[ - index, tile_size * tile.y : tile_size * (tile.y + 1), tile_size * tile.x : tile_size * (tile.x + 1) + index, + tile_size * tile.y : tile_size * (tile.y + 1), + tile_size * tile.x : tile_size * (tile.x + 1), ] if any(dim > 4000 for dim in data.shape): @@ -164,7 +170,11 @@ def _to_rgba( # noqa: C901 if nodata_lower: mask_nodata = data <= nodata_lower if nodata_upper: - mask_nodata = (mask_nodata | (data >= nodata_upper)) if mask_nodata is not None else (data >= nodata_upper) + mask_nodata = ( + (mask_nodata | (data >= nodata_upper)) + if mask_nodata is not None + else (data >= nodata_upper) + ) if min_value is None: min_value = np.nanmin(data) @@ -189,17 +199,27 @@ def _to_rgba( # noqa: C901 result[mask_le_min] = 1 del mask_ge_max, mask_le_min - final = red[result] + (green[result] << 8) + (blue[result] << 16) + (a[result] << 24) + final = ( + red[result] + + (green[result] << 8) + + (blue[result] << 16) + + (a[result] << 24) + ) return final @staticmethod def test_store(path: str): store = zarr.storage.MemoryStore(root="hazard.zarr") root = zarr.open(store=store, mode="w") - x, y = np.meshgrid((np.arange(1000) - 500.0) / 500.0, (np.arange(1000) - 500.0) / 500.0) + x, y = np.meshgrid( + (np.arange(1000) - 500.0) / 500.0, (np.arange(1000) - 500.0) / 500.0 + ) im = np.exp(-(x**2 + y**2)) z = root.create_dataset( # type: ignore - path, shape=(1, im.shape[0], im.shape[1]), chunks=(1, im.shape[0], im.shape[1]), dtype="f4" + path, + shape=(1, im.shape[0], im.shape[1]), + chunks=(1, im.shape[0], im.shape[1]), + dtype="f4", ) z[0, :, :] = im return store diff --git a/src/physrisk/data/inventory.py b/src/physrisk/data/inventory.py index f4ab6246..d766d0fb 100644 --- a/src/physrisk/data/inventory.py +++ b/src/physrisk/data/inventory.py @@ -2,11 +2,10 @@ import hashlib import importlib.resources import json -import logging from collections import defaultdict from typing import DefaultDict, Dict, Iterable, List, Tuple -from pydantic import TypeAdapter, parse_obj_as +from pydantic import TypeAdapter import physrisk.data.colormap_provider as colormap_provider import physrisk.data.static.hazard @@ -27,10 +26,14 @@ def __init__(self, hazard_resources: Iterable[HazardResource]): hazard_resources (Iterable[HazardResource]): list of resources """ self.resources: Dict[str, HazardResource] = {} - self.resources_by_type_id: DefaultDict[Tuple[str, str], List[HazardResource]] = defaultdict(list) + self.resources_by_type_id: DefaultDict[ + Tuple[str, str], List[HazardResource] + ] = defaultdict(list) for resource in hazard_resources: self.resources[resource.key()] = resource - self.resources_by_type_id[(resource.hazard_type, resource.indicator_id)].append(resource) + self.resources_by_type_id[ + (resource.hazard_type, resource.indicator_id) + ].append(resource) def json_ordered(self): sorted_resources = sorted(self.resources_by_type_id.items()) @@ -49,7 +52,9 @@ class EmbeddedInventory(Inventory): """ def __init__(self): - with importlib.resources.open_text(physrisk.data.static.hazard, "inventory.json") as f: + with importlib.resources.open_text( + physrisk.data.static.hazard, "inventory.json" + ) as f: models = TypeAdapter(HazardModels).validate_python(json.load(f)).resources expanded_models = expand(models) super().__init__(expanded_models) @@ -95,7 +100,9 @@ def expand(resources: List[HazardResource]) -> List[HazardResource]: scenario.periods = [] for year in scenario.years: name_format = model.map.path - path = name_format.format(scenario=scenario.id, year=year, return_period=1000) + path = name_format.format( + scenario=scenario.id, year=year, return_period=1000 + ) id = alphanumeric(path)[0:6] scenario.periods.append(Period(year=year, map_id=id)) # if a period was specified explicitly, we check that hash is the same: a build-in check diff --git a/src/physrisk/data/inventory_reader.py b/src/physrisk/data/inventory_reader.py index c566f0c0..09b07207 100644 --- a/src/physrisk/data/inventory_reader.py +++ b/src/physrisk/data/inventory_reader.py @@ -48,7 +48,9 @@ def read(self, path: str) -> List[HazardResource]: if not self._fs.exists(self._full_path(path)): return [] json_str = self.read_json(path) - models = TypeAdapter(HazardModels).validate_python(json.loads(json_str)).resources + models = ( + TypeAdapter(HazardModels).validate_python(json.loads(json_str)).resources + ) return models def read_description_markdown(self, paths: List[str]) -> Dict[str, str]: diff --git a/src/physrisk/data/pregenerated_hazard_model.py b/src/physrisk/data/pregenerated_hazard_model.py index d813394b..abbb33bd 100644 --- a/src/physrisk/data/pregenerated_hazard_model.py +++ b/src/physrisk/data/pregenerated_hazard_model.py @@ -55,13 +55,16 @@ def _get_hazard_events( # noqa: C901 # can change max_workers=1 for debugging with concurrent.futures.ThreadPoolExecutor(max_workers=8) as executor: futures = [ - executor.submit(self._get_hazard_events_batch, batches[key], responses) for key in batches.keys() + executor.submit(self._get_hazard_events_batch, batches[key], responses) + for key in batches.keys() ] concurrent.futures.wait(futures) return responses def _get_hazard_events_batch( - self, batch: List[HazardDataRequest], responses: MutableMapping[HazardDataRequest, HazardDataResponse] + self, + batch: List[HazardDataRequest], + responses: MutableMapping[HazardDataRequest, HazardDataResponse], ): try: hazard_type, indicator_id, scenario, year, hint, buffer = ( @@ -75,7 +78,9 @@ def _get_hazard_events_batch( longitudes = [req.longitude for req in batch] latitudes = [req.latitude for req in batch] if indicator_data(hazard_type, indicator_id) == IndicatorData.EVENT: - intensities, return_periods, units, path = self.hazard_data_providers[hazard_type].get_data( + intensities, return_periods, units, path = self.hazard_data_providers[ + hazard_type + ].get_data( longitudes, latitudes, indicator_id=indicator_id, @@ -87,21 +92,38 @@ def _get_hazard_events_batch( for i, req in enumerate(batch): valid = ~np.isnan(intensities[i, :]) - valid_periods, valid_intensities = return_periods[valid], intensities[i, :][valid] + valid_periods, valid_intensities = ( + return_periods[valid], + intensities[i, :][valid], + ) if len(valid_periods) == 0: - valid_periods, valid_intensities = np.array([100]), np.array([0]) - responses[req] = HazardEventDataResponse(valid_periods, valid_intensities, units, path) + valid_periods, valid_intensities = ( + np.array([100]), + np.array([0]), + ) + responses[req] = HazardEventDataResponse( + valid_periods, valid_intensities, units, path + ) else: # type: ignore - parameters, defns, units, path = self.hazard_data_providers[hazard_type].get_data( - longitudes, latitudes, indicator_id=indicator_id, scenario=scenario, year=year, hint=hint + parameters, defns, units, path = self.hazard_data_providers[ + hazard_type + ].get_data( + longitudes, + latitudes, + indicator_id=indicator_id, + scenario=scenario, + year=year, + hint=hint, ) for i, req in enumerate(batch): valid = ~np.isnan(parameters[i, :]) - responses[req] = HazardParameterDataResponse(parameters[i, :][valid], defns[valid], units, path) + responses[req] = HazardParameterDataResponse( + parameters[i, :][valid], defns[valid], units, path + ) except Exception as err: # e.g. the requested data is unavailable - for i, req in enumerate(batch): + for _i, req in enumerate(batch): responses[req] = HazardDataFailedResponse(err) return @@ -120,7 +142,9 @@ def __init__( super().__init__( { - t: HazardDataProvider(sp, zarr_reader=zarr_reader, interpolation=interpolation) + t: HazardDataProvider( + sp, zarr_reader=zarr_reader, interpolation=interpolation + ) for t, sp in source_paths.items() } ) diff --git a/src/physrisk/data/static/world.py b/src/physrisk/data/static/world.py index 84acec03..b82f6e7e 100644 --- a/src/physrisk/data/static/world.py +++ b/src/physrisk/data/static/world.py @@ -24,13 +24,18 @@ def get_countries_json(): countries = [ Country(continent=continent, country=country, country_iso_a3=code) - for (continent, country, code) in zip(world["continent"], world["name"], world["iso_a3"]) + for (continent, country, code) in zip( + world["continent"], world["name"], world["iso_a3"] + ) ] return json.dumps(Countries(items=countries).dict(), sort_keys=True, indent=4) -def get_countries_and_continents(longitudes: Union[List[float], np.ndarray], latitudes: Union[List[float], np.ndarray]): +def get_countries_and_continents( + longitudes: Union[List[float], np.ndarray], + latitudes: Union[List[float], np.ndarray], +): """Only for use when on-boarding; look up country and continent (e.g. for use in vulnerability models) by latitude and longitude.""" @@ -42,7 +47,9 @@ def get_countries_and_continents(longitudes: Union[List[float], np.ndarray], lat # consider using map here https://gadm.org/download_world.html world = gpd.read_file(gpd.datasets.get_path("naturalearth_lowres")) # type: ignore - gdf = gpd.GeoDataFrame(crs=world.crs, geometry=gpd.points_from_xy(longitudes, latitudes)) + gdf = gpd.GeoDataFrame( + crs=world.crs, geometry=gpd.points_from_xy(longitudes, latitudes) + ) result = gpd.sjoin(gdf, world, how="left") return list(result["name"]), list(result["continent"]) diff --git a/src/physrisk/data/zarr_reader.py b/src/physrisk/data/zarr_reader.py index c34bb088..580275be 100644 --- a/src/physrisk/data/zarr_reader.py +++ b/src/physrisk/data/zarr_reader.py @@ -45,7 +45,9 @@ def __init__( if store is None: # if no store is provided, attempt to connect to an S3 bucket if get_env is None: - raise TypeError("if no store specified, get_env is required to provide credentials") + raise TypeError( + "if no store specified, get_env is required to provide credentials" + ) store = ZarrReader.create_s3_zarr_store(get_env) @@ -54,12 +56,16 @@ def __init__( pass def all_data(self, set_id: str): - path = self._path_provider(set_id) if self._path_provider is not None else set_id + path = ( + self._path_provider(set_id) if self._path_provider is not None else set_id + ) z = self._root[path] # e.g. inundation/wri/v2/ return z @classmethod - def create_s3_zarr_store(cls, get_env: Callable[[str, Optional[str]], str] = get_env): + def create_s3_zarr_store( + cls, get_env: Callable[[str, Optional[str]], str] = get_env + ): access_key = get_env(cls.__access_key, "") secret_key = get_env(cls.__secret_key, "") s3_bucket = get_env(cls.__S3_bucket, "physrisk-hazard-indicators") @@ -90,7 +96,9 @@ def get_curves(self, set_id, longitudes, latitudes, interpolation="floor"): # assume that calls to this are large, not chatty if len(longitudes) != len(latitudes): raise ValueError("length of longitudes and latitudes not equal") - path = self._path_provider(set_id) if self._path_provider is not None else set_id + path = ( + self._path_provider(set_id) if self._path_provider is not None else set_id + ) z = self._root[path] # e.g. inundation/wri/v2/ # OSC-specific attributes contain transform and return periods @@ -102,7 +110,11 @@ def get_curves(self, set_id, longitudes, latitudes, interpolation="floor"): # in the case of acute risks, index_values will contain the return periods index_values = self.get_index_values(z) image_coords = self._get_coordinates( - longitudes, latitudes, crs, transform, pixel_is_area=interpolation != "floor" + longitudes, + latitudes, + crs, + transform, + pixel_is_area=interpolation != "floor", ) if interpolation == "floor": @@ -113,14 +125,22 @@ def get_curves(self, set_id, longitudes, latitudes, interpolation="floor"): ix = np.repeat(image_coords[0, :], len(index_values)) data = z.get_coordinate_selection((iz, iy, ix)) # type: ignore - return data.reshape([len(longitudes), len(index_values)]), np.array(index_values), units + return ( + data.reshape([len(longitudes), len(index_values)]), + np.array(index_values), + units, + ) elif interpolation in ["linear", "max", "min"]: - res = ZarrReader._linear_interp_frac_coordinates(z, image_coords, index_values, interpolation=interpolation) + res = ZarrReader._linear_interp_frac_coordinates( + z, image_coords, index_values, interpolation=interpolation + ) return res, np.array(index_values), units else: - raise ValueError("interpolation must have value 'floor', 'linear', 'max' or 'min") + raise ValueError( + "interpolation must have value 'floor', 'linear', 'max' or 'min" + ) def get_index_values(self, z: zarr.Array): index_values = z.attrs.get("index_values", [0]) @@ -142,7 +162,9 @@ def get_max_curves(self, set_id, shapes, interpolation="floor"): return_periods: return periods in years. units: units. """ - path = self._path_provider(set_id) if self._path_provider is not None else set_id + path = ( + self._path_provider(set_id) if self._path_provider is not None else set_id + ) z = self._root[path] # e.g. inundation/wri/v2/ # in the case of acute risks, index_values will contain the return periods @@ -153,32 +175,52 @@ def get_max_curves(self, set_id, shapes, interpolation="floor"): units: str = z.attrs.get("units", "default") matrix = np.array(~transform).reshape(3, 3).transpose()[:, :-1].reshape(6) - transformed_shapes = [affinity.affine_transform(shape, matrix) for shape in shapes] + transformed_shapes = [ + affinity.affine_transform(shape, matrix) for shape in shapes + ] pixel_offset = 0.5 if interpolation != "floor" else 0.0 multipoints = [ MultiPoint( [ (x - pixel_offset, y - pixel_offset) - for x in range(int(np.floor(shape.bounds[0])), int(np.ceil(shape.bounds[2])) + 1) - for y in range(int(np.floor(shape.bounds[1])), int(np.ceil(shape.bounds[3])) + 1) + for x in range( + int(np.floor(shape.bounds[0])), + int(np.ceil(shape.bounds[2])) + 1, + ) + for y in range( + int(np.floor(shape.bounds[1])), + int(np.ceil(shape.bounds[3])) + 1, + ) ] ).intersection(shape) for shape in transformed_shapes ] multipoints = [ ( - Point(0.5 * (shape.bounds[0] + shape.bounds[2]), 0.5 * (shape.bounds[1] + shape.bounds[3])) + Point( + 0.5 * (shape.bounds[0] + shape.bounds[2]), + 0.5 * (shape.bounds[1] + shape.bounds[3]), + ) if multipoint.is_empty else multipoint ) for shape, multipoint in zip(transformed_shapes, multipoints) ] - multipoints = [MultiPoint([(point.x, point.y)]) if isinstance(point, Point) else point for point in multipoints] + multipoints = [ + MultiPoint([(point.x, point.y)]) if isinstance(point, Point) else point + for point in multipoints + ] if interpolation == "floor": image_coords = np.floor( - np.array([[point.x, point.y] for multipoint in multipoints for point in multipoint.geoms]).transpose() + np.array( + [ + [point.x, point.y] + for multipoint in multipoints + for point in multipoint.geoms + ] + ).transpose() ).astype(int) image_coords[0, :] %= z.shape[2] @@ -196,10 +238,16 @@ def get_max_curves(self, set_id, shapes, interpolation="floor"): if isinstance(transformed_shape, Point) else MultiPoint(transformed_shape.exterior.coords) ) - for transformed_shape, multipoint in zip(transformed_shapes, multipoints) + for transformed_shape, multipoint in zip( + transformed_shapes, multipoints + ) ] image_coords = np.array( - [[point.x, point.y] for multipoint in multipoints for point in multipoint.geoms] + [ + [point.x, point.y] + for multipoint in multipoints + for point in multipoint.geoms + ] ).transpose() curves = ZarrReader._linear_interp_frac_coordinates( @@ -207,9 +255,13 @@ def get_max_curves(self, set_id, shapes, interpolation="floor"): ) else: - raise ValueError("interpolation must have value 'floor', 'linear', 'max' or 'min") + raise ValueError( + "interpolation must have value 'floor', 'linear', 'max' or 'min" + ) - numbers_of_points_per_shape = [len(multipoint.geoms) for multipoint in multipoints] + numbers_of_points_per_shape = [ + len(multipoint.geoms) for multipoint in multipoints + ] numbers_of_points_per_shape_cumulated = np.cumsum(numbers_of_points_per_shape) curves_max = np.array( [ @@ -222,7 +274,15 @@ def get_max_curves(self, set_id, shapes, interpolation="floor"): return curves_max, np.array(index_values), units - def get_max_curves_on_grid(self, set_id, longitudes, latitudes, interpolation="floor", delta_km=1.0, n_grid=5): + def get_max_curves_on_grid( + self, + set_id, + longitudes, + latitudes, + interpolation="floor", + delta_km=1.0, + n_grid=5, + ): """Get maximal intensity curve for a grid around a given latitude and longitude coordinate pair. It is almost equivalent to: self.get_max_curves @@ -265,7 +325,8 @@ def get_max_curves_on_grid(self, set_id, longitudes, latitudes, interpolation="f np.array(latitudes).reshape(n_data, 1, 1), (len(latitudes), n_grid, n_grid) ) lons_grid_baseline = np.broadcast_to( - np.array(longitudes).reshape(n_data, 1, 1), (len(longitudes), n_grid, n_grid) + np.array(longitudes).reshape(n_data, 1, 1), + (len(longitudes), n_grid, n_grid), ) lats_grid_offsets = delta_deg * grid.reshape((1, n_grid, 1)) lons_grid_offsets = ( @@ -276,24 +337,37 @@ def get_max_curves_on_grid(self, set_id, longitudes, latitudes, interpolation="f lats_grid = lats_grid_baseline + lats_grid_offsets lons_grid = lons_grid_baseline + lons_grid_offsets curves, return_periods, _ = self.get_curves( - set_id, lons_grid.reshape(-1), lats_grid.reshape(-1), interpolation=interpolation + set_id, + lons_grid.reshape(-1), + lats_grid.reshape(-1), + interpolation=interpolation, + ) + curves_max = np.nanmax( + curves.reshape((n_data, n_grid * n_grid, len(return_periods))), axis=1 ) - curves_max = np.nanmax(curves.reshape((n_data, n_grid * n_grid, len(return_periods))), axis=1) return curves_max, return_periods @staticmethod - def _linear_interp_frac_coordinates(z, image_coords, return_periods, interpolation="linear"): + def _linear_interp_frac_coordinates( + z, image_coords, return_periods, interpolation="linear" + ): """Return linear interpolated data from fractional row and column coordinates.""" icx = np.floor(image_coords[0, :]).astype(int)[..., None] # note periodic boundary condition ix = np.concatenate( - [icx % z.shape[2], icx % z.shape[2], (icx + 1) % z.shape[2], (icx + 1) % z.shape[2]], axis=1 - )[..., None].repeat( - len(return_periods), axis=2 - ) # points, 4, return_periods + [ + icx % z.shape[2], + icx % z.shape[2], + (icx + 1) % z.shape[2], + (icx + 1) % z.shape[2], + ], + axis=1, + )[..., None].repeat(len(return_periods), axis=2) # points, 4, return_periods icy = np.floor(image_coords[1, :]).astype(int)[..., None] - iy = np.concatenate([icy, icy + 1, icy, icy + 1], axis=1)[..., None].repeat(len(return_periods), axis=2) + iy = np.concatenate([icy, icy + 1, icy, icy + 1], axis=1)[..., None].repeat( + len(return_periods), axis=2 + ) iz = ( np.arange(len(return_periods), dtype=int)[None, ...] @@ -324,7 +398,10 @@ def _linear_interp_frac_coordinates(z, image_coords, return_periods, interpolati data = np.nan_to_num(data, 0) w_good = w * mask w_good_sum = np.transpose( - np.sum(w_good, axis=1).reshape(tuple([1]) + np.sum(w_good, axis=1).shape), axes=(1, 0, 2) + np.sum(w_good, axis=1).reshape( + tuple([1]) + np.sum(w_good, axis=1).shape + ), + axes=(1, 0, 2), ) w_used = np.divide(w_good, np.where(w_good_sum == 0.0, np.nan, w_good_sum)) return np.nan_to_num(np.sum(w_used * data, axis=1), nan=nan_output_value) @@ -332,7 +409,9 @@ def _linear_interp_frac_coordinates(z, image_coords, return_periods, interpolati elif interpolation == "max": data = np.where(np.isnan(data), -np.inf, data) return np.nan_to_num( - np.maximum.reduce([data[:, 0, :], data[:, 1, :], data[:, 2, :], data[:, 3, :]]), + np.maximum.reduce( + [data[:, 0, :], data[:, 1, :], data[:, 2, :], data[:, 3, :]] + ), nan=nan_output_value, neginf=nan_output_value, ) @@ -340,7 +419,9 @@ def _linear_interp_frac_coordinates(z, image_coords, return_periods, interpolati elif interpolation == "min": data = np.where(np.isnan(data), np.inf, data) return np.nan_to_num( - np.minimum.reduce([data[:, 0, :], data[:, 1, :], data[:, 2, :], data[:, 3, :]]), + np.minimum.reduce( + [data[:, 0, :], data[:, 1, :], data[:, 2, :], data[:, 3, :]] + ), nan=nan_output_value, posinf=nan_output_value, ) @@ -349,7 +430,9 @@ def _linear_interp_frac_coordinates(z, image_coords, return_periods, interpolati raise ValueError("interpolation must have value 'linear', 'max' or 'min") @staticmethod - def _get_coordinates(longitudes, latitudes, crs: str, transform: Affine, pixel_is_area: bool): + def _get_coordinates( + longitudes, latitudes, crs: str, transform: Affine, pixel_is_area: bool + ): if crs.lower() != "epsg:4236": transproj = Transformer.from_crs("epsg:4326", crs, always_xy=True) x, y = transproj.transform(longitudes, latitudes) diff --git a/src/physrisk/hazard_models/core_hazards.py b/src/physrisk/hazard_models/core_hazards.py index f576aef2..f8d16ec0 100644 --- a/src/physrisk/hazard_models/core_hazards.py +++ b/src/physrisk/hazard_models/core_hazards.py @@ -4,7 +4,12 @@ from physrisk.data.hazard_data_provider import HazardDataHint, SourcePath from physrisk.data.inventory import EmbeddedInventory, Inventory from physrisk.kernel import hazards -from physrisk.kernel.hazards import ChronicHeat, CoastalInundation, RiverineInundation, Wind +from physrisk.kernel.hazards import ( + ChronicHeat, + CoastalInundation, + RiverineInundation, + Wind, +) class ResourceSubset: @@ -31,7 +36,9 @@ def with_model_gcm(self, gcm: str): return ResourceSubset(r for r in self.resources if r.indicator_model_gcm == gcm) def with_model_id(self, model_id: str): - return ResourceSubset(r for r in self.resources if r.indicator_model_id == model_id) + return ResourceSubset( + r for r in self.resources if r.indicator_model_id == model_id + ) class ResourceSelector(Protocol): @@ -40,7 +47,12 @@ class ResourceSelector(Protocol): all matches. The selection rule depends on scenario and year.""" def __call__( - self, *, candidates: ResourceSubset, scenario: str, year: int, hint: Optional[HazardDataHint] = None + self, + *, + candidates: ResourceSubset, + scenario: str, + year: int, + hint: Optional[HazardDataHint] = None, ) -> HazardResource: ... @@ -59,81 +71,138 @@ def __init__(self, inventory: Inventory): self._selectors: Dict[ResourceSelectorKey, ResourceSelector] = {} def source_paths(self) -> Dict[type, SourcePath]: - all_hazard_types = list(set(htype for ((htype, _), _) in self._inventory.resources_by_type_id.items())) + all_hazard_types = list( + set( + htype + for ((htype, _), _) in self._inventory.resources_by_type_id.items() + ) + ) source_paths: Dict[type, SourcePath] = {} for hazard_type in all_hazard_types: - source_paths[hazards.hazard_class(hazard_type)] = self._get_resource_source_path( - hazard_type, + source_paths[hazards.hazard_class(hazard_type)] = ( + self._get_resource_source_path( + hazard_type, + ) ) return source_paths - def add_selector(self, hazard_type: type, indicator_id: str, selector: ResourceSelector): + def add_selector( + self, hazard_type: type, indicator_id: str, selector: ResourceSelector + ): self._selectors[ResourceSelectorKey(hazard_type, indicator_id)] = selector def _get_resource_source_path(self, hazard_type: str): - def _get_source_path(*, indicator_id: str, scenario: str, year: int, hint: Optional[HazardDataHint] = None): + def _get_source_path( + *, + indicator_id: str, + scenario: str, + year: int, + hint: Optional[HazardDataHint] = None, + ): # all matching resources in the inventory selector = self._selectors.get( - ResourceSelectorKey(hazard_type=hazards.hazard_class(hazard_type), indicator_id=indicator_id), + ResourceSelectorKey( + hazard_type=hazards.hazard_class(hazard_type), + indicator_id=indicator_id, + ), self._no_selector, ) - resources = self._inventory.resources_by_type_id[(hazard_type, indicator_id)] + resources = self._inventory.resources_by_type_id[ + (hazard_type, indicator_id) + ] if len(resources) == 0: raise RuntimeError( - f"unable to find any resources for hazard {hazard_type} " f"and indicator ID {indicator_id}" + f"unable to find any resources for hazard {hazard_type} " + f"and indicator ID {indicator_id}" ) candidates = ResourceSubset(resources) try: if hint is not None: resource = candidates.match(hint) else: - resource = selector(candidates=candidates, scenario=scenario, year=year, hint=hint) + resource = selector( + candidates=candidates, scenario=scenario, year=year, hint=hint + ) except Exception: raise RuntimeError( - f"unable to select unique resource for hazard {hazard_type} " f"and indicator ID {indicator_id}" + f"unable to select unique resource for hazard {hazard_type} " + f"and indicator ID {indicator_id}" ) proxy_scenario = ( cmip6_scenario_to_rcp(scenario) - if resource.scenarios[0].id.startswith("rcp") or resource.scenarios[-1].id.startswith("rcp") + if resource.scenarios[0].id.startswith("rcp") + or resource.scenarios[-1].id.startswith("rcp") else scenario ) if scenario == "historical": - scenarios = next(iter(s for s in resource.scenarios if s.id == "historical"), None) + scenarios = next( + iter(s for s in resource.scenarios if s.id == "historical"), None + ) if scenarios is None: - scenarios = next(s for s in sorted(resource.scenarios, key=lambda s: next(y for y in s.years))) + scenarios = next( + s + for s in sorted( + resource.scenarios, key=lambda s: next(y for y in s.years) + ) + ) proxy_scenario = scenarios.id year = next(s for s in scenarios.years) - return resource.path.format(id=indicator_id, scenario=proxy_scenario, year=year) + return resource.path.format( + id=indicator_id, scenario=proxy_scenario, year=year + ) return _get_source_path @staticmethod - def _no_selector(candidates: ResourceSubset, scenario: str, year: int, hint: Optional[HazardDataHint] = None): + def _no_selector( + candidates: ResourceSubset, + scenario: str, + year: int, + hint: Optional[HazardDataHint] = None, + ): return candidates.first() class CoreInventorySourcePaths(InventorySourcePaths): def __init__(self, inventory: Inventory): super().__init__(inventory) - for indicator_id in ["mean_work_loss/low", "mean_work_loss/medium", "mean_work_loss/high"]: + for indicator_id in [ + "mean_work_loss/low", + "mean_work_loss/medium", + "mean_work_loss/high", + ]: self.add_selector(ChronicHeat, indicator_id, self._select_chronic_heat) - self.add_selector(ChronicHeat, "mean/degree/days/above/32c", self._select_chronic_heat) - self.add_selector(RiverineInundation, "flood_depth", self._select_riverine_inundation) - self.add_selector(CoastalInundation, "flood_depth", self._select_coastal_inundation) + self.add_selector( + ChronicHeat, "mean/degree/days/above/32c", self._select_chronic_heat + ) + self.add_selector( + RiverineInundation, "flood_depth", self._select_riverine_inundation + ) + self.add_selector( + CoastalInundation, "flood_depth", self._select_coastal_inundation + ) self.add_selector(Wind, "max_speed", self._select_wind) def resources_with(self, *, hazard_type: type, indicator_id: str): - return ResourceSubset(self._inventory.resources_by_type_id[(hazard_type.__name__, indicator_id)]) + return ResourceSubset( + self._inventory.resources_by_type_id[(hazard_type.__name__, indicator_id)] + ) @staticmethod def _select_chronic_heat( - candidates: ResourceSubset, scenario: str, year: int, hint: Optional[HazardDataHint] = None + candidates: ResourceSubset, + scenario: str, + year: int, + hint: Optional[HazardDataHint] = None, ): return candidates.with_model_gcm("ACCESS-CM2").first() @staticmethod def _select_coastal_inundation( - candidates: ResourceSubset, scenario: str, year: int, hint: Optional[HazardDataHint] = None + candidates: ResourceSubset, + scenario: str, + year: int, + hint: Optional[HazardDataHint] = None, ): return ( candidates.with_model_id("nosub").first() @@ -143,7 +212,10 @@ def _select_coastal_inundation( @staticmethod def _select_riverine_inundation( - candidates: ResourceSubset, scenario: str, year: int, hint: Optional[HazardDataHint] = None + candidates: ResourceSubset, + scenario: str, + year: int, + hint: Optional[HazardDataHint] = None, ): return ( candidates.with_model_gcm("historical").first() @@ -152,7 +224,12 @@ def _select_riverine_inundation( ) @staticmethod - def _select_wind(candidates: ResourceSubset, scenario: str, year: int, hint: Optional[HazardDataHint] = None): + def _select_wind( + candidates: ResourceSubset, + scenario: str, + year: int, + hint: Optional[HazardDataHint] = None, + ): return candidates.prefer_group_id("iris_osc").first() diff --git a/src/physrisk/kernel/assets.py b/src/physrisk/kernel/assets.py index 9ff8a3fc..6d5bb30c 100644 --- a/src/physrisk/kernel/assets.py +++ b/src/physrisk/kernel/assets.py @@ -40,7 +40,9 @@ class TurbineKind(Enum): class Asset: - def __init__(self, latitude: float, longitude: float, id: Optional[str] = None, **kwargs): + def __init__( + self, latitude: float, longitude: float, id: Optional[str] = None, **kwargs + ): self.latitude = latitude self.longitude = longitude self.id = id @@ -93,7 +95,14 @@ def __init__( capacity: Optional[float] = None, **kwargs, ): - super().__init__(latitude, longitude, type=type, location=location, capacity=capacity, **kwargs) + super().__init__( + latitude, + longitude, + type=type, + location=location, + capacity=capacity, + **kwargs, + ) self.turbine: Optional[TurbineKind] = None self.cooling: Optional[CoolingKind] = None @@ -117,7 +126,9 @@ def get_inundation_protection_return_period(self): class RealEstateAsset(Asset): - def __init__(self, latitude: float, longitude: float, *, location: str, type: str, **kwargs): + def __init__( + self, latitude: float, longitude: float, *, location: str, type: str, **kwargs + ): super().__init__(latitude, longitude, **kwargs) self.location = location self.type = type @@ -125,7 +136,13 @@ def __init__(self, latitude: float, longitude: float, *, location: str, type: st class ManufacturingAsset(Asset): def __init__( - self, latitude: float, longitude: float, *, location: Optional[str] = None, type: Optional[str] = None, **kwargs + self, + latitude: float, + longitude: float, + *, + location: Optional[str] = None, + type: Optional[str] = None, + **kwargs, ): super().__init__(latitude, longitude, **kwargs) self.location = location @@ -133,7 +150,15 @@ def __init__( class IndustrialActivity(Asset): - def __init__(self, latitude: float, longitude: float, *, location: Optional[str] = None, type: str, **kwargs): + def __init__( + self, + latitude: float, + longitude: float, + *, + location: Optional[str] = None, + type: str, + **kwargs, + ): super().__init__(latitude, longitude, **kwargs) self.location = location self.type = type diff --git a/src/physrisk/kernel/calculation.py b/src/physrisk/kernel/calculation.py index f5ce57e8..4ba2df9f 100644 --- a/src/physrisk/kernel/calculation.py +++ b/src/physrisk/kernel/calculation.py @@ -26,7 +26,13 @@ ThermalPowerGenerationWaterTemperatureModel, ) -from .assets import IndustrialActivity, PowerGeneratingAsset, RealEstateAsset, TestAsset, ThermalPowerGeneratingAsset +from .assets import ( + IndustrialActivity, + PowerGeneratingAsset, + RealEstateAsset, + TestAsset, + ThermalPowerGeneratingAsset, +) from .hazard_model import HazardModel from .vulnerability_model import VulnerabilityModelBase @@ -46,10 +52,16 @@ def get_default_vulnerability_models() -> Dict[type, Sequence[VulnerabilityModel RealEstateRiverineInundationModel(), GenericTropicalCycloneModel(), PlaceholderVulnerabilityModel("fire_probability", Fire, ImpactType.damage), - PlaceholderVulnerabilityModel("days/above/35c", ChronicHeat, ImpactType.damage), + PlaceholderVulnerabilityModel( + "days/above/35c", ChronicHeat, ImpactType.damage + ), PlaceholderVulnerabilityModel("days/above/5cm", Hail, ImpactType.damage), - PlaceholderVulnerabilityModel("months/spei3m/below/-2", Drought, ImpactType.damage), - PlaceholderVulnerabilityModel("max/daily/water_equivalent", Precipitation, ImpactType.damage), + PlaceholderVulnerabilityModel( + "months/spei3m/below/-2", Drought, ImpactType.damage + ), + PlaceholderVulnerabilityModel( + "max/daily/water_equivalent", Precipitation, ImpactType.damage + ), ], PowerGeneratingAsset: [pgam.InundationModel()], RealEstateAsset: [ diff --git a/src/physrisk/kernel/curve.py b/src/physrisk/kernel/curve.py index 08e3d03e..bbeb86f3 100644 --- a/src/physrisk/kernel/curve.py +++ b/src/physrisk/kernel/curve.py @@ -95,7 +95,9 @@ def process_bin_edges_for_graph(bin_edges, range_fraction=0.05): if j >= len(bin_edges): delta = r * range_fraction / (j - i - 1) else: - delta = min(r * range_fraction, 0.25 * (bin_edges[j] - bin_edges[i])) / (j - i - 1) + delta = min(r * range_fraction, 0.25 * (bin_edges[j] - bin_edges[i])) / ( + j - i - 1 + ) offset = delta for k in range(i + 1, j): new_edges[k] = new_edges[k] + offset @@ -119,7 +121,11 @@ class ExceedanceCurve: __slots__ = ["probs", "values"] - def __init__(self, probs: Union[List[float], np.ndarray], values: Union[List[float], np.ndarray]): + def __init__( + self, + probs: Union[List[float], np.ndarray], + values: Union[List[float], np.ndarray], + ): """Create a new asset event distribution. Args: probs: exceedance probabilities (must be sorted and decreasing). @@ -163,11 +169,17 @@ def get_probability_bins(self, include_last: bool = False): value_bins = self.values[:] probs = self.probs[:-1] - self.probs[1:] # type: ignore if include_last or len(self.values) == 1: - value_bins = np.append(value_bins, value_bins[-1]) # last bin has zero width + value_bins = np.append( + value_bins, value_bins[-1] + ) # last bin has zero width probs = np.append(probs, self.probs[-1]) return value_bins, probs def get_samples(self, uniforms): """Return value, v, for each probability p in uniforms such that p is the probability that the random variable < v.""" - return np.where(uniforms > (1.0 - self.probs[0]), np.interp(uniforms, 1.0 - self.probs, self.values), 0.0) + return np.where( + uniforms > (1.0 - self.probs[0]), + np.interp(uniforms, 1.0 - self.probs, self.values), + 0.0, + ) diff --git a/src/physrisk/kernel/events.py b/src/physrisk/kernel/events.py index d010b938..8bf76c97 100644 --- a/src/physrisk/kernel/events.py +++ b/src/physrisk/kernel/events.py @@ -6,14 +6,18 @@ # @njit(cache=True) -def calculate_cumulative_probs(bins_lower: np.ndarray, bins_upper: np.ndarray, probs: np.ndarray): +def calculate_cumulative_probs( + bins_lower: np.ndarray, bins_upper: np.ndarray, probs: np.ndarray +): # note: in some circumstances we could exclude the two extreme points and rely on flat extrapolation # this implementation retains points for clarity, sacrificing some performance assert bins_lower.size == bins_upper.size assert probs.shape[1] == bins_lower.size nb_bins = bins_lower.size # aka M: number of bins n = probs.shape[0] - nb_points = bins_lower.size * 2 - np.count_nonzero(bins_lower[1:] == bins_upper[:-1]) + nb_points = bins_lower.size * 2 - np.count_nonzero( + bins_lower[1:] == bins_upper[:-1] + ) cum_prob = np.zeros(n) values = np.zeros(nb_points) cum_probs = np.zeros(shape=(n, nb_points)) @@ -40,7 +44,9 @@ def calculate_cumulative_probs(bins_lower: np.ndarray, bins_upper: np.ndarray, p @njit(cache=True) -def sample_from_cumulative_probs(values: np.ndarray, cum_probs: np.ndarray, uniforms: np.ndarray): +def sample_from_cumulative_probs( + values: np.ndarray, cum_probs: np.ndarray, uniforms: np.ndarray +): n = cum_probs.shape[0] nb_samples = uniforms.shape[1] assert uniforms.shape[0] == n @@ -57,7 +63,9 @@ def inv_cumulative_marginal_probs(self, cum_probs: np.ndarray): ... class EmpiricalMultivariateDistribution(MultivariateDistribution): """Stores an N dimensional empirical probability density function.""" - def __init__(self, bins_lower: np.ndarray, bins_upper: np.ndarray, probs: np.ndarray): + def __init__( + self, bins_lower: np.ndarray, bins_upper: np.ndarray, probs: np.ndarray + ): """N marginal probability distributions are each represented as a set of bins of uniform probability density. @@ -69,11 +77,17 @@ def __init__(self, bins_lower: np.ndarray, bins_upper: np.ndarray, probs: np.nda Raises: ValueError: _description_ """ - if bins_lower.ndim > 1 or bins_upper.ndim > 1 or bins_lower.size != bins_upper.size: + if ( + bins_lower.ndim > 1 + or bins_upper.ndim > 1 + or bins_lower.size != bins_upper.size + ): raise ValueError("bin upper and lower bounds must be 1-D and same size.") if probs.ndim != 2 or probs.shape[1] != bins_upper.size: raise ValueError("probabilities must be (N, M).") - if np.any(bins_lower[1:] < bins_lower[:-1]) or np.any(bins_upper[1:] < bins_upper[:-1]): + if np.any(bins_lower[1:] < bins_lower[:-1]) or np.any( + bins_upper[1:] < bins_upper[:-1] + ): raise ValueError("bins must be non-decreasing.") if np.any(bins_upper[1:] < bins_lower[:-1]): raise ValueError("bins may not overlap.") @@ -91,13 +105,19 @@ def inv_cumulative_marginal_probs(self, cum_probs: np.ndarray): cum_probs (np.ndarray): Cumulative probabilities (N, P), P being number of samples. axis (int): Specifies the axis of the N events. """ - values_dist, cum_probs_dist = calculate_cumulative_probs(self.bins_lower, self.bins_upper, self.probs) + values_dist, cum_probs_dist = calculate_cumulative_probs( + self.bins_lower, self.bins_upper, self.probs + ) return sample_from_cumulative_probs(values_dist, cum_probs_dist, cum_probs) -def event_samples(impacts_bins: np.ndarray, probs: List[np.ndarray], nb_events: int, nb_samples: int): +def event_samples( + impacts_bins: np.ndarray, probs: List[np.ndarray], nb_events: int, nb_samples: int +): if any([p.size != 1 and p.size != nb_events for p in probs]): - raise ValueError(f"probabilities must be scalar or vector or length {nb_events}.") + raise ValueError( + f"probabilities must be scalar or vector or length {nb_events}." + ) return event_samples_numba(impacts_bins, probs, nb_events, nb_samples) @@ -117,7 +137,9 @@ def find(elements: np.ndarray, value): @njit(cache=True) -def event_samples_numba(impacts_bins: np.ndarray, probs: List[np.ndarray], nb_events: int, nb_samples: int): +def event_samples_numba( + impacts_bins: np.ndarray, probs: List[np.ndarray], nb_events: int, nb_samples: int +): samples = np.zeros(shape=(nb_samples, nb_events)) np.random.seed(111) cum_probs = np.zeros(len(probs)) diff --git a/src/physrisk/kernel/exposure.py b/src/physrisk/kernel/exposure.py index b2e16a2c..c1ac327a 100644 --- a/src/physrisk/kernel/exposure.py +++ b/src/physrisk/kernel/exposure.py @@ -16,7 +16,14 @@ HazardModel, HazardParameterDataResponse, ) -from physrisk.kernel.hazards import ChronicHeat, CombinedInundation, Drought, Fire, Hail, Wind +from physrisk.kernel.hazards import ( + ChronicHeat, + CombinedInundation, + Drought, + Fire, + Hail, + Wind, +) from physrisk.kernel.impact import _request_consolidated from physrisk.kernel.vulnerability_model import DataRequester from physrisk.utils.helpers import get_iterable @@ -56,7 +63,9 @@ class JupterExposureMeasure(ExposureMeasure): def __init__(self): self.exposure_bins = self.get_exposure_bins() - def get_data_requests(self, asset: Asset, *, scenario: str, year: int) -> Iterable[HazardDataRequest]: + def get_data_requests( + self, asset: Asset, *, scenario: str, year: int + ) -> Iterable[HazardDataRequest]: return [ HazardDataRequest( hazard_type, @@ -66,7 +75,9 @@ def get_data_requests(self, asset: Asset, *, scenario: str, year: int) -> Iterab year=year, indicator_id=indicator_id, # select specific model for wind for consistency with thresholds - hint=HazardDataHint(path="wind/jupiter/v1/max_1min_{scenario}_{year}") if hazard_type == Wind else None, + hint=HazardDataHint(path="wind/jupiter/v1/max_1min_{scenario}_{year}") + if hazard_type == Wind + else None, ) for (hazard_type, indicator_id) in self.exposure_bins.keys() ] @@ -157,10 +168,16 @@ def bounds_to_lookup(self, bounds: Iterable[Bounds]): def calculate_exposures( - assets: List[Asset], hazard_model: HazardModel, exposure_measure: ExposureMeasure, scenario: str, year: int + assets: List[Asset], + hazard_model: HazardModel, + exposure_measure: ExposureMeasure, + scenario: str, + year: int, ) -> Dict[Asset, AssetExposureResult]: requester_assets: Dict[DataRequester, List[Asset]] = {exposure_measure: assets} - asset_requests, responses = _request_consolidated(hazard_model, requester_assets, scenario, year) + asset_requests, responses = _request_consolidated( + hazard_model, requester_assets, scenario, year + ) logging.info( "Applying exposure measure {0} to {1} assets of type {2}".format( @@ -170,8 +187,12 @@ def calculate_exposures( result: Dict[Asset, AssetExposureResult] = {} for asset in assets: - requests = asset_requests[(exposure_measure, asset)] # (ordered) requests for a given asset + requests = asset_requests[ + (exposure_measure, asset) + ] # (ordered) requests for a given asset hazard_data = [responses[req] for req in get_iterable(requests)] - result[asset] = AssetExposureResult(hazard_categories=exposure_measure.get_exposures(asset, hazard_data)) + result[asset] = AssetExposureResult( + hazard_categories=exposure_measure.get_exposures(asset, hazard_data) + ) return result diff --git a/src/physrisk/kernel/financial_model.py b/src/physrisk/kernel/financial_model.py index 64dd25b5..bb1494b4 100644 --- a/src/physrisk/kernel/financial_model.py +++ b/src/physrisk/kernel/financial_model.py @@ -14,7 +14,9 @@ def get_asset_value(self, asset: Asset, currency: str) -> float: ... @abstractmethod - def get_asset_aggregate_cashflows(self, asset: Asset, start: datetime, end: datetime, currency: str) -> float: + def get_asset_aggregate_cashflows( + self, asset: Asset, start: datetime, end: datetime, currency: str + ) -> float: """Return the expected sum of the cashflows generated by the Asset between start and end, in specified currency.""" ... @@ -27,7 +29,9 @@ def damage_to_loss(self, asset: Asset, impact: np.ndarray, currency: str): ... @abstractmethod - def disruption_to_loss(self, asset: Asset, impact: np.ndarray, year: int, currency: str): + def disruption_to_loss( + self, asset: Asset, impact: np.ndarray, year: int, currency: str + ): """Convert the fractional annual disruption of the specified asset to a financial loss.""" ... @@ -41,7 +45,9 @@ def __init__(self, data_provider: FinancialDataProvider): def damage_to_loss(self, asset: Asset, impact: np.ndarray, currency: str): return self.data_provider.get_asset_value(asset, currency) * impact - def disruption_to_loss(self, asset: Asset, impact: np.ndarray, year: int, currency: str): + def disruption_to_loss( + self, asset: Asset, impact: np.ndarray, year: int, currency: str + ): return ( self.data_provider.get_asset_aggregate_cashflows( asset, datetime(year, 1, 1), datetime(year, 12, 31), currency @@ -57,7 +63,13 @@ def __init__(self, financial_models: Dict[type, FinancialModelBase]): self.financial_models = financial_models def damage_to_loss(self, asset: Asset, impact: np.ndarray, currency: str): - return self.financial_models[type(asset)].damage_to_loss(asset, impact, currency) + return self.financial_models[type(asset)].damage_to_loss( + asset, impact, currency + ) - def disruption_to_loss(self, asset: Asset, impact: np.ndarray, year: int, currency: str): - return self.financial_models[type(asset)].disruption_to_loss(asset, impact, year, currency) + def disruption_to_loss( + self, asset: Asset, impact: np.ndarray, year: int, currency: str + ): + return self.financial_models[type(asset)].disruption_to_loss( + asset, impact, year, currency + ) diff --git a/src/physrisk/kernel/hazard_model.py b/src/physrisk/kernel/hazard_model.py index 40bbc958..48c66b97 100644 --- a/src/physrisk/kernel/hazard_model.py +++ b/src/physrisk/kernel/hazard_model.py @@ -73,7 +73,11 @@ class HazardEventDataResponse(HazardDataResponse): """Response to HazardDataRequest for acute hazards.""" def __init__( - self, return_periods: np.ndarray, intensities: np.ndarray, units: str = "default", path: str = "unknown" + self, + return_periods: np.ndarray, + intensities: np.ndarray, + units: str = "default", + path: str = "unknown", ): """Create HazardEventDataResponse. @@ -127,7 +131,9 @@ def parameter(self) -> float: class HazardModelFactory(Protocol): - def hazard_model(self, interpolation: str = "floor", provider_max_requests: Dict[str, int] = {}): + def hazard_model( + self, interpolation: str = "floor", provider_max_requests: Dict[str, int] = {} + ): """Create a HazardModel instance based on a number of options. Args: @@ -144,7 +150,9 @@ class HazardModel(ABC): EventDataResponses.""" @abstractmethod - def get_hazard_events(self, requests: List[HazardDataRequest]) -> Mapping[HazardDataRequest, HazardDataResponse]: + def get_hazard_events( + self, requests: List[HazardDataRequest] + ) -> Mapping[HazardDataRequest, HazardDataResponse]: """Process the hazard data requests and return responses.""" ... @@ -161,7 +169,9 @@ class CompositeHazardModel(HazardModel): def __init__(self, hazard_models: Dict[type, HazardModel]): self.hazard_models = hazard_models - def get_hazard_events(self, requests: List[HazardDataRequest]) -> Mapping[HazardDataRequest, HazardDataResponse]: + def get_hazard_events( + self, requests: List[HazardDataRequest] + ) -> Mapping[HazardDataRequest, HazardDataResponse]: requests_by_event_type = defaultdict(list) for request in requests: diff --git a/src/physrisk/kernel/hazards.py b/src/physrisk/kernel/hazards.py index 5056928e..4c8ed8ae 100644 --- a/src/physrisk/kernel/hazards.py +++ b/src/physrisk/kernel/hazards.py @@ -25,7 +25,11 @@ def hazard_kind(hazard_type: Type[Hazard]): def indicator_data(hazard_type: Type[Hazard], indicator_id: str): - default = IndicatorData.EVENT if hazard_type.kind == HazardKind.ACUTE else IndicatorData.PARAMETERS + default = ( + IndicatorData.EVENT + if hazard_type.kind == HazardKind.ACUTE + else IndicatorData.PARAMETERS + ) return hazard_type.indicator_data.get(indicator_id, default) @@ -36,7 +40,10 @@ class ChronicHeat(Hazard): class Inundation(Hazard): kind = HazardKind.ACUTE - indicator_data = {"flood_depth": IndicatorData.EVENT, "flood_sop": IndicatorData.PARAMETERS} + indicator_data = { + "flood_depth": IndicatorData.EVENT, + "flood_sop": IndicatorData.PARAMETERS, + } pass @@ -107,7 +114,9 @@ class Subsidence(Hazard): def all_hazards(): return [ - obj for _, obj in inspect.getmembers(sys.modules[__name__]) if inspect.isclass(obj) and issubclass(obj, Hazard) + obj + for _, obj in inspect.getmembers(sys.modules[__name__]) + if inspect.isclass(obj) and issubclass(obj, Hazard) ] diff --git a/src/physrisk/kernel/impact.py b/src/physrisk/kernel/impact.py index 0a03a336..91db4faa 100644 --- a/src/physrisk/kernel/impact.py +++ b/src/physrisk/kernel/impact.py @@ -5,7 +5,12 @@ from physrisk.kernel.assets import Asset from physrisk.kernel.hazard_event_distrib import HazardEventDistrib -from physrisk.kernel.hazard_model import HazardDataFailedResponse, HazardDataRequest, HazardDataResponse, HazardModel +from physrisk.kernel.hazard_model import ( + HazardDataFailedResponse, + HazardDataRequest, + HazardDataResponse, + HazardModel, +) from physrisk.kernel.impact_distrib import ImpactDistrib from physrisk.kernel.vulnerability_distrib import VulnerabilityDistrib from physrisk.kernel.vulnerability_model import ( @@ -41,7 +46,9 @@ class AssetImpactResult: impact: ImpactDistrib vulnerability: Optional[VulnerabilityDistrib] = None event: Optional[HazardEventDistrib] = None - hazard_data: Optional[Iterable[HazardDataResponse]] = None # optional detailed results for drill-down + hazard_data: Optional[Iterable[HazardDataResponse]] = ( + None # optional detailed results for drill-down + ) def calculate_impacts( # noqa: C901 @@ -65,7 +72,9 @@ def calculate_impacts( # noqa: C901 model_assets[mapping].append(asset) results: Dict[ImpactKey, List[AssetImpactResult]] = {} - asset_requests, responses = _request_consolidated(hazard_model, model_assets, scenario, year) + asset_requests, responses = _request_consolidated( + hazard_model, model_assets, scenario, year + ) logging.info("Calculating impacts") for model, assets in model_assets.items(): @@ -78,20 +87,52 @@ def calculate_impacts( # noqa: C901 requests = asset_requests[(model, asset)] hazard_data = [responses[req] for req in get_iterable(requests)] assert isinstance(model, VulnerabilityModelBase) - if ImpactKey(asset=asset, hazard_type=model.hazard_type, scenario=scenario, key_year=year) not in results: - results[ImpactKey(asset=asset, hazard_type=model.hazard_type, scenario=scenario, key_year=year)] = [] + if ( + ImpactKey( + asset=asset, + hazard_type=model.hazard_type, + scenario=scenario, + key_year=year, + ) + not in results + ): + results[ + ImpactKey( + asset=asset, + hazard_type=model.hazard_type, + scenario=scenario, + key_year=year, + ) + ] = [] if any(isinstance(hd, HazardDataFailedResponse) for hd in hazard_data): continue try: if isinstance(model, VulnerabilityModelAcuteBase): impact, vul, event = model.get_impact_details(asset, hazard_data) results[ - ImpactKey(asset=asset, hazard_type=model.hazard_type, scenario=scenario, key_year=year) - ].append(AssetImpactResult(impact, vulnerability=vul, event=event, hazard_data=hazard_data)) + ImpactKey( + asset=asset, + hazard_type=model.hazard_type, + scenario=scenario, + key_year=year, + ) + ].append( + AssetImpactResult( + impact, + vulnerability=vul, + event=event, + hazard_data=hazard_data, + ) + ) elif isinstance(model, VulnerabilityModelBase): impact = model.get_impact(asset, hazard_data) results[ - ImpactKey(asset=asset, hazard_type=model.hazard_type, scenario=scenario, key_year=year) + ImpactKey( + asset=asset, + hazard_type=model.hazard_type, + scenario=scenario, + key_year=year, + ) ].append(AssetImpactResult(impact, hazard_data=hazard_data)) except Exception as e: logger.exception(e) @@ -99,21 +140,31 @@ def calculate_impacts( # noqa: C901 def _request_consolidated( - hazard_model: HazardModel, requester_assets: Dict[DataRequester, List[Asset]], scenario: str, year: int + hazard_model: HazardModel, + requester_assets: Dict[DataRequester, List[Asset]], + scenario: str, + year: int, ): """As an important performance optimization, data requests are consolidated for all requesters (e.g. vulnerability model) because different requesters may query the same hazard data sets note that key for a single request is (requester, asset). """ # the list of requests for each requester and asset - asset_requests: Dict[Tuple[DataRequester, Asset], Union[HazardDataRequest, Iterable[HazardDataRequest]]] = {} + asset_requests: Dict[ + Tuple[DataRequester, Asset], + Union[HazardDataRequest, Iterable[HazardDataRequest]], + ] = {} logging.info("Generating hazard data requests for requesters") for requester, assets in requester_assets.items(): for asset in assets: - asset_requests[(requester, asset)] = requester.get_data_requests(asset, scenario=scenario, year=year) + asset_requests[(requester, asset)] = requester.get_data_requests( + asset, scenario=scenario, year=year + ) logging.info("Retrieving hazard data") - flattened_requests = [req for requests in asset_requests.values() for req in get_iterable(requests)] + flattened_requests = [ + req for requests in asset_requests.values() for req in get_iterable(requests) + ] responses = hazard_model.get_hazard_events(flattened_requests) return asset_requests, responses diff --git a/src/physrisk/kernel/impact_distrib.py b/src/physrisk/kernel/impact_distrib.py index b16072d3..d6627124 100644 --- a/src/physrisk/kernel/impact_distrib.py +++ b/src/physrisk/kernel/impact_distrib.py @@ -41,7 +41,9 @@ def impact_bins_explicit(self): return zip(self.__impact_bins[0:-1], self.__impact_bins[1:]) def mean_impact(self): - return np.sum((self.__impact_bins[:-1] + self.__impact_bins[1:]) * self.__prob / 2) + return np.sum( + (self.__impact_bins[:-1] + self.__impact_bins[1:]) * self.__prob / 2 + ) def stddev_impact(self): mean = self.mean_impact() @@ -56,10 +58,16 @@ def above_mean_stddev_impact(self): return 0.0 if len(above_mean_bins) == len(self.__prob): return self.stddev_impact() - above_mean_probs = self.__prob[-len(above_mean_bins) :] / np.sum(self.__prob[-len(above_mean_bins) :]) + above_mean_probs = self.__prob[-len(above_mean_bins) :] / np.sum( + self.__prob[-len(above_mean_bins) :] + ) above_mean_mean = np.sum(above_mean_bins * above_mean_probs / 2) return np.sqrt( - np.sum(above_mean_probs * (above_mean_bins - above_mean_mean) * (above_mean_bins - above_mean_mean)) + np.sum( + above_mean_probs + * (above_mean_bins - above_mean_mean) + * (above_mean_bins - above_mean_mean) + ) ) def to_exceedance_curve(self): diff --git a/src/physrisk/kernel/risk.py b/src/physrisk/kernel/risk.py index d2ee115d..dfcd1d4c 100644 --- a/src/physrisk/kernel/risk.py +++ b/src/physrisk/kernel/risk.py @@ -1,6 +1,17 @@ import concurrent.futures from dataclasses import dataclass -from typing import Dict, List, NamedTuple, Optional, Protocol, Sequence, Set, Tuple, Type, Union +from typing import ( + Dict, + List, + NamedTuple, + Optional, + Protocol, + Sequence, + Set, + Tuple, + Type, + Union, +) from physrisk.api.v1.impact_req_resp import Category, ScoreBasedRiskMeasureDefinition from physrisk.kernel.assets import Asset @@ -25,26 +36,36 @@ class RiskModel: """Base class for a risk model (i.e. a calculation of risk that makes use of hazard and vulnerability models).""" - def __init__(self, hazard_model: HazardModel, vulnerability_models: VulnerabilityModels): + def __init__( + self, hazard_model: HazardModel, vulnerability_models: VulnerabilityModels + ): self._hazard_model = hazard_model self._vulnerability_models = vulnerability_models - def calculate_risk_measures(self, assets: Sequence[Asset], prosp_scens: Sequence[str], years: Sequence[int]): ... + def calculate_risk_measures( + self, assets: Sequence[Asset], prosp_scens: Sequence[str], years: Sequence[int] + ): ... def _calculate_all_impacts( - self, assets: Sequence[Asset], prosp_scens: Sequence[str], years: Sequence[int], include_histo: bool = False + self, + assets: Sequence[Asset], + prosp_scens: Sequence[str], + years: Sequence[int], + include_histo: bool = False, ): # ensure "historical" is present, e.g. needed for risk measures - scenarios = set(["historical"] + list(prosp_scens)) if include_histo else prosp_scens + scenarios = ( + set(["historical"] + list(prosp_scens)) if include_histo else prosp_scens + ) impact_results: Dict[ImpactKey, List[AssetImpactResult]] = {} # in case of multiple calculation, run on separate threads with concurrent.futures.ThreadPoolExecutor(max_workers=8) as executor: # with concurrent.futures.ProcessPoolExecutor(max_workers=8) as executor: tagged_futures = { - executor.submit(self._calculate_single_impact, assets, scenario, year): BatchId( - scenario, None if scenario == "historical" else year - ) + executor.submit( + self._calculate_single_impact, assets, scenario, year + ): BatchId(scenario, None if scenario == "historical" else year) for scenario in scenarios for year in years } @@ -66,9 +87,17 @@ def _calculate_all_impacts( print("%r generated an exception: %s" % (tag, exc)) return impact_results - def _calculate_single_impact(self, assets: Sequence[Asset], scenario: str, year: int): + def _calculate_single_impact( + self, assets: Sequence[Asset], scenario: str, year: int + ): """Calculate impacts for a single scenario and year.""" - return calculate_impacts(assets, self._hazard_model, self._vulnerability_models, scenario=scenario, year=year) + return calculate_impacts( + assets, + self._hazard_model, + self._vulnerability_models, + scenario=scenario, + year=year, + ) class MeasureKey(NamedTuple): @@ -87,10 +116,15 @@ class Measure: class RiskMeasureCalculator(Protocol): def calc_measure( - self, hazard_type: Type[Hazard], base_impact: AssetImpactResult, impact: AssetImpactResult + self, + hazard_type: Type[Hazard], + base_impact: AssetImpactResult, + impact: AssetImpactResult, ) -> Optional[Measure]: ... - def get_definition(self, hazard_type: Type[Hazard]) -> ScoreBasedRiskMeasureDefinition: ... + def get_definition( + self, hazard_type: Type[Hazard] + ) -> ScoreBasedRiskMeasureDefinition: ... def supported_hazards(self) -> Set[type]: ... @@ -118,19 +152,26 @@ def __init__( super().__init__(hazard_model, vulnerability_models) self._measure_calculators = measure_calculators - def calculate_impacts(self, assets: Sequence[Asset], prosp_scens: Sequence[str], years: Sequence[int]): + def calculate_impacts( + self, assets: Sequence[Asset], prosp_scens: Sequence[str], years: Sequence[int] + ): impacts = self._calculate_all_impacts(assets, prosp_scens, years) return impacts def populate_measure_definitions( self, assets: Sequence[Asset] - ) -> Tuple[Dict[Type[Hazard], List[str]], Dict[ScoreBasedRiskMeasureDefinition, str]]: + ) -> Tuple[ + Dict[Type[Hazard], List[str]], Dict[ScoreBasedRiskMeasureDefinition, str] + ]: hazards = all_hazards() # the identifiers of the score-based risk measures used for each asset, for each hazard type measure_ids_for_hazard: Dict[Type[Hazard], List[str]] = {} # one calcs_by_asset = [ - self._measure_calculators.get(type(asset), self._measure_calculators.get(Asset, None)) for asset in assets + self._measure_calculators.get( + type(asset), self._measure_calculators.get(Asset, None) + ) + for asset in assets ] # match to specific asset and if no match then use the generic calculator assigned to Asset used_calcs = {c for c in calcs_by_asset if c is not None} @@ -141,14 +182,18 @@ def populate_measure_definitions( set( item for item in ( - cal.get_definition(hazard_type=hazard_type) for hazard_type in hazards for cal in used_calcs + cal.get_definition(hazard_type=hazard_type) + for hazard_type in hazards + for cal in used_calcs ) if item is not None ) ) } - def get_measure_id(measure_calc: Union[RiskMeasureCalculator, None], hazard_type: type): + def get_measure_id( + measure_calc: Union[RiskMeasureCalculator, None], hazard_type: type + ): if measure_calc is None: return "na" measure = measure_calc.get_definition(hazard_type=hazard_type) @@ -159,8 +204,12 @@ def get_measure_id(measure_calc: Union[RiskMeasureCalculator, None], hazard_type measure_ids_for_hazard[hazard_type] = measure_ids return measure_ids_for_hazard, measure_id_lookup - def calculate_risk_measures(self, assets: Sequence[Asset], prosp_scens: Sequence[str], years: Sequence[int]): - impacts = self._calculate_all_impacts(assets, prosp_scens, years, include_histo=True) + def calculate_risk_measures( + self, assets: Sequence[Asset], prosp_scens: Sequence[str], years: Sequence[int] + ): + impacts = self._calculate_all_impacts( + assets, prosp_scens, years, include_histo=True + ) measures: Dict[MeasureKey, Measure] = {} for asset in assets: @@ -171,17 +220,35 @@ def calculate_risk_measures(self, assets: Sequence[Asset], prosp_scens: Sequence for year in years: for hazard_type in measure_calc.supported_hazards(): base_impacts = impacts.get( - ImpactKey(asset=asset, hazard_type=hazard_type, scenario="historical", key_year=None) + ImpactKey( + asset=asset, + hazard_type=hazard_type, + scenario="historical", + key_year=None, + ) ) prosp_impacts = impacts.get( - ImpactKey(asset=asset, hazard_type=hazard_type, scenario=prosp_scen, key_year=year) + ImpactKey( + asset=asset, + hazard_type=hazard_type, + scenario=prosp_scen, + key_year=year, + ) ) risk_inds = [ - measure_calc.calc_measure(hazard_type, base_impact, prosp_impact) - for base_impact, prosp_impact in zip(base_impacts, prosp_impacts) + measure_calc.calc_measure( + hazard_type, base_impact, prosp_impact + ) + for base_impact, prosp_impact in zip( + base_impacts, prosp_impacts + ) + ] + risk_ind = [ + risk_ind for risk_ind in risk_inds if risk_ind is not None ] - risk_ind = [risk_ind for risk_ind in risk_inds if risk_ind is not None] if 0 < len(risk_ind): # TODO: Aggregate measures instead of picking the first value. - measures[MeasureKey(asset, prosp_scen, year, hazard_type)] = risk_ind[0] + measures[ + MeasureKey(asset, prosp_scen, year, hazard_type) + ] = risk_ind[0] return impacts, measures diff --git a/src/physrisk/kernel/vulnerability_model.py b/src/physrisk/kernel/vulnerability_model.py index eecd23d8..6886d076 100644 --- a/src/physrisk/kernel/vulnerability_model.py +++ b/src/physrisk/kernel/vulnerability_model.py @@ -1,7 +1,17 @@ import importlib.resources import json from abc import ABC, abstractmethod -from typing import Dict, Iterable, List, Optional, Protocol, Sequence, Tuple, Type, Union +from typing import ( + Dict, + Iterable, + List, + Optional, + Protocol, + Sequence, + Tuple, + Type, + Union, +) import numpy as np from scipy import stats @@ -43,7 +53,9 @@ def decorator_events(func): def get_vulnerability_curves_from_resource(id: str) -> VulnerabilityCurves: - with importlib.resources.open_text(physrisk.data.static.vulnerability, id + ".json") as f: + with importlib.resources.open_text( + physrisk.data.static.vulnerability, id + ".json" + ) as f: curve_set = VulnerabilityCurves(**json.load(f)) return curve_set @@ -72,7 +84,9 @@ def get_data_requests( class EventBased(Protocol): - def impact_samples(self, asset: Asset, data_responses: Iterable[HazardDataResponse]) -> np.ndarray: + def impact_samples( + self, asset: Asset, data_responses: Iterable[HazardDataResponse] + ) -> np.ndarray: # event-based models generate impact samples based on events received by the hazard model # the events may be in the form of an array of severities in the form of return periods. ... @@ -95,11 +109,15 @@ def get_data_requests( ... @abstractmethod - def get_impact(self, asset: Asset, event_data: List[HazardDataResponse]) -> ImpactDistrib: ... + def get_impact( + self, asset: Asset, event_data: List[HazardDataResponse] + ) -> ImpactDistrib: ... class VulnerabilityModels(Protocol): - def vuln_model_for_asset_of_type(self, type: Type[Asset]) -> Sequence[VulnerabilityModelBase]: + def vuln_model_for_asset_of_type( + self, type: Type[Asset] + ) -> Sequence[VulnerabilityModelBase]: """Returns for a given asset type the vulnerability models for each hazard required. Returns: @@ -134,7 +152,9 @@ class VulnerabilityModelAcuteBase(VulnerabilityModelBase): """ def __init__(self, indicator_id: str, hazard_type: type, impact_type: ImpactType): - super().__init__(indicator_id=indicator_id, hazard_type=hazard_type, impact_type=impact_type) + super().__init__( + indicator_id=indicator_id, hazard_type=hazard_type, impact_type=impact_type + ) @abstractmethod def get_distributions( @@ -178,7 +198,9 @@ def get_impact_details( def _check_event_type(self): if self.hazard_type not in self._event_types: - raise NotImplementedError(f"model does not support events of type {self.hazard_type.__name__}") + raise NotImplementedError( + f"model does not support events of type {self.hazard_type.__name__}" + ) class VulnerabilityModel(VulnerabilityModelAcuteBase): @@ -217,7 +239,9 @@ def get_distributions( (event_data,) = event_data_responses assert isinstance(event_data, HazardEventDataResponse) - intensity_curve = ExceedanceCurve(1.0 / event_data.return_periods, event_data.intensities) + intensity_curve = ExceedanceCurve( + 1.0 / event_data.return_periods, event_data.intensities + ) intensity_bin_edges, probs = intensity_curve.get_probability_bins() intensity_bin_centres = (intensity_bin_edges[1:] + intensity_bin_edges[:-1]) / 2 @@ -226,27 +250,41 @@ def get_distributions( intensity_bin_edges, self.impact_bin_edges, # np.eye(8, 11) - self.get_impact_curve(intensity_bin_centres, asset).to_prob_matrix(self.impact_bin_edges), + self.get_impact_curve(intensity_bin_centres, asset).to_prob_matrix( + self.impact_bin_edges + ), ) - event = HazardEventDistrib(self.hazard_type, intensity_bin_edges, probs, [event_data.path]) + event = HazardEventDistrib( + self.hazard_type, intensity_bin_edges, probs, [event_data.path] + ) return vul, event @abstractmethod - def get_impact_curve(self, intensity_bin_centres: np.ndarray, asset: Asset) -> VulnMatrixProvider: + def get_impact_curve( + self, intensity_bin_centres: np.ndarray, asset: Asset + ) -> VulnMatrixProvider: """Defines a VulnMatrixProvider. The VulnMatrixProvider returns probabilities of specified impact bins for the intensity bin centres.""" ... class CurveBasedVulnerabilityModel(VulnerabilityModel): - def get_impact_curve(self, intensity_bin_centres: np.ndarray, asset: Asset) -> VulnMatrixProvider: + def get_impact_curve( + self, intensity_bin_centres: np.ndarray, asset: Asset + ) -> VulnMatrixProvider: curve: VulnerabilityCurve = self.get_vulnerability_curve(asset) - impact_means = np.interp(intensity_bin_centres, curve.intensity, curve.impact_mean) - impact_stddevs = np.interp(intensity_bin_centres, curve.intensity, curve.impact_std) + impact_means = np.interp( + intensity_bin_centres, curve.intensity, curve.impact_mean + ) + impact_stddevs = np.interp( + intensity_bin_centres, curve.intensity, curve.impact_std + ) return VulnMatrixProvider( intensity_bin_centres, - impact_cdfs=[checked_beta_distrib(m, s) for m, s in zip(impact_means, impact_stddevs)], + impact_cdfs=[ + checked_beta_distrib(m, s) for m, s in zip(impact_means, impact_stddevs) + ], ) @abstractmethod @@ -279,7 +317,9 @@ def __init__( indicator_id (str): ID of the hazard indicator to which this applies. Defaults to "". buffer (Optional[int]): Delimitation of the area for the hazard data in metres (within [0,1000]). """ - super().__init__(indicator_id=indicator_id, hazard_type=hazard_type, impact_type=impact_type) + super().__init__( + indicator_id=indicator_id, hazard_type=hazard_type, impact_type=impact_type + ) self.damage_curve_intensities = damage_curve_intensities self.damage_curve_impacts = damage_curve_impacts self.buffer = buffer @@ -303,17 +343,28 @@ def get_distributions( (event_data,) = event_data_responses assert isinstance(event_data, HazardEventDataResponse) - intensity_curve = ExceedanceCurve(1.0 / event_data.return_periods, event_data.intensities) + intensity_curve = ExceedanceCurve( + 1.0 / event_data.return_periods, event_data.intensities + ) intensity_bin_edges, probs = intensity_curve.get_probability_bins() # look up the impact bin edges - impact_bins_edges = np.interp(intensity_bin_edges, self.damage_curve_intensities, self.damage_curve_impacts) + impact_bins_edges = np.interp( + intensity_bin_edges, + self.damage_curve_intensities, + self.damage_curve_impacts, + ) # the vulnerability distribution probabilities are an identity matrix: # we assume that if the intensity falls within a certain bin then the impacts *will* fall within the # bin where the edges are obtained by applying the damage curve to the intensity bin edges. vul = VulnerabilityDistrib( - type(self.hazard_type), intensity_bin_edges, impact_bins_edges, np.eye(len(impact_bins_edges) - 1) + type(self.hazard_type), + intensity_bin_edges, + impact_bins_edges, + np.eye(len(impact_bins_edges) - 1), + ) + event = HazardEventDistrib( + type(self.hazard_type), intensity_bin_edges, probs, [event_data.path] ) - event = HazardEventDistrib(type(self.hazard_type), intensity_bin_edges, probs, [event_data.path]) return vul, event diff --git a/src/physrisk/requests.py b/src/physrisk/requests.py index 854cf790..0b50e6ad 100644 --- a/src/physrisk/requests.py +++ b/src/physrisk/requests.py @@ -8,7 +8,12 @@ import physrisk.data.static.example_portfolios from physrisk.api.v1.common import Distribution, ExceedanceCurve, VulnerabilityDistrib -from physrisk.api.v1.exposure_req_resp import AssetExposure, AssetExposureRequest, AssetExposureResponse, Exposure +from physrisk.api.v1.exposure_req_resp import ( + AssetExposure, + AssetExposureRequest, + AssetExposureResponse, + Exposure, +) from physrisk.api.v1.hazard_image import HazardImageRequest from physrisk.data.hazard_data_provider import HazardDataHint from physrisk.data.inventory import expand @@ -19,7 +24,13 @@ from physrisk.kernel.hazards import all_hazards from physrisk.kernel.impact import AssetImpactResult, ImpactKey # , ImpactKey from physrisk.kernel.impact_distrib import EmptyImpactDistrib -from physrisk.kernel.risk import AssetLevelRiskModel, Measure, MeasureKey, RiskMeasureCalculator, RiskMeasuresFactory +from physrisk.kernel.risk import ( + AssetLevelRiskModel, + Measure, + MeasureKey, + RiskMeasureCalculator, + RiskMeasuresFactory, +) from physrisk.kernel.vulnerability_model import ( DictBasedVulnerabilityModels, VulnerabilityModels, @@ -60,7 +71,11 @@ from .kernel import calculation as calc from .kernel.hazard_model import HazardDataRequest as hmHazardDataRequest from .kernel.hazard_model import HazardEventDataResponse as hmHazardEventDataResponse -from .kernel.hazard_model import HazardModel, HazardModelFactory, HazardParameterDataResponse +from .kernel.hazard_model import ( + HazardModel, + HazardModelFactory, + HazardParameterDataResponse, +) Colormaps = Dict[str, Any] @@ -90,30 +105,47 @@ def get(self, *, request_id, request_dict): if request_id == "get_hazard_data": request = HazardDataRequest(**request_dict) hazard_model = self.hazard_model_factory.hazard_model( - interpolation=request.interpolation, provider_max_requests=request.provider_max_requests + interpolation=request.interpolation, + provider_max_requests=request.provider_max_requests, ) - return json.dumps(_get_hazard_data(request, hazard_model=hazard_model).model_dump()) # , allow_nan=False) + return json.dumps( + _get_hazard_data(request, hazard_model=hazard_model).model_dump() + ) # , allow_nan=False) elif request_id == "get_hazard_data_availability": request = HazardAvailabilityRequest(**request_dict) - return json.dumps(_get_hazard_data_availability(request, self.inventory, self.colormaps).model_dump()) + return json.dumps( + _get_hazard_data_availability( + request, self.inventory, self.colormaps + ).model_dump() + ) elif request_id == "get_hazard_data_description": request = HazardDescriptionRequest(**request_dict) return json.dumps(_get_hazard_data_description(request).dict()) elif request_id == "get_asset_exposure": request = AssetExposureRequest(**request_dict) hazard_model = self.hazard_model_factory.hazard_model( - interpolation=request.calc_settings.hazard_interp, provider_max_requests=request.provider_max_requests + interpolation=request.calc_settings.hazard_interp, + provider_max_requests=request.provider_max_requests, + ) + return json.dumps( + _get_asset_exposures(request, hazard_model).model_dump( + exclude_none=True + ) ) - return json.dumps(_get_asset_exposures(request, hazard_model).model_dump(exclude_none=True)) elif request_id == "get_asset_impact": request = AssetImpactRequest(**request_dict) hazard_model = self.hazard_model_factory.hazard_model( - interpolation=request.calc_settings.hazard_interp, provider_max_requests=request.provider_max_requests + interpolation=request.calc_settings.hazard_interp, + provider_max_requests=request.provider_max_requests, + ) + vulnerability_models = ( + self.vulnerability_models_factory.vulnerability_models() ) - vulnerability_models = self.vulnerability_models_factory.vulnerability_models() measure_calculators = self.measures_factory.calculators(request.use_case_id) return dumps( - _get_asset_impacts(request, hazard_model, vulnerability_models, measure_calculators).model_dump() + _get_asset_impacts( + request, hazard_model, vulnerability_models, measure_calculators + ).model_dump() ) elif request_id == "get_example_portfolios": return dumps(_get_example_portfolios()) @@ -124,7 +156,9 @@ def get_image(self, *, request_dict): inventory = self.inventory zarr_reader = self.zarr_reader request = HazardImageRequest(**request_dict) - if not _read_permitted(request.group_ids, inventory.resources[request.resource]): + if not _read_permitted( + request.group_ids, inventory.resources[request.resource] + ): raise PermissionError() model = inventory.resources[request.resource] len(PosixPath(model.map.path).parts) @@ -133,14 +167,24 @@ def get_image(self, *, request_dict): if len(PosixPath(model.map.path).parts) == 1 else model.map.path ).format(scenario=request.scenario_id, year=request.year) - colormap = request.colormap if request.colormap is not None else model.map.colormap.name + colormap = ( + request.colormap + if request.colormap is not None + else model.map.colormap.name + ) creator = ImageCreator(zarr_reader) # store=ImageCreator.test_store(path)) return creator.convert( - path, colormap=colormap, tile=request.tile, min_value=request.min_value, max_value=request.max_value + path, + colormap=colormap, + tile=request.tile, + min_value=request.min_value, + max_value=request.max_value, ) -def _create_inventory(reader: Optional[InventoryReader] = None, sources: Optional[List[str]] = None): +def _create_inventory( + reader: Optional[InventoryReader] = None, sources: Optional[List[str]] = None +): resources: List[HazardResource] = [] colormaps: Dict[str, Dict[str, Any]] = {} request_sources = ["embedded"] if sources is None else [s.lower() for s in sources] @@ -185,14 +229,18 @@ def _read_permitted(group_ids: List[str], resource: HazardResource): return ("osc" in group_ids) or resource.group_id == "public" -def _get_hazard_data_availability(request: HazardAvailabilityRequest, inventory: Inventory, colormaps: dict): +def _get_hazard_data_availability( + request: HazardAvailabilityRequest, inventory: Inventory, colormaps: dict +): response = HazardAvailabilityResponse( models=list(inventory.resources.values()), colormaps=colormaps ) # type: ignore return response -def _get_hazard_data_description(request: HazardDescriptionRequest, reader: InventoryReader): +def _get_hazard_data_description( + request: HazardDescriptionRequest, reader: InventoryReader +): descriptions = reader.read_description_markdown(request.paths) return HazardDescriptionResponse(descriptions=descriptions) @@ -207,21 +255,33 @@ def _get_hazard_data(request: HazardDataRequest, hazard_model: HazardModel): # get hazard event types: event_types = Hazard.__subclasses__() event_dict = dict((et.__name__, et) for et in event_types) - event_dict.update((est.__name__, est) for et in event_types for est in et.__subclasses__()) + event_dict.update( + (est.__name__, est) for et in event_types for est in et.__subclasses__() + ) # flatten list to let event processor decide how to group item_requests = [] all_requests = [] for item in request.items: hazard_type = ( - item.hazard_type if item.hazard_type is not None else item.event_type if item.event_type is not None else "" + item.hazard_type + if item.hazard_type is not None + else item.event_type + if item.event_type is not None + else "" ) event_type = event_dict[hazard_type] hint = None if item.path is None else HazardDataHint(path=item.path) data_requests = [ hmHazardDataRequest( - event_type, lon, lat, indicator_id=item.indicator_id, scenario=item.scenario, year=item.year, hint=hint + event_type, + lon, + lat, + indicator_id=item.indicator_id, + scenario=item.scenario, + year=item.year, + hint=hint, ) for (lon, lat) in zip(item.longitudes, item.latitudes) ] @@ -255,7 +315,12 @@ def _get_hazard_data(request: HazardDataRequest, hazard_model: HazardModel): return_periods=[], ) if isinstance(resp, HazardParameterDataResponse) - else IntensityCurve(intensities=[], index_values=[], index_name="", return_periods=[]) + else IntensityCurve( + intensities=[], + index_values=[], + index_name="", + return_periods=[], + ) ) ) for resp in resps @@ -278,7 +343,9 @@ def create_assets(api_assets: Assets, assets: Optional[List[Asset]] = None): # """Create list of Asset objects from the Assets API object:""" if assets is not None: if len(api_assets.items) != 0: - raise ValueError("Cannot provide asset items in the request while specifying an explicit asset list") + raise ValueError( + "Cannot provide asset items in the request while specifying an explicit asset list" + ) return assets else: module = import_module("physrisk.kernel.assets") @@ -302,11 +369,15 @@ def create_assets(api_assets: Assets, assets: Optional[List[Asset]] = None): # def _get_asset_exposures( - request: AssetExposureRequest, hazard_model: HazardModel, assets: Optional[List[Asset]] = None + request: AssetExposureRequest, + hazard_model: HazardModel, + assets: Optional[List[Asset]] = None, ): _assets = create_assets(request.assets, assets) measure = JupterExposureMeasure() - results = calculate_exposures(_assets, hazard_model, measure, scenario="ssp585", year=2030) + results = calculate_exposures( + _assets, hazard_model, measure, scenario="ssp585", year=2030 + ) return AssetExposureResponse( items=[ AssetExposure( @@ -337,23 +408,43 @@ def _get_asset_impacts( # based on asset_class: _assets = create_assets(request.assets, assets) measure_calculators = ( - calc.get_default_risk_measure_calculators() if measure_calculators is None else measure_calculators + calc.get_default_risk_measure_calculators() + if measure_calculators is None + else measure_calculators + ) + risk_model = AssetLevelRiskModel( + hazard_model, vulnerability_models, measure_calculators ) - risk_model = AssetLevelRiskModel(hazard_model, vulnerability_models, measure_calculators) - scenarios = [request.scenario] if request.scenarios is None or len(request.scenarios) == 0 else request.scenarios - years = [request.year] if request.years is None or len(request.years) == 0 else request.years + scenarios = ( + [request.scenario] + if request.scenarios is None or len(request.scenarios) == 0 + else request.scenarios + ) + years = ( + [request.year] + if request.years is None or len(request.years) == 0 + else request.years + ) risk_measures = None if request.include_measures: - impacts, measures = risk_model.calculate_risk_measures(_assets, scenarios, years) - measure_ids_for_asset, definitions = risk_model.populate_measure_definitions(_assets) + impacts, measures = risk_model.calculate_risk_measures( + _assets, scenarios, years + ) + measure_ids_for_asset, definitions = risk_model.populate_measure_definitions( + _assets + ) # create object for API: - risk_measures = _create_risk_measures(measures, measure_ids_for_asset, definitions, _assets, scenarios, years) + risk_measures = _create_risk_measures( + measures, measure_ids_for_asset, definitions, _assets, scenarios, years + ) elif request.include_asset_level: impacts = risk_model.calculate_impacts(_assets, scenarios, years) if request.include_asset_level: - asset_impacts = compile_asset_impacts(impacts, _assets, request.include_calc_details) + asset_impacts = compile_asset_impacts( + impacts, _assets, request.include_calc_details + ) else: asset_impacts = None @@ -361,7 +452,9 @@ def _get_asset_impacts( def compile_asset_impacts( - impacts: Dict[ImpactKey, List[AssetImpactResult]], assets: List[Asset], include_calc_details: bool + impacts: Dict[ImpactKey, List[AssetImpactResult]], + assets: List[Asset], + include_calc_details: bool, ): """Convert (internal) list of AssetImpactResult objects to a list of AssetLevelImpact objects ready for serialization. @@ -390,10 +483,12 @@ def compile_asset_impacts( ) calc_details = AcuteHazardCalculationDetails( hazard_exceedance=ExceedanceCurve( - values=hazard_exceedance.values, exceed_probabilities=hazard_exceedance.probs + values=hazard_exceedance.values, + exceed_probabilities=hazard_exceedance.probs, ), hazard_distribution=Distribution( - bin_edges=v.event.intensity_bin_edges, probabilities=v.event.prob + bin_edges=v.event.intensity_bin_edges, + probabilities=v.event.prob, ), vulnerability_distribution=vulnerability_distribution, hazard_path=v.impact.path, @@ -405,21 +500,31 @@ def compile_asset_impacts( continue impact_exceedance = v.impact.to_exceedance_curve() - key = APIImpactKey(hazard_type=k.hazard_type.__name__, scenario_id=k.scenario, year=str(k.key_year)) + key = APIImpactKey( + hazard_type=k.hazard_type.__name__, + scenario_id=k.scenario, + year=str(k.key_year), + ) hazard_impacts = AssetSingleImpact( key=key, impact_type=v.impact.impact_type.name, impact_exceedance=ExceedanceCurve( - values=impact_exceedance.values, exceed_probabilities=impact_exceedance.probs + values=impact_exceedance.values, + exceed_probabilities=impact_exceedance.probs, + ), + impact_distribution=Distribution( + bin_edges=v.impact.impact_bins, probabilities=v.impact.prob ), - impact_distribution=Distribution(bin_edges=v.impact.impact_bins, probabilities=v.impact.prob), impact_mean=v.impact.mean_impact(), impact_std_deviation=v.impact.stddev_impact(), calc_details=None if v.event is None else calc_details, ) ordered_impacts[k.asset].append(hazard_impacts) # note that this does rely on ordering of dictionary (post 3.6) - return [AssetLevelImpact(asset_id=k.id if k.id is not None else "", impacts=v) for k, v in ordered_impacts.items()] + return [ + AssetLevelImpact(asset_id=k.id if k.id is not None else "", impacts=v) + for k, v in ordered_impacts.items() + ] def _create_risk_measures( @@ -454,24 +559,39 @@ def _create_risk_measures( for year in years: # we calculate and tag results for each scenario, year and hazard score_key = RiskMeasureKey( - hazard_type=hazard_type.__name__, scenario_id=scenario_id, year=str(year), measure_id=measure_set_id + hazard_type=hazard_type.__name__, + scenario_id=scenario_id, + year=str(year), + measure_id=measure_set_id, ) scores = [-1] * len(assets) # measures_0 = [float("nan")] * len(assets) measures_0 = [nan_value] * len(assets) for i, asset in enumerate(assets): # look up result using the MeasureKey: - measure_key = MeasureKey(asset=asset, prosp_scen=scenario_id, year=year, hazard_type=hazard_type) + measure_key = MeasureKey( + asset=asset, + prosp_scen=scenario_id, + year=year, + hazard_type=hazard_type, + ) measure = measures.get(measure_key, None) if measure is not None: scores[i] = measure.score measures_0[i] = measure.measure_0 measures_for_assets.append( - RiskMeasuresForAssets(key=score_key, scores=scores, measures_0=measures_0, measures_1=None) + RiskMeasuresForAssets( + key=score_key, + scores=scores, + measures_0=measures_0, + measures_1=None, + ) ) score_based_measure_set_defn = ScoreBasedRiskMeasureSetDefinition( measure_set_id=measure_set_id, - asset_measure_ids_for_hazard={k.__name__: v for k, v in measure_ids_for_asset.items()}, + asset_measure_ids_for_hazard={ + k.__name__: v for k, v in measure_ids_for_asset.items() + }, score_definitions={v: k for (k, v) in definitions.items()}, ) return RiskMeasures( @@ -488,7 +608,9 @@ def _get_example_portfolios() -> List[Assets]: for file in importlib.resources.contents(physrisk.data.static.example_portfolios): if not str(file).endswith(".json"): continue - with importlib.resources.open_text(physrisk.data.static.example_portfolios, file) as f: + with importlib.resources.open_text( + physrisk.data.static.example_portfolios, file + ) as f: portfolio = Assets(**json.load(f)) portfolios.append(portfolio) return portfolios diff --git a/src/physrisk/risk_models/generic_risk_model.py b/src/physrisk/risk_models/generic_risk_model.py index 33f2baef..c4c78c79 100644 --- a/src/physrisk/risk_models/generic_risk_model.py +++ b/src/physrisk/risk_models/generic_risk_model.py @@ -12,8 +12,19 @@ RiskScoreValue, ScoreBasedRiskMeasureDefinition, ) -from physrisk.kernel.hazard_model import HazardEventDataResponse, HazardParameterDataResponse -from physrisk.kernel.hazards import ChronicHeat, Drought, Fire, Hail, Hazard, Precipitation, Wind +from physrisk.kernel.hazard_model import ( + HazardEventDataResponse, + HazardParameterDataResponse, +) +from physrisk.kernel.hazards import ( + ChronicHeat, + Drought, + Fire, + Hail, + Hazard, + Precipitation, + Wind, +) from physrisk.kernel.impact import AssetImpactResult from physrisk.kernel.risk import Measure, RiskMeasureCalculator @@ -109,7 +120,9 @@ def __init__(self): } self._definition_lookup[Wind] = ScoreBasedRiskMeasureDefinition( hazard_types=[Wind.__name__], - values=self._definition_values(self.wind_bounds, self.wind_label_description), + values=self._definition_values( + self.wind_bounds, self.wind_label_description + ), underlying_measures=[ RiskMeasureDefinition( measure_id="measure_wind", @@ -124,7 +137,9 @@ def __init__(self): ) self._definition_lookup[Hail] = ScoreBasedRiskMeasureDefinition( hazard_types=[Hail.__name__], - values=self._definition_values(self.hail_bounds, self.hail_label_description), + values=self._definition_values( + self.hail_bounds, self.hail_label_description + ), underlying_measures=[ RiskMeasureDefinition( measure_id="measure_hail", @@ -139,7 +154,9 @@ def __init__(self): ) self._definition_lookup[Drought] = ScoreBasedRiskMeasureDefinition( hazard_types=[Drought.__name__], - values=self._definition_values(self.drought_bounds, self.drought_label_description), + values=self._definition_values( + self.drought_bounds, self.drought_label_description + ), underlying_measures=[ RiskMeasureDefinition( measure_id="measure_drought", @@ -163,7 +180,9 @@ def __init__(self): ) self._definition_lookup[Fire] = ScoreBasedRiskMeasureDefinition( hazard_types=[Fire.__name__], - values=self._definition_values(self.fire_bounds, self.fire_label_description), + values=self._definition_values( + self.fire_bounds, self.fire_label_description + ), underlying_measures=[ RiskMeasureDefinition( measure_id="measure_fire", @@ -184,7 +203,9 @@ def __init__(self): ) self._definition_lookup[Precipitation] = ScoreBasedRiskMeasureDefinition( hazard_types=[Precipitation.__name__], - values=self._definition_values(self.precipitation_bounds, self.precipitation_label_description), + values=self._definition_values( + self.precipitation_bounds, self.precipitation_label_description + ), underlying_measures=[ RiskMeasureDefinition( measure_id="measure_precipitation", @@ -199,7 +220,9 @@ def __init__(self): ) self._definition_lookup[ChronicHeat] = ScoreBasedRiskMeasureDefinition( hazard_types=[ChronicHeat.__name__], - values=self._definition_values(self.heat_bounds, self.heat_label_description), + values=self._definition_values( + self.heat_bounds, self.heat_label_description + ), underlying_measures=[ RiskMeasureDefinition( measure_id="measure_chronicHeat", @@ -235,7 +258,9 @@ def _definition_values( return risk_score_values def wind_label_description(self, bounds: HazardIndicatorBounds): - label = f"Max wind speed between {bounds.lower} and {bounds.upper} {bounds.units}" + label = ( + f"Max wind speed between {bounds.lower} and {bounds.upper} {bounds.units}" + ) description = f"Max sustained wind speed between {bounds.lower} and {bounds.upper} {bounds.units}" return label, description @@ -277,7 +302,10 @@ def heat_label_description(self, bounds: HazardIndicatorBounds): return label, description def calc_measure( - self, hazard_type: Type[Hazard], base_impact_res: AssetImpactResult, impact_res: AssetImpactResult + self, + hazard_type: Type[Hazard], + base_impact_res: AssetImpactResult, + impact_res: AssetImpactResult, ) -> Measure: # in general we want to use the impact distribution, but in certain circumstances we can use # the underlying hazard data some care is needed given that vulnerability models are interchangeable @@ -295,17 +323,23 @@ def calc_measure( param = resp.parameter elif isinstance(resp, HazardEventDataResponse): return_period = bounds[0].indicator_return - param = float(np.interp(return_period, resp.return_periods, resp.intensities)) + param = float( + np.interp(return_period, resp.return_periods, resp.intensities) + ) if resp.units != "default": param = ureg.convert(param, resp.units, bounds[0].units) if math.isnan(param): return Measure( - score=Category.NODATA, measure_0=float(param), definition=self.get_definition(hazard_type) + score=Category.NODATA, + measure_0=float(param), + definition=self.get_definition(hazard_type), ) else: index = np.searchsorted(lower_bounds, param, side="right") - 1 return Measure( - score=categories[index], measure_0=float(param), definition=self.get_definition(hazard_type) + score=categories[index], + measure_0=float(param), + definition=self.get_definition(hazard_type), ) else: raise NotImplementedError("impact distribution case not implemented yet") @@ -314,4 +348,6 @@ def get_definition(self, hazard_type: Type[Hazard]): return self._definition_lookup.get(hazard_type, None) def supported_hazards(self) -> Set[type]: - return set([Wind, Fire, Hail, ChronicHeat, Drought, Precipitation]) # RiverineInundation, CoastalInundation, + return set( + [Wind, Fire, Hail, ChronicHeat, Drought, Precipitation] + ) # RiverineInundation, CoastalInundation, diff --git a/src/physrisk/risk_models/loss_model.py b/src/physrisk/risk_models/loss_model.py index 335f7b39..337d596b 100644 --- a/src/physrisk/risk_models/loss_model.py +++ b/src/physrisk/risk_models/loss_model.py @@ -6,11 +6,17 @@ from physrisk.kernel.impact_distrib import ImpactDistrib, ImpactType from ..kernel.assets import Asset -from ..kernel.calculation import get_default_hazard_model, get_default_vulnerability_models +from ..kernel.calculation import ( + get_default_hazard_model, + get_default_vulnerability_models, +) from ..kernel.financial_model import FinancialModelBase from ..kernel.hazard_model import HazardModel from ..kernel.impact import calculate_impacts -from ..kernel.vulnerability_model import DictBasedVulnerabilityModels, VulnerabilityModels +from ..kernel.vulnerability_model import ( + DictBasedVulnerabilityModels, + VulnerabilityModels, +) class Aggregator(ABC): @@ -29,7 +35,9 @@ def __init__( hazard_model: Optional[HazardModel] = None, vulnerability_models: Optional[VulnerabilityModels] = None, ): - self.hazard_model = get_default_hazard_model() if hazard_model is None else hazard_model + self.hazard_model = ( + get_default_hazard_model() if hazard_model is None else hazard_model + ) self.vulnerability_models = ( DictBasedVulnerabilityModels(get_default_vulnerability_models()) if vulnerability_models is None @@ -54,7 +62,13 @@ def get_financial_impacts( aggregation_pools: Dict[str, np.ndarray] = {} - results = calculate_impacts(assets, self.hazard_model, self.vulnerability_models, scenario=scenario, year=year) + results = calculate_impacts( + assets, + self.hazard_model, + self.vulnerability_models, + scenario=scenario, + year=year, + ) # the impacts in the results are either fractional damage or a fractional disruption rg = np.random.Generator(np.random.MT19937(seed=111)) @@ -71,9 +85,13 @@ def get_financial_impacts( impact_samples = self.uncorrelated_samples(impact, sims, rg) if impact.impact_type == ImpactType.damage: - loss = financial_model.damage_to_loss(impact_key.asset, impact_samples, currency) + loss = financial_model.damage_to_loss( + impact_key.asset, impact_samples, currency + ) else: # impact.impact_type == ImpactType.disruption: - loss = financial_model.disruption_to_loss(impact_key.asset, impact_samples, year, currency) + loss = financial_model.disruption_to_loss( + impact_key.asset, impact_samples, year, currency + ) for key in keys: if key not in aggregation_pools: @@ -91,5 +109,7 @@ def get_financial_impacts( return measures - def uncorrelated_samples(self, impact: ImpactDistrib, samples: int, generator: np.random.Generator) -> np.ndarray: + def uncorrelated_samples( + self, impact: ImpactDistrib, samples: int, generator: np.random.Generator + ) -> np.ndarray: return impact.to_exceedance_curve().get_samples(generator.uniform(size=samples)) diff --git a/src/physrisk/risk_models/risk_models.py b/src/physrisk/risk_models/risk_models.py index 01723ad4..edb1531d 100644 --- a/src/physrisk/risk_models/risk_models.py +++ b/src/physrisk/risk_models/risk_models.py @@ -7,7 +7,13 @@ RiskScoreValue, ScoreBasedRiskMeasureDefinition, ) -from physrisk.kernel.hazards import ChronicHeat, CoastalInundation, Hazard, RiverineInundation, Wind +from physrisk.kernel.hazards import ( + ChronicHeat, + CoastalInundation, + Hazard, + RiverineInundation, + Wind, +) from physrisk.kernel.impact import AssetImpactResult from physrisk.kernel.impact_distrib import EmptyImpactDistrib, ImpactDistrib from physrisk.kernel.risk import Measure, RiskMeasureCalculator @@ -24,7 +30,9 @@ class RealEstateToyRiskMeasures(RiskMeasureCalculator): def __init__(self): self.model_summary = {"*Toy* model for real estate risk assessment."} - self.return_period = 100.0 # criteria based on 1 in 100-year flood or cyclone events + self.return_period = ( + 100.0 # criteria based on 1 in 100-year flood or cyclone events + ) self.measure_thresholds_acute = { Threshold.ABS_HIGH: 0.1, # fraction Threshold.ABS_LOW: 0.03, # fraction @@ -47,7 +55,11 @@ def __init__(self): def _definition_acute(self): definition = ScoreBasedRiskMeasureDefinition( - hazard_types=[RiverineInundation.__name__, CoastalInundation.__name__, Wind.__name__], + hazard_types=[ + RiverineInundation.__name__, + CoastalInundation.__name__, + Wind.__name__, + ], values=self._definition_values(self._acute_description), underlying_measures=[ RiskMeasureDefinition( @@ -101,7 +113,9 @@ def _definition_values(self, description: Callable[[Category], str]): label="No material impact.", description=description(Category.LOW), ), - RiskScoreValue(value=Category.NODATA, label="No data.", description="No data."), + RiskScoreValue( + value=Category.NODATA, label="No data.", description="No data." + ), ] def _acute_description(self, category: Category): @@ -161,16 +175,27 @@ def _cooling_description(self, category: Category): return description def calc_measure( - self, hazard_type: Type[Hazard], base_impact_res: AssetImpactResult, impact_res: AssetImpactResult + self, + hazard_type: Type[Hazard], + base_impact_res: AssetImpactResult, + impact_res: AssetImpactResult, ) -> Optional[Measure]: - if isinstance(base_impact_res.impact, EmptyImpactDistrib) or isinstance(impact_res.impact, EmptyImpactDistrib): + if isinstance(base_impact_res.impact, EmptyImpactDistrib) or isinstance( + impact_res.impact, EmptyImpactDistrib + ): return None if hazard_type == ChronicHeat: - return self.calc_measure_cooling(hazard_type, base_impact_res.impact, impact_res.impact) + return self.calc_measure_cooling( + hazard_type, base_impact_res.impact, impact_res.impact + ) else: - return self.calc_measure_acute(hazard_type, base_impact_res.impact, impact_res.impact) + return self.calc_measure_acute( + hazard_type, base_impact_res.impact, impact_res.impact + ) - def calc_measure_acute(self, hazard_type: type, base_impact: ImpactDistrib, impact: ImpactDistrib) -> Measure: + def calc_measure_acute( + self, hazard_type: type, base_impact: ImpactDistrib, impact: ImpactDistrib + ) -> Measure: return_period = 100.0 # criterion based on 1 in 100-year flood events histo_loss = base_impact.to_exceedance_curve().get_value(1.0 / return_period) future_loss = impact.to_exceedance_curve().get_value(1.0 / return_period) @@ -193,9 +218,15 @@ def calc_measure_acute(self, hazard_type: type, base_impact: ImpactDistrib, impa score = Category.MEDIUM else: score = Category.LOW - return Measure(score=score, measure_0=future_loss, definition=self.get_definition(hazard_type)) + return Measure( + score=score, + measure_0=future_loss, + definition=self.get_definition(hazard_type), + ) - def calc_measure_cooling(self, hazard_type: type, base_impact: ImpactDistrib, impact: ImpactDistrib) -> Measure: + def calc_measure_cooling( + self, hazard_type: type, base_impact: ImpactDistrib, impact: ImpactDistrib + ) -> Measure: histo_cooling = base_impact.mean_impact() future_cooling = impact.mean_impact() cooling_change = (future_cooling - histo_cooling) / histo_cooling @@ -217,7 +248,11 @@ def calc_measure_cooling(self, hazard_type: type, base_impact: ImpactDistrib, im score = Category.MEDIUM else: score = Category.LOW - return Measure(score=score, measure_0=future_cooling, definition=self.get_definition(hazard_type)) + return Measure( + score=score, + measure_0=future_cooling, + definition=self.get_definition(hazard_type), + ) def get_definition(self, hazard_type: type): return self._definition_lookup.get(hazard_type, None) diff --git a/src/physrisk/vulnerability_models/chronic_heat_models.py b/src/physrisk/vulnerability_models/chronic_heat_models.py index eb5a65f9..e922b7cc 100644 --- a/src/physrisk/vulnerability_models/chronic_heat_models.py +++ b/src/physrisk/vulnerability_models/chronic_heat_models.py @@ -4,7 +4,11 @@ from scipy.stats import norm from physrisk.kernel.assets import Asset, IndustrialActivity -from physrisk.kernel.hazard_model import HazardDataRequest, HazardDataResponse, HazardParameterDataResponse +from physrisk.kernel.hazard_model import ( + HazardDataRequest, + HazardDataResponse, + HazardParameterDataResponse, +) from physrisk.kernel.hazards import ChronicHeat from physrisk.kernel.impact_distrib import ImpactDistrib, ImpactType from physrisk.kernel.vulnerability_model import VulnerabilityModelBase @@ -18,10 +22,16 @@ class ChronicHeatGZNModel(VulnerabilityModelBase): def __init__(self, indicator_id: str = "mean_degree_days/above/32c", delta=True): super().__init__( - indicator_id=indicator_id, hazard_type=ChronicHeat, impact_type=ImpactType.disruption + indicator_id=indicator_id, + hazard_type=ChronicHeat, + impact_type=ImpactType.disruption, ) # opportunity to give a model hint, but blank here - self.time_lost_per_degree_day = 4.671 # This comes from the paper converted to celsius - self.time_lost_per_degree_day_se = 2.2302 # This comes from the paper converted to celsius + self.time_lost_per_degree_day = ( + 4.671 # This comes from the paper converted to celsius + ) + self.time_lost_per_degree_day_se = ( + 2.2302 # This comes from the paper converted to celsius + ) self.total_labour_hours = 107460 self.delta = delta @@ -59,7 +69,9 @@ def get_data_requests( ), ] - def get_impact(self, asset: Asset, data_responses: List[HazardDataResponse]) -> ImpactDistrib: + def get_impact( + self, asset: Asset, data_responses: List[HazardDataResponse] + ) -> ImpactDistrib: """Calcaulate impact (disruption) of asset based on the hazard data returned. Args: @@ -75,14 +87,25 @@ def get_impact(self, asset: Asset, data_responses: List[HazardDataResponse]) -> assert isinstance(baseline_dd_above_mean, HazardParameterDataResponse) assert isinstance(scenario_dd_above_mean, HazardParameterDataResponse) - delta_dd_above_mean: float = scenario_dd_above_mean.parameter - baseline_dd_above_mean.parameter * self.delta + delta_dd_above_mean: float = ( + scenario_dd_above_mean.parameter + - baseline_dd_above_mean.parameter * self.delta + ) hours_worked = self.total_labour_hours - fraction_loss_mean = (delta_dd_above_mean * self.time_lost_per_degree_day) / hours_worked - fraction_loss_std = (delta_dd_above_mean * self.time_lost_per_degree_day_se) / hours_worked + fraction_loss_mean = ( + delta_dd_above_mean * self.time_lost_per_degree_day + ) / hours_worked + fraction_loss_std = ( + delta_dd_above_mean * self.time_lost_per_degree_day_se + ) / hours_worked return get_impact_distrib( - fraction_loss_mean, fraction_loss_std, ChronicHeat, [scenario_dd_above_mean.path], ImpactType.disruption + fraction_loss_mean, + fraction_loss_std, + ChronicHeat, + [scenario_dd_above_mean.path], + ImpactType.disruption, ) @@ -92,10 +115,16 @@ class ChronicHeatWBGTGZNModel(ChronicHeatGZNModel): results based on applying both GZN and WBGT.""" def __init__(self, indicator_id: str = "mean_work_loss_high"): - super().__init__(indicator_id=indicator_id) # opportunity to give a model hint, but blank here + super().__init__( + indicator_id=indicator_id + ) # opportunity to give a model hint, but blank here def work_type_mapping(self): - return {"low": ["low", "medium"], "medium": ["medium", "low", "high"], "high": ["high", "medium"]} + return { + "low": ["low", "medium"], + "medium": ["medium", "low", "high"], + "high": ["high", "medium"], + } def get_data_requests( self, asset: Asset, *, scenario: str, year: int @@ -162,57 +191,83 @@ def get_data_requests( ), ] + wbgt_data_requests - def get_impact(self, asset: Asset, data_responses: List[HazardDataResponse]) -> ImpactDistrib: + def get_impact( + self, asset: Asset, data_responses: List[HazardDataResponse] + ) -> ImpactDistrib: """ Function to return the impact distribution of the wbgt model. """ assert isinstance(asset, IndustrialActivity) - wbgt_responses = [cast(HazardParameterDataResponse, r) for r in data_responses[2:]] + wbgt_responses = [ + cast(HazardParameterDataResponse, r) for r in data_responses[2:] + ] - hazard_paths = [cast(HazardParameterDataResponse, r).path for r in data_responses] + hazard_paths = [ + cast(HazardParameterDataResponse, r).path for r in data_responses + ] baseline_dd_above_mean = cast(HazardParameterDataResponse, data_responses[0]) scenario_dd_above_mean = cast(HazardParameterDataResponse, data_responses[1]) hours_worked = self.total_labour_hours - fraction_loss_mean_base_gzn = (baseline_dd_above_mean.parameter * self.time_lost_per_degree_day) / hours_worked + fraction_loss_mean_base_gzn = ( + baseline_dd_above_mean.parameter * self.time_lost_per_degree_day + ) / hours_worked fraction_loss_mean_scenario_gzn = ( scenario_dd_above_mean.parameter * self.time_lost_per_degree_day ) / hours_worked - fraction_loss_std_base = (baseline_dd_above_mean.parameter * self.time_lost_per_degree_day_se) / hours_worked + fraction_loss_std_base = ( + baseline_dd_above_mean.parameter * self.time_lost_per_degree_day_se + ) / hours_worked fraction_loss_std_scenario = ( scenario_dd_above_mean.parameter * self.time_lost_per_degree_day_se ) / hours_worked - baseline_work_ability = (1 - fraction_loss_mean_base_gzn) * (1 - wbgt_responses[0].parameter) - scenario_work_ability = (1 - fraction_loss_mean_scenario_gzn) * (1 - wbgt_responses[1].parameter) + baseline_work_ability = (1 - fraction_loss_mean_base_gzn) * ( + 1 - wbgt_responses[0].parameter + ) + scenario_work_ability = (1 - fraction_loss_mean_scenario_gzn) * ( + 1 - wbgt_responses[1].parameter + ) # Getting the parameters required for the uniform distribution. if asset.type in ["low", "high"]: a_historical = ( - wbgt_responses[0].parameter - abs((wbgt_responses[2].parameter - wbgt_responses[0].parameter)) / 2 + wbgt_responses[0].parameter + - abs((wbgt_responses[2].parameter - wbgt_responses[0].parameter)) / 2 ) b_historical = ( - wbgt_responses[0].parameter + abs((wbgt_responses[2].parameter - wbgt_responses[0].parameter)) / 2 + wbgt_responses[0].parameter + + abs((wbgt_responses[2].parameter - wbgt_responses[0].parameter)) / 2 ) a_scenario = ( - wbgt_responses[1].parameter - abs((wbgt_responses[3].parameter - wbgt_responses[1].parameter)) / 2 + wbgt_responses[1].parameter + - abs((wbgt_responses[3].parameter - wbgt_responses[1].parameter)) / 2 ) b_scenario = ( - wbgt_responses[1].parameter + abs((wbgt_responses[3].parameter - wbgt_responses[1].parameter)) / 2 + wbgt_responses[1].parameter + + abs((wbgt_responses[3].parameter - wbgt_responses[1].parameter)) / 2 ) elif asset.type == "medium": - a_historical = wbgt_responses[0].parameter - (wbgt_responses[2].parameter - wbgt_responses[0].parameter) / 2 - b_historical = wbgt_responses[0].parameter + (wbgt_responses[4].parameter - wbgt_responses[0].parameter) / 2 + a_historical = ( + wbgt_responses[0].parameter + - (wbgt_responses[2].parameter - wbgt_responses[0].parameter) / 2 + ) + b_historical = ( + wbgt_responses[0].parameter + + (wbgt_responses[4].parameter - wbgt_responses[0].parameter) / 2 + ) a_scenario = ( - wbgt_responses[1].parameter - abs((wbgt_responses[3].parameter - wbgt_responses[1].parameter)) / 2 + wbgt_responses[1].parameter + - abs((wbgt_responses[3].parameter - wbgt_responses[1].parameter)) / 2 ) b_scenario = ( - wbgt_responses[1].parameter + abs((wbgt_responses[5].parameter - wbgt_responses[1].parameter)) / 2 + wbgt_responses[1].parameter + + abs((wbgt_responses[5].parameter - wbgt_responses[1].parameter)) / 2 ) # Estimation of the variance @@ -236,7 +291,13 @@ def get_impact(self, asset: Asset, data_responses: List[HazardDataResponse]) -> total_work_loss_delta: float = baseline_work_ability - scenario_work_ability - return get_impact_distrib(total_work_loss_delta, std_delta, self.hazard_type, hazard_paths, self.impact_type) + return get_impact_distrib( + total_work_loss_delta, + std_delta, + self.hazard_type, + hazard_paths, + self.impact_type, + ) def two_variable_joint_variance(ex, varx, ey, vary): @@ -275,9 +336,11 @@ def get_impact_distrib( ] ) - probs_cumulative = np.vectorize(lambda x: norm.cdf(x, loc=fraction_loss_mean, scale=max(1e-12, fraction_loss_std)))( - impact_bins - ) + probs_cumulative = np.vectorize( + lambda x: norm.cdf( + x, loc=fraction_loss_mean, scale=max(1e-12, fraction_loss_std) + ) + )(impact_bins) probs_cumulative[-1] = np.maximum(probs_cumulative[-1], 1.0) probs = np.diff(probs_cumulative) diff --git a/src/physrisk/vulnerability_models/example_models.py b/src/physrisk/vulnerability_models/example_models.py index 4817ed75..79c5f45b 100644 --- a/src/physrisk/vulnerability_models/example_models.py +++ b/src/physrisk/vulnerability_models/example_models.py @@ -8,15 +8,25 @@ from ..kernel.impact_distrib import EmptyImpactDistrib, ImpactDistrib, ImpactType from ..kernel.vulnerability_matrix_provider import VulnMatrixProvider -from ..kernel.vulnerability_model import VulnerabilityModel, VulnerabilityModelBase, checked_beta_distrib +from ..kernel.vulnerability_model import ( + VulnerabilityModel, + VulnerabilityModelBase, + checked_beta_distrib, +) class ExampleCdfBasedVulnerabilityModel(VulnerabilityModel): def __init__(self, *, indicator_id: str, hazard_type: type): self.intensities = np.array([0, 0.01, 0.5, 1.0, 1.5, 2, 3, 4, 5, 6]) - self.impact_means = np.array([0, 0.2, 0.44, 0.58, 0.68, 0.78, 0.85, 0.92, 0.96, 1.0]) - self.impact_stddevs = np.array([0, 0.17, 0.14, 0.14, 0.17, 0.14, 0.13, 0.10, 0.06, 0]) - impact_bin_edges = np.array([0, 0.01, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]) + self.impact_means = np.array( + [0, 0.2, 0.44, 0.58, 0.68, 0.78, 0.85, 0.92, 0.96, 1.0] + ) + self.impact_stddevs = np.array( + [0, 0.17, 0.14, 0.14, 0.17, 0.14, 0.13, 0.10, 0.06, 0] + ) + impact_bin_edges = np.array( + [0, 0.01, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0] + ) super().__init__( indicator_id=indicator_id, hazard_type=hazard_type, @@ -29,7 +39,10 @@ def get_impact_curve(self, intensities, asset): impact_means = np.interp(intensities, self.intensities, self.impact_means) impact_stddevs = np.interp(intensities, self.intensities, self.impact_stddevs) return VulnMatrixProvider( - intensities, impact_cdfs=[checked_beta_distrib(m, s) for m, s in zip(impact_means, impact_stddevs)] + intensities, + impact_cdfs=[ + checked_beta_distrib(m, s) for m, s in zip(impact_means, impact_stddevs) + ], ) @@ -39,7 +52,9 @@ class PlaceholderVulnerabilityModel(VulnerabilityModelBase): """ def __init__(self, indicator_id: str, hazard_type: type, impact_type: ImpactType): - super().__init__(indicator_id=indicator_id, hazard_type=hazard_type, impact_type=impact_type) + super().__init__( + indicator_id=indicator_id, hazard_type=hazard_type, impact_type=impact_type + ) def get_data_requests(self, asset: Asset, *, scenario: str, year: int): return HazardDataRequest( @@ -51,7 +66,9 @@ def get_data_requests(self, asset: Asset, *, scenario: str, year: int): indicator_id=self.indicator_id, ) - def get_impact(self, asset: Asset, data_responses: Iterable[HazardDataResponse]) -> ImpactDistrib: + def get_impact( + self, asset: Asset, data_responses: Iterable[HazardDataResponse] + ) -> ImpactDistrib: return EmptyImpactDistrib() @@ -80,7 +97,9 @@ def __init__(self): impact_type=ImpactType.damage, ) - def get_impact(self, asset: Asset, data_responses: Iterable[HazardDataResponse]) -> ImpactDistrib: + def get_impact( + self, asset: Asset, data_responses: Iterable[HazardDataResponse] + ) -> ImpactDistrib: # params = data_responses # assert isinstance(data_response, HazardParameterDataResponse) # TODO add model with stochastic impact diff --git a/src/physrisk/vulnerability_models/power_generating_asset_models.py b/src/physrisk/vulnerability_models/power_generating_asset_models.py index 30931933..f3fdedc8 100644 --- a/src/physrisk/vulnerability_models/power_generating_asset_models.py +++ b/src/physrisk/vulnerability_models/power_generating_asset_models.py @@ -5,11 +5,19 @@ from ..kernel.assets import Asset, PowerGeneratingAsset from ..kernel.curve import ExceedanceCurve from ..kernel.hazard_event_distrib import HazardEventDistrib -from ..kernel.hazard_model import HazardDataRequest, HazardDataResponse, HazardEventDataResponse +from ..kernel.hazard_model import ( + HazardDataRequest, + HazardDataResponse, + HazardEventDataResponse, +) from ..kernel.hazards import Drought, RiverineInundation from ..kernel.impact_distrib import ImpactType from ..kernel.vulnerability_distrib import VulnerabilityDistrib -from ..kernel.vulnerability_model import DeterministicVulnerabilityModel, VulnerabilityModelAcuteBase, applies_to_assets +from ..kernel.vulnerability_model import ( + DeterministicVulnerabilityModel, + VulnerabilityModelAcuteBase, + applies_to_assets, +) @applies_to_assets([PowerGeneratingAsset]) @@ -19,7 +27,11 @@ def __init__(self, indicator_id="flood_depth"): self.__curve_depth = np.array([0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 1]) self.__curve_impact = np.array([0, 1, 2, 7, 14, 30, 60, 180, 365]) self.__indicator_id = indicator_id - super().__init__(indicator_id=indicator_id, hazard_type=RiverineInundation, impact_type=ImpactType.disruption) + super().__init__( + indicator_id=indicator_id, + hazard_type=RiverineInundation, + impact_type=ImpactType.disruption, + ) pass def get_data_requests( @@ -48,7 +60,9 @@ def get_data_requests( return histo, future - def get_distributions(self, asset: Asset, event_data_responses: Iterable[HazardDataResponse]): + def get_distributions( + self, asset: Asset, event_data_responses: Iterable[HazardDataResponse] + ): """Return distributions for asset, based on hazard event date: VulnerabilityDistrib and HazardEventDistrib.""" @@ -66,14 +80,18 @@ def get_distributions(self, asset: Asset, event_data_responses: Iterable[HazardD depth_bins, probs = curve_future.get_probability_bins() - impact_bins = np.interp(depth_bins, self.__curve_depth, self.__curve_impact) / 365.0 + impact_bins = ( + np.interp(depth_bins, self.__curve_depth, self.__curve_impact) / 365.0 + ) # keep all bins, but make use of vulnerability matrix to apply protection level # for improved performance we could truncate (and treat identify matrix as a special case) # but this general version allows model uncertainties to be added probs_protected = np.where(depth_bins[1:] <= protection_depth, 0.0, 1.0) - vul = VulnerabilityDistrib(RiverineInundation, depth_bins, impact_bins, np.diag(probs_protected)) + vul = VulnerabilityDistrib( + RiverineInundation, depth_bins, impact_bins, np.diag(probs_protected) + ) event = HazardEventDistrib(RiverineInundation, depth_bins, probs, [future.path]) return vul, event diff --git a/src/physrisk/vulnerability_models/real_estate_models.py b/src/physrisk/vulnerability_models/real_estate_models.py index efa13528..e41eeafc 100644 --- a/src/physrisk/vulnerability_models/real_estate_models.py +++ b/src/physrisk/vulnerability_models/real_estate_models.py @@ -5,12 +5,22 @@ from physrisk.api.v1.common import VulnerabilityCurve, VulnerabilityCurves from physrisk.kernel.assets import Asset, RealEstateAsset -from physrisk.kernel.hazard_model import HazardDataRequest, HazardDataResponse, HazardParameterDataResponse +from physrisk.kernel.hazard_model import ( + HazardDataRequest, + HazardDataResponse, + HazardParameterDataResponse, +) from physrisk.kernel.impact_distrib import ImpactDistrib, ImpactType from physrisk.kernel.vulnerability_matrix_provider import VulnMatrixProvider from physrisk.kernel.vulnerability_model import VulnerabilityModel -from ..kernel.hazards import ChronicHeat, CoastalInundation, PluvialInundation, RiverineInundation, Wind +from ..kernel.hazards import ( + ChronicHeat, + CoastalInundation, + PluvialInundation, + RiverineInundation, + Wind, +) from ..kernel.vulnerability_model import ( DeterministicVulnerabilityModel, VulnerabilityModelBase, @@ -21,7 +31,9 @@ class RealEstateInundationModel(VulnerabilityModel): - _default_impact_bin_edges = np.array([0, 0.01, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]) + _default_impact_bin_edges = np.array( + [0, 0.01, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0] + ) _default_resource = "EU JRC global flood depth-damage functions" def __init__( @@ -42,10 +54,14 @@ def __init__( impact_bin_edges: specifies the impact (fractional damage/disruption bins). """ - curve_set: VulnerabilityCurves = get_vulnerability_curves_from_resource(resource) + curve_set: VulnerabilityCurves = get_vulnerability_curves_from_resource( + resource + ) # for this model, key for looking up curves is (location, asset_type), e.g. ('Asia', 'Building/Industrial') - self.vulnerability_curves = dict(((c.location, c.asset_type), c) for c in curve_set.items) + self.vulnerability_curves = dict( + ((c.location, c.asset_type), c) for c in curve_set.items + ) self.vuln_curves_by_type = defaultdict(list) self.proxy_curves: Dict[Tuple[str, str], VulnerabilityCurve] = {} for item in curve_set.items: @@ -55,9 +71,10 @@ def __init__( impact_type = ( ImpactType.damage if len(self.vulnerability_curves) == 0 - else [ImpactType[self.vulnerability_curves[key].impact_type.lower()] for key in self.vulnerability_curves][ - 0 - ] + else [ + ImpactType[self.vulnerability_curves[key].impact_type.lower()] + for key in self.vulnerability_curves + ][0] ) super().__init__( indicator_id=indicator_id, @@ -79,23 +96,39 @@ def get_impact_curve(self, intensity_bin_centres: np.ndarray, asset: Asset): self.proxy_curves[key] = self.closest_curve_of_type(curve, asset) std_curve = self.proxy_curves[key] - impact_means = np.interp(intensity_bin_centres, curve.intensity, curve.impact_mean) - impact_stddevs = np.interp(intensity_bin_centres, std_curve.intensity, std_curve.impact_std) + impact_means = np.interp( + intensity_bin_centres, curve.intensity, curve.impact_mean + ) + impact_stddevs = np.interp( + intensity_bin_centres, std_curve.intensity, std_curve.impact_std + ) return VulnMatrixProvider( intensity_bin_centres, - impact_cdfs=[checked_beta_distrib(m, s) for m, s in zip(impact_means, impact_stddevs)], + impact_cdfs=[ + checked_beta_distrib(m, s) for m, s in zip(impact_means, impact_stddevs) + ], ) def closest_curve_of_type(self, curve: VulnerabilityCurve, asset: RealEstateAsset): # we return the standard deviations of the damage curve most similar to the asset location - candidate_set = list(cand for cand in self.vuln_curves_by_type[asset.type] if (len(cand.impact_std) > 0)) + candidate_set = list( + cand + for cand in self.vuln_curves_by_type[asset.type] + if (len(cand.impact_std) > 0) + ) sum_square_diff = (self.sum_square_diff(curve, cand) for cand in candidate_set) lowest = np.argmin(np.array(list(sum_square_diff))) return candidate_set[lowest] def sum_square_diff(self, curve1: VulnerabilityCurve, curve2: VulnerabilityCurve): - return np.sum((curve1.impact_mean - np.interp(curve1.intensity, curve2.intensity, curve2.impact_mean)) ** 2) + return np.sum( + ( + curve1.impact_mean + - np.interp(curve1.intensity, curve2.intensity, curve2.impact_mean) + ) + ** 2 + ) @applies_to_events([CoastalInundation]) @@ -224,12 +257,22 @@ def get_data_requests(self, asset: Asset, *, scenario: str, year: int): indicator_id=self.indicator_id, ) - def get_impact(self, asset: Asset, data_responses: List[HazardDataResponse]) -> ImpactDistrib: + def get_impact( + self, asset: Asset, data_responses: List[HazardDataResponse] + ) -> ImpactDistrib: (data,) = data_responses assert isinstance(data, HazardParameterDataResponse) # we interpolate the specific threshold from the different values - deg_days = float(np.interp(self.threshold_temp_c, data.param_defns, data.parameters)) # [0] - heat_transfer = deg_days * self._default_transfer_coeff * 24 / 1000 # kWh of heat removed from asset - annual_electricity = heat_transfer / self._default_cooling_cop # kWh of electricity required for heat removal + deg_days = float( + np.interp(self.threshold_temp_c, data.param_defns, data.parameters) + ) # [0] + heat_transfer = ( + deg_days * self._default_transfer_coeff * 24 / 1000 + ) # kWh of heat removed from asset + annual_electricity = ( + heat_transfer / self._default_cooling_cop + ) # kWh of electricity required for heat removal # this is non-probabilistic model: probability of 1 of electricity use - return ImpactDistrib(ChronicHeat, [annual_electricity, annual_electricity], [1], [data.path]) + return ImpactDistrib( + ChronicHeat, [annual_electricity, annual_electricity], [1], [data.path] + ) diff --git a/src/physrisk/vulnerability_models/thermal_power_generation_models.py b/src/physrisk/vulnerability_models/thermal_power_generation_models.py index 5df51c30..a424ef10 100644 --- a/src/physrisk/vulnerability_models/thermal_power_generation_models.py +++ b/src/physrisk/vulnerability_models/thermal_power_generation_models.py @@ -7,7 +7,10 @@ from physrisk.api.v1.common import VulnerabilityCurve, VulnerabilityCurves from physrisk.kernel.assets import Asset, ThermalPowerGeneratingAsset, TurbineKind from physrisk.kernel.impact_distrib import EmptyImpactDistrib, ImpactDistrib, ImpactType -from physrisk.kernel.vulnerability_model import DeterministicVulnerabilityModel, VulnerabilityModelBase +from physrisk.kernel.vulnerability_model import ( + DeterministicVulnerabilityModel, + VulnerabilityModelBase, +) from ..kernel.curve import ExceedanceCurve from ..kernel.hazard_event_distrib import HazardEventDistrib @@ -27,7 +30,11 @@ WaterTemperature, ) from ..kernel.vulnerability_distrib import VulnerabilityDistrib -from ..kernel.vulnerability_model import applies_to_assets, applies_to_events, get_vulnerability_curves_from_resource +from ..kernel.vulnerability_model import ( + applies_to_assets, + applies_to_events, + get_vulnerability_curves_from_resource, +) class ThermalPowerGenerationInundationModel(DeterministicVulnerabilityModel): @@ -38,7 +45,12 @@ class ThermalPowerGenerationInundationModel(DeterministicVulnerabilityModel): _default_buffer = 1000 def __init__( - self, *, hazard_type: type, indicator_id: str, resource: str = _default_resource, buffer: int = _default_buffer + self, + *, + hazard_type: type, + indicator_id: str, + resource: str = _default_resource, + buffer: int = _default_buffer, ): """ Inundation vulnerability model for thermal power generation. @@ -51,22 +63,29 @@ def __init__( buffer (int): delimitation of the area for the hazard data expressed in metres (within [0,1000]). """ - curve_set: VulnerabilityCurves = get_vulnerability_curves_from_resource(resource) + curve_set: VulnerabilityCurves = get_vulnerability_curves_from_resource( + resource + ) # for this model, key for looking up curves is asset_type, e.g. 'Steam/Recirculating' self.vulnerability_curves = dict( - (c.asset_type, c) for c in curve_set.items if c.event_type == hazard_type.__base__.__name__ # type:ignore + (c.asset_type, c) + for c in curve_set.items + if c.event_type == hazard_type.__base__.__name__ # type:ignore ) self.vuln_curves_by_type = defaultdict(list) for key in self.vulnerability_curves: - self.vuln_curves_by_type[TurbineKind[key.split("/")[0]]].append(self.vulnerability_curves[key]) + self.vuln_curves_by_type[TurbineKind[key.split("/")[0]]].append( + self.vulnerability_curves[key] + ) impact_type = ( ImpactType.disruption if len(self.vulnerability_curves) == 0 - else [ImpactType[self.vulnerability_curves[key].impact_type.lower()] for key in self.vulnerability_curves][ - 0 - ] + else [ + ImpactType[self.vulnerability_curves[key].impact_type.lower()] + for key in self.vulnerability_curves + ][0] ) # global circulation parameter 'model' is a hint; can be overriden by hazard model @@ -113,24 +132,34 @@ def get_distributions( assert isinstance(response_scenario, HazardEventDataResponse) assert isinstance(response_baseline, HazardEventDataResponse) - baseline_curve = ExceedanceCurve(1.0 / response_baseline.return_periods, response_baseline.intensities) + baseline_curve = ExceedanceCurve( + 1.0 / response_baseline.return_periods, response_baseline.intensities + ) protection_depth = ( 0.0 if len(response_baseline.intensities) == 0 - else baseline_curve.get_value(1.0 / asset.get_inundation_protection_return_period()) + else baseline_curve.get_value( + 1.0 / asset.get_inundation_protection_return_period() + ) ) - intensity_curve = ExceedanceCurve(1.0 / response_scenario.return_periods, response_scenario.intensities) + intensity_curve = ExceedanceCurve( + 1.0 / response_scenario.return_periods, response_scenario.intensities + ) if 0 < len(intensity_curve.values): if intensity_curve.values[0] < protection_depth: if protection_depth < intensity_curve.values[-1]: intensity_curve = intensity_curve.add_value_point(protection_depth) - intensity_bins, probability_bins = intensity_curve.get_probability_bins(include_last=True) + intensity_bins, probability_bins = intensity_curve.get_probability_bins( + include_last=True + ) curves: List[VulnerabilityCurve] = [] if asset.turbine is None: - curves = [self.vulnerability_curves[key] for key in self.vulnerability_curves] + curves = [ + self.vulnerability_curves[key] for key in self.vulnerability_curves + ] elif asset.cooling is not None: key = "/".join([asset.turbine.name, asset.cooling.name]) if key in self.vulnerability_curves: @@ -141,7 +170,13 @@ def get_distributions( if 0 < len(curves): impact_bins = [ ( - np.max([np.interp(intensity, curve.intensity, curve.impact_mean) for curve in curves]) / 365.0 + np.max( + [ + np.interp(intensity, curve.intensity, curve.impact_mean) + for curve in curves + ] + ) + / 365.0 if protection_depth < intensity else 0.0 ) @@ -150,14 +185,20 @@ def get_distributions( else: impact_bins = [0.0 for _ in intensity_bins] - vul = VulnerabilityDistrib(self.hazard_type, intensity_bins, impact_bins, np.eye(len(probability_bins))) - event = HazardEventDistrib(self.hazard_type, intensity_bins, probability_bins, [response_scenario.path]) + vul = VulnerabilityDistrib( + self.hazard_type, intensity_bins, impact_bins, np.eye(len(probability_bins)) + ) + event = HazardEventDistrib( + self.hazard_type, intensity_bins, probability_bins, [response_scenario.path] + ) return vul, event @applies_to_events([CoastalInundation]) @applies_to_assets([ThermalPowerGeneratingAsset]) -class ThermalPowerGenerationCoastalInundationModel(ThermalPowerGenerationInundationModel): +class ThermalPowerGenerationCoastalInundationModel( + ThermalPowerGenerationInundationModel +): def __init__( self, *, @@ -165,12 +206,16 @@ def __init__( resource: str = ThermalPowerGenerationInundationModel._default_resource, ): # by default include subsidence and 95% sea-level rise - super().__init__(hazard_type=CoastalInundation, indicator_id=indicator_id, resource=resource) + super().__init__( + hazard_type=CoastalInundation, indicator_id=indicator_id, resource=resource + ) @applies_to_events([RiverineInundation]) @applies_to_assets([ThermalPowerGeneratingAsset]) -class ThermalPowerGenerationRiverineInundationModel(ThermalPowerGenerationInundationModel): +class ThermalPowerGenerationRiverineInundationModel( + ThermalPowerGenerationInundationModel +): def __init__( self, *, @@ -178,7 +223,9 @@ def __init__( resource: str = ThermalPowerGenerationInundationModel._default_resource, ): # by default request HazardModel to use "MIROC-ESM-CHEM" GCM - super().__init__(hazard_type=RiverineInundation, indicator_id=indicator_id, resource=resource) + super().__init__( + hazard_type=RiverineInundation, indicator_id=indicator_id, resource=resource + ) @applies_to_events([Drought]) @@ -203,28 +250,37 @@ def __init__( """ hazard_type = Drought - curve_set: VulnerabilityCurves = get_vulnerability_curves_from_resource(resource) + curve_set: VulnerabilityCurves = get_vulnerability_curves_from_resource( + resource + ) # for this model, key for looking up curves is asset_type, e.g. 'Steam/Recirculating' self.vulnerability_curves = dict( - (c.asset_type, c) for c in curve_set.items if c.event_type == hazard_type.__name__ + (c.asset_type, c) + for c in curve_set.items + if c.event_type == hazard_type.__name__ ) self.vuln_curves_by_type = defaultdict(list) for key in self.vulnerability_curves: - self.vuln_curves_by_type[TurbineKind[key.split("/")[0]]].append(self.vulnerability_curves[key]) + self.vuln_curves_by_type[TurbineKind[key.split("/")[0]]].append( + self.vulnerability_curves[key] + ) impact_type = ( ImpactType.disruption if len(self.vulnerability_curves) == 0 - else [ImpactType[self.vulnerability_curves[key].impact_type.lower()] for key in self.vulnerability_curves][ - 0 - ] + else [ + ImpactType[self.vulnerability_curves[key].impact_type.lower()] + for key in self.vulnerability_curves + ][0] ) # global circulation parameter 'model' is a hint; can be overriden by hazard model super().__init__( - indicator_id="months/spei3m/below/-2" if impact_based_on_a_single_point else "months/spei12m/below/index", + indicator_id="months/spei3m/below/-2" + if impact_based_on_a_single_point + else "months/spei12m/below/index", hazard_type=hazard_type, impact_type=impact_type, ) @@ -241,24 +297,34 @@ def get_data_requests( indicator_id=self.indicator_id, ) - def get_impact(self, asset: Asset, data_responses: List[HazardDataResponse]) -> ImpactDistrib: + def get_impact( + self, asset: Asset, data_responses: List[HazardDataResponse] + ) -> ImpactDistrib: assert isinstance(asset, ThermalPowerGeneratingAsset) - hazard_paths = [cast(HazardParameterDataResponse, r).path for r in data_responses] + hazard_paths = [ + cast(HazardParameterDataResponse, r).path for r in data_responses + ] # The unit being number of months per year, we divide by 12 to express the result as a year fraction. - intensities = np.array(cast(HazardParameterDataResponse, data_responses[0]).parameters / 12.0) + intensities = np.array( + cast(HazardParameterDataResponse, data_responses[0]).parameters / 12.0 + ) if len(intensities) == 1: thresholds = np.array([-2.0]) # hard-coded probability_bins = intensities else: - thresholds = np.array(cast(HazardParameterDataResponse, data_responses[0]).param_defns) + thresholds = np.array( + cast(HazardParameterDataResponse, data_responses[0]).param_defns + ) probability_bins = intensities[:-1] - intensities[1:] probability_bins = np.append(probability_bins, intensities[-1]) curves: List[VulnerabilityCurve] = [] if asset.turbine is None: - curves = [self.vulnerability_curves[key] for key in self.vulnerability_curves] + curves = [ + self.vulnerability_curves[key] for key in self.vulnerability_curves + ] elif asset.cooling is not None: key = "/".join([asset.turbine.name, asset.cooling.name]) if key in self.vulnerability_curves: @@ -271,13 +337,21 @@ def get_impact(self, asset: Asset, data_responses: List[HazardDataResponse]) -> impact = 0.0 denominator = norm.cdf(thresholds[0]) for curve in curves: - cdf = np.array([min(norm.cdf(threshold) / denominator, 1.0) for threshold in curve.intensity]) + cdf = np.array( + [ + min(norm.cdf(threshold) / denominator, 1.0) + for threshold in curve.intensity + ] + ) impact = max( impact, curve.impact_mean[-1] * cdf[-1] + np.sum( (cdf[:-1] - cdf[1:]) - * (np.array(curve.impact_mean[:-1]) + np.array(curve.impact_mean[1:])) + * ( + np.array(curve.impact_mean[:-1]) + + np.array(curve.impact_mean[1:]) + ) / 2 ), ) @@ -286,7 +360,14 @@ def get_impact(self, asset: Asset, data_responses: List[HazardDataResponse]) -> impact_bins = np.array( [ np.max( - [np.interp(threshold, curve.intensity[::-1], curve.impact_mean[::-1]) for curve in curves] + [ + np.interp( + threshold, + curve.intensity[::-1], + curve.impact_mean[::-1], + ) + for curve in curves + ] ) for threshold in thresholds ] @@ -296,7 +377,13 @@ def get_impact(self, asset: Asset, data_responses: List[HazardDataResponse]) -> impact_bins = np.append(impact_bins, impact_bins[-1]) - impact_distrib = ImpactDistrib(self.hazard_type, impact_bins, probability_bins, hazard_paths, self.impact_type) + impact_distrib = ImpactDistrib( + self.hazard_type, + impact_bins, + probability_bins, + hazard_paths, + self.impact_type, + ) return impact_distrib @@ -307,7 +394,12 @@ class ThermalPowerGenerationAirTemperatureModel(VulnerabilityModelBase): _default_resource = "WRI thermal power plant physical climate vulnerability factors" _default_temperatures = [25.0, 30.0, 35.0, 40.0, 45.0, 50.0, 55.0] - def __init__(self, *, resource: str = _default_resource, temperatures: List[float] = _default_temperatures): + def __init__( + self, + *, + resource: str = _default_resource, + temperatures: List[float] = _default_temperatures, + ): """ Air temperature vulnerability model for thermal power generation. @@ -315,26 +407,39 @@ def __init__(self, *, resource: str = _default_resource, temperatures: List[floa resource (str): embedded resource identifier used to infer vulnerability table. temperatures (list[Float]): thresholds of the "days with average temperature above". """ - curve_set: VulnerabilityCurves = get_vulnerability_curves_from_resource(resource) + curve_set: VulnerabilityCurves = get_vulnerability_curves_from_resource( + resource + ) # for this model, key for looking up curves is asset_type, e.g. 'Steam/Recirculating' - self.vulnerability_curves = dict((c.asset_type, c) for c in curve_set.items if c.event_type == "AirTemperature") + self.vulnerability_curves = dict( + (c.asset_type, c) + for c in curve_set.items + if c.event_type == "AirTemperature" + ) self.vuln_curves_by_type = defaultdict(list) for key in self.vulnerability_curves: - self.vuln_curves_by_type[TurbineKind[key.split("/")[0]]].append(self.vulnerability_curves[key]) + self.vuln_curves_by_type[TurbineKind[key.split("/")[0]]].append( + self.vulnerability_curves[key] + ) impact_type = ( ImpactType.disruption if len(self.vulnerability_curves) == 0 - else [ImpactType[self.vulnerability_curves[key].impact_type.lower()] for key in self.vulnerability_curves][ - 0 - ] + else [ + ImpactType[self.vulnerability_curves[key].impact_type.lower()] + for key in self.vulnerability_curves + ][0] ) self.temperatures = temperatures # global circulation parameter 'model' is a hint; can be overriden by hazard model - super().__init__(indicator_id="days_tas/above/{temp_c}c", hazard_type=AirTemperature, impact_type=impact_type) + super().__init__( + indicator_id="days_tas/above/{temp_c}c", + hazard_type=AirTemperature, + impact_type=impact_type, + ) def get_data_requests( self, asset: Asset, *, scenario: str, year: int @@ -364,12 +469,16 @@ def get_data_requests( ) return data_request - def get_impact(self, asset: Asset, data_responses: List[HazardDataResponse]) -> ImpactDistrib: + def get_impact( + self, asset: Asset, data_responses: List[HazardDataResponse] + ) -> ImpactDistrib: assert isinstance(asset, ThermalPowerGeneratingAsset) assert 2 * len(self.temperatures) == len(data_responses) - hazard_paths = [cast(HazardParameterDataResponse, r).path for r in data_responses] + hazard_paths = [ + cast(HazardParameterDataResponse, r).path for r in data_responses + ] # The unit being number of days per year, we divide by 365 to express the result as a year fraction. baseline = [ @@ -394,7 +503,9 @@ def get_impact(self, asset: Asset, data_responses: List[HazardDataResponse]) -> curves: List[VulnerabilityCurve] = [] if asset.turbine is None: - curves = [self.vulnerability_curves[key] for key in self.vulnerability_curves] + curves = [ + self.vulnerability_curves[key] for key in self.vulnerability_curves + ] elif asset.cooling is not None: key = "/".join([asset.turbine.name, asset.cooling.name]) if key in self.vulnerability_curves: @@ -413,7 +524,11 @@ def get_impact(self, asset: Asset, data_responses: List[HazardDataResponse]) -> if temperature < design_air_temperature else np.max( [ - np.interp(temperature - design_air_temperature, curve.intensity, curve.impact_mean) + np.interp( + temperature - design_air_temperature, + curve.intensity, + curve.impact_mean, + ) for curve in curves ] ) @@ -427,7 +542,13 @@ def get_impact(self, asset: Asset, data_responses: List[HazardDataResponse]) -> impact_bins = np.append(impact_bins, impact_bins[-1]) - impact_distrib = ImpactDistrib(self.hazard_type, impact_bins, probability_bins, hazard_paths, self.impact_type) + impact_distrib = ImpactDistrib( + self.hazard_type, + impact_bins, + probability_bins, + hazard_paths, + self.impact_type, + ) return impact_distrib @@ -438,7 +559,12 @@ class ThermalPowerGenerationWaterTemperatureModel(VulnerabilityModelBase): _default_resource = "WRI thermal power plant physical climate vulnerability factors" _default_correlation = 0.5 - def __init__(self, *, resource: str = _default_resource, correlation: float = _default_correlation): + def __init__( + self, + *, + resource: str = _default_resource, + correlation: float = _default_correlation, + ): """ Water temperature vulnerability model for thermal power generation. @@ -447,33 +573,47 @@ def __init__(self, *, resource: str = _default_resource, correlation: float = _d correlation (float): correlation specifying the Gaussian copula which joins the marginal distributions of water temperature and WBGT. """ - curve_set: VulnerabilityCurves = get_vulnerability_curves_from_resource(resource) + curve_set: VulnerabilityCurves = get_vulnerability_curves_from_resource( + resource + ) self.gaussian_copula = multivariate_normal( - mean=np.array([0.0, 0.0]), cov=np.array([[1.0, correlation], [correlation, 1.0]]) + mean=np.array([0.0, 0.0]), + cov=np.array([[1.0, correlation], [correlation, 1.0]]), ) # for this model, key for looking up curves is asset_type, e.g. 'Steam/Recirculating' self.vulnerability_curves = dict( - (c.asset_type, c) for c in curve_set.items if c.event_type == "WaterTemperature" + (c.asset_type, c) + for c in curve_set.items + if c.event_type == "WaterTemperature" ) self.vuln_curves_by_type = defaultdict(list) for key in self.vulnerability_curves: - self.vuln_curves_by_type[TurbineKind[key.split("/")[0]]].append(self.vulnerability_curves[key]) + self.vuln_curves_by_type[TurbineKind[key.split("/")[0]]].append( + self.vulnerability_curves[key] + ) self.regulatory_discharge_curves = dict( - (c.asset_type, c) for c in curve_set.items if c.event_type == "RegulatoryDischargeWaterLimit" + (c.asset_type, c) + for c in curve_set.items + if c.event_type == "RegulatoryDischargeWaterLimit" ) impact_type = ( ImpactType.disruption if len(self.vulnerability_curves) == 0 - else [ImpactType[self.vulnerability_curves[key].impact_type.lower()] for key in self.vulnerability_curves][ - 0 - ] + else [ + ImpactType[self.vulnerability_curves[key].impact_type.lower()] + for key in self.vulnerability_curves + ][0] ) # global circulation parameter 'model' is a hint; can be overriden by hazard model - super().__init__(indicator_id="weeks_water_temp_above", hazard_type=WaterTemperature, impact_type=impact_type) + super().__init__( + indicator_id="weeks_water_temp_above", + hazard_type=WaterTemperature, + impact_type=impact_type, + ) def get_data_requests( self, asset: Asset, *, scenario: str, year: int @@ -521,11 +661,15 @@ def get_data_requests( ) return data_request - def get_impact(self, asset: Asset, data_responses: List[HazardDataResponse]) -> ImpactDistrib: + def get_impact( + self, asset: Asset, data_responses: List[HazardDataResponse] + ) -> ImpactDistrib: assert isinstance(asset, ThermalPowerGeneratingAsset) assert len(data_responses) == 4 - hazard_paths = [cast(HazardParameterDataResponse, r).path for r in data_responses] + hazard_paths = [ + cast(HazardParameterDataResponse, r).path for r in data_responses + ] # Water temperature below which the power plant does not experience generation losses. design_intake_water_temperature = cast( @@ -533,7 +677,9 @@ def get_impact(self, asset: Asset, data_responses: List[HazardDataResponse]) -> np.interp( 0.9, # The unit being number of weeks per year, we divide by 52 to express the result as a year fraction. - 1.0 - cast(HazardParameterDataResponse, data_responses[1]).parameters / 52.0, + 1.0 + - cast(HazardParameterDataResponse, data_responses[1]).parameters + / 52.0, cast(HazardParameterDataResponse, data_responses[1]).param_defns, ), ) @@ -551,7 +697,9 @@ def get_impact(self, asset: Asset, data_responses: List[HazardDataResponse]) -> np.interp( 0.99, # The unit being number of days per year, we divide by 365 to express the result as a year fraction. - 1.0 - cast(HazardParameterDataResponse, data_responses[3]).parameters / 365.0, + 1.0 + - cast(HazardParameterDataResponse, data_responses[3]).parameters + / 365.0, cast(HazardParameterDataResponse, data_responses[3]).param_defns, ), ) @@ -566,20 +714,37 @@ def get_impact(self, asset: Asset, data_responses: List[HazardDataResponse]) -> / 365.0, ) - intake_water_temperatures = cast(HazardParameterDataResponse, data_responses[0]).param_defns - intake_water_temperature_intensities = cast(HazardParameterDataResponse, data_responses[0]).parameters / 52.0 + intake_water_temperatures = cast( + HazardParameterDataResponse, data_responses[0] + ).param_defns + intake_water_temperature_intensities = ( + cast(HazardParameterDataResponse, data_responses[0]).parameters / 52.0 + ) - if impact_scale_for_recirculating_steam_unit <= 0.0 or 1.0 <= impact_scale_for_recirculating_steam_unit: - intake_water_temperature_intensities_for_recirculating_steam_unit = intake_water_temperature_intensities + if ( + impact_scale_for_recirculating_steam_unit <= 0.0 + or 1.0 <= impact_scale_for_recirculating_steam_unit + ): + intake_water_temperature_intensities_for_recirculating_steam_unit = ( + intake_water_temperature_intensities + ) else: - gaussian_threshold: float = norm.ppf(impact_scale_for_recirculating_steam_unit) + gaussian_threshold: float = norm.ppf( + impact_scale_for_recirculating_steam_unit + ) intake_water_temperature_intensities_for_recirculating_steam_unit = np.array( [ ( max(0.0, min(1.0, intake_water_temperature_intensity)) - if intake_water_temperature_intensity <= 0.0 or 1.0 <= intake_water_temperature_intensity + if intake_water_temperature_intensity <= 0.0 + or 1.0 <= intake_water_temperature_intensity else self.gaussian_copula.cdf( - np.array([norm.ppf(intake_water_temperature_intensity), gaussian_threshold]) + np.array( + [ + norm.ppf(intake_water_temperature_intensity), + gaussian_threshold, + ] + ) ) / impact_scale_for_recirculating_steam_unit ) @@ -588,24 +753,30 @@ def get_impact(self, asset: Asset, data_responses: List[HazardDataResponse]) -> ) intake_water_temperature_probability_bins = ( - intake_water_temperature_intensities[:-1] - intake_water_temperature_intensities[1:] + intake_water_temperature_intensities[:-1] + - intake_water_temperature_intensities[1:] ) intake_water_temperature_probability_bins = np.append( - intake_water_temperature_probability_bins, intake_water_temperature_intensities[-1] + intake_water_temperature_probability_bins, + intake_water_temperature_intensities[-1], ) intake_water_temperature_probability_bins_for_recirculating_steam_unit = ( intake_water_temperature_intensities_for_recirculating_steam_unit[:-1] - intake_water_temperature_intensities_for_recirculating_steam_unit[1:] ) - intake_water_temperature_probability_bins_for_recirculating_steam_unit = np.append( - intake_water_temperature_probability_bins_for_recirculating_steam_unit, - intake_water_temperature_intensities_for_recirculating_steam_unit[-1], + intake_water_temperature_probability_bins_for_recirculating_steam_unit = ( + np.append( + intake_water_temperature_probability_bins_for_recirculating_steam_unit, + intake_water_temperature_intensities_for_recirculating_steam_unit[-1], + ) ) curves: List[VulnerabilityCurve] = [] if asset.turbine is None: - curves = [self.vulnerability_curves[key] for key in self.vulnerability_curves] + curves = [ + self.vulnerability_curves[key] for key in self.vulnerability_curves + ] elif asset.cooling is not None: key = "/".join([asset.turbine.name, asset.cooling.name]) if key in self.vulnerability_curves: @@ -641,7 +812,9 @@ def get_impact(self, asset: Asset, data_responses: List[HazardDataResponse]) -> ] ) if curve.asset_type in self.regulatory_discharge_curves: - regulatory_discharge_curve = self.regulatory_discharge_curves[curve.asset_type] + regulatory_discharge_curve = self.regulatory_discharge_curves[ + curve.asset_type + ] impact_bins = np.array( [ max( @@ -655,22 +828,37 @@ def get_impact(self, asset: Asset, data_responses: List[HazardDataResponse]) -> ), ), ) - for impact, intake_water_temperature in zip(impact_bins, intake_water_temperatures) + for impact, intake_water_temperature in zip( + impact_bins, intake_water_temperatures + ) ] ) impact_bins = np.append(impact_bins, impact_bins[-1]) impact_distrib_by_curve.append( - ImpactDistrib(self.hazard_type, impact_bins, probability_bins, hazard_paths, self.impact_type) + ImpactDistrib( + self.hazard_type, + impact_bins, + probability_bins, + hazard_paths, + self.impact_type, + ) ) if 0 < len(impact_distrib_by_curve): - impact_distrib = sorted(impact_distrib_by_curve, key=lambda x: x.mean_impact())[-1] + impact_distrib = sorted( + impact_distrib_by_curve, key=lambda x: x.mean_impact() + )[-1] else: impact_distrib = ImpactDistrib( self.hazard_type, - [0.0 for _ in range(0, len(intake_water_temperature_probability_bins) + 1)], + [ + 0.0 + for _ in range( + 0, len(intake_water_temperature_probability_bins) + 1 + ) + ], intake_water_temperature_probability_bins, hazard_paths, self.impact_type, @@ -692,24 +880,33 @@ def __init__(self, *, resource: str = _default_resource): Args: resource (str): embedded resource identifier used to infer vulnerability table. """ - curve_set: VulnerabilityCurves = get_vulnerability_curves_from_resource(resource) + curve_set: VulnerabilityCurves = get_vulnerability_curves_from_resource( + resource + ) # for this model, key for looking up curves is asset_type, e.g. 'Steam/Recirculating' - self.vulnerability_curves = dict((c.asset_type, c) for c in curve_set.items if c.event_type == "WaterStress") + self.vulnerability_curves = dict( + (c.asset_type, c) for c in curve_set.items if c.event_type == "WaterStress" + ) self.vuln_curves_by_type = defaultdict(list) for key in self.vulnerability_curves: - self.vuln_curves_by_type[TurbineKind[key.split("/")[0]]].append(self.vulnerability_curves[key]) + self.vuln_curves_by_type[TurbineKind[key.split("/")[0]]].append( + self.vulnerability_curves[key] + ) impact_type = ( ImpactType.disruption if len(self.vulnerability_curves) == 0 - else [ImpactType[self.vulnerability_curves[key].impact_type.lower()] for key in self.vulnerability_curves][ - 0 - ] + else [ + ImpactType[self.vulnerability_curves[key].impact_type.lower()] + for key in self.vulnerability_curves + ][0] ) # global circulation parameter 'model' is a hint; can be overriden by hazard model - super().__init__(indicator_id="water_stress", hazard_type=WaterRisk, impact_type=impact_type) + super().__init__( + indicator_id="water_stress", hazard_type=WaterRisk, impact_type=impact_type + ) def get_data_requests( self, asset: Asset, *, scenario: str, year: int @@ -747,11 +944,15 @@ def get_data_requests( ) return data_request - def get_impact(self, asset: Asset, data_responses: List[HazardDataResponse]) -> ImpactDistrib: + def get_impact( + self, asset: Asset, data_responses: List[HazardDataResponse] + ) -> ImpactDistrib: assert isinstance(asset, ThermalPowerGeneratingAsset) assert len(data_responses) == 3 - hazard_paths = [cast(HazardParameterDataResponse, r).path for r in data_responses] + hazard_paths = [ + cast(HazardParameterDataResponse, r).path for r in data_responses + ] if ( len(cast(HazardParameterDataResponse, data_responses[0]).parameters) == 0 @@ -762,19 +963,31 @@ def get_impact(self, asset: Asset, data_responses: List[HazardDataResponse]) -> # We (naively) assume that water stress follows a shifted uniform distribution: water_stress - 0.5 + U(0,1): probability_water_stress_above_40pct = max( - 0.0, min(1.0, 0.1 + cast(HazardParameterDataResponse, data_responses[0]).parameter) + 0.0, + min( + 1.0, + 0.1 + cast(HazardParameterDataResponse, data_responses[0]).parameter, + ), ) - baseline_water_supply = cast(HazardParameterDataResponse, data_responses[2]).parameter + baseline_water_supply = cast( + HazardParameterDataResponse, data_responses[2] + ).parameter supply_reduction_rate = ( 0.0 if baseline_water_supply == 0.0 - else (cast(HazardParameterDataResponse, data_responses[1]).parameter / baseline_water_supply - 1.0) + else ( + cast(HazardParameterDataResponse, data_responses[1]).parameter + / baseline_water_supply + - 1.0 + ) ) curves: List[VulnerabilityCurve] = [] if asset.turbine is None: - curves = [self.vulnerability_curves[key] for key in self.vulnerability_curves] + curves = [ + self.vulnerability_curves[key] for key in self.vulnerability_curves + ] elif asset.cooling is not None: key = "/".join([asset.turbine.name, asset.cooling.name]) if key in self.vulnerability_curves: @@ -783,12 +996,23 @@ def get_impact(self, asset: Asset, data_responses: List[HazardDataResponse]) -> curves = self.vuln_curves_by_type[asset.turbine] impact = ( - np.max([np.interp(-supply_reduction_rate, curve.intensity, curve.impact_mean) for curve in curves]) + np.max( + [ + np.interp( + -supply_reduction_rate, curve.intensity, curve.impact_mean + ) + for curve in curves + ] + ) if 0 < len(curves) else 0.0 ) impact_distrib = ImpactDistrib( - self.hazard_type, [impact, impact], [probability_water_stress_above_40pct], hazard_paths, self.impact_type + self.hazard_type, + [impact, impact], + [probability_water_stress_above_40pct], + hazard_paths, + self.impact_type, ) return impact_distrib diff --git a/tests/api/container_test.py b/tests/api/container_test.py index 819b6434..c0875190 100644 --- a/tests/api/container_test.py +++ b/tests/api/container_test.py @@ -1,7 +1,6 @@ import fsspec.implementations.local as local from dependency_injector import containers, providers -from physrisk.container import Container from physrisk.data.inventory_reader import InventoryReader from ..data.hazard_model_store_test import mock_hazard_model_store_heat @@ -12,6 +11,8 @@ class TestContainer(containers.DeclarativeContainer): config = providers.Configuration(default={"zarr_sources": ["embedded"]}) - inventory_reader = providers.Singleton(lambda: InventoryReader(fs=local.LocalFileSystem(), base_path="")) + inventory_reader = providers.Singleton( + lambda: InventoryReader(fs=local.LocalFileSystem(), base_path="") + ) zarr_store = providers.Singleton(lambda: mock_hazard_model_store_heat([0], [0])) diff --git a/tests/api/data_requests_test.py b/tests/api/data_requests_test.py index 59dd2854..6076751a 100644 --- a/tests/api/data_requests_test.py +++ b/tests/api/data_requests_test.py @@ -40,13 +40,20 @@ def test_hazard_data_description(self): # test that validation passes: container = Container() requester = container.requester - _ = requester.get(request_id="get_hazard_data_description", request_dict={"paths": ["test_path.md"]}) + _ = requester.get( + request_id="get_hazard_data_description", + request_dict={"paths": ["test_path.md"]}, + ) def test_generic_source_path(self): inventory = EmbeddedInventory() source_paths = get_default_source_paths(inventory) - result_heat = source_paths[ChronicHeat](indicator_id="mean_degree_days/above/32c", scenario="rcp8p5", year=2050) - result_flood = source_paths[RiverineInundation](indicator_id="flood_depth", scenario="rcp8p5", year=2050) + result_heat = source_paths[ChronicHeat]( + indicator_id="mean_degree_days/above/32c", scenario="rcp8p5", year=2050 + ) + result_flood = source_paths[RiverineInundation]( + indicator_id="flood_depth", scenario="rcp8p5", year=2050 + ) result_flood_hist = source_paths[RiverineInundation]( indicator_id="flood_depth", scenario="historical", year=2080 ) @@ -54,13 +61,24 @@ def test_generic_source_path(self): indicator_id="mean_degree_days/above/32c", scenario="rcp8p5", year=2050, - hint=HazardDataHint(path="chronic_heat/osc/v2/mean_degree_days_v2_above_32c_CMCC-ESM2_{scenario}_{year}"), + hint=HazardDataHint( + path="chronic_heat/osc/v2/mean_degree_days_v2_above_32c_CMCC-ESM2_{scenario}_{year}" + ), ) - assert result_heat == "chronic_heat/osc/v2/mean_degree_days_v2_above_32c_ACCESS-CM2_rcp8p5_2050" + assert ( + result_heat + == "chronic_heat/osc/v2/mean_degree_days_v2_above_32c_ACCESS-CM2_rcp8p5_2050" + ) assert result_flood == "inundation/wri/v2/inunriver_rcp8p5_MIROC-ESM-CHEM_2050" - assert result_flood_hist == "inundation/wri/v2/inunriver_historical_000000000WATCH_1980" - assert result_heat_hint == "chronic_heat/osc/v2/mean_degree_days_v2_above_32c_CMCC-ESM2_rcp8p5_2050" + assert ( + result_flood_hist + == "inundation/wri/v2/inunriver_historical_000000000WATCH_1980" + ) + assert ( + result_heat_hint + == "chronic_heat/osc/v2/mean_degree_days_v2_above_32c_CMCC-ESM2_rcp8p5_2050" + ) def test_zarr_reading(self): request_dict = { @@ -68,7 +86,9 @@ def test_zarr_reading(self): { "request_item_id": "test_inundation", "event_type": "RiverineInundation", - "longitudes": TestData.longitudes[0:3], # coords['longitudes'][0:100], + "longitudes": TestData.longitudes[ + 0:3 + ], # coords['longitudes'][0:100], "latitudes": TestData.latitudes[0:3], # coords['latitudes'][0:100], "year": 2080, "scenario": "rcp8p5", @@ -84,14 +104,22 @@ def test_zarr_reading(self): result = requests._get_hazard_data( request, - ZarrHazardModel(source_paths=get_default_source_paths(EmbeddedInventory()), reader=ZarrReader(store=store)), + ZarrHazardModel( + source_paths=get_default_source_paths(EmbeddedInventory()), + reader=ZarrReader(store=store), + ), ) - numpy.testing.assert_array_almost_equal_nulp(result.items[0].intensity_curve_set[0].intensities, np.zeros((9))) numpy.testing.assert_array_almost_equal_nulp( - result.items[0].intensity_curve_set[1].intensities, np.linspace(0.1, 1.0, 9, dtype="f4") + result.items[0].intensity_curve_set[0].intensities, np.zeros((9)) + ) + numpy.testing.assert_array_almost_equal_nulp( + result.items[0].intensity_curve_set[1].intensities, + np.linspace(0.1, 1.0, 9, dtype="f4"), + ) + numpy.testing.assert_array_almost_equal_nulp( + result.items[0].intensity_curve_set[2].intensities, np.zeros((9)) ) - numpy.testing.assert_array_almost_equal_nulp(result.items[0].intensity_curve_set[2].intensities, np.zeros((9))) def test_zarr_reading_chronic(self): request_dict = { @@ -100,7 +128,9 @@ def test_zarr_reading_chronic(self): { "request_item_id": "test_inundation", "event_type": "ChronicHeat", - "longitudes": TestData.longitudes[0:3], # coords['longitudes'][0:100], + "longitudes": TestData.longitudes[ + 0:3 + ], # coords['longitudes'][0:100], "latitudes": TestData.latitudes[0:3], # coords['latitudes'][0:100], "year": 2050, "scenario": "ssp585", @@ -115,9 +145,12 @@ def test_zarr_reading_chronic(self): source_paths = get_default_source_paths(EmbeddedInventory()) result = requests._get_hazard_data( - request, ZarrHazardModel(source_paths=source_paths, reader=ZarrReader(store)) + request, + ZarrHazardModel(source_paths=source_paths, reader=ZarrReader(store)), + ) + numpy.testing.assert_array_almost_equal_nulp( + result.items[0].intensity_curve_set[0].intensities[0], 600.0 ) - numpy.testing.assert_array_almost_equal_nulp(result.items[0].intensity_curve_set[0].intensities[0], 600.0) # request_with_hint = request.copy() # request_with_hint.items[0].path = "chronic_heat/osc/v2/mean_degree_days_v2_above_32c_CMCC-ESM2_rcp8p5_2050" @@ -166,13 +199,21 @@ def test_zarr_reading_live(self): ], } - response_floor = requester.get(request_id="get_hazard_data", request_dict=request1) + response_floor = requester.get( + request_id="get_hazard_data", request_dict=request1 + ) request1["interpolation"] = "linear" # type: ignore - response_linear = requester.get(request_id="get_hazard_data", request_dict=request1) + response_linear = requester.get( + request_id="get_hazard_data", request_dict=request1 + ) print(response_linear) - floor = json.loads(response_floor)["items"][0]["intensity_curve_set"][5]["intensities"] - linear = json.loads(response_linear)["items"][0]["intensity_curve_set"][5]["intensities"] + floor = json.loads(response_floor)["items"][0]["intensity_curve_set"][5][ + "intensities" + ] + linear = json.loads(response_linear)["items"][0]["intensity_curve_set"][5][ + "intensities" + ] print(floor) print(linear) diff --git a/tests/api/impact_requests_test.py b/tests/api/impact_requests_test.py index 4a6e86c4..3e0ad05f 100644 --- a/tests/api/impact_requests_test.py +++ b/tests/api/impact_requests_test.py @@ -12,7 +12,11 @@ from physrisk.data.pregenerated_hazard_model import ZarrHazardModel from physrisk.data.zarr_reader import ZarrReader from physrisk.hazard_models.core_hazards import get_default_source_paths -from physrisk.kernel.assets import PowerGeneratingAsset, RealEstateAsset, ThermalPowerGeneratingAsset +from physrisk.kernel.assets import ( + PowerGeneratingAsset, + RealEstateAsset, + ThermalPowerGeneratingAsset, +) from physrisk.kernel.vulnerability_model import DictBasedVulnerabilityModels from physrisk.vulnerability_models.power_generating_asset_models import InundationModel from physrisk.vulnerability_models.real_estate_models import ( @@ -115,14 +119,21 @@ def test_impact_request(self): request = requests.AssetImpactRequest(**request_dict) # type: ignore - curve = np.array([0.0596, 0.333, 0.505, 0.715, 0.864, 1.003, 1.149, 1.163, 1.163]) - store = mock_hazard_model_store_inundation(TestData.longitudes, TestData.latitudes, curve) + curve = np.array( + [0.0596, 0.333, 0.505, 0.715, 0.864, 1.003, 1.149, 1.163, 1.163] + ) + store = mock_hazard_model_store_inundation( + TestData.longitudes, TestData.latitudes, curve + ) source_paths = get_default_source_paths(EmbeddedInventory()) vulnerability_models = DictBasedVulnerabilityModels( { PowerGeneratingAsset: [InundationModel()], - RealEstateAsset: [RealEstateCoastalInundationModel(), RealEstateRiverineInundationModel()], + RealEstateAsset: [ + RealEstateCoastalInundationModel(), + RealEstateRiverineInundationModel(), + ], } ) @@ -132,7 +143,9 @@ def test_impact_request(self): vulnerability_models=vulnerability_models, ) - self.assertEqual(response.asset_impacts[0].impacts[0].hazard_type, "CoastalInundation") + self.assertEqual( + response.asset_impacts[0].impacts[0].hazard_type, "CoastalInundation" + ) def test_risk_model_impact_request(self): """Tests the risk model functionality of the impact request.""" @@ -167,14 +180,21 @@ def test_risk_model_impact_request(self): request = requests.AssetImpactRequest(**request_dict) # type: ignore - curve = np.array([0.0596, 0.333, 0.505, 0.715, 0.864, 1.003, 1.149, 1.163, 1.163]) - store = mock_hazard_model_store_inundation(TestData.longitudes, TestData.latitudes, curve) + curve = np.array( + [0.0596, 0.333, 0.505, 0.715, 0.864, 1.003, 1.149, 1.163, 1.163] + ) + store = mock_hazard_model_store_inundation( + TestData.longitudes, TestData.latitudes, curve + ) source_paths = get_default_source_paths(EmbeddedInventory()) vulnerability_models = DictBasedVulnerabilityModels( { PowerGeneratingAsset: [InundationModel()], - RealEstateAsset: [RealEstateCoastalInundationModel(), RealEstateRiverineInundationModel()], + RealEstateAsset: [ + RealEstateCoastalInundationModel(), + RealEstateRiverineInundationModel(), + ], } ) @@ -184,10 +204,11 @@ def test_risk_model_impact_request(self): vulnerability_models=vulnerability_models, ) - self.assertEqual(response.asset_impacts[0].impacts[0].hazard_type, "CoastalInundation") + self.assertEqual( + response.asset_impacts[0].impacts[0].hazard_type, "CoastalInundation" + ) def test_thermal_power_generation(self): - latitudes = np.array([32.6017]) longitudes = np.array([-87.7811]) @@ -501,7 +522,23 @@ def test_thermal_power_generation(self): ) # Add mock water temperature data: - return_periods = [5, 7.5, 10, 12.5, 15, 17.5, 20, 22.5, 25, 27.5, 30, 32.5, 35, 37.5, 40] + return_periods = [ + 5, + 7.5, + 10, + 12.5, + 15, + 17.5, + 20, + 22.5, + 25, + 27.5, + 30, + 32.5, + 35, + 37.5, + 40, + ] shape, t = shape_transform_21600_43200(return_periods=return_periods) add_curves( root, @@ -537,7 +574,25 @@ def test_thermal_power_generation(self): latitudes, "chronic_heat/nluu/v2/weeks_water_temp_above_GFDL_rcp8p5_2050", shape, - np.array([51.85, 51.5, 50.25, 46.75, 41.95, 35.35, 29.4, 24.55, 20.15, 13.85, 6.75, 3.5, 1.3, 0.25, 0.1]), + np.array( + [ + 51.85, + 51.5, + 50.25, + 46.75, + 41.95, + 35.35, + 29.4, + 24.55, + 20.15, + 13.85, + 6.75, + 3.5, + 1.3, + 0.25, + 0.1, + ] + ), return_periods, t, ) @@ -577,7 +632,20 @@ def test_thermal_power_generation(self): "chronic_heat/osc/v2/days_wbgt_above_ACCESS-CM2_historical_2005", shape, np.array( - [361.95273, 342.51804, 278.8146, 213.5123, 157.4511, 101.78238, 12.6897545, 0.0, 0.0, 0.0, 0.0, 0.0] + [ + 361.95273, + 342.51804, + 278.8146, + 213.5123, + 157.4511, + 101.78238, + 12.6897545, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + ] ), return_periods, t, @@ -604,49 +672,95 @@ def test_thermal_power_generation(self): ) # Air Temperature - self.assertAlmostEqual(response.asset_impacts[0].impacts[0].impact_mean, 0.0075618606988512764) - self.assertAlmostEqual(response.asset_impacts[1].impacts[0].impact_mean, 0.0075618606988512764) - self.assertAlmostEqual(response.asset_impacts[2].impacts[0].impact_mean, 0.0025192163596997963) - self.assertAlmostEqual(response.asset_impacts[3].impacts[0].impact_mean, 0.0025192163596997963) + self.assertAlmostEqual( + response.asset_impacts[0].impacts[0].impact_mean, 0.0075618606988512764 + ) + self.assertAlmostEqual( + response.asset_impacts[1].impacts[0].impact_mean, 0.0075618606988512764 + ) + self.assertAlmostEqual( + response.asset_impacts[2].impacts[0].impact_mean, 0.0025192163596997963 + ) + self.assertAlmostEqual( + response.asset_impacts[3].impacts[0].impact_mean, 0.0025192163596997963 + ) self.assertAlmostEqual(response.asset_impacts[4].impacts[0].impact_mean, 0.0) self.assertAlmostEqual(response.asset_impacts[5].impacts[0].impact_mean, 0.0) # Drought - self.assertAlmostEqual(response.asset_impacts[0].impacts[1].impact_mean, 0.0008230079663917424) + self.assertAlmostEqual( + response.asset_impacts[0].impacts[1].impact_mean, 0.0008230079663917424 + ) self.assertAlmostEqual(response.asset_impacts[1].impacts[1].impact_mean, 0.0) - self.assertAlmostEqual(response.asset_impacts[2].impacts[1].impact_mean, 0.0008230079663917424) + self.assertAlmostEqual( + response.asset_impacts[2].impacts[1].impact_mean, 0.0008230079663917424 + ) self.assertAlmostEqual(response.asset_impacts[3].impacts[1].impact_mean, 0.0) - self.assertAlmostEqual(response.asset_impacts[4].impacts[1].impact_mean, 0.0008230079663917424) - self.assertAlmostEqual(response.asset_impacts[5].impacts[1].impact_mean, 0.0008230079663917424) + self.assertAlmostEqual( + response.asset_impacts[4].impacts[1].impact_mean, 0.0008230079663917424 + ) + self.assertAlmostEqual( + response.asset_impacts[5].impacts[1].impact_mean, 0.0008230079663917424 + ) # Riverine Inundation - self.assertAlmostEqual(response.asset_impacts[0].impacts[2].impact_mean, 0.0046864436945997625) - self.assertAlmostEqual(response.asset_impacts[1].impacts[2].impact_mean, 0.0046864436945997625) - self.assertAlmostEqual(response.asset_impacts[2].impacts[2].impact_mean, 0.0046864436945997625) - self.assertAlmostEqual(response.asset_impacts[3].impacts[2].impact_mean, 0.0046864436945997625) - self.assertAlmostEqual(response.asset_impacts[4].impacts[2].impact_mean, 0.0046864436945997625) - self.assertAlmostEqual(response.asset_impacts[5].impacts[2].impact_mean, 0.0046864436945997625) + self.assertAlmostEqual( + response.asset_impacts[0].impacts[2].impact_mean, 0.0046864436945997625 + ) + self.assertAlmostEqual( + response.asset_impacts[1].impacts[2].impact_mean, 0.0046864436945997625 + ) + self.assertAlmostEqual( + response.asset_impacts[2].impacts[2].impact_mean, 0.0046864436945997625 + ) + self.assertAlmostEqual( + response.asset_impacts[3].impacts[2].impact_mean, 0.0046864436945997625 + ) + self.assertAlmostEqual( + response.asset_impacts[4].impacts[2].impact_mean, 0.0046864436945997625 + ) + self.assertAlmostEqual( + response.asset_impacts[5].impacts[2].impact_mean, 0.0046864436945997625 + ) # Water Stress - self.assertAlmostEqual(response.asset_impacts[0].impacts[3].impact_mean, 0.010181435900296947) + self.assertAlmostEqual( + response.asset_impacts[0].impacts[3].impact_mean, 0.010181435900296947 + ) self.assertAlmostEqual(response.asset_impacts[1].impacts[3].impact_mean, 0.0) - self.assertAlmostEqual(response.asset_impacts[2].impacts[3].impact_mean, 0.010181435900296947) + self.assertAlmostEqual( + response.asset_impacts[2].impacts[3].impact_mean, 0.010181435900296947 + ) self.assertAlmostEqual(response.asset_impacts[3].impacts[3].impact_mean, 0.0) - self.assertAlmostEqual(response.asset_impacts[4].impacts[3].impact_mean, 0.010181435900296947) - self.assertAlmostEqual(response.asset_impacts[5].impacts[3].impact_mean, 0.010181435900296947) + self.assertAlmostEqual( + response.asset_impacts[4].impacts[3].impact_mean, 0.010181435900296947 + ) + self.assertAlmostEqual( + response.asset_impacts[5].impacts[3].impact_mean, 0.010181435900296947 + ) # Water Temperature - self.assertAlmostEqual(response.asset_impacts[0].impacts[4].impact_mean, 0.1448076958069578) + self.assertAlmostEqual( + response.asset_impacts[0].impacts[4].impact_mean, 0.1448076958069578 + ) self.assertAlmostEqual(response.asset_impacts[1].impacts[4].impact_mean, 0.0) - self.assertAlmostEqual(response.asset_impacts[2].impacts[4].impact_mean, 0.1448076958069578) + self.assertAlmostEqual( + response.asset_impacts[2].impacts[4].impact_mean, 0.1448076958069578 + ) self.assertAlmostEqual(response.asset_impacts[3].impacts[4].impact_mean, 0.0) - self.assertAlmostEqual(response.asset_impacts[4].impacts[4].impact_mean, 0.1448076958069578) - self.assertAlmostEqual(response.asset_impacts[5].impacts[4].impact_mean, 0.005896707722257193) + self.assertAlmostEqual( + response.asset_impacts[4].impacts[4].impact_mean, 0.1448076958069578 + ) + self.assertAlmostEqual( + response.asset_impacts[5].impacts[4].impact_mean, 0.005896707722257193 + ) vulnerability_models = DictBasedVulnerabilityModels( { ThermalPowerGeneratingAsset: [ - ThermalPowerGenerationDroughtModel(impact_based_on_a_single_point=True), + ThermalPowerGenerationDroughtModel( + impact_based_on_a_single_point=True + ), ] } ) @@ -659,12 +773,20 @@ def test_thermal_power_generation(self): ) # Drought (Jupiter) - self.assertAlmostEqual(response.asset_impacts[0].impacts[0].impact_mean, 0.0005859470850072303) + self.assertAlmostEqual( + response.asset_impacts[0].impacts[0].impact_mean, 0.0005859470850072303 + ) self.assertAlmostEqual(response.asset_impacts[1].impacts[0].impact_mean, 0.0) - self.assertAlmostEqual(response.asset_impacts[2].impacts[0].impact_mean, 0.0005859470850072303) + self.assertAlmostEqual( + response.asset_impacts[2].impacts[0].impact_mean, 0.0005859470850072303 + ) self.assertAlmostEqual(response.asset_impacts[3].impacts[0].impact_mean, 0.0) - self.assertAlmostEqual(response.asset_impacts[4].impacts[0].impact_mean, 0.0005859470850072303) - self.assertAlmostEqual(response.asset_impacts[5].impacts[0].impact_mean, 0.0005859470850072303) + self.assertAlmostEqual( + response.asset_impacts[4].impacts[0].impact_mean, 0.0005859470850072303 + ) + self.assertAlmostEqual( + response.asset_impacts[5].impacts[0].impact_mean, 0.0005859470850072303 + ) @unittest.skip("example, not test") def test_example_portfolios(self): @@ -679,7 +801,9 @@ def test_example_portfolios(self): } container = Container() requester = container.requester() - response = requester.get(request_id="get_asset_impact", request_dict=request_dict) + response = requester.get( + request_id="get_asset_impact", request_dict=request_dict + ) with open("out.json", "w") as f: f.write(response) assert response is not None @@ -709,10 +833,21 @@ def test_example_portfolios_risk_measures(self): } container = Container() requester = container.requester() - response = requester.get(request_id="get_asset_impact", request_dict=request_dict) + response = requester.get( + request_id="get_asset_impact", request_dict=request_dict + ) risk_measures_dict = json.loads(response)["risk_measures"] - helper = RiskMeasuresHelper(TypeAdapter(RiskMeasures).validate_python(risk_measures_dict)) - for hazard_type in ["RiverineInundation", "CoastalInundation", "ChronicHeat", "Wind"]: - scores, measure_values, measure_defns = helper.get_measure(hazard_type, "ssp585", 2050) + helper = RiskMeasuresHelper( + TypeAdapter(RiskMeasures).validate_python(risk_measures_dict) + ) + for hazard_type in [ + "RiverineInundation", + "CoastalInundation", + "ChronicHeat", + "Wind", + ]: + scores, measure_values, measure_defns = helper.get_measure( + hazard_type, "ssp585", 2050 + ) label, description = helper.get_score_details(scores[0], measure_defns[0]) print(label) diff --git a/tests/data/events_retrieval_test.py b/tests/data/events_retrieval_test.py index 6c5f216a..9f5f02dd 100644 --- a/tests/data/events_retrieval_test.py +++ b/tests/data/events_retrieval_test.py @@ -9,7 +9,11 @@ from fsspec.implementations.memory import MemoryFileSystem from shapely import Polygon -from physrisk.api.v1.hazard_data import HazardAvailabilityRequest, HazardResource, Scenario +from physrisk.api.v1.hazard_data import ( + HazardAvailabilityRequest, + HazardResource, + Scenario, +) from physrisk.data.inventory import EmbeddedInventory, Inventory from physrisk.data.inventory_reader import InventoryReader from physrisk.data.pregenerated_hazard_model import ZarrHazardModel @@ -20,7 +24,10 @@ # from pathlib import PurePosixPath from ..base_test import TestWithCredentials -from ..data.hazard_model_store_test import ZarrStoreMocker, mock_hazard_model_store_inundation +from ..data.hazard_model_store_test import ( + ZarrStoreMocker, + mock_hazard_model_store_inundation, +) class TestEventRetrieval(TestWithCredentials): @@ -30,14 +37,19 @@ def test_inventory_change(self): embedded = EmbeddedInventory() resources1 = embedded.to_resources() inventory = Inventory(resources1).json_ordered() - with open(os.path.join(os.path.dirname(os.path.abspath(__file__)), "inventory.json"), "w") as f: + with open( + os.path.join(os.path.dirname(os.path.abspath(__file__)), "inventory.json"), + "w", + ) as f: f.write(inventory) def test_hazard_data_availability_summary(self): # check validation passes calling in service-like way inventory = EmbeddedInventory() response = _get_hazard_data_availability( - HazardAvailabilityRequest(sources=["embedded"]), inventory, inventory.colormaps() + HazardAvailabilityRequest(sources=["embedded"]), + inventory, + inventory.colormaps(), ) # , "hazard_test"]) assert len(response.models) > 0 # rely on Pydantic validation for test @@ -91,8 +103,12 @@ def test_zarr_bilinear(self): image_coords_surr = np.stack( [ - np.concatenate([np.floor(x), np.floor(x) + 1, np.floor(x), np.floor(x) + 1]), - np.concatenate([np.floor(y), np.floor(y), np.floor(y) + 1, np.floor(y) + 1]), + np.concatenate( + [np.floor(x), np.floor(x) + 1, np.floor(x), np.floor(x) + 1] + ), + np.concatenate( + [np.floor(y), np.floor(y), np.floor(y) + 1, np.floor(y) + 1] + ), ] ) values_surr = ZarrReader._linear_interp_frac_coordinates( @@ -152,7 +168,17 @@ def test_zarr_geomax_on_grid(self): lons_ = np.array([3.92783]) lats_ = np.array([50.882394]) curve = np.array( - [0.00, 0.06997928, 0.2679602, 0.51508933, 0.69842442, 0.88040525, 1.11911115, 1.29562478, 1.47200677] + [ + 0.00, + 0.06997928, + 0.2679602, + 0.51508933, + 0.69842442, + 0.88040525, + 1.11911115, + 1.29562478, + 1.47200677, + ] ) set_id = r"inundation/wri/v2\\inunriver_rcp8p5_MIROC-ESM-CHEM_2080" interpolation = "linear" @@ -164,17 +190,44 @@ def test_zarr_geomax_on_grid(self): lons_ = np.array([3.92916667, 3.925] + list(lons_)) lats_ = np.array([50.87916667, 50.88333333] + list(lats_)) curves_max_candidate, _ = zarrreader_.get_max_curves_on_grid( - set_id, lons_, lats_, interpolation=interpolation, delta_km=delta_km, n_grid=n_grid + set_id, + lons_, + lats_, + interpolation=interpolation, + delta_km=delta_km, + n_grid=n_grid, ) curves_max_expected = np.array( [ curve, - [0.0, 0.02272942, 0.08703404, 0.16730212, 0.22684974, 0.28595751, 0.3634897, 0.42082168, 0.47811095], - [0.0, 0.0432026, 0.16542863, 0.31799695, 0.43118118, 0.54352937, 0.69089751, 0.7998704, 0.90876211], + [ + 0.0, + 0.02272942, + 0.08703404, + 0.16730212, + 0.22684974, + 0.28595751, + 0.3634897, + 0.42082168, + 0.47811095, + ], + [ + 0.0, + 0.0432026, + 0.16542863, + 0.31799695, + 0.43118118, + 0.54352937, + 0.69089751, + 0.7998704, + 0.90876211, + ], ] ) - numpy.testing.assert_allclose(curves_max_candidate, curves_max_expected, rtol=1e-6) + numpy.testing.assert_allclose( + curves_max_candidate, curves_max_expected, rtol=1e-6 + ) def test_zarr_geomax(self): longitudes = np.array([3.926]) @@ -199,11 +252,19 @@ def test_zarr_geomax(self): zarr_reader = ZarrReader(store) curves_max_expected = np.array([[0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4]]) - curves_max_candidate, _, _ = zarr_reader.get_max_curves(set_id, shapes, interpolation="floor") - numpy.testing.assert_allclose(curves_max_candidate, curves_max_expected, rtol=1e-6) + curves_max_candidate, _, _ = zarr_reader.get_max_curves( + set_id, shapes, interpolation="floor" + ) + numpy.testing.assert_allclose( + curves_max_candidate, curves_max_expected, rtol=1e-6 + ) - curves_max_candidate, _, _ = zarr_reader.get_max_curves(set_id, shapes, interpolation="linear") - numpy.testing.assert_allclose(curves_max_candidate, curves_max_expected / 4, rtol=1e-6) + curves_max_candidate, _, _ = zarr_reader.get_max_curves( + set_id, shapes, interpolation="linear" + ) + numpy.testing.assert_allclose( + curves_max_candidate, curves_max_expected / 4, rtol=1e-6 + ) def test_reproject(self): """Test adding data in a non-ESPG-4326 coordinate reference system. Check that attribute @@ -222,9 +283,22 @@ def test_reproject(self): [1.0, 2.0, 3.0], ) - source_paths = {RiverineInundation: lambda indicator_id, scenario, year, hint: "test"} + source_paths = { + RiverineInundation: lambda indicator_id, scenario, year, hint: "test" + } hazard_model = ZarrHazardModel(source_paths=source_paths, store=mocker.store) response = hazard_model.get_hazard_events( - [HazardDataRequest(RiverineInundation, lons[0], lats[0], indicator_id="", scenario="", year=2050)] + [ + HazardDataRequest( + RiverineInundation, + lons[0], + lats[0], + indicator_id="", + scenario="", + year=2050, + ) + ] + ) + numpy.testing.assert_equal( + next(iter(response.values())).intensities, [1.0, 2.0, 3.0] ) - numpy.testing.assert_equal(next(iter(response.values())).intensities, [1.0, 2.0, 3.0]) diff --git a/tests/data/hazard_model_store_test.py b/tests/data/hazard_model_store_test.py index ff0db588..293a84a9 100644 --- a/tests/data/hazard_model_store_test.py +++ b/tests/data/hazard_model_store_test.py @@ -49,11 +49,26 @@ def add_curves_global( units: str = "default", ): crs = "epsg:4326" - crs, shape, trans = self._crs_shape_transform_global(return_periods=return_periods, width=width, height=height) - self._add_curves(array_path, longitudes, latitudes, crs, shape, trans, return_periods, intensities, units=units) + crs, shape, trans = self._crs_shape_transform_global( + return_periods=return_periods, width=width, height=height + ) + self._add_curves( + array_path, + longitudes, + latitudes, + crs, + shape, + trans, + return_periods, + intensities, + units=units, + ) def _crs_shape_transform_global( - self, width: int = 43200, height: int = 21600, return_periods: Union[List[float], npt.NDArray] = [0.0] + self, + width: int = 43200, + height: int = 21600, + return_periods: Union[List[float], npt.NDArray] = [0.0], ): return self._crs_shape_transform(width, height, return_periods) @@ -70,7 +85,10 @@ def _add_curves( units: str = "default", ): z = self._root.create_dataset( # type: ignore - array_path, shape=(shape[0], shape[1], shape[2]), chunks=(shape[0], 1000, 1000), dtype="f4" + array_path, + shape=(shape[0], shape[1], shape[2]), + chunks=(shape[0], 1000, 1000), + dtype="f4", ) z.attrs["transform_mat3x3"] = trans z.attrs["index_values"] = return_periods @@ -96,13 +114,20 @@ def _add_curves( for j in range(len(x)): z[:, image_coords[1, j], image_coords[0, j]] = intensities - def _crs_shape_transform(self, width: int, height: int, return_periods: Union[List[float], npt.NDArray] = [0.0]): + def _crs_shape_transform( + self, + width: int, + height: int, + return_periods: Union[List[float], npt.NDArray] = [0.0], + ): t = [360.0 / width, 0.0, -180.0, 0.0, -180.0 / height, 90.0, 0.0, 0.0, 1.0] return "epsg:4326", (len(return_periods), height, width), t def shape_transform_21600_43200( - width: int = 43200, height: int = 21600, return_periods: Union[List[float], npt.NDArray] = [0.0] + width: int = 43200, + height: int = 21600, + return_periods: Union[List[float], npt.NDArray] = [0.0], ): t = [360.0 / width, 0.0, -180.0, 0.0, -180.0 / height, 90.0, 0.0, 0.0, 1.0] return (len(return_periods), height, width), t @@ -124,7 +149,10 @@ def add_curves( trans: List[float], ): z = root.create_dataset( # type: ignore - array_path, shape=(shape[0], shape[1], shape[2]), chunks=(shape[0], 1000, 1000), dtype="f4" + array_path, + shape=(shape[0], shape[1], shape[2]), + chunks=(shape[0], 1000, 1000), + dtype="f4", ) z.attrs["transform_mat3x3"] = trans z.attrs["index_values"] = return_periods @@ -146,13 +174,28 @@ def get_mock_hazard_model_store_single_curve(): is applied at all locations.""" return_periods = inundation_return_periods() - t = [0.008333333333333333, 0.0, -180.0, 0.0, -0.008333333333333333, 90.0, 0.0, 0.0, 1.0] + t = [ + 0.008333333333333333, + 0.0, + -180.0, + 0.0, + -0.008333333333333333, + 90.0, + 0.0, + 0.0, + 1.0, + ] shape = (len(return_periods), 21600, 43200) store = zarr.storage.MemoryStore(root="hazard.zarr") root = zarr.open(store=store, mode="w") - array_path = get_source_path_wri_riverine_inundation(model="MIROC-ESM-CHEM", scenario="rcp8p5", year=2080) + array_path = get_source_path_wri_riverine_inundation( + model="MIROC-ESM-CHEM", scenario="rcp8p5", year=2080 + ) z = root.create_dataset( # type: ignore - array_path, shape=(shape[0], shape[1], shape[2]), chunks=(shape[0], 1000, 1000), dtype="f4" + array_path, + shape=(shape[0], shape[1], shape[2]), + chunks=(shape[0], 1000, 1000), + dtype="f4", ) z.attrs["transform_mat3x3"] = t z.attrs["index_values"] = return_periods @@ -173,15 +216,21 @@ def get_mock_hazard_model_store_single_curve(): def mock_hazard_model_store_heat(longitudes, latitudes): - return mock_hazard_model_store_for_parameter_sets(longitudes, latitudes, degree_day_heat_parameter_set()) + return mock_hazard_model_store_for_parameter_sets( + longitudes, latitudes, degree_day_heat_parameter_set() + ) def mock_hazard_model_store_heat_wbgt(longitudes, latitudes): - return mock_hazard_model_store_for_parameter_sets(longitudes, latitudes, wbgt_gzn_joint_parameter_set()) + return mock_hazard_model_store_for_parameter_sets( + longitudes, latitudes, wbgt_gzn_joint_parameter_set() + ) def mock_hazard_model_store_inundation(longitudes, latitudes, curve): - return mock_hazard_model_store_single_curve_for_paths(longitudes, latitudes, curve, inundation_paths) + return mock_hazard_model_store_single_curve_for_paths( + longitudes, latitudes, curve, inundation_paths + ) def mock_hazard_model_store_for_parameter_sets(longitudes, latitudes, path_parameters): @@ -191,12 +240,24 @@ def mock_hazard_model_store_for_parameter_sets(longitudes, latitudes, path_param return_periods = None shape = (1, 21600, 43200) - t = [0.008333333333333333, 0.0, -180.0, 0.0, -0.008333333333333333, 90.0, 0.0, 0.0, 1.0] + t = [ + 0.008333333333333333, + 0.0, + -180.0, + 0.0, + -0.008333333333333333, + 90.0, + 0.0, + 0.0, + 1.0, + ] store = zarr.storage.MemoryStore(root="hazard.zarr") root = zarr.open(store=store, mode="w") for path, parameter in path_parameters.items(): - add_curves(root, longitudes, latitudes, path, shape, parameter, return_periods, t) + add_curves( + root, longitudes, latitudes, path, shape, parameter, return_periods, t + ) return store @@ -205,9 +266,15 @@ def mock_hazard_model_store_single_curve_for_paths(longitudes, latitudes, curve, """Create a MemoryStore for creation of Zarr hazard model to be used with unit tests, with the specified longitudes and latitudes set to the curve supplied.""" - return_periods = [0.0] if len(curve) == 1 else [2.0, 5.0, 10.0, 25.0, 50.0, 100.0, 250.0, 500.0, 1000.0] + return_periods = ( + [0.0] + if len(curve) == 1 + else [2.0, 5.0, 10.0, 25.0, 50.0, 100.0, 250.0, 500.0, 1000.0] + ) if len(curve) != len(return_periods): - raise ValueError(f"curve must be single value or of length {len(return_periods)}") + raise ValueError( + f"curve must be single value or of length {len(return_periods)}" + ) shape, t = shape_transform_21600_43200(return_periods=return_periods) store, root = zarr_memory_store() @@ -222,11 +289,23 @@ def inundation_return_periods(): return [2.0, 5.0, 10.0, 25.0, 50.0, 100.0, 250.0, 500.0, 1000.0] -def mock_hazard_model_store_path_curves(longitudes, latitudes, path_curves: Dict[str, np.ndarray]): +def mock_hazard_model_store_path_curves( + longitudes, latitudes, path_curves: Dict[str, np.ndarray] +): """Create a MemoryStore for creation of Zarr hazard model to be used with unit tests, with the specified longitudes and latitudes set to the curve supplied.""" - t = [0.008333333333333333, 0.0, -180.0, 0.0, -0.008333333333333333, 90.0, 0.0, 0.0, 1.0] + t = [ + 0.008333333333333333, + 0.0, + -180.0, + 0.0, + -0.008333333333333333, + 90.0, + 0.0, + 0.0, + 1.0, + ] store = zarr.storage.MemoryStore(root="hazard.zarr") root = zarr.open(store=store, mode="w") @@ -237,7 +316,9 @@ def mock_hazard_model_store_path_curves(longitudes, latitudes, path_curves: Dict else: return_periods = [2.0, 5.0, 10.0, 25.0, 50.0, 100.0, 250.0, 500.0, 1000.0] if len(curve) != len(return_periods): - raise ValueError(f"curve must be single value or of length {len(return_periods)}") + raise ValueError( + f"curve must be single value or of length {len(return_periods)}" + ) shape = (len(return_periods), 21600, 43200) add_curves(root, longitudes, latitudes, path, shape, curve, return_periods, t) @@ -254,7 +335,9 @@ def degree_day_heat_parameter_set(): ("mean_degree_days/above/32c/ACCESS-CM2", "ssp585", 2050), ("mean_degree_days/below/32c/ACCESS-CM2", "ssp585", 2050), ]: - paths.append(get_source_path_osc_chronic_heat(model=model, scenario=scenario, year=year)) + paths.append( + get_source_path_osc_chronic_heat(model=model, scenario=scenario, year=year) + ) parameters = [300, 300, 600, -200] return dict(zip(paths, parameters)) @@ -268,28 +351,43 @@ def wbgt_gzn_joint_parameter_set(): ("mean_degree_days/above/32c/ACCESS-CM2", "ssp585", 2050), ("mean_degree_days/below/32c/ACCESS-CM2", "ssp585", 2050), ]: - paths.append(get_source_path_osc_chronic_heat(model=model, scenario=scenario, year=year)) + paths.append( + get_source_path_osc_chronic_heat(model=model, scenario=scenario, year=year) + ) for model, scenario, year in [ ("mean_work_loss/high/ACCESS-CM2", "historical", 2005), # 2005 ("mean_work_loss/medium/ACCESS-CM2", "historical", 2005), ("mean_work_loss/high/ACCESS-CM2", "ssp585", 2050), ("mean_work_loss/medium/ACCESS-CM2", "ssp585", 2050), ]: - paths.append(get_source_path_osc_chronic_heat(model=model, scenario=scenario, year=year)) + paths.append( + get_source_path_osc_chronic_heat(model=model, scenario=scenario, year=year) + ) parameters = [300, 300, 600, -200, 0.05, 0.003, 0.11, 0.013] return dict(zip(paths, parameters)) def inundation_paths(): paths = [] - for model, scenario, year in [("MIROC-ESM-CHEM", "rcp8p5", 2080), ("000000000WATCH", "historical", 1980)]: - paths.append(get_source_path_wri_riverine_inundation(model=model, scenario=scenario, year=year)) + for model, scenario, year in [ + ("MIROC-ESM-CHEM", "rcp8p5", 2080), + ("000000000WATCH", "historical", 1980), + ]: + paths.append( + get_source_path_wri_riverine_inundation( + model=model, scenario=scenario, year=year + ) + ) for model, scenario, year in [ ("wtsub/95", "rcp8p5", "2080"), ("wtsub", "historical", "hist"), ("nosub", "historical", "hist"), ]: - paths.append(get_source_path_wri_coastal_inundation(model=model, scenario=scenario, year=year)) + paths.append( + get_source_path_wri_coastal_inundation( + model=model, scenario=scenario, year=year + ) + ) return paths @@ -308,16 +406,22 @@ def get_source_path_wri_coastal_inundation(*, model: str, scenario: str, year: s model_components = model.split("/") sub = model_components[0] if sub not in _subsidence_set: - raise ValueError("expected model input of the form {subsidence/percentile}, e.g. wtsub/95, nosub/5, wtsub/50") + raise ValueError( + "expected model input of the form {subsidence/percentile}, e.g. wtsub/95, nosub/5, wtsub/50" + ) perc = "95" if len(model_components) == 1 else model_components[1] return os.path.join( - _wri_inundation_prefix(), f"inun{type}_{cmip6_scenario_to_rcp(scenario)}_{sub}_{year}_{_percentiles_map[perc]}" + _wri_inundation_prefix(), + f"inun{type}_{cmip6_scenario_to_rcp(scenario)}_{sub}_{year}_{_percentiles_map[perc]}", ) def get_source_path_wri_riverine_inundation(*, model: str, scenario: str, year: int): type = "river" - return os.path.join(_wri_inundation_prefix(), f"inun{type}_{cmip6_scenario_to_rcp(scenario)}_{model}_{year}") + return os.path.join( + _wri_inundation_prefix(), + f"inun{type}_{cmip6_scenario_to_rcp(scenario)}_{model}_{year}", + ) def _osc_chronic_heat_prefix(): @@ -331,12 +435,20 @@ def get_source_path_osc_chronic_heat(*, model: str, scenario: str, year: int): assert levels[0] in ["above", "below"] # above or below assert levels[1] in ["18c", "32c"] # threshold temperature assert levels[2] in ["ACCESS-CM2"] # gcms - return _osc_chronic_heat_prefix() + "/" + f"{type}_v2_{levels[0]}_{levels[1]}_{levels[2]}_{scenario}_{year}" + return ( + _osc_chronic_heat_prefix() + + "/" + + f"{type}_v2_{levels[0]}_{levels[1]}_{levels[2]}_{scenario}_{year}" + ) elif type == "mean_work_loss": assert levels[0] in ["low", "medium", "high"] # work intensity assert levels[1] in ["ACCESS-CM2"] # gcms - return _osc_chronic_heat_prefix() + "/" + f"{type}_{levels[0]}_{levels[1]}_{scenario}_{year}" + return ( + _osc_chronic_heat_prefix() + + "/" + + f"{type}_{levels[0]}_{levels[1]}_{scenario}_{year}" + ) else: raise ValueError("valid types are {valid_types}") diff --git a/tests/data/static_data_test.py b/tests/data/static_data_test.py index af4a4e9f..cc9f57c7 100644 --- a/tests/data/static_data_test.py +++ b/tests/data/static_data_test.py @@ -1,14 +1,22 @@ import unittest -from physrisk.data.static.world import World, get_countries_and_continents, get_countries_json +from physrisk.data.static.world import ( + World, + get_countries_and_continents, + get_countries_json, +) from ..data.hazard_model_store_test import TestData class TestStaticDate(unittest.TestCase): - @unittest.skip("example that requires geopandas (consider adding for running tests only)") + @unittest.skip( + "example that requires geopandas (consider adding for running tests only)" + ) def test_get_countries_and_continents(self): - countries, continents = get_countries_and_continents(TestData.longitudes, TestData.latitudes) + countries, continents = get_countries_and_continents( + TestData.longitudes, TestData.latitudes + ) self.assertEqual(countries[0:3], ["Afghanistan", "Afghanistan", "Albania"]) @unittest.skip("not really a test; just showing how world.json was generated") diff --git a/tests/kernel/asset_impact_test.py b/tests/kernel/asset_impact_test.py index 6fb87441..5726f12a 100644 --- a/tests/kernel/asset_impact_test.py +++ b/tests/kernel/asset_impact_test.py @@ -1,4 +1,4 @@ -""" Test asset impact calculations.""" +"""Test asset impact calculations.""" import unittest @@ -24,10 +24,22 @@ def test_impact_curve(self): """Testing the generation of an asset when only an impact curve (e.g. damage curve is available)""" # exceedance curve - return_periods = np.array([2.0, 5.0, 10.0, 25.0, 50.0, 100.0, 250.0, 500.0, 1000.0]) + return_periods = np.array( + [2.0, 5.0, 10.0, 25.0, 50.0, 100.0, 250.0, 500.0, 1000.0] + ) exceed_probs = 1.0 / return_periods depths = np.array( - [0.059601218, 0.33267087, 0.50511575, 0.71471703, 0.8641244, 1.0032823, 1.1491022, 1.1634114, 1.1634114] + [ + 0.059601218, + 0.33267087, + 0.50511575, + 0.71471703, + 0.8641244, + 1.0032823, + 1.1491022, + 1.1634114, + 1.1634114, + ] ) curve = ExceedanceCurve(exceed_probs, depths) @@ -52,9 +64,21 @@ def test_impact_curve(self): self.assertAlmostEqual(mean, 4.8453897) def test_protection_level(self): - return_periods = np.array([2.0, 5.0, 10.0, 25.0, 50.0, 100.0, 250.0, 500.0, 1000.0]) + return_periods = np.array( + [2.0, 5.0, 10.0, 25.0, 50.0, 100.0, 250.0, 500.0, 1000.0] + ) base_depth = np.array( - [0.0, 0.22372675, 0.3654859, 0.5393629, 0.6642473, 0.78564394, 0.9406518, 1.0539534, 1.1634114] + [ + 0.0, + 0.22372675, + 0.3654859, + 0.5393629, + 0.6642473, + 0.78564394, + 0.9406518, + 1.0539534, + 1.1634114, + ] ) # future_depth = np.array( # [0.059601218, 0.33267087, 0.50511575, 0.71471703, 0.8641244, 1.0032823, 1.1491022, 1.1634114, 1.1634114] @@ -63,16 +87,30 @@ def test_protection_level(self): exceed_probs = 1.0 / return_periods protection_return_period = 250.0 # protection level of 250 years - protection_depth = np.interp(1.0 / protection_return_period, exceed_probs[::-1], base_depth[::-1]) + protection_depth = np.interp( + 1.0 / protection_return_period, exceed_probs[::-1], base_depth[::-1] + ) self.assertAlmostEqual(protection_depth, 0.9406518) # type: ignore def test_single_asset_impact(self): # exceedance curve - return_periods = np.array([2.0, 5.0, 10.0, 25.0, 50.0, 100.0, 250.0, 500.0, 1000.0]) + return_periods = np.array( + [2.0, 5.0, 10.0, 25.0, 50.0, 100.0, 250.0, 500.0, 1000.0] + ) exceed_probs = 1.0 / return_periods depths = np.array( - [0.059601218, 0.33267087, 0.50511575, 0.71471703, 0.8641244, 1.0032823, 1.1491022, 1.1634114, 1.1634114] + [ + 0.059601218, + 0.33267087, + 0.50511575, + 0.71471703, + 0.8641244, + 1.0032823, + 1.1491022, + 1.1634114, + 1.1634114, + ] ) curve = ExceedanceCurve(exceed_probs, depths) @@ -94,7 +132,9 @@ def test_single_asset_impact(self): type(RiverineInundation), depth_bins, impact_bins, np.diag(probs_w_cutoff) ) # np.eye(n_bins, n_bins)) hazard_paths = ["unknown"] - event = HazardEventDistrib(type(RiverineInundation), depth_bins, probs, hazard_paths) # type: ignore + event = HazardEventDistrib( + type(RiverineInundation), depth_bins, probs, hazard_paths + ) # type: ignore impact_prob = vul.prob_matrix.T @ event.prob impact = ImpactDistrib(vul.event_type, vul.impact_bins, impact_prob, event.path) @@ -110,9 +150,15 @@ def test_performance_hazardlookup(self): start = time.time() - assets = [RealEstateAsset(latitude=0, longitude=0, location="", type="") for _ in range(10000)] + assets = [ + RealEstateAsset(latitude=0, longitude=0, location="", type="") + for _ in range(10000) + ] - vulnerability_models = [RealEstateCoastalInundationModel(), RealEstateRiverineInundationModel()] + vulnerability_models = [ + RealEstateCoastalInundationModel(), + RealEstateRiverineInundationModel(), + ] time_assets = time.time() - start print(f"Time for asset generation {time_assets}s ") @@ -124,7 +170,14 @@ def test_performance_hazardlookup(self): for v in vulnerability_models: for a in assets: asset_requests[(v, a)] = [ - HazardDataRequest(RiverineInundation, 0, 0, indicator_id="", scenario="", year=2030) + HazardDataRequest( + RiverineInundation, + 0, + 0, + indicator_id="", + scenario="", + year=2030, + ) ] time_requests = time.time() - start diff --git a/tests/kernel/chronic_asset_impact_test.py b/tests/kernel/chronic_asset_impact_test.py index 018a8f60..1f5e74b6 100644 --- a/tests/kernel/chronic_asset_impact_test.py +++ b/tests/kernel/chronic_asset_impact_test.py @@ -1,5 +1,5 @@ import unittest -from typing import Iterable, List, Union, cast +from typing import Iterable, List, Union import numpy as np from scipy.stats import norm @@ -7,11 +7,18 @@ from physrisk.data.pregenerated_hazard_model import ZarrHazardModel from physrisk.hazard_models.core_hazards import get_default_source_paths from physrisk.kernel.assets import Asset, IndustrialActivity -from physrisk.kernel.hazard_model import HazardDataRequest, HazardDataResponse, HazardParameterDataResponse +from physrisk.kernel.hazard_model import ( + HazardDataRequest, + HazardDataResponse, + HazardParameterDataResponse, +) from physrisk.kernel.hazards import ChronicHeat from physrisk.kernel.impact import calculate_impacts from physrisk.kernel.impact_distrib import ImpactDistrib, ImpactType -from physrisk.kernel.vulnerability_model import DictBasedVulnerabilityModels, VulnerabilityModelBase +from physrisk.kernel.vulnerability_model import ( + DictBasedVulnerabilityModels, + VulnerabilityModelBase, +) from physrisk.vulnerability_models.chronic_heat_models import ChronicHeatGZNModel from ..data.hazard_model_store_test import TestData, mock_hazard_model_store_heat @@ -24,14 +31,24 @@ class ExampleChronicHeatModel(VulnerabilityModelBase): https://www.sphinx-doc.org/en/master/usage/extensions/example_google.html """ - def __init__(self, indicator_id: str = "mean_degree_days_above_32c", delta: bool = True): + def __init__( + self, indicator_id: str = "mean_degree_days_above_32c", delta: bool = True + ): super().__init__( - indicator_id=indicator_id, hazard_type=ChronicHeat, impact_type=ImpactType.disruption + indicator_id=indicator_id, + hazard_type=ChronicHeat, + impact_type=ImpactType.disruption, ) # opportunity to give a model hint, but blank here - self.time_lost_per_degree_day = 4.671 # This comes from the paper converted to celsius - self.time_lost_per_degree_day_se = 2.2302 # This comes from the paper converted to celsius - self.total_labour_hours = 107460 # OECD Average hours worked within the USA in one year. + self.time_lost_per_degree_day = ( + 4.671 # This comes from the paper converted to celsius + ) + self.time_lost_per_degree_day_se = ( + 2.2302 # This comes from the paper converted to celsius + ) + self.total_labour_hours = ( + 107460 # OECD Average hours worked within the USA in one year. + ) self.delta = delta def get_data_requests( @@ -71,7 +88,9 @@ def get_data_requests( ), ] - def get_impact(self, asset: Asset, data_responses: List[HazardDataResponse]) -> ImpactDistrib: + def get_impact( + self, asset: Asset, data_responses: List[HazardDataResponse] + ) -> ImpactDistrib: """Calcaulate impact (disruption) of asset based on the hazard data returned. Args: @@ -94,13 +113,24 @@ def get_impact(self, asset: Asset, data_responses: List[HazardDataResponse]) -> # use hazard data requests via: # Allow for either a delta approach or a level estimate. - delta_dd_above_mean: float = scenario_dd_above_mean.parameter - baseline_dd_above_mean.parameter * self.delta + delta_dd_above_mean: float = ( + scenario_dd_above_mean.parameter + - baseline_dd_above_mean.parameter * self.delta + ) hours_worked = self.total_labour_hours - fraction_loss_mean = (delta_dd_above_mean * self.time_lost_per_degree_day) / hours_worked - fraction_loss_std = (delta_dd_above_mean * self.time_lost_per_degree_day_se) / hours_worked + fraction_loss_mean = ( + delta_dd_above_mean * self.time_lost_per_degree_day + ) / hours_worked + fraction_loss_std = ( + delta_dd_above_mean * self.time_lost_per_degree_day_se + ) / hours_worked return get_impact_distrib( - fraction_loss_mean, fraction_loss_std, ChronicHeat, [scenario_dd_above_mean.path], ImpactType.disruption + fraction_loss_mean, + fraction_loss_std, + ChronicHeat, + [scenario_dd_above_mean.path], + ImpactType.disruption, ) @@ -133,9 +163,11 @@ def get_impact_distrib( ] ) - probs_cumulative = np.vectorize(lambda x: norm.cdf(x, loc=fraction_loss_mean, scale=max(1e-12, fraction_loss_std)))( - impact_bins - ) + probs_cumulative = np.vectorize( + lambda x: norm.cdf( + x, loc=fraction_loss_mean, scale=max(1e-12, fraction_loss_std) + ) + )(impact_bins) probs_cumulative[-1] = np.maximum(probs_cumulative[-1], 1.0) probs = np.diff(probs_cumulative) @@ -159,20 +191,26 @@ def test_chronic_vulnerability_model(self): """Testing the generation of an asset when only an impact curve (e.g. damage curve is available)""" store = mock_hazard_model_store_heat(TestData.longitudes, TestData.latitudes) - hazard_model = ZarrHazardModel(source_paths=get_default_source_paths(), store=store) + hazard_model = ZarrHazardModel( + source_paths=get_default_source_paths(), store=store + ) # to run a live calculation, we omit the store parameter scenario = "ssp585" year = 2050 - vulnerability_models = DictBasedVulnerabilityModels({IndustrialActivity: [ChronicHeatGZNModel()]}) + vulnerability_models = DictBasedVulnerabilityModels( + {IndustrialActivity: [ChronicHeatGZNModel()]} + ) assets = [ IndustrialActivity(lat, lon, type="Construction") for lon, lat in zip(TestData.longitudes, TestData.latitudes) ][:1] - results = calculate_impacts(assets, hazard_model, vulnerability_models, scenario=scenario, year=year) + results = calculate_impacts( + assets, hazard_model, vulnerability_models, scenario=scenario, year=year + ) value_test = list(results.values())[0][0].impact.mean_impact() value_test = list(results.values())[0][0].impact.prob diff --git a/tests/kernel/curves_test.py b/tests/kernel/curves_test.py index f8c22a19..bf488772 100644 --- a/tests/kernel/curves_test.py +++ b/tests/kernel/curves_test.py @@ -1,4 +1,4 @@ -""" Test asset impact calculations.""" +"""Test asset impact calculations.""" import unittest @@ -11,7 +11,9 @@ class TestAssetImpact(unittest.TestCase): """Tests asset impact calculations.""" def test_return_period_data(self): - return_periods = np.array([2.0, 5.0, 10.0, 25.0, 50.0, 100.0, 250.0, 500.0, 1000.0]) + return_periods = np.array( + [2.0, 5.0, 10.0, 25.0, 50.0, 100.0, 250.0, 500.0, 1000.0] + ) depths = np.array([0.059, 0.33, 0.51, 0.71, 0.86, 1.00, 1.15, 1.16, 1.16]) # say we need to add an extra depth point because the damage below that point is zero diff --git a/tests/kernel/exposure_test.py b/tests/kernel/exposure_test.py index 8c726a5d..72dad183 100644 --- a/tests/kernel/exposure_test.py +++ b/tests/kernel/exposure_test.py @@ -4,7 +4,10 @@ import numpy as np import physrisk.api.v1.common -from physrisk.api.v1.exposure_req_resp import AssetExposureRequest, AssetExposureResponse +from physrisk.api.v1.exposure_req_resp import ( + AssetExposureRequest, + AssetExposureResponse, +) from physrisk.container import ZarrHazardModelFactory from physrisk.data.inventory import EmbeddedInventory from physrisk.data.inventory_reader import InventoryReader @@ -12,8 +15,19 @@ from physrisk.hazard_models.core_hazards import get_default_source_paths from physrisk.kernel.assets import Asset from physrisk.kernel.calculation import DefaultMeasuresFactory -from physrisk.kernel.exposure import Category, JupterExposureMeasure, calculate_exposures -from physrisk.kernel.hazards import ChronicHeat, CombinedInundation, Drought, Fire, Hail, Wind +from physrisk.kernel.exposure import ( + Category, + JupterExposureMeasure, + calculate_exposures, +) +from physrisk.kernel.hazards import ( + ChronicHeat, + CombinedInundation, + Drought, + Fire, + Hail, + Wind, +) from physrisk.requests import Requester from ..base_test import TestWithCredentials @@ -35,12 +49,16 @@ def test_jupiter_exposure_service(self): ) assets_api = physrisk.api.v1.common.Assets( items=[ - physrisk.api.v1.common.Asset(asset_class="Asset", latitude=a.latitude, longitude=a.longitude) + physrisk.api.v1.common.Asset( + asset_class="Asset", latitude=a.latitude, longitude=a.longitude + ) for a in assets[0:1] ] ) request = AssetExposureRequest(assets=assets_api, scenario="ssp585", year=2050) - response = requester.get(request_id="get_asset_exposure", request_dict=request.model_dump()) + response = requester.get( + request_id="get_asset_exposure", request_dict=request.model_dump() + ) result = AssetExposureResponse(**json.loads(response)).items[0] expected = dict((k.__name__, v) for (k, v) in expected.items()) for key in result.exposures.keys(): @@ -52,7 +70,11 @@ def test_jupiter_exposure(self): measure = JupterExposureMeasure() results = calculate_exposures( - [asset], hazard_model_factory.hazard_model(), measure, scenario="ssp585", year=2030 + [asset], + hazard_model_factory.hazard_model(), + measure, + scenario="ssp585", + year=2030, ) categories = results[asset].hazard_categories for k, v in expected.items(): @@ -84,11 +106,19 @@ def _get_components(self): } def path_curves(): - return dict((r.path.format(scenario="ssp585", year=2030), v) for (r, v) in zip(resources, values)) - - assets = [Asset(lat, lon) for (lat, lon) in zip(TestData.latitudes, TestData.longitudes)] + return dict( + (r.path.format(scenario="ssp585", year=2030), v) + for (r, v) in zip(resources, values) + ) + + assets = [ + Asset(lat, lon) + for (lat, lon) in zip(TestData.latitudes, TestData.longitudes) + ] - store = mock_hazard_model_store_path_curves(TestData.longitudes, TestData.latitudes, path_curves()) + store = mock_hazard_model_store_path_curves( + TestData.longitudes, TestData.latitudes, path_curves() + ) hazard_model_factory = ZarrHazardModelFactory( source_paths=get_default_source_paths(EmbeddedInventory()), store=store diff --git a/tests/kernel/financial_model_test.py b/tests/kernel/financial_model_test.py index 9914eb24..e06a3a49 100644 --- a/tests/kernel/financial_model_test.py +++ b/tests/kernel/financial_model_test.py @@ -16,7 +16,9 @@ class MockFinancialDataProvider(FinancialDataProvider): def get_asset_value(self, asset: Asset, currency: str) -> float: return 1000 - def get_asset_aggregate_cashflows(self, asset: Asset, start: datetime, end: datetime, currency: str) -> float: + def get_asset_aggregate_cashflows( + self, asset: Asset, start: datetime, end: datetime, currency: str + ) -> float: return 1000 @@ -25,28 +27,51 @@ class TestAssetImpact(unittest.TestCase): def test_financial_model(self): curve = np.array( - [0.059601218, 0.33267087, 0.50511575, 0.71471703, 0.8641244, 1.0032823, 1.1491022, 1.1634114, 1.1634114] + [ + 0.059601218, + 0.33267087, + 0.50511575, + 0.71471703, + 0.8641244, + 1.0032823, + 1.1491022, + 1.1634114, + 1.1634114, + ] + ) + store = mock_hazard_model_store_inundation( + TestData.longitudes, TestData.latitudes, curve ) - store = mock_hazard_model_store_inundation(TestData.longitudes, TestData.latitudes, curve) # we need to define # 1) The hazard models # 2) The vulnerability models # 3) The financial models - hazard_model = ZarrHazardModel(source_paths=get_default_source_paths(), store=store) + hazard_model = ZarrHazardModel( + source_paths=get_default_source_paths(), store=store + ) model = LossModel(hazard_model=hazard_model) data_provider = MockFinancialDataProvider() financial_model = FinancialModel(data_provider) - assets = [PowerGeneratingAsset(lat, lon) for lon, lat in zip(TestData.longitudes, TestData.latitudes)] + assets = [ + PowerGeneratingAsset(lat, lon) + for lon, lat in zip(TestData.longitudes, TestData.latitudes) + ] measures = model.get_financial_impacts( - assets, financial_model=financial_model, scenario="ssp585", year=2080, sims=100000 + assets, + financial_model=financial_model, + scenario="ssp585", + year=2080, + sims=100000, ) np.testing.assert_array_almost_equal_nulp( measures["RiverineInundation"]["percentile_values"], - np.array([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1000.0, 1000.0, 1000.0, 2000.0]), + np.array( + [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1000.0, 1000.0, 1000.0, 2000.0] + ), ) diff --git a/tests/kernel/hazard_models_test.py b/tests/kernel/hazard_models_test.py index 4e2802ab..7950b9a1 100644 --- a/tests/kernel/hazard_models_test.py +++ b/tests/kernel/hazard_models_test.py @@ -46,22 +46,40 @@ def __init__(self, points: Sequence[SinglePointData]): 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 + 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): - return PointsKey(latitude=round(latitude, 3), longitude=round(longitude, 3), scenario=scenario, year=year) + 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]: + def get_hazard_events( + self, requests: List[HazardDataRequest] + ) -> Mapping[HazardDataRequest, HazardDataResponse]: response: Dict[HazardDataRequest, HazardDataResponse] = {} for request in requests: - point = self.points[self._get_key(request.latitude, request.longitude, request.scenario, request.year)] + 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 + 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( + np.array(point.chronic_heat_intensity) ) - elif request.hazard_type == ChronicHeat and request.indicator_id == "days/above/35c": - response[request] = HazardParameterDataResponse(np.array(point.chronic_heat_intensity)) # etc return response @@ -89,8 +107,12 @@ def test_using_point_based_hazard_model(): ) hazard_model = PointBasedHazardModel([point]) - vulnerability_models = DictBasedVulnerabilityModels({RealEstateAsset: [GenericTropicalCycloneModel()]}) - results = calculate_impacts(assets, hazard_model, vulnerability_models, scenario=scenario, year=year) + vulnerability_models = DictBasedVulnerabilityModels( + {RealEstateAsset: [GenericTropicalCycloneModel()]} + ) + results = calculate_impacts( + assets, hazard_model, vulnerability_models, scenario=scenario, year=year + ) impact_distrib = results[(assets[0], Wind, scenario, year)][0].impact mean_impact = impact_distrib.mean_impact() np.testing.assert_almost_equal(mean_impact, 0.009909858317497338) diff --git a/tests/kernel/image_creation_test.py b/tests/kernel/image_creation_test.py index 1e2638e7..975bea9c 100644 --- a/tests/kernel/image_creation_test.py +++ b/tests/kernel/image_creation_test.py @@ -20,7 +20,10 @@ def test_image_creation(self): im = np.array([[1.2, 0.8], [0.5, 0.4]]) z = root.create_dataset( # type: ignore - path, shape=(1, im.shape[0], im.shape[1]), chunks=(1, im.shape[0], im.shape[1]), dtype="f4" + path, + shape=(1, im.shape[0], im.shape[1]), + chunks=(1, im.shape[0], im.shape[1]), + dtype="f4", ) z[0, :, :] = im converter = ImageCreator(reader=ZarrReader(store)) @@ -31,8 +34,15 @@ def get_colors(index: int): result = converter._to_rgba(im, get_colors) # Max should be 255, min should be 1. Other values span the 253 elements from 2 to 254. - expected = np.array([[255, 2 + (0.8 - 0.4) * 253 / (1.2 - 0.4)], [2 + (0.5 - 0.4) * 253 / (1.2 - 0.4), 1]]) - converter.convert(path, colormap="test") # check no error running through mocked example. + expected = np.array( + [ + [255, 2 + (0.8 - 0.4) * 253 / (1.2 - 0.4)], + [2 + (0.5 - 0.4) * 253 / (1.2 - 0.4), 1], + ] + ) + converter.convert( + path, colormap="test" + ) # check no error running through mocked example. np.testing.assert_equal(result, expected.astype(np.uint8)) @unittest.skip("just example") @@ -41,6 +51,8 @@ def test_write_file(self): # useful for testing image generation test_output_dir = "{set me}" test_path = "wildfire/jupiter/v1/wildfire_probability_ssp585_2050_map" - store = zarr.DirectoryStore(os.path.join(test_output_dir, "hazard_test", "hazard.zarr")) + store = zarr.DirectoryStore( + os.path.join(test_output_dir, "hazard_test", "hazard.zarr") + ) creator = ImageCreator(ZarrReader(store)) creator.to_file(os.path.join(test_output_dir, "test.png"), test_path) diff --git a/tests/kernel/live_services_test.py b/tests/kernel/live_services_test.py index 4750d0d6..9c2e7900 100644 --- a/tests/kernel/live_services_test.py +++ b/tests/kernel/live_services_test.py @@ -10,7 +10,13 @@ def test_live_exposure(): request = { "assets": { "items": [ - {"asset_class": "Asset", "type": None, "location": None, "latitude": 34.556, "longitude": 69.4787} + { + "asset_class": "Asset", + "type": None, + "location": None, + "latitude": 34.556, + "longitude": 69.4787, + } ] }, "calc_settings": {"hazard_interp": "floor"}, diff --git a/tests/models/example_models_test.py b/tests/models/example_models_test.py index 7d00c4c1..3026865c 100644 --- a/tests/models/example_models_test.py +++ b/tests/models/example_models_test.py @@ -11,17 +11,31 @@ from physrisk.kernel.impact import calculate_impacts from physrisk.kernel.impact_distrib import ImpactType from physrisk.kernel.vulnerability_matrix_provider import VulnMatrixProvider -from physrisk.kernel.vulnerability_model import DictBasedVulnerabilityModels, VulnerabilityModel -from physrisk.vulnerability_models.example_models import ExampleCdfBasedVulnerabilityModel -from tests.data.hazard_model_store_test import TestData, mock_hazard_model_store_inundation +from physrisk.kernel.vulnerability_model import ( + DictBasedVulnerabilityModels, + VulnerabilityModel, +) +from physrisk.vulnerability_models.example_models import ( + ExampleCdfBasedVulnerabilityModel, +) +from tests.data.hazard_model_store_test import ( + TestData, + mock_hazard_model_store_inundation, +) class ExampleRealEstateInundationModel(VulnerabilityModel): def __init__(self): self.intensities = np.array([0, 0.01, 0.5, 1.0, 1.5, 2, 3, 4, 5, 6]) - self.impact_means = np.array([0, 0.2, 0.44, 0.58, 0.68, 0.78, 0.85, 0.92, 0.96, 1.0]) - self.impact_stddevs = np.array([0, 0.17, 0.14, 0.14, 0.17, 0.14, 0.13, 0.10, 0.06, 0]) - impact_bin_edges = np.array([0, 0.01, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]) + self.impact_means = np.array( + [0, 0.2, 0.44, 0.58, 0.68, 0.78, 0.85, 0.92, 0.96, 1.0] + ) + self.impact_stddevs = np.array( + [0, 0.17, 0.14, 0.14, 0.17, 0.14, 0.13, 0.10, 0.06, 0] + ) + impact_bin_edges = np.array( + [0, 0.01, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0] + ) super().__init__( indicator_id="flood_depth", hazard_type=RiverineInundation, @@ -34,7 +48,10 @@ def get_impact_curve(self, intensities, asset): impact_means = np.interp(intensities, self.intensities, self.impact_means) impact_stddevs = np.interp(intensities, self.intensities, self.impact_stddevs) return VulnMatrixProvider( - intensities, impact_cdfs=[checked_beta_distrib(m, s) for m, s in zip(impact_means, impact_stddevs)] + intensities, + impact_cdfs=[ + checked_beta_distrib(m, s) for m, s in zip(impact_means, impact_stddevs) + ], ) @@ -60,14 +77,28 @@ def beta_distrib(mean, std): class TestExampleModels(unittest.TestCase): def test_pdf_based_vulnerability_model(self): - model = ExampleCdfBasedVulnerabilityModel(indicator_id="", hazard_type=Inundation) + model = ExampleCdfBasedVulnerabilityModel( + indicator_id="", hazard_type=Inundation + ) latitude, longitude = 45.268405, 19.885738 asset = Asset(latitude, longitude) - return_periods = np.array([2.0, 5.0, 10.0, 25.0, 50.0, 100.0, 250.0, 500.0, 1000.0]) + return_periods = np.array( + [2.0, 5.0, 10.0, 25.0, 50.0, 100.0, 250.0, 500.0, 1000.0] + ) intensities = np.array( - [0.059601218, 0.33267087, 0.50511575, 0.71471703, 0.8641244, 1.0032823, 1.1491022, 1.1634114, 1.1634114] + [ + 0.059601218, + 0.33267087, + 0.50511575, + 0.71471703, + 0.8641244, + 1.0032823, + 1.1491022, + 1.1634114, + 1.1634114, + ] ) mock_response = HazardEventDataResponse(return_periods, intensities) @@ -76,23 +107,44 @@ def test_pdf_based_vulnerability_model(self): def test_user_supplied_model(self): curve = np.array( - [0.059601218, 0.33267087, 0.50511575, 0.71471703, 0.8641244, 1.0032823, 1.1491022, 1.1634114, 1.1634114] + [ + 0.059601218, + 0.33267087, + 0.50511575, + 0.71471703, + 0.8641244, + 1.0032823, + 1.1491022, + 1.1634114, + 1.1634114, + ] + ) + store = mock_hazard_model_store_inundation( + TestData.longitudes, TestData.latitudes, curve + ) + hazard_model = ZarrHazardModel( + source_paths=get_default_source_paths(), store=store ) - store = mock_hazard_model_store_inundation(TestData.longitudes, TestData.latitudes, curve) - hazard_model = ZarrHazardModel(source_paths=get_default_source_paths(), store=store) scenario = "rcp8p5" year = 2080 - vulnerability_models = DictBasedVulnerabilityModels({RealEstateAsset: [ExampleRealEstateInundationModel()]}) + vulnerability_models = DictBasedVulnerabilityModels( + {RealEstateAsset: [ExampleRealEstateInundationModel()]} + ) assets = [ RealEstateAsset(lat, lon, location="Asia", type="Building/Industrial") for lon, lat in zip(TestData.longitudes, TestData.latitudes) ] - results = calculate_impacts(assets, hazard_model, vulnerability_models, scenario=scenario, year=year) + results = calculate_impacts( + assets, hazard_model, vulnerability_models, scenario=scenario, year=year + ) self.assertAlmostEqual( - results[assets[0], RiverineInundation, scenario, year][0].impact.to_exceedance_curve().probs[0], 0.499 + results[assets[0], RiverineInundation, scenario, year][0] + .impact.to_exceedance_curve() + .probs[0], + 0.499, ) diff --git a/tests/models/power_generating_asset_models_test.py b/tests/models/power_generating_asset_models_test.py index 56996daf..5478686c 100644 --- a/tests/models/power_generating_asset_models_test.py +++ b/tests/models/power_generating_asset_models_test.py @@ -1,4 +1,4 @@ -""" Test asset impact calculations.""" +"""Test asset impact calculations.""" import os import unittest @@ -9,7 +9,11 @@ import physrisk.api.v1.common import physrisk.data.static.world as wd from physrisk.kernel import Asset, PowerGeneratingAsset, calculation -from physrisk.kernel.assets import IndustrialActivity, RealEstateAsset, ThermalPowerGeneratingAsset +from physrisk.kernel.assets import ( + IndustrialActivity, + RealEstateAsset, + ThermalPowerGeneratingAsset, +) from physrisk.kernel.hazard_model import HazardEventDataResponse from physrisk.kernel.impact import calculate_impacts from physrisk.kernel.impact_distrib import EmptyImpactDistrib @@ -26,12 +30,34 @@ class TestPowerGeneratingAssetModels(TestWithCredentials): def test_inundation(self): # exceedance curve - return_periods = np.array([2.0, 5.0, 10.0, 25.0, 50.0, 100.0, 250.0, 500.0, 1000.0]) + return_periods = np.array( + [2.0, 5.0, 10.0, 25.0, 50.0, 100.0, 250.0, 500.0, 1000.0] + ) base_depth = np.array( - [0.0, 0.22372675, 0.3654859, 0.5393629, 0.6642473, 0.78564394, 0.9406518, 1.0539534, 1.1634114] + [ + 0.0, + 0.22372675, + 0.3654859, + 0.5393629, + 0.6642473, + 0.78564394, + 0.9406518, + 1.0539534, + 1.1634114, + ] ) future_depth = np.array( - [0.059601218, 0.33267087, 0.50511575, 0.71471703, 0.8641244, 1.0032823, 1.1491022, 1.1634114, 1.1634114] + [ + 0.059601218, + 0.33267087, + 0.50511575, + 0.71471703, + 0.8641244, + 1.0032823, + 1.1491022, + 1.1634114, + 1.1634114, + ] ) # we mock the response of the data request @@ -66,12 +92,18 @@ def test_create_synthetic_portfolios_and_test(self): primary_fuel = np.array(filtered["primary_fuel"]) generation = np.array(filtered["estimated_generation_gwh"]) - _, continents = wd.get_countries_and_continents(latitudes=latitudes, longitudes=longitudes) + _, continents = wd.get_countries_and_continents( + latitudes=latitudes, longitudes=longitudes + ) # Power generating assets that are of interest assets = [ - PowerGeneratingAsset(lat, lon, generation=gen, location=continent, type=prim_fuel) - for lon, lat, gen, prim_fuel, continent in zip(longitudes, latitudes, generation, primary_fuel, continents) + PowerGeneratingAsset( + lat, lon, generation=gen, location=continent, type=prim_fuel + ) + for lon, lat, gen, prim_fuel, continent in zip( + longitudes, latitudes, generation, primary_fuel, continents + ) ] detailed_results = calculate_impacts(assets, scenario="ssp585", year=2030) keys = list(detailed_results.keys()) @@ -79,7 +111,10 @@ def test_create_synthetic_portfolios_and_test(self): means = np.array([detailed_results[key].impact.mean_impact() for key in keys]) interesting = [k for (k, m) in zip(keys, means) if m > 0] assets_out = self.api_assets(item[0] for item in interesting[0:10]) - with open(os.path.join(cache_folder, "assets_example_power_generating_small.json"), "w") as f: + with open( + os.path.join(cache_folder, "assets_example_power_generating_small.json"), + "w", + ) as f: f.write(assets_out.model_dump_json(indent=4)) # Synthetic portfolio; industrial activity at different locations @@ -90,10 +125,18 @@ def test_create_synthetic_portfolios_and_test(self): assets = [assets[i] for i in [0, 100, 200, 300, 400, 500, 600, 700, 800, 900]] detailed_results = calculate_impacts(assets, scenario="ssp585", year=2030) keys = list(detailed_results.keys()) - means = np.array([detailed_results[key][0].impact.mean_impact() for key in detailed_results.keys()]) + means = np.array( + [ + detailed_results[key][0].impact.mean_impact() + for key in detailed_results.keys() + ] + ) interesting = [k for (k, m) in zip(keys, means) if m > 0] assets_out = self.api_assets(item[0] for item in interesting[0:10]) - with open(os.path.join(cache_folder, "assets_example_industrial_activity_small.json"), "w") as f: + with open( + os.path.join(cache_folder, "assets_example_industrial_activity_small.json"), + "w", + ) as f: f.write(assets_out.model_dump_json(indent=4)) # Synthetic portfolio; real estate assets at different locations @@ -104,10 +147,17 @@ def test_create_synthetic_portfolios_and_test(self): ] detailed_results = calculate_impacts(assets, scenario="ssp585", year=2030) keys = list(detailed_results.keys()) - means = np.array([detailed_results[key][0].impact.mean_impact() for key in detailed_results.keys()]) + means = np.array( + [ + detailed_results[key][0].impact.mean_impact() + for key in detailed_results.keys() + ] + ) interesting = [k for (k, m) in zip(keys, means) if m > 0] assets_out = self.api_assets(item[0] for item in interesting[0:10]) - with open(os.path.join(cache_folder, "assets_example_real_estate_small.json"), "w") as f: + with open( + os.path.join(cache_folder, "assets_example_real_estate_small.json"), "w" + ) as f: f.write(assets_out.model_dump_json(indent=4)) self.assertAlmostEqual(1, 1) @@ -116,25 +166,38 @@ def test_thermal_power_generation_portfolio(self): cache_folder = os.environ.get("CREDENTIAL_DOTENV_DIR", os.getcwd()) asset_list = pd.read_csv(os.path.join(cache_folder, "wri-all.csv")) - filtered = asset_list.loc[asset_list["primary_fuel"].isin(["Coal", "Gas", "Nuclear", "Oil"])] + filtered = asset_list.loc[ + asset_list["primary_fuel"].isin(["Coal", "Gas", "Nuclear", "Oil"]) + ] filtered = filtered[-60 < filtered["latitude"]] longitudes = np.array(filtered["longitude"]) latitudes = np.array(filtered["latitude"]) primary_fuels = np.array( - [primary_fuel.replace(" and ", "And").replace(" ", "") for primary_fuel in filtered["primary_fuel"]] + [ + primary_fuel.replace(" and ", "And").replace(" ", "") + for primary_fuel in filtered["primary_fuel"] + ] ) # Capacity describes a maximum electric power rate. # Generation describes the actual electricity output of the plant over a period of time. capacities = np.array(filtered["capacity_mw"]) - _, continents = wd.get_countries_and_continents(latitudes=latitudes, longitudes=longitudes) + _, continents = wd.get_countries_and_continents( + latitudes=latitudes, longitudes=longitudes + ) # Power generating assets that are of interest assets = [ - ThermalPowerGeneratingAsset(latitude, longitude, type=primary_fuel, location=continent, capacity=capacity) + ThermalPowerGeneratingAsset( + latitude, + longitude, + type=primary_fuel, + location=continent, + capacity=capacity, + ) for latitude, longitude, capacity, primary_fuel, continent in zip( latitudes, longitudes, @@ -148,9 +211,13 @@ def test_thermal_power_generation_portfolio(self): year = 2030 hazard_model = calculation.get_default_hazard_model() - vulnerability_models = DictBasedVulnerabilityModels(calculation.get_default_vulnerability_models()) + vulnerability_models = DictBasedVulnerabilityModels( + calculation.get_default_vulnerability_models() + ) - results = calculate_impacts(assets, hazard_model, vulnerability_models, scenario=scenario, year=year) + results = calculate_impacts( + assets, hazard_model, vulnerability_models, scenario=scenario, year=year + ) out = [ { "asset": type(result.asset).__name__, @@ -169,7 +236,14 @@ def test_thermal_power_generation_portfolio(self): for result, key in zip(results, results.keys()) ] pd.DataFrame.from_dict(out).to_csv( - os.path.join(cache_folder, "thermal_power_generation_example_" + scenario + "_" + str(year) + ".csv") + os.path.join( + cache_folder, + "thermal_power_generation_example_" + + scenario + + "_" + + str(year) + + ".csv", + ) ) self.assertAlmostEqual(1, 1) diff --git a/tests/models/real_estate_models_test.py b/tests/models/real_estate_models_test.py index 92e7283e..1d1df679 100644 --- a/tests/models/real_estate_models_test.py +++ b/tests/models/real_estate_models_test.py @@ -1,4 +1,4 @@ -""" Test asset impact calculations.""" +"""Test asset impact calculations.""" import unittest @@ -22,9 +22,15 @@ class TestRealEstateModels(unittest.TestCase): """Tests RealEstateInundationModel.""" def test_real_estate_model_details(self): - curve = np.array([0.0596, 0.333, 0.505, 0.715, 0.864, 1.003, 1.149, 1.163, 1.163]) - store = mock_hazard_model_store_inundation(TestData.longitudes, TestData.latitudes, curve) - hazard_model = ZarrHazardModel(source_paths=get_default_source_paths(), store=store) + curve = np.array( + [0.0596, 0.333, 0.505, 0.715, 0.864, 1.003, 1.149, 1.163, 1.163] + ) + store = mock_hazard_model_store_inundation( + TestData.longitudes, TestData.latitudes, curve + ) + hazard_model = ZarrHazardModel( + source_paths=get_default_source_paths(), store=store + ) # location="Europe", type="Buildings/Residential" assets = [ @@ -35,14 +41,20 @@ def test_real_estate_model_details(self): scenario = "rcp8p5" year = 2080 - vulnerability_models = DictBasedVulnerabilityModels({RealEstateAsset: [RealEstateRiverineInundationModel()]}) + vulnerability_models = DictBasedVulnerabilityModels( + {RealEstateAsset: [RealEstateRiverineInundationModel()]} + ) - results = calculate_impacts(assets, hazard_model, vulnerability_models, scenario=scenario, year=year) + results = calculate_impacts( + assets, hazard_model, vulnerability_models, scenario=scenario, year=year + ) - hazard_bin_edges = results[ImpactKey(assets[0], RiverineInundation, scenario, year)][ - 0 - ].event.intensity_bin_edges - hazard_bin_probs = results[ImpactKey(assets[0], RiverineInundation, scenario, year)][0].event.prob + hazard_bin_edges = results[ + ImpactKey(assets[0], RiverineInundation, scenario, year) + ][0].event.intensity_bin_edges + hazard_bin_probs = results[ + ImpactKey(assets[0], RiverineInundation, scenario, year) + ][0].event.prob # check one: # the probability of inundation greater than 0.505m in a year is 1/10.0 @@ -52,15 +64,17 @@ def test_real_estate_model_details(self): np.testing.assert_almost_equal(hazard_bin_probs[1], 0.1) # check that intensity bin edges for vulnerability matrix are same as for hazard - vulnerability_intensity_bin_edges = results[ImpactKey(assets[0], RiverineInundation, scenario, year)][ - 0 - ].vulnerability.intensity_bins - np.testing.assert_almost_equal(vulnerability_intensity_bin_edges, hazard_bin_edges) + vulnerability_intensity_bin_edges = results[ + ImpactKey(assets[0], RiverineInundation, scenario, year) + ][0].vulnerability.intensity_bins + np.testing.assert_almost_equal( + vulnerability_intensity_bin_edges, hazard_bin_edges + ) # check the impact distribution the matrix is size [len(intensity_bins) - 1, len(impact_bins) - 1] - cond_probs = results[ImpactKey(assets[0], RiverineInundation, scenario, year)][0].vulnerability.prob_matrix[ - 1, : - ] + cond_probs = results[ImpactKey(assets[0], RiverineInundation, scenario, year)][ + 0 + ].vulnerability.prob_matrix[1, :] # check conditional prob for inundation intensity 0.333..0.505 mean, std = np.mean(cond_probs), np.std(cond_probs) np.testing.assert_almost_equal(cond_probs.sum(), 1) @@ -69,13 +83,17 @@ def test_real_estate_model_details(self): # probability that impact occurs between impact bin edge 1 and impact bin edge 2 prob_impact = np.dot( hazard_bin_probs, - results[ImpactKey(assets[0], RiverineInundation, scenario, year)][0].vulnerability.prob_matrix[:, 1], + results[ImpactKey(assets[0], RiverineInundation, scenario, year)][ + 0 + ].vulnerability.prob_matrix[:, 1], ) np.testing.assert_almost_equal(prob_impact, 0.19350789547968042) # no check with pre-calculated values for others: np.testing.assert_allclose( - results[ImpactKey(assets[0], RiverineInundation, scenario, year)][0].impact.prob, + results[ImpactKey(assets[0], RiverineInundation, scenario, year)][ + 0 + ].impact.prob, np.array( [ 0.02815762, @@ -97,24 +115,36 @@ def test_real_estate_model_details(self): def test_coastal_real_estate_model(self): curve = np.array([0.223, 0.267, 0.29, 0.332, 0.359, 0.386, 0.422, 0.449, 0.476]) - store = mock_hazard_model_store_inundation(TestData.coastal_longitudes, TestData.coastal_latitudes, curve) - hazard_model = ZarrHazardModel(source_paths=get_default_source_paths(), store=store) + store = mock_hazard_model_store_inundation( + TestData.coastal_longitudes, TestData.coastal_latitudes, curve + ) + hazard_model = ZarrHazardModel( + source_paths=get_default_source_paths(), store=store + ) # location="Europe", type="Buildings/Residential" assets = [ RealEstateAsset(lat, lon, location="Asia", type="Buildings/Industrial") - for lon, lat in zip(TestData.coastal_longitudes[0:1], TestData.coastal_latitudes[0:1]) + for lon, lat in zip( + TestData.coastal_longitudes[0:1], TestData.coastal_latitudes[0:1] + ) ] scenario = "rcp8p5" year = 2080 - vulnerability_models = DictBasedVulnerabilityModels({RealEstateAsset: [RealEstateCoastalInundationModel()]}) + vulnerability_models = DictBasedVulnerabilityModels( + {RealEstateAsset: [RealEstateCoastalInundationModel()]} + ) - results = calculate_impacts(assets, hazard_model, vulnerability_models, scenario=scenario, year=year) + results = calculate_impacts( + assets, hazard_model, vulnerability_models, scenario=scenario, year=year + ) np.testing.assert_allclose( - results[ImpactKey(assets[0], CoastalInundation, scenario, year)][0].impact.prob, + results[ImpactKey(assets[0], CoastalInundation, scenario, year)][ + 0 + ].impact.prob, np.array( [ 2.78081230e-02, @@ -135,14 +165,30 @@ def test_coastal_real_estate_model(self): def test_commercial_real_estate_model_details(self): curve = np.array( - [2.8302893e-06, 0.09990284, 0.21215445, 0.531271, 0.7655724, 0.99438345, 1.2871761, 1.502281, 1.7134278] + [ + 2.8302893e-06, + 0.09990284, + 0.21215445, + 0.531271, + 0.7655724, + 0.99438345, + 1.2871761, + 1.502281, + 1.7134278, + ] + ) + store = mock_hazard_model_store_inundation( + TestData.longitudes, TestData.latitudes, curve + ) + hazard_model = ZarrHazardModel( + source_paths=get_default_source_paths(), store=store ) - store = mock_hazard_model_store_inundation(TestData.longitudes, TestData.latitudes, curve) - hazard_model = ZarrHazardModel(source_paths=get_default_source_paths(), store=store) # location="South America", type="Buildings/Commercial" assets = [ - RealEstateAsset(lat, lon, location="South America", type="Buildings/Commercial") + RealEstateAsset( + lat, lon, location="South America", type="Buildings/Commercial" + ) for lon, lat in zip(TestData.longitudes[-4:-3], TestData.latitudes[-4:-3]) ] @@ -175,30 +221,38 @@ def test_commercial_real_estate_model_details(self): } ) - results = calculate_impacts(assets, hazard_model, vulnerability_models, scenario=scenario, year=year) + results = calculate_impacts( + assets, hazard_model, vulnerability_models, scenario=scenario, year=year + ) - hazard_bin_edges = results[ImpactKey(assets[0], RiverineInundation, scenario, year)][ - 0 - ].event.intensity_bin_edges - hazard_bin_probs = results[ImpactKey(assets[0], RiverineInundation, scenario, year)][0].event.prob + hazard_bin_edges = results[ + ImpactKey(assets[0], RiverineInundation, scenario, year) + ][0].event.intensity_bin_edges + hazard_bin_probs = results[ + ImpactKey(assets[0], RiverineInundation, scenario, year) + ][0].event.prob # check one: # the probability of inundation greater than 0.531271m in a year is 1/25 # the probability of inundation greater than 0.21215445m in a year is 1/10 # therefore the probability of an inundation between 0.21215445 and 0.531271 in a year is 1/10 - 1/25 - np.testing.assert_almost_equal(hazard_bin_edges[2:4], np.array([0.21215445, 0.531271])) + np.testing.assert_almost_equal( + hazard_bin_edges[2:4], np.array([0.21215445, 0.531271]) + ) np.testing.assert_almost_equal(hazard_bin_probs[2], 0.06) # check that intensity bin edges for vulnerability matrix are same as for hazard - vulnerability_intensity_bin_edges = results[ImpactKey(assets[0], RiverineInundation, scenario, year)][ - 0 - ].vulnerability.intensity_bins - np.testing.assert_almost_equal(vulnerability_intensity_bin_edges, hazard_bin_edges) + vulnerability_intensity_bin_edges = results[ + ImpactKey(assets[0], RiverineInundation, scenario, year) + ][0].vulnerability.intensity_bins + np.testing.assert_almost_equal( + vulnerability_intensity_bin_edges, hazard_bin_edges + ) # check the impact distribution the matrix is size [len(intensity_bins) - 1, len(impact_bins) - 1] - cond_probs = results[ImpactKey(assets[0], RiverineInundation, scenario, year)][0].vulnerability.prob_matrix[ - 2, : - ] + cond_probs = results[ImpactKey(assets[0], RiverineInundation, scenario, year)][ + 0 + ].vulnerability.prob_matrix[2, :] # check conditional prob for inundation intensity at 0.371712725m mean, std = np.mean(cond_probs), np.std(cond_probs) np.testing.assert_almost_equal(cond_probs.sum(), 1) @@ -207,13 +261,17 @@ def test_commercial_real_estate_model_details(self): # probability that impact occurs between impact bin edge 2 and impact bin edge 3 prob_impact = np.dot( hazard_bin_probs, - results[ImpactKey(assets[0], RiverineInundation, scenario, year)][0].vulnerability.prob_matrix[:, 2], + results[ImpactKey(assets[0], RiverineInundation, scenario, year)][ + 0 + ].vulnerability.prob_matrix[:, 2], ) np.testing.assert_almost_equal(prob_impact, 0.10040196672295522) # no check with pre-calculated values for others: np.testing.assert_allclose( - results[ImpactKey(assets[0], RiverineInundation, scenario, year)][0].impact.prob, + results[ImpactKey(assets[0], RiverineInundation, scenario, year)][ + 0 + ].impact.prob, np.array( [ 2.009085e-07, diff --git a/tests/models/wbgt_model_test.py b/tests/models/wbgt_model_test.py index 990ab20b..ea4f5195 100644 --- a/tests/models/wbgt_model_test.py +++ b/tests/models/wbgt_model_test.py @@ -6,12 +6,19 @@ from physrisk.data.pregenerated_hazard_model import ZarrHazardModel from physrisk.hazard_models.core_hazards import get_default_source_paths from physrisk.kernel.assets import Asset, IndustrialActivity -from physrisk.kernel.hazard_model import HazardDataRequest, HazardDataResponse, HazardParameterDataResponse +from physrisk.kernel.hazard_model import ( + HazardDataRequest, + HazardDataResponse, + HazardParameterDataResponse, +) from physrisk.kernel.hazards import ChronicHeat from physrisk.kernel.impact import calculate_impacts from physrisk.kernel.impact_distrib import ImpactDistrib, ImpactType from physrisk.kernel.vulnerability_model import DictBasedVulnerabilityModels -from physrisk.vulnerability_models.chronic_heat_models import ChronicHeatGZNModel, get_impact_distrib +from physrisk.vulnerability_models.chronic_heat_models import ( + ChronicHeatGZNModel, + get_impact_distrib, +) from ..data.hazard_model_store_test import TestData, mock_hazard_model_store_heat_wbgt @@ -22,10 +29,16 @@ class ExampleWBGTGZNJointModel(ChronicHeatGZNModel): results based on applying both GZN and WBGT""" def __init__(self, indicator_id: str = "mean_work_loss_high"): - super().__init__(indicator_id, ChronicHeat) # opportunity to give a model hint, but blank here + super().__init__( + indicator_id, ChronicHeat + ) # opportunity to give a model hint, but blank here def work_type_mapping(self): - return {"low": ["low", "medium"], "medium": ["medium", "low", "high"], "high": ["high", "medium"]} + return { + "low": ["low", "medium"], + "medium": ["medium", "low", "high"], + "high": ["high", "medium"], + } def get_data_requests( self, asset: Asset, *, scenario: str, year: int @@ -92,57 +105,83 @@ def get_data_requests( ), ] + wbgt_data_requests - def get_impact(self, asset: Asset, data_responses: List[HazardDataResponse]) -> ImpactDistrib: + def get_impact( + self, asset: Asset, data_responses: List[HazardDataResponse] + ) -> ImpactDistrib: """ Function to return the impact distribution of the wbgt model. """ assert isinstance(asset, IndustrialActivity) - wbgt_responses = [cast(HazardParameterDataResponse, r) for r in data_responses[2:]] + wbgt_responses = [ + cast(HazardParameterDataResponse, r) for r in data_responses[2:] + ] - hazard_paths = [cast(HazardParameterDataResponse, r).path for r in data_responses] + hazard_paths = [ + cast(HazardParameterDataResponse, r).path for r in data_responses + ] baseline_dd_above_mean = cast(HazardParameterDataResponse, data_responses[0]) scenario_dd_above_mean = cast(HazardParameterDataResponse, data_responses[1]) hours_worked = self.total_labour_hours - fraction_loss_mean_base_gzn = (baseline_dd_above_mean.parameter * self.time_lost_per_degree_day) / hours_worked + fraction_loss_mean_base_gzn = ( + baseline_dd_above_mean.parameter * self.time_lost_per_degree_day + ) / hours_worked fraction_loss_mean_scenario_gzn = ( scenario_dd_above_mean.parameter * self.time_lost_per_degree_day ) / hours_worked - fraction_loss_std_base = (baseline_dd_above_mean.parameter * self.time_lost_per_degree_day_se) / hours_worked + fraction_loss_std_base = ( + baseline_dd_above_mean.parameter * self.time_lost_per_degree_day_se + ) / hours_worked fraction_loss_std_scenario = ( scenario_dd_above_mean.parameter * self.time_lost_per_degree_day_se ) / hours_worked - baseline_work_ability = (1 - fraction_loss_mean_base_gzn) * (1 - wbgt_responses[0].parameter) - scenario_work_ability = (1 - fraction_loss_mean_scenario_gzn) * (1 - wbgt_responses[1].parameter) + baseline_work_ability = (1 - fraction_loss_mean_base_gzn) * ( + 1 - wbgt_responses[0].parameter + ) + scenario_work_ability = (1 - fraction_loss_mean_scenario_gzn) * ( + 1 - wbgt_responses[1].parameter + ) # Getting the parameters required for the uniform distribution. if asset.type in ["low", "high"]: a_historical = ( - wbgt_responses[0].parameter - abs((wbgt_responses[2].parameter - wbgt_responses[0].parameter)) / 2 + wbgt_responses[0].parameter + - abs((wbgt_responses[2].parameter - wbgt_responses[0].parameter)) / 2 ) b_historical = ( - wbgt_responses[0].parameter + abs((wbgt_responses[2].parameter - wbgt_responses[0].parameter)) / 2 + wbgt_responses[0].parameter + + abs((wbgt_responses[2].parameter - wbgt_responses[0].parameter)) / 2 ) a_scenario = ( - wbgt_responses[1].parameter - abs((wbgt_responses[3].parameter - wbgt_responses[1].parameter)) / 2 + wbgt_responses[1].parameter + - abs((wbgt_responses[3].parameter - wbgt_responses[1].parameter)) / 2 ) b_scenario = ( - wbgt_responses[1].parameter + abs((wbgt_responses[3].parameter - wbgt_responses[1].parameter)) / 2 + wbgt_responses[1].parameter + + abs((wbgt_responses[3].parameter - wbgt_responses[1].parameter)) / 2 ) elif asset.type == "medium": - a_historical = wbgt_responses[0].parameter - (wbgt_responses[2].parameter - wbgt_responses[0].parameter) / 2 - b_historical = wbgt_responses[0].parameter + (wbgt_responses[4].parameter - wbgt_responses[0].parameter) / 2 + a_historical = ( + wbgt_responses[0].parameter + - (wbgt_responses[2].parameter - wbgt_responses[0].parameter) / 2 + ) + b_historical = ( + wbgt_responses[0].parameter + + (wbgt_responses[4].parameter - wbgt_responses[0].parameter) / 2 + ) a_scenario = ( - wbgt_responses[1].parameter - abs((wbgt_responses[3].parameter - wbgt_responses[1].parameter)) / 2 + wbgt_responses[1].parameter + - abs((wbgt_responses[3].parameter - wbgt_responses[1].parameter)) / 2 ) b_scenario = ( - wbgt_responses[1].parameter + abs((wbgt_responses[5].parameter - wbgt_responses[1].parameter)) / 2 + wbgt_responses[1].parameter + + abs((wbgt_responses[5].parameter - wbgt_responses[1].parameter)) / 2 ) # Estimation of the variance @@ -166,7 +205,13 @@ def get_impact(self, asset: Asset, data_responses: List[HazardDataResponse]) -> total_work_loss_delta: float = baseline_work_ability - scenario_work_ability - return get_impact_distrib(total_work_loss_delta, std_delta, ChronicHeat, hazard_paths, ImpactType.disruption) + return get_impact_distrib( + total_work_loss_delta, + std_delta, + ChronicHeat, + hazard_paths, + ImpactType.disruption, + ) def two_variable_joint_variance(ex, varx, ey, vary): @@ -180,19 +225,28 @@ class TestChronicAssetImpact(unittest.TestCase): """Tests the impact on an asset of a chronic hazard model.""" def test_wbgt_vulnerability(self): - store = mock_hazard_model_store_heat_wbgt(TestData.longitudes, TestData.latitudes) - hazard_model = ZarrHazardModel(source_paths=get_default_source_paths(), store=store) + store = mock_hazard_model_store_heat_wbgt( + TestData.longitudes, TestData.latitudes + ) + hazard_model = ZarrHazardModel( + source_paths=get_default_source_paths(), store=store + ) # 'chronic_heat/osc/v2/mean_work_loss_high_ACCESS-CM2_historical_2005' scenario = "ssp585" year = 2050 - vulnerability_models = DictBasedVulnerabilityModels({IndustrialActivity: [ExampleWBGTGZNJointModel()]}) + vulnerability_models = DictBasedVulnerabilityModels( + {IndustrialActivity: [ExampleWBGTGZNJointModel()]} + ) assets = [ - IndustrialActivity(lat, lon, type="high") for lon, lat in zip(TestData.longitudes, TestData.latitudes) + IndustrialActivity(lat, lon, type="high") + for lon, lat in zip(TestData.longitudes, TestData.latitudes) ][:1] - results = calculate_impacts(assets, hazard_model, vulnerability_models, scenario=scenario, year=year) + results = calculate_impacts( + assets, hazard_model, vulnerability_models, scenario=scenario, year=year + ) value_test = list(results.values())[0][0].impact.prob diff --git a/tests/models/wind_models_test.py b/tests/models/wind_models_test.py index 349e2c0c..7b7e64b9 100644 --- a/tests/models/wind_models_test.py +++ b/tests/models/wind_models_test.py @@ -2,7 +2,10 @@ import tests.data.hazard_model_store_test as hms from physrisk.data.pregenerated_hazard_model import ZarrHazardModel -from physrisk.hazard_models.core_hazards import ResourceSubset, get_default_source_path_provider +from physrisk.hazard_models.core_hazards import ( + ResourceSubset, + get_default_source_path_provider, +) from physrisk.kernel.assets import RealEstateAsset from physrisk.kernel.hazards import Wind from physrisk.kernel.impact import calculate_impacts @@ -20,14 +23,25 @@ def test_wind_real_estate_model(): 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 ]) # noqa # fmt: on 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) + path = f"wind/iris/v1/max_speed_{scenario}_{year}".format( + scenario=scenario, year=year + ) hms.add_curves( - root, hms.TestData.longitudes, hms.TestData.latitudes, path, shape, intensity, return_periods, transform + root, + hms.TestData.longitudes, + hms.TestData.latitudes, + path, + shape, + intensity, + return_periods, + transform, ) provider = get_default_source_path_provider() - def select_iris_osc(candidates: ResourceSubset, scenario: str, year: int, hint=None): + def select_iris_osc( + candidates: ResourceSubset, scenario: str, year: int, hint=None + ): return candidates.with_group_id("iris_osc").first() # specify use of IRIS (OSC contribution) @@ -38,13 +52,19 @@ def select_iris_osc(candidates: ResourceSubset, scenario: str, year: int, hint=N RealEstateAsset(lat, lon, location="Asia", type="Buildings/Industrial") for lon, lat in zip(hms.TestData.longitudes[0:1], hms.TestData.latitudes[0:1]) ] - vulnerability_models = DictBasedVulnerabilityModels({RealEstateAsset: [GenericTropicalCycloneModel()]}) - results = calculate_impacts(assets, hazard_model, vulnerability_models, scenario=scenario, year=year) + vulnerability_models = DictBasedVulnerabilityModels( + {RealEstateAsset: [GenericTropicalCycloneModel()]} + ) + results = calculate_impacts( + assets, hazard_model, vulnerability_models, scenario=scenario, year=year + ) # check calculation cum_probs = 1.0 / np.array(return_periods) probs = cum_probs[:-1] - cum_probs[1:] model = GenericTropicalCycloneModel() - edges = np.interp(intensity, model.damage_curve_intensities, model.damage_curve_impacts) + edges = np.interp( + intensity, model.damage_curve_intensities, model.damage_curve_impacts + ) centres = (edges[1:] + edges[:-1]) / 2 mean_check = np.sum(probs * centres) diff --git a/tests/models/wind_turbine_models_test.py b/tests/models/wind_turbine_models_test.py index c2da4c10..6c01e371 100644 --- a/tests/models/wind_turbine_models_test.py +++ b/tests/models/wind_turbine_models_test.py @@ -15,7 +15,9 @@ class SupportsEventImpact(typing.Protocol[TAsset]): - def get_impact(self, asset: TAsset, event_data: HazardEventDataResponse) -> MultivariateDistribution: + def get_impact( + self, asset: TAsset, event_data: HazardEventDataResponse + ) -> MultivariateDistribution: pass @@ -95,9 +97,15 @@ def test_cumulative_probs(self): bins_upper = np.array([2.0, 3.0, 5.0, 5.5, 6.0]) probs = np.array([[0.1, 0.2, 0.3, 0.2, 0.2], [0.1, 0.1, 0.1, 0.1, 0.6]]) values, cum_probs = calculate_cumulative_probs(bins_lower, bins_upper, probs) - np.testing.assert_almost_equal(values, [2.0, 2.0, 3.0, 3.0, 4.0, 5.0, 5.5, 6.0, 6.0]) - np.testing.assert_almost_equal(cum_probs[0, :], [0, 0.1, 0.1, 0.3, 0.3, 0.6, 0.8, 0.8, 1.0]) - np.testing.assert_almost_equal(cum_probs[1, :], [0, 0.1, 0.1, 0.2, 0.2, 0.3, 0.4, 0.4, 1.0]) + np.testing.assert_almost_equal( + values, [2.0, 2.0, 3.0, 3.0, 4.0, 5.0, 5.5, 6.0, 6.0] + ) + np.testing.assert_almost_equal( + cum_probs[0, :], [0, 0.1, 0.1, 0.3, 0.3, 0.6, 0.8, 0.8, 1.0] + ) + np.testing.assert_almost_equal( + cum_probs[1, :], [0, 0.1, 0.1, 0.2, 0.2, 0.3, 0.4, 0.4, 1.0] + ) def test_sampling(self): """Test sampling from probability distributions comprising two lumped probabilities.""" @@ -127,7 +135,7 @@ def test_performance(self): gen = np.random.Generator(np.random.MT19937(111)) uniforms = gen.random(size=(nb_events, nb_samples)) samples = pdf.inv_cumulative_marginal_probs(uniforms) - for i in range(1000): + for _i in range(1000): uniforms = gen.random(size=(nb_events, nb_samples)) samples = pdf.inv_cumulative_marginal_probs(uniforms) print(samples) @@ -144,8 +152,12 @@ def test_rotor_damage_event_based(self): rng = np.random.Generator(np.random.MT19937(111)) asset1, asset2 = WindTurbine(), WindTurbine() nb_events = 20 - response_asset1 = HazardEventDataResponse(np.array([1.0]), np.array(rng.weibull(a=4, size=[1, nb_events]))) - response_asset2 = HazardEventDataResponse(np.array([1.0]), np.array(rng.weibull(a=4, size=[1, nb_events]))) + response_asset1 = HazardEventDataResponse( + np.array([1.0]), np.array(rng.weibull(a=4, size=[1, nb_events])) + ) + response_asset2 = HazardEventDataResponse( + np.array([1.0]), np.array(rng.weibull(a=4, size=[1, nb_events])) + ) turbine_model = WindTurbineModel() @@ -160,5 +172,9 @@ def test_rotor_damage_event_based(self): # we can then combine samples and calculate measures... # for now just sanity-check that we get the approx. 0.3 of total loss in events from placeholder. - np.testing.assert_almost_equal(np.count_nonzero(samples_asset1 == 1.0) / samples_asset1.size, 0.31) - np.testing.assert_almost_equal(np.count_nonzero(samples_asset2 == 1.0) / samples_asset2.size, 0.31) + np.testing.assert_almost_equal( + np.count_nonzero(samples_asset1 == 1.0) / samples_asset1.size, 0.31 + ) + np.testing.assert_almost_equal( + np.count_nonzero(samples_asset2 == 1.0) / samples_asset2.size, 0.31 + ) diff --git a/tests/risk_models/risk_models_test.py b/tests/risk_models/risk_models_test.py index 4a5b4da2..f7f624d0 100644 --- a/tests/risk_models/risk_models_test.py +++ b/tests/risk_models/risk_models_test.py @@ -1,17 +1,17 @@ -""" Test asset impact calculations.""" +"""Test asset impact calculations.""" from typing import Dict, Sequence -import fsspec.implementations.local as local import numpy as np from dependency_injector import providers -from physrisk import requests -from physrisk.api.v1.impact_req_resp import AssetImpactResponse, RiskMeasureKey, RiskMeasuresHelper +from physrisk.api.v1.impact_req_resp import ( + AssetImpactResponse, + RiskMeasureKey, + RiskMeasuresHelper, +) from physrisk.container import Container -from physrisk.data.inventory_reader import InventoryReader from physrisk.data.pregenerated_hazard_model import ZarrHazardModel -from physrisk.data.zarr_reader import ZarrReader from physrisk.hazard_models.core_hazards import get_default_source_paths from physrisk.kernel.assets import Asset, RealEstateAsset from physrisk.kernel.calculation import get_default_vulnerability_models @@ -31,10 +31,13 @@ from physrisk.requests import _create_risk_measures from physrisk.risk_models.generic_risk_model import GenericScoreBasedRiskMeasures from physrisk.risk_models.risk_models import RealEstateToyRiskMeasures -from tests.api.container_test import TestContainer from ..base_test import TestWithCredentials -from ..data.hazard_model_store_test import TestData, ZarrStoreMocker, inundation_return_periods +from ..data.hazard_model_store_test import ( + TestData, + ZarrStoreMocker, + inundation_return_periods, +) class TestRiskModels(TestWithCredentials): @@ -51,16 +54,22 @@ def test_risk_indicator_model(self): {RealEstateAsset: RealEstateToyRiskMeasures()}, ) measure_ids_for_asset, definitions = model.populate_measure_definitions(assets) - _, measures = model.calculate_risk_measures(assets, prosp_scens=scenarios, years=years) + _, measures = model.calculate_risk_measures( + assets, prosp_scens=scenarios, years=years + ) # how to get a score using the MeasureKey - measure = measures[MeasureKey(assets[0], scenarios[0], years[0], RiverineInundation)] + measure = measures[ + MeasureKey(assets[0], scenarios[0], years[0], RiverineInundation) + ] score = measure.score measure_0 = measure.measure_0 np.testing.assert_allclose([measure_0], [0.89306593179]) # packing up the risk measures, e.g. for JSON transmission: - risk_measures = _create_risk_measures(measures, measure_ids_for_asset, definitions, assets, scenarios, years) + risk_measures = _create_risk_measures( + measures, measure_ids_for_asset, definitions, assets, scenarios, years + ) # we still have a key, but no asset: key = RiskMeasureKey( hazard_type="RiverineInundation", @@ -75,13 +84,20 @@ def test_risk_indicator_model(self): assert measure_0 == measure_0_2 helper = RiskMeasuresHelper(risk_measures) - asset_scores, measures, definitions = helper.get_measure("CoastalInundation", scenarios[0], years[0]) + asset_scores, measures, definitions = helper.get_measure( + "CoastalInundation", scenarios[0], years[0] + ) label, description = helper.get_score_details(asset_scores[0], definitions[0]) assert asset_scores[0] == 4 def _create_assets(self): assets = [ - RealEstateAsset(TestData.latitudes[0], TestData.longitudes[0], location="Asia", type="Buildings/Industrial") + RealEstateAsset( + TestData.latitudes[0], + TestData.longitudes[0], + location="Asia", + type="Buildings/Industrial", + ) for i in range(2) ] return assets @@ -109,40 +125,70 @@ def _create_hazard_model(self, scenarios, years): source_paths = get_default_source_paths() def sp_riverine(scenario, year): - return source_paths[RiverineInundation](indicator_id="flood_depth", scenario=scenario, year=year) + return source_paths[RiverineInundation]( + indicator_id="flood_depth", scenario=scenario, year=year + ) def sp_coastal(scenario, year): - return source_paths[CoastalInundation](indicator_id="flood_depth", scenario=scenario, year=year) + return source_paths[CoastalInundation]( + indicator_id="flood_depth", scenario=scenario, year=year + ) def sp_wind(scenario, year): - return source_paths[Wind](indicator_id="max_speed", scenario=scenario, year=year) + return source_paths[Wind]( + indicator_id="max_speed", scenario=scenario, year=year + ) def sp_heat(scenario, year): - return source_paths[ChronicHeat](indicator_id="days/above/35c", scenario=scenario, year=year) + return source_paths[ChronicHeat]( + indicator_id="days/above/35c", scenario=scenario, year=year + ) def sp_fire(scenario, year): - return source_paths[Fire](indicator_id="fire_probability", scenario=scenario, year=year) + return source_paths[Fire]( + indicator_id="fire_probability", scenario=scenario, year=year + ) def sp_hail(scenario, year): - return source_paths[Hail](indicator_id="days/above/5cm", scenario=scenario, year=year) + return source_paths[Hail]( + indicator_id="days/above/5cm", scenario=scenario, year=year + ) def sp_drought(scenario, year): - return source_paths[Drought](indicator_id="months/spei3m/below/-2", scenario=scenario, year=year) + return source_paths[Drought]( + indicator_id="months/spei3m/below/-2", scenario=scenario, year=year + ) def sp_precipitation(scenario, year): - return source_paths[Precipitation](indicator_id="max/daily/water_equivalent", scenario=scenario, year=year) + return source_paths[Precipitation]( + indicator_id="max/daily/water_equivalent", scenario=scenario, year=year + ) mocker = ZarrStoreMocker() return_periods = inundation_return_periods() - flood_histo_curve = np.array([0.0596, 0.333, 0.505, 0.715, 0.864, 1.003, 1.149, 1.163, 1.163]) - flood_projected_curve = np.array([0.0596, 0.333, 0.605, 0.915, 1.164, 1.503, 1.649, 1.763, 1.963]) + flood_histo_curve = np.array( + [0.0596, 0.333, 0.505, 0.715, 0.864, 1.003, 1.149, 1.163, 1.163] + ) + flood_projected_curve = np.array( + [0.0596, 0.333, 0.605, 0.915, 1.164, 1.503, 1.649, 1.763, 1.963] + ) for path in [sp_riverine("historical", 1980), sp_coastal("historical", 1980)]: - mocker.add_curves_global(path, TestData.longitudes, TestData.latitudes, return_periods, flood_histo_curve) + mocker.add_curves_global( + path, + TestData.longitudes, + TestData.latitudes, + return_periods, + flood_histo_curve, + ) for path in [sp_riverine("rcp8p5", 2050), sp_coastal("rcp8p5", 2050)]: mocker.add_curves_global( - path, TestData.longitudes, TestData.latitudes, return_periods, flood_projected_curve + path, + TestData.longitudes, + TestData.latitudes, + return_periods, + flood_projected_curve, ) mocker.add_curves_global( @@ -232,7 +278,9 @@ def sp_precipitation(scenario, year): [70], ) - return ZarrHazardModel(source_paths=get_default_source_paths(), store=mocker.store) + return ZarrHazardModel( + source_paths=get_default_source_paths(), store=mocker.store + ) def test_via_requests(self): scenarios = ["ssp585"] @@ -256,11 +304,19 @@ def test_via_requests(self): container = Container() class TestHazardModelFactory(HazardModelFactory): - def hazard_model(self, interpolation: str = "floor", provider_max_requests: Dict[str, int] = ...): + def hazard_model( + self, + interpolation: str = "floor", + provider_max_requests: Dict[str, int] = {}, + ): return hazard_model - container.override_providers(hazard_model_factory=providers.Factory(TestHazardModelFactory)) - container.override_providers(config=providers.Configuration(default={"zarr_sources": ["embedded"]})) + container.override_providers( + hazard_model_factory=providers.Factory(TestHazardModelFactory) + ) + container.override_providers( + config=providers.Configuration(default={"zarr_sources": ["embedded"]}) + ) container.override_providers(inventory_reader=None) container.override_providers(zarr_reader=None) @@ -274,7 +330,9 @@ def hazard_model(self, interpolation: str = "floor", provider_max_requests: Dict # vulnerability_models=DictBasedVulnerabilityModels(get_default_vulnerability_models()), # ) res = next( - ma for ma in response.risk_measures.measures_for_assets if ma.key.hazard_type == "RiverineInundation" + ma + for ma in response.risk_measures.measures_for_assets + if ma.key.hazard_type == "RiverineInundation" ) np.testing.assert_allclose(res.measures_0, [0.89306593179, 0.89306593179]) # json_str = json.dumps(response.model_dump(), cls=NumpyArrayEncoder) @@ -283,7 +341,9 @@ def test_generic_model(self): scenarios = ["rcp8p5"] years = [2050] - assets = [Asset(TestData.latitudes[0], TestData.longitudes[0]) for i in range(2)] + assets = [ + Asset(TestData.latitudes[0], TestData.longitudes[0]) for i in range(2) + ] hazard_model = self._create_hazard_model(scenarios, years) model = AssetLevelRiskModel( @@ -292,7 +352,10 @@ def test_generic_model(self): {Asset: GenericScoreBasedRiskMeasures()}, ) measure_ids_for_asset, definitions = model.populate_measure_definitions(assets) - _, measures = model.calculate_risk_measures(assets, prosp_scens=scenarios, years=years) + _, measures = model.calculate_risk_measures( + assets, prosp_scens=scenarios, years=years + ) np.testing.assert_approx_equal( - measures[MeasureKey(assets[0], scenarios[0], years[0], Wind)].measure_0, 214.01549835205077 + measures[MeasureKey(assets[0], scenarios[0], years[0], Wind)].measure_0, + 214.01549835205077, ) diff --git a/tox.ini b/tox.ini index 9ccadae3..e1df6a6b 100644 --- a/tox.ini +++ b/tox.ini @@ -4,32 +4,23 @@ envlist = py38, static [testenv] deps = pytest + pandas>=2.0.3 + dependency-injector>=4.41.0 + geopandas<1.0,>=0.13.2 commands = pytest {posargs} [testenv:static] deps = mypy - isort - black - flake8 commands = mypy --install-types --non-interactive src - isort --check . - black --check . - flake8 src [testenv:cov] usedevelop = True deps = pytest-cov -commands = pytest --cov-report=html {posargs} + pandas>=2.0.3 + dependency-injector>=4.41.0 + geopandas<1.0,>=0.13.2 -[flake8] -count = True -max-line-length = 120 -max-complexity = 10 -# Allow __init__ files to have unused imports. -per-file-ignores = __init__.py:F401 -extend-ignore = - # Allow spacing before colon (to favor Black). - E203 +commands = pytest --cov-report=html {posargs}