valid.py 7.63 KB
Newer Older
1
2
3
4
"""
Defines a set of input validation methods to check physically correct or
consistent quantities
"""
5
6
7
import numpy as np
from copy import deepcopy
from geopandas import GeoDataFrame
8
from typing import Tuple, Optional, Union, Dict
9
from openquake.hazardlib.geo.nodalplane import NodalPlane
10
11
12
13
14
15
16
17
18
from openquake.hazardlib.gsim import get_available_gsims
from shakyground2.magnitude_scaling_relations import (
    BaseScalingRelation,
    PEERScalingRelation,
)


# OpenQuake GMPE List
GSIM_LIST = get_available_gsims()
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96


def longitude(lon: float) -> float:
    """
    Verify that the longitude is between -180 and 180 degrees
    """
    if lon < -180.0 or lon > 180.0:
        raise ValueError("Longitude %.4f is not between -180 and 180" % lon)
    return lon


def latitude(lat: float) -> float:
    """
    Verify that the latitude is between -90 and 90 degrees
    """
    if lat < -90.0 or lat > 90.0:
        raise ValueError("Latitude %.4f is not between -90 and 90" % lat)
    return lat


def positive_float(value: float, key: str) -> float:
    """
    Verify the value is a positive float (or zero) and raise error otherwise
    """
    if value < 0.0:
        raise ValueError("%s must be positive float or zero, %.6f input" % (key, value))
    return value


def strike(value: Optional[float]) -> float:
    """
    Verify that the strike is within the range 0 to 360 degrees, or else return None if
    unspecified
    """
    if value is not None and (value < 0 or value >= 360.0):
        raise ValueError("Strike %.2f not in the range 0 to 360 degrees" % value)
    return value


def dip(value: Optional[float]) -> float:
    """
    Verify that the dip is within the range 0 to 90 degrees, or else return None if unspecified
    """
    if value is not None and (value < 0 or value > 90.0):
        raise ValueError("Dip %.2f not in the range 0 to 90 degrees" % value)
    return value


def rake(value: Optional[float]) -> float:
    """
    Verify that the rake is within the range -180 to 180 degrees, according to the Aki &
    Richards (1980) convention, or else return None if unspecified
    """
    if value is not None and (value < -180.0 or value > 180.0):
        raise ValueError("Rake %.2f not in the range -180 to 180 degrees" % value)
    return value


def mechanism(
    istrike: float, idip: float, irake: float
) -> Union[Tuple[Optional[float], Optional[float], Optional[float]], NodalPlane]:
    """
    Verifies that a valid focal mechanism is defined. A valid focal mechanism requires a
    strike, dip and rake value, which are limited to the range (0, 360), (0, 90) and
    (-180, 180) respectively. The Aki & Richards (1980) definition of rake is applied.

    Note that this only checks validity of the mechanism in terms of input value,
    not in terms of the complete physical properties of the mechanism
    """
    mechanism = (strike(istrike), dip(idip), rake(irake))
    if None in mechanism:
        # Not a valid mechanism, return tuple
        return (istrike, idip, irake)
    else:
        # Is a valid mechanism, so return openquake.hazardlib.nodalplane.NodalPlane object
        return NodalPlane(*mechanism)


97
98
99
100
101
102
103
104
105
106
107
108
def focal_mechanism(focal_mech: Optional[Dict]) -> Dict:
    """
    A focal mechanism is represented by two orthogonal planes, each plane described by a
    strike, dip and rake
    """
    if not focal_mech:
        return focal_mech
    assert "nodal_plane_1" in focal_mech, "Focal mechanism missing nodal plane 1"
    assert "nodal_plane_2" in focal_mech, "Focal mechanism missing nodal plane 2"
    focal_mechanism = {}
    for plane in ["nodal_plane_1", "nodal_plane_2"]:
        focal_mechanism[plane] = mechanism(
109
110
111
            focal_mech[plane]["strike"],
            focal_mech[plane]["dip"],
            focal_mech[plane]["rake"],
112
113
114
115
        )
    return focal_mechanism


116
117
118
119
120
121
122
123
124
def seismogenic_thickness(
    upper_seismo_depth: float, lower_seismo_depth: float
) -> Tuple[float, float]:
    """
    Verifies that a valid seismogenic thickness of the crust is defined
    """
    usd = positive_float(upper_seismo_depth, "upper seismogenic depth")
    lsd = positive_float(lower_seismo_depth, "lower seismogenic depth")
    if lsd < usd:
Marius Kriegerowski's avatar
linting    
Marius Kriegerowski committed
125
126
127
128
        raise ValueError(
            "Lower seismogenic depth %.2f km shallower than upper seismogenic "
            "depth %.2f km" % (lsd, usd)
        )
129
130
131
132
133
134
135
136
137
138
    return usd, lsd


def hypocenter_position(hypo_pos: Tuple[float, float]) -> Tuple[float, float]:
    """
    Verifies that a hypocenter position is valid within the range [0, 1] for both the
    along-strike and down-dip cases
    """
    along_strike, down_dip = hypo_pos
    if along_strike < 0.0 or along_strike > 1.0:
Marius Kriegerowski's avatar
linting    
Marius Kriegerowski committed
139
140
141
        raise ValueError(
            "Along strike position %.3f should be in the range 0 to 1" % along_strike
        )
142
143
144
    if down_dip < 0.0 or down_dip > 1.0:
        raise ValueError("Down dip position %.3f should be in the range 0 to 1" % down_dip)
    return along_strike, down_dip
145
146
147
148
149
150
151
152
153
154
155


def scaling_relation(msr: Optional[BaseScalingRelation]):
    """
    Verifies that the magnitude scaling relation is one supported by the
    software, or return a default is none is provided
    """
    if not msr:
        # If no relation is defined then use the default
        return PEERScalingRelation()
    if not isinstance(msr, BaseScalingRelation):
Marius Kriegerowski's avatar
linting    
Marius Kriegerowski committed
156
157
158
        raise TypeError(
            "Magnitude Scaling Relation %s not instance of BaseScalingRelation" % str(msr)
        )
159
    return msr
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210


def regionalization_mapping(mapping: Dict) -> Dict:
    """
    Velidates a ground motion mapping to parse the ground motion model strings to instances
    of the ground motion models. Also checks the weights sum correctly to 1.0
    """
    new_mapping = {}
    for key in mapping:
        new_mapping[key] = []
        # Verify that weights sum to 1.0
        weight_sum = sum([gmm["weight"] for gmm in mapping[key]])
        weight_check = np.isclose(weight_sum, 1.0)
        assert (
            weight_check
        ), "Ground motion model weights for region %s do not sum to 1 " "(sum = %.6f)" % (
            key,
            weight_sum,
        )
        for gmm in deepcopy(mapping[key]):
            gmm_id = gmm.pop("id")
            gmm_weight = gmm.pop("weight")
            gmm_name = gmm.pop("model")
            new_mapping[key].append((gmm_id, GSIM_LIST[gmm_name](**gmm), gmm_weight))
    return new_mapping


def regionalization(regions: GeoDataFrame, mapping: Dict) -> Tuple[GeoDataFrame, Dict]:
    """
    A regionalisation is represented by a geometry set (as a geopandas.GeoDataFrame) and a
    corresponding dictionary to map the regions in the geometry set to a set of ground motion
    models (as strings of the OpenQuake input names) and respective weights. Function verifies
    that the region file has the necessary information and that a mapping for each region is
    present. Returns the region set and the mapping with instantiated ground motion model.
    """
    if not regions.crs:
        # If no coordinate reference system is defined then assume WGS84
        regions.crs = {"init": "epsg:4326"}
    if str(regions.crs) != "+init=epsg:4326 +type=crs":
        regions = regions.to_crs({"init": "epsg:4326"})
    # Verify that the region set has the necessary columns
    for col in ["REGION", "UPPER DEPTH", "LOWER DEPTH", "geometry"]:
        if col not in regions.columns:
            raise IOError("Regionalization has missing attribute %s" % col)
    # Verify that every region in the regionalization has a corresponding mapping
    for region in regions.REGION.unique():
        if region not in mapping:
            raise IOError(
                "Region %s has no corresponding ground motion model in mapping" % region
            )
    return regions, regionalization_mapping(mapping)