IMPACT: set the default maximum_distance based on trt and mag by interpolation

openquake/commands/mosaic.py

@@ -302,7 +302,7 @@ def impact(exposure_hdf5='', *,
            rupfname=FAMOUS,
            stations='',
            mosaic_model='',
-           maximum_distance='300',
+           maximum_distance='',
            maximum_distance_stations='',
            asset_hazard_distance='15',
            number_of_ground_motion_fields='10'):

openquake/commonlib/oqvalidation.py

@@ -977,7 +977,7 @@
         'BC': {'display_name': 'BC - Soft rock', 'vs30': 760},
         'C': {'display_name': 'C - Very dense sand or hard clay', 'vs30': 530},
         'CD': {'display_name':
-               'CD - Dense sand or ery stiff clay', 'vs30': 365},
+               'CD - Dense sand or very stiff clay', 'vs30': 365},
         'D': {'display_name':
               'D - Medium dense sand or stiff clay', 'vs30': 260},
         'DE': {'display_name':
@@ -1370,7 +1370,8 @@ def _set_hazard_imtls(self, names_vals):
                     'Each IMT must have the same number of levels, instead '
                     'you have %s' % dic)
         elif 'intensity_measure_types' in names_vals:
-            self.hazard_imtls = dict.fromkeys(self.intensity_measure_types, [0])
+            self.hazard_imtls = dict.fromkeys(
+                self.intensity_measure_types, [0])
             delattr(self, 'intensity_measure_types')
 
     def _normalize_minimum_intensity(self):
@@ -1426,6 +1427,11 @@ def __init__(self, **names_vals):
         if hasattr(self, 'maximum_distance'):
             # can be missing in post-calculations
             self.maximum_distance.cut(self.minimum_magnitude)
+            logging.info(
+                f'Tectonic region type:'
+                f' {self.tectonic_region_type}; '
+                f'Maximum source-to-site distance:'
+                f' {self.maximum_distance[self.tectonic_region_type]}')
 
         self.check_hazard()
         self.check_gsim_lt()
@@ -1590,7 +1596,8 @@ def check_hazard(self):
                 self.raise_invalid('poes_disagg or iml_disagg must be set')
             elif self.poes_disagg and self.iml_disagg:
                 self.raise_invalid(
-                    'iml_disagg and poes_disagg cannot be set at the same time')
+                    'iml_disagg and poes_disagg cannot be set'
+                    ' at the same time')
             if not self.disagg_bin_edges:
                 for k in ('mag_bin_width', 'distance_bin_width',
                           'coordinate_bin_width', 'num_epsilon_bins'):
@@ -1620,7 +1627,7 @@ def set_loss_types(self):
             with datastore.read(self.hazard_calculation_id) as ds:
                 self._parent = ds['oqparam']
             if not self.total_losses:
-                self.total_losses = self._parent.total_losses            
+                self.total_losses = self._parent.total_losses
         else:
             self._parent = None
         # set all_cost_types
@@ -1699,21 +1706,22 @@ def check_gsims(self, gsims):
         has_sites = self.sites or 'site_model' in self.inputs
         if not has_sites:
             return
-        
+
         # Check if using IMT-dependent weights per GMM so can raise error
         # if not supported and non-zero weight has been assigned to IMT
         weight = {
-            (gmm, from_string(imt).name): 0.0 for gmm in gsims for imt
-              in self.imtls}
-        imt_weighted = {gmm: False for gmm in gsims} # Req to track if weighted
-
+            (gmm, from_string(imt).name): 0.0
+            for gmm in gsims for imt in self.imtls}
+        imt_weighted = {gmm: False
+                        for gmm in gsims}  # Req to track if weighted
+
         if "gsim_logic_tree" in self.inputs:
             branches = GsimLogicTree(self.inputs["gsim_logic_tree"]).branches
             for branch in branches:
                 if branch.gsim not in imt_weighted:
                     continue
                 gimts = [cls.__name__ for cls in
-                        branch.gsim.DEFINED_FOR_INTENSITY_MEASURE_TYPES]
+                         branch.gsim.DEFINED_FOR_INTENSITY_MEASURE_TYPES]
                 if isinstance(branch.weight, ImtWeight):
                     imt_weighted[branch.gsim] = True
                     for imt in self.imtls:
@@ -1725,25 +1733,25 @@ def check_gsims(self, gsims):
         for gsim in gsims:
             params = getattr(gsim, 'params', {})
             ok_imts = {cls.__name__ for cls in gsim.
-                    DEFINED_FOR_INTENSITY_MEASURE_TYPES}
+                       DEFINED_FOR_INTENSITY_MEASURE_TYPES}
             if 'conditional_gmpe' in params:
                 ok_imts |= set(params['conditional_gmpe'])
             elif hasattr(gsim, 'gmpe'):
                 ok_imts |= {cls.__name__ for cls in gsim.gmpe.
-                        DEFINED_FOR_INTENSITY_MEASURE_TYPES}
+                            DEFINED_FOR_INTENSITY_MEASURE_TYPES}
             if ok_imts:
                 invalid_imts = imts - ok_imts
                 # Don't collect IMT if unsupported buts has zero wt
                 bad_wt_imts = [imt for imt in invalid_imts if
-                            weight[(gsim, imt)] > 0]
+                               weight[(gsim, imt)] > 0]
                 if (invalid_imts and not bad_wt_imts and not
-                    imt_weighted[gsim]
-                    ):
+                        imt_weighted[gsim]):
                     raise ValueError(
                         'The IMT %s is not accepted by the GSIM %s' % (
                             ', '.join(invalid_imts), gsim))
-                if bad_wt_imts: # If IMT-dependent and wt greater than
-                                # zero for unsupported IMT then raise
+                if bad_wt_imts:
+                    # If IMT-dependent and wt greater than
+                    # zero for unsupported IMT then raise
                     raise ValueError(
                         'Non-zero weights assigned to unsupported IMT(s) '
                         '(%s) for the GSIM %s' % (', '.join(bad_wt_imts), gsim)
@@ -1765,6 +1773,7 @@ def check_gsims(self, gsims):
                             raise ValueError(
                                 'Please set a value for %r, this is required '
                                 'by the GSIM %s' % (param_name, gsim))
+
     @property
     def tses(self):
         """
@@ -2673,7 +2682,7 @@ def to_ini(key, val):
     elif key in ('reqv_ignore_sources', 'poes', 'quantiles', 'disagg_outputs',
                  'source_id', 'source_nodes', 'soil_intensities'):
         return f"{key} = {' '.join(map(str, val))}"
-    elif key  == 'cache':
+    elif key == 'cache':
         # parameter not affecting the numbers
         return ''
     else:

openquake/engine/impact.py

@@ -60,7 +60,7 @@ def main_web(allparams, jobctxs,
 def main_cmd(usgs_id, rupture_file=None,
              callback=trivial_callback, *,
              shakemap_version=None, time_event='day',
-             maximum_distance='300', mosaic_model=None, trt=None,
+             maximum_distance='', mosaic_model=None, trt=None,
              truncation_level='3',
              number_of_ground_motion_fields='10', asset_hazard_distance='15',
              ses_seed='42', local_timestamp='',

openquake/hazardlib/calc/filters.py

@@ -39,6 +39,132 @@
 MAX_DISTANCE = 2000  # km, ultra big distance used if there is no filter
 trt_smr = operator.attrgetter('trt_smr')
 
+
+# The default value of the maximum_distance parameter is based on the
+# tectonic region type and magnitude by interpolation using the following
+# (mag, rrup) points, rounding the interpolated distance to the nearest 5 km.
+# These are calibrated approximately to the distance beyond which
+# the μ + 1σ  PGA ≤ 0.06g for vs30 = 360m/s, i.e., assets beyond
+# this distance are highly unlikely to incur damage and can safely
+# be excluded from the analysis. There might be the rare earthquakes
+# where the shaking is high even beyond these distances
+# (eg. 1985 M8.0 Mexico City earthquake) but in such cases the user
+# can always override the default value.
+MAGDIST_BY_TRT = {
+    'ACTIVE_SHALLOW_CRUST': [
+        (5.0,  30), (5.5,  40), (6.0,  50), (6.5,  75),
+        (7.0, 110), (7.5, 150), (8.0, 200),
+    ],
+    'STABLE_CONTINENTAL': [
+        (5.0,  40), (5.5,  50), (6.0,  75), (6.5, 125),
+        (7.0, 200), (7.5, 300), (8.0, 400),
+    ],
+    'SUBDUCTION_INTERFACE': [
+        (6.0,  60), (6.5,  90), (7.0, 120), (7.5, 170),
+        (8.0, 225), (8.5, 270), (9.0, 305), (9.5, 340),
+    ],
+    'SUBDUCTION_INTRASLAB': [
+        (6.0, 100), (6.5, 150), (7.0, 200), (7.5, 250), (8.0, 300),
+    ],
+    'VOLCANIC': [(4.5, 10), (5.0, 15)]
+}
+
+
+# Mapping from every known mosaic TRT string -> canonical TRT string
+MOSAIC_TRT_TO_CANONICAL_TRT = {
+    # Active Shallow Crust
+    "Active Shallow Crust ALS":           "ACTIVE_SHALLOW_CRUST",
+    "Active Shallow Crust CAM":           "ACTIVE_SHALLOW_CRUST",
+    "Active Shallow Crust CAN":           "ACTIVE_SHALLOW_CRUST",
+    "Active Shallow Crust CCA":           "ACTIVE_SHALLOW_CRUST",
+    "Active Shallow Crust CHN":           "ACTIVE_SHALLOW_CRUST",
+    "Active Shallow Crust EMME":          "ACTIVE_SHALLOW_CRUST",
+    "Active Shallow Crust JPN":           "ACTIVE_SHALLOW_CRUST",
+    "Active Shallow Crust MEX":           "ACTIVE_SHALLOW_CRUST",
+    "Active Shallow Crust NWA-NEA":       "ACTIVE_SHALLOW_CRUST",
+    "Active Shallow Crust PHL":           "ACTIVE_SHALLOW_CRUST",
+    "Active Shallow Crust SAM":           "ACTIVE_SHALLOW_CRUST",
+    "Active Shallow Crust SEA":           "ACTIVE_SHALLOW_CRUST",
+    "Active Shallow Crust TEM":           "ACTIVE_SHALLOW_CRUST",
+    "Active Shallow Crust TWN":           "ACTIVE_SHALLOW_CRUST",
+    "Active Shallow Crust USA":           "ACTIVE_SHALLOW_CRUST",
+    "Active Shallow Crust":               "ACTIVE_SHALLOW_CRUST",
+    "active shallow crust":               "ACTIVE_SHALLOW_CRUST",
+    "Himalayan Thrust":                   "ACTIVE_SHALLOW_CRUST",
+    "Iceland Atlantic Active Region":     "ACTIVE_SHALLOW_CRUST",
+    "Oceanic Crust without Volcanism":    "ACTIVE_SHALLOW_CRUST",
+    "Philippine Active Shallow Crust":    "ACTIVE_SHALLOW_CRUST",
+    "Shallow Default":                    "ACTIVE_SHALLOW_CRUST",
+    "Spreading and Transforms":           "ACTIVE_SHALLOW_CRUST",
+    "TECTONIC_REGION_2":                  "ACTIVE_SHALLOW_CRUST",
+    "TECTONIC_REGION_3":                  "ACTIVE_SHALLOW_CRUST",
+    # Stable Continental
+    "Active-Stable Shallow Crust":        "STABLE_CONTINENTAL",
+    "Craton":                             "STABLE_CONTINENTAL",
+    "Cratonic Crust":                     "STABLE_CONTINENTAL",
+    "cratonic":                           "STABLE_CONTINENTAL",
+    "Cratonic":                           "STABLE_CONTINENTAL",
+    "intraplate margin lower":            "STABLE_CONTINENTAL",
+    "intraplate margin upper":            "STABLE_CONTINENTAL",
+    "Non_cratonic":                       "STABLE_CONTINENTAL",
+    "SCR Non Craton":                     "STABLE_CONTINENTAL",
+    "SSC KOR":                            "STABLE_CONTINENTAL",
+    "SSC":                                "STABLE_CONTINENTAL",
+    "Stable Continental Crust":           "STABLE_CONTINENTAL",
+    "Stable Continental Region":          "STABLE_CONTINENTAL",
+    "Stable Continental Crust NWA-NEA":   "STABLE_CONTINENTAL",
+    "Stable Shallow Crust EMME":          "STABLE_CONTINENTAL",
+    "Stable Shallow Crust MEX":           "STABLE_CONTINENTAL",
+    "Stable Shallow Crust NWA-NEA":       "STABLE_CONTINENTAL",
+    "Stable Shallow Crust SAM":           "STABLE_CONTINENTAL",
+    "Stable Shallow Crust USA":           "STABLE_CONTINENTAL",
+    "Stable Shallow Crust":               "STABLE_CONTINENTAL",
+    "stable shallow crust":               "STABLE_CONTINENTAL",
+    "TECTONIC_REGION_1":                  "STABLE_CONTINENTAL",
+    "TECTONIC_REGION_4":                  "STABLE_CONTINENTAL",
+    "Tectonic_Type_A":                    "STABLE_CONTINENTAL",
+    "Tectonic_Type_B":                    "STABLE_CONTINENTAL",
+    "Tectonic_Type_C":                    "STABLE_CONTINENTAL",
+    "Tectonic_Type_D":                    "STABLE_CONTINENTAL",
+    "Tectonic_Type_E":                    "STABLE_CONTINENTAL",
+    # Subduction Interface
+    "Interface subduction zone":                       "SUBDUCTION_INTERFACE",
+    "Philippine Subduction":                           "SUBDUCTION_INTERFACE",
+    "Subduction Interface - North East Correction":    "SUBDUCTION_INTERFACE",
+    "Subduction Interface - South West Correction":    "SUBDUCTION_INTERFACE",
+    "Subduction Interface CAM":                        "SUBDUCTION_INTERFACE",
+    "Subduction Interface CAN":                        "SUBDUCTION_INTERFACE",
+    "Subduction Interface CCA":                        "SUBDUCTION_INTERFACE",
+    "Subduction Interface EMME":                       "SUBDUCTION_INTERFACE",
+    "Subduction Interface SAM":                        "SUBDUCTION_INTERFACE",
+    "Subduction Interface":                            "SUBDUCTION_INTERFACE",
+    "subduction interface":                            "SUBDUCTION_INTERFACE",
+    # Subduction IntraSlab
+    "Subduction IntraSlab":                            "SUBDUCTION_INTRASLAB",
+    "Deep Crust 1":                                    "SUBDUCTION_INTRASLAB",
+    "Deep Seismicity EMME":                            "SUBDUCTION_INTRASLAB",
+    "Deep Seismicity":                                 "SUBDUCTION_INTRASLAB",
+    "Intraslab subduction zone":                       "SUBDUCTION_INTRASLAB",
+    "Non-Subduction Deep":                             "SUBDUCTION_INTRASLAB",
+    "Subduction Inslab EMME":                          "SUBDUCTION_INTRASLAB",
+    "Subduction Inslab":                               "SUBDUCTION_INTRASLAB",
+    "Subduction IntraSlab - North East Correction":    "SUBDUCTION_INTRASLAB",
+    "Subduction IntraSlab - South West Correction":    "SUBDUCTION_INTRASLAB",
+    "Subduction IntraSlab CAM":                        "SUBDUCTION_INTRASLAB",
+    "Subduction IntraSlab CCA":                        "SUBDUCTION_INTRASLAB",
+    "Subduction IntraSlab SAM":                        "SUBDUCTION_INTRASLAB",
+    "Subduction Intraslab":                            "SUBDUCTION_INTRASLAB",
+    "subduction intraslab":                            "SUBDUCTION_INTRASLAB",
+    "Subduction IntraSlab30":                          "SUBDUCTION_INTRASLAB",
+    "Subduction IntraSlab55":                          "SUBDUCTION_INTRASLAB",
+    "Subduction":                                      "SUBDUCTION_INTRASLAB",
+    "DeepSeismicity":                                  "SUBDUCTION_INTRASLAB",
+    # Volcanic
+    "Oceanic Crust with Volcanism":       "VOLCANIC",
+    "Volcanic":                           "VOLCANIC",
+}
+
+
 class FilteredAway(Exception):
     pass
 
@@ -289,6 +415,12 @@ def get_dist_bins(self, trt, nbins=51):
         return .01 + numpy.arange(nbins) * self(trt) / (nbins - 1)
 
 
+DEFAULT_INTEGRATION_DISTANCE = IntegrationDistance(
+        (mosaic_trt, MAGDIST_BY_TRT[canonical_trt])
+        for mosaic_trt, canonical_trt in MOSAIC_TRT_TO_CANONICAL_TRT.items()
+)
+
+
 def split_source(src):
     """
     :param src: a splittable (or not splittable) source
@@ -324,7 +456,7 @@ def split_source(src):
         if has_grp_id:
             split.grp_id = src.grp_id
         offset += split.num_ruptures
-        #split.nsites = src.nsites
+        # split.nsites = src.nsites
     return splits
 
 
@@ -467,6 +599,7 @@ def __repr__(self):
         return '<%s mag=%.1f dist=%.0f>' % (self.__class__.__name__,
                                             self.rup.mag, self.dist)
 
+
 class SourceFilter(object):
     """
     Filter objects have a .filter method yielding filtered sources
@@ -613,7 +746,7 @@ def get_close(self, secparams):
         distr = distance.cdist(xyz, spherical_to_cartesian(tr0, tr1))
         dists = numpy.min([distl, distr], axis=0)  # shape (N, S)
         maxdist = self.integration_distance.y[-1]
-        return (dists <= maxdist).sum(axis=0) # shape (N, S) => S
+        return (dists <= maxdist).sum(axis=0)  # shape (N, S) => S
 
     def __getitem__(self, slc):
         if slc.start is None and slc.stop is None:

openquake/hazardlib/shakemap/validate.py

@@ -29,6 +29,7 @@
 from openquake.hazardlib.geo.utils import SiteAssociationError
 from openquake.hazardlib.scalerel import get_available_magnitude_scalerel
 from openquake.hazardlib.nrml import validators as nrml_validators
+from openquake.hazardlib.calc.filters import DEFAULT_INTEGRATION_DISTANCE
 
 MOSAIC_DIR = config.directory.mosaic_dir or os.path.dirname(mosaic.__file__)
 
@@ -46,6 +47,7 @@ class ImpactParam:
     rupture_file: str = None
     station_data_file: str = None
     mmi_file: str = None
+    maximum_distance: float = None
     maximum_distance_stations: float = None
     mosaic_model: str = None
     trt: str = None
@@ -99,6 +101,8 @@ def get_oqparams(self, usgs_id, mosaic_models, trts, use_shakemap):
             params['gsim'] = '[FromFile]'
         if self.local_timestamp is not None:
             params['local_timestamp'] = self.local_timestamp
+        if self.maximum_distance is not None:
+            params['maximum_distance'] = str(self.maximum_distance)
         if self.maximum_distance_stations is not None:
             params['maximum_distance_stations'] = str(
                 self.maximum_distance_stations)
@@ -110,6 +114,9 @@ def get_oqparams(self, usgs_id, mosaic_models, trts, use_shakemap):
                 params['description'] = (
                     f'{rupdic["usgs_id"]}: M {rupdic["mag"]}'
                     f' ({rupdic["lat"]}, {rupdic["lon"]})')
+        if self.maximum_distance is None:
+            md = DEFAULT_INTEGRATION_DISTANCE[self.trt]
+            params['maximum_distance'] = str({self.trt: md})
         oq = readinput.get_oqparam(params)
         # NB: fake h5 to cache `get_site_model` and avoid multiple associations
         _sitecol, assetcol, _discarded, _exp = readinput.get_sitecol_assetcol(
@@ -198,7 +205,7 @@ def get_oqparams(self, usgs_id, mosaic_models, trts, use_shakemap):
     'time_event': 'day',
     'dip': '90',
     'strike': '0',
-    'maximum_distance': '300',
+    'maximum_distance': '',
     'truncation_level': '3',
     'number_of_ground_motion_fields': '100',
     'asset_hazard_distance': '15',
@@ -218,8 +225,8 @@ def get_oqparams(self, usgs_id, mosaic_models, trts, use_shakemap):
     'use_shakemap_from_usgs': 'Use ShakeMap from the USGS',
     'use_pnt_rup_from_usgs': 'Use point rupture from the USGS',
     'build_rup_from_usgs': 'Build rupture from USGS nodal plane solutions',
-    'use_shakemap_fault_rup_from_usgs': (
-        'Use ShakeMap fault rupture from the USGS'),
+    'use_shakemap_fault_rup_from_usgs':
+        'Use ShakeMap fault rupture from the USGS',
     'use_finite_fault_model_from_usgs': 'Use finite fault model from the USGS',
     'provide_rup': 'Provide earthquake rupture in OpenQuake NRML format',
     'provide_rup_params': 'Provide earthquake rupture parameters',
@@ -270,7 +277,8 @@ def _validate(POST):
         try:
             value = validation_func(POST.get(field))
         except Exception as exc:
-            blankable = ['dip', 'strike', 'maximum_distance_stations',
+            blankable = ['dip', 'strike',
+                         'maximum_distance', 'maximum_distance_stations',
                          'local_timestamp']
             if field in blankable and POST.get(field) == '':
                 if field in inputdic:

openquake/hazardlib/valid.py

@@ -725,7 +725,7 @@ def positivefloatorsentinel(value):
     :returns: positive float or -999 (sentinel)
     """
     f = float(not_empty(value))
-    if f < 0 and f!= -999:
+    if f < 0 and f != -999:
         msg = 'float %s < 0 or not equal to -999' % f
         raise ValueError(msg)
     return f