diff --git a/openquake/commands/mosaic.py b/openquake/commands/mosaic.py index 47f2eb5a770..bae618a1894 100644 --- a/openquake/commands/mosaic.py +++ b/openquake/commands/mosaic.py @@ -300,7 +300,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'): diff --git a/openquake/commonlib/oqvalidation.py b/openquake/commonlib/oqvalidation.py index 66cfc7cf341..afe2b1de558 100644 --- a/openquake/commonlib/oqvalidation.py +++ b/openquake/commonlib/oqvalidation.py @@ -971,7 +971,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': @@ -1351,7 +1351,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): @@ -1401,6 +1402,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() @@ -1564,7 +1570,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'): @@ -1594,7 +1601,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 @@ -1666,21 +1673,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: @@ -1692,25 +1700,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) @@ -1732,6 +1740,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): """ @@ -2621,7 +2630,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: diff --git a/openquake/engine/impact.py b/openquake/engine/impact.py index 4ca8752e061..cee6138b828 100644 --- a/openquake/engine/impact.py +++ b/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='', diff --git a/openquake/hazardlib/calc/filters.py b/openquake/hazardlib/calc/filters.py index 66bada58933..7479512b8d0 100644 --- a/openquake/hazardlib/calc/filters.py +++ b/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: diff --git a/openquake/hazardlib/shakemap/validate.py b/openquake/hazardlib/shakemap/validate.py index 95dcf8d1c83..33afb48e5a8 100644 --- a/openquake/hazardlib/shakemap/validate.py +++ b/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__) @@ -37,7 +38,6 @@ class AristotleParam: rupture_dict: dict time_event: str - maximum_distance: float truncation_level: float number_of_ground_motion_fields: int asset_hazard_distance: float @@ -46,6 +46,7 @@ class AristotleParam: 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 @@ -83,7 +84,6 @@ def get_oqparams(self, usgs_id, mosaic_models, trts, use_shakemap): calculation_mode='scenario_risk', rupture_dict=str(rupdic), time_event=self.time_event, - maximum_distance=str(self.maximum_distance), mosaic_model=self.mosaic_model, tectonic_region_type=self.trt, truncation_level=str(self.truncation_level), @@ -99,6 +99,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 +112,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 +203,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,14 +223,16 @@ 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', } -msr_choices = [msr.__class__.__name__ for msr in get_available_magnitude_scalerel()] +msr_choices = [msr.__class__.__name__ + for msr in get_available_magnitude_scalerel()] validators = { 'approach': valid.Choice(*IMPACT_APPROACHES), @@ -258,16 +265,18 @@ def _validate(POST): validation_errs = {} invalid_inputs = [] params = {} - inputdic = dict(approach=None, usgs_id=None, lon=None, lat=None, dep=None, - mag=None, msr=None, aspect_ratio=None, rake=None, dip=None, - strike=None, description=None) + inputdic = dict( + approach=None, usgs_id=None, lon=None, lat=None, dep=None, + mag=None, msr=None, aspect_ratio=None, rake=None, dip=None, + strike=None, description=None) for field, validation_func in validators.items(): if field not in POST: continue 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: @@ -364,7 +373,7 @@ def impact_validate(POST, user, rupture_file=None, station_data_file=None, mosaic_models = readinput.get_close_mosaic_models( rupdic['lon'], rupdic['lat'], 5) except ValueError as exc: - # e.g. '(-139.0, 35.0) is farther than 5 deg from any mosaic model!' + # e.g. '(-139.0, 35.0) is farther than 5 deg from any mosaic model' err = {"status": "failed", "error_msg": str(exc)} return rup, rupdic, params, err for mosaic_model in mosaic_models: diff --git a/openquake/hazardlib/valid.py b/openquake/hazardlib/valid.py index 8c7ea8f26f9..c5d9d1d17cf 100644 --- a/openquake/hazardlib/valid.py +++ b/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