Merge branch 'gem:master' into gsim_CB2003_World_BZP

.github/workflows/engine_pr_test.yml

@@ -70,6 +70,7 @@ jobs:
           ruff check --preview .
           pytest calculators -k 'risk or damage or bcr' -x --doctest-modules --disable-warnings --color=yes --durations=5
           pytest -xs --doctest-modules --disable-warnings --color=yes --durations=8 sep hmtk risklib commonlib baselib hazardlib
+          pytest --doctest-modules ../doc/contributing/*.rst
 
   server_public_mode:
     runs-on: ${{ matrix.os }}

doc/contributing/developing-with-the-engine.rst

@@ -262,7 +262,8 @@ Then the hazard curve can be computed as follows::
 	>>> sitecol = readinput.get_site_collection(oq)
 	>>> gsims = readinput.get_gsim_lt(oq).values['*']
 	>>> calc_hazard_curve(sitecol, src, gsims, oq)
-	[[0.00507997]]>
+	array([[0.00508004]], dtype=float32)
+
 
 Working with GMPEs directly: the ContextMaker
 ---------------------------------------------
@@ -345,6 +346,61 @@ the result provided by ``calc_hazard_curve(sitecol, src, gsims, oq)`` in the sec
 If you want to know exactly how ``get_pmap`` works you are invited to look at the source code in
 ``openquake.hazardlib.contexts``.
 
+Generating ground motion fields from a rupture
+----------------------------------------------
+
+The easiest way to create a finite size rupture (a.k.a. planar rupture)
+is to use the factory function `get_planar`:
+
+>>> from openquake.hazardlib.source.rupture import get_planar
+
+The function requires in input a site and a magnitude scaling relationship,
+so first you have to build such objects:
+
+>>> [site] = sitecol  # since there is a single site
+>>> from openquake.hazardlib.scalerel import WC1994
+>>> msr = WC1994()  # magnitude scaling relationship
+>>> mag = 6.
+>>> rup = get_planar(site, msr, mag, aratio=1., strike=11., dip=38.,
+...                  rake=55., trt=cmaker.trt)
+
+If you want to generate the GMF produced by a rupture (i.e. to emulate
+a scenario calculation) you need to supplement the number of
+occurrences of the rupture and a random seed, i.e. you need to convert the
+hazardlib rupture into an EBRupture:
+
+>>> from openquake.hazardlib.source.rupture import EBRupture
+>>> ebr = EBRupture(rup, n_occ=2, seed=42)
+
+Then you can use the GmfComputer class to perform the calculation:
+
+>>> from openquake.hazardlib.calc.gmf import GmfComputer
+>>> gc = GmfComputer(ebr, sitecol, cmaker)
+>>> gc.compute_all()  # returns a DataFrame
+      gmv_0  eid  sid  rlz
+0  0.660239    0    0    0
+1  0.301583    1    0    0
+
+`gmv_0` is the value of the ground motion field for the first IMT (i.e PGA in this
+case), `eid` the event ID, `sid` the site ID (there is a single site in this case)
+and `rlz` the realization index.
+In scenario calculations there is a realization for each GSIM and in this case
+there is a single GSIM, so rlz=0. The total number of events is the number
+of realizations times the number of occurrences and therefore in this case
+the event ID (`eid`) can only have the values 0 or 1.
+
+It is also possible to perform calculations with point-like ruptures
+(i.e. ignoring the finite-size effects):
+
+>>> from openquake.hazardlib.source.rupture import PointRupture
+>>> occ_rate = None  # not used in the GmfComputer
+>>> rup =  PointRupture(mag, cmaker.trt, site.location, occ_rate, cmaker.tom)
+>>> ebr = EBRupture(rup, n_occ=2, seed=42)
+>>> GmfComputer(ebr, sitecol, cmaker).compute_all()
+      gmv_0  eid  sid  rlz
+0  0.541180    0    0    0
+1  0.247199    1    0    0
+
 Working with verification tables
 --------------------------------
 

openquake/baselib/zeromq.py

@@ -96,7 +96,7 @@ class Socket(object):
     """
     # NB: the timeout has to be large since starting a workerpool can be
     # slow due to numba compiling everything, so you have to wait
-    def __init__(self, end_point, socket_type, mode, timeout=30):
+    def __init__(self, end_point, socket_type, mode, timeout=60):
         assert socket_type in (zmq.REP, zmq.REQ, zmq.PULL, zmq.PUSH)
         assert mode in ('bind', 'connect'), mode
         if mode == 'bind':

openquake/commonlib/oqvalidation.py

@@ -1472,12 +1472,12 @@ def set_loss_types(self):
         else:
             self._parent = None
         # set all_cost_types
-        # rt has the form 'vulnerability/structural', 'fragility/...', ...
-        costtypes = set(rt.rsplit('/')[1] for rt in self.risk_files)
+        # rt has the form 'earthquake/vulnerability/structural', ...
+        costtypes = set(rt.split('/')[2] for rt in self.risk_files)
         if not costtypes and self.hazard_calculation_id:
             try:
                 self._risk_files = rfs = get_risk_files(self._parent.inputs)
-                costtypes = set(rt.rsplit('/')[1] for rt in rfs)
+                costtypes = set(rt.split('/')[2] for rt in rfs)
             except OSError:  # FileNotFound for wrong hazard_calculation_id
                 pass
         self.all_cost_types = sorted(costtypes)  # including occupants

openquake/hazardlib/contexts.py

@@ -199,13 +199,18 @@ class Oq(object):
     mea_tau_phi = False
     split_sources = True
     use_rates = False
+    af = None
 
     def __init__(self, **hparams):
         vars(self).update(hparams)
 
     @property
     def min_iml(self):
-        return numpy.array([1E-10 for imt in self.imtls])
+        try:
+            imtls = self.imtls
+        except AttributeError:
+            imtls = self.hazard_imtls
+        return numpy.array([1E-10 for imt in imtls])
 
     def get_reqv(self):
         if 'reqv' not in self.inputs:
@@ -531,7 +536,10 @@ class ContextMaker(object):
 
     def __init__(self, trt, gsims, oq, monitor=Monitor(), extraparams=()):
         self.trt = trt
-        self.gsims = gsims
+        if isinstance(gsims, dict):
+            self.gsims = gsims
+        else:
+            self.gsims = {gsim: U32([i]) for i, gsim in enumerate(gsims)}
         # NB: the gid array can be overridden later on
         self.gid = numpy.arange(len(gsims), dtype=numpy.uint16)
         if isinstance(oq, dict):
@@ -595,16 +603,15 @@ def _init1(self, param):
         self.disagg_bin_edges = param.get('disagg_bin_edges', {})
         self.ps_grid_spacing = param.get('ps_grid_spacing')
         self.split_sources = self.oq.split_sources
-
-    def _init2(self, param, extraparams):
         for gsim in self.gsims:
             if hasattr(gsim, 'set_tables'):
                 if len(self.mags) == 0 and not is_modifiable(gsim):
                     raise ValueError(
                         'You must supply a list of magnitudes as 2-digit '
                         'strings, like mags=["6.00", "6.10", "6.20"]')
                 gsim.set_tables(self.mags, self.imtls)
-        self.effect = param.get('effect')
+
+    def _init2(self, param, extraparams):
         for req in self.REQUIRES:
             reqset = set()
             for gsim in self.gsims:
@@ -1221,7 +1228,9 @@ def get_4MN(self, ctxs, gsim):
     def get_att_curves(self, site, msr, mag, aratio=1., strike=0.,
                        dip=45., rake=-90):
         """
-        :returns: 4 attenuation curves mu, sig, tau, phi
+        :returns:
+            4 attenuation curves mea, sig, tau, phi
+            (up to 500 km from the site at steps of 5 km)
         """
         from openquake.hazardlib.source import rupture
         rup = rupture.get_planar(

openquake/hazardlib/source/point.py

@@ -305,8 +305,8 @@ def _gen_ruptures(self, shift_hypo=False, step=1, iruptures=False):
                     for d, (drate, cdep) in hdd_[::step]:
                         rate = mrate * nrate * drate
                         yield PointRupture(
-                            mag, np.rake, self.tectonic_region_type,
-                            Point(clon, clat, cdep), np.strike, np.dip, rate,
+                            mag, self.tectonic_region_type,
+                            Point(clon, clat, cdep), rate,
                             self.temporal_occurrence_model,
                             self.lower_seismogenic_depth)
 

openquake/hazardlib/source/rupture.py

@@ -583,19 +583,15 @@ def get_cdppvalue(self, target, buf=1.0, delta=0.01, space=2.):
 class PointSurface:
     """
     A fake surface used in PointRuptures.
-    The parameters `hypocenter` and `strike` are determined by
-    collapsing the corresponding parameters in the original PointSource.
     """
-    def __init__(self, hypocenter, strike, dip):
+    def __init__(self, hypocenter):
         self.hypocenter = hypocenter
-        self.strike = strike
-        self.dip = dip
 
     def get_strike(self):
-        return self.strike
+        return 0
 
     def get_dip(self):
-        return self.dip
+        return 0
 
     def get_top_edge_depth(self):
         return self.hypocenter.depth
@@ -628,17 +624,17 @@ class PointRupture(ParametricProbabilisticRupture):
     A rupture coming from a far away PointSource, so that the finite
     size effects can be neglected.
     """
-    def __init__(self, mag, rake, tectonic_region_type, hypocenter, strike,
-                 dip, occurrence_rate, temporal_occurrence_model, zbot):
+    def __init__(self, mag, tectonic_region_type, hypocenter,
+                 occurrence_rate, temporal_occurrence_model, zbot=0):
+        self.mag = mag
         self.tectonic_region_type = tectonic_region_type
         self.hypocenter = hypocenter
-        self.mag = mag
-        self.strike = strike
-        self.rake = rake
-        self.dip = dip
         self.occurrence_rate = occurrence_rate
         self.temporal_occurrence_model = temporal_occurrence_model
-        self.surface = PointSurface(hypocenter, strike, dip)
+        self.surface = PointSurface(hypocenter)
+        self.rake = 0
+        self.dip = 0
+        self.strike = 0
         self.zbot = zbot  # bottom edge depth, used in Campbell-Bozorgnia
 
 
@@ -719,7 +715,7 @@ class ExportedRupture(object):
     Simplified Rupture class with attributes rupid, events_by_ses, indices
     and others, used in export.
 
-    :param rupid: rupture.seed ID
+    :param rupid: rupture ID
     :param events_by_ses: dictionary ses_idx -> event records
     :param indices: site indices
     """
@@ -744,13 +740,14 @@ class EBRupture(object):
     """
     seed = 'NA'  # set by the engine
 
-    def __init__(self, rupture, source_id, trt_smr, n_occ=1, id=None, e0=0):
+    def __init__(self, rupture, source_id=0, trt_smr=0, n_occ=1, id=0, e0=0, seed=42):
         self.rupture = rupture
         self.source_id = source_id
         self.trt_smr = trt_smr
         self.n_occ = n_occ
         self.id = source_id * TWO30 + id
         self.e0 = e0
+        self.seed = seed
 
     @property
     def tectonic_region_type(self):

openquake/hazardlib/tests/calc/_conditioned_gmfs_test_data.py

@@ -46,11 +46,8 @@
 
 RUP = PointRupture(
     mag=MAG,
-    rake=0,
     tectonic_region_type=const.TRT.ACTIVE_SHALLOW_CRUST,
     hypocenter=Point(LON, LAT, DEP),
-    strike=0,
-    dip=90,
     occurrence_rate=1,
     temporal_occurrence_model=PoissonTOM(1.0),
     zbot=DEP,

openquake/hazardlib/tests/gsim/mgmpe/modifiable_gmpe_test.py

@@ -31,11 +31,11 @@
 class ModifiableGMPETest(unittest.TestCase):
 
     def test_AlAtik2015Sigma(self):
+        gmpe = valid.gsim('YenierAtkinson2015BSSA')
         params1 = {"tau_model": "global", "ergodic": False}
         params2 = {"tau_model": "cena", "ergodic": True}
         params3 = {}
-        gsims = [ModifiableGMPE(gmpe={'YenierAtkinson2015BSSA': {}},
-                                sigma_model_alatik2015=params)
+        gsims = [valid.modified_gsim(gmpe, sigma_model_alatik2015=params)
                  for params in [params1, params2, params3]]
         cmaker = simple_cmaker(gsims, ['PGA'])
         ctx = cmaker.new_ctx(4)
@@ -56,9 +56,10 @@ def test_AlAtik2015Sigma(self):
         aae(tau[2], 0.36855)
 
         # now test with_betw_ratio
-        cmaker.gsims[0] = ModifiableGMPE(
+        gsims[0] = ModifiableGMPE(
             gmpe={'Campbell2003': {}},
             add_between_within_stds={'with_betw_ratio': 0.6})
+        cmaker = simple_cmaker(gsims, ['PGA'])
         _mea, _sig, tau, phi = cmaker.get_mean_stds([ctx])  # (G,M,N)
         aae(tau[0], 0.44075136)
         aae(phi[0], 0.26445082)
@@ -70,9 +71,9 @@ def test_AlAtik2015Sigma(self):
 
     def test_AkkarEtAlRjb2014(self):
         # check mean and stds
-        gsims = [ModifiableGMPE(gmpe={'AkkarEtAlRjb2014': {}},
-                                set_between_epsilon={'epsilon_tau': 0.5}),
-                 valid.gsim('AkkarEtAlRjb2014')]
+        gmm = valid.gsim('AkkarEtAlRjb2014')
+        gsims = [valid.modified_gsim(
+            gmm, set_between_epsilon={'epsilon_tau': 0.5}), gmm]
         cmaker = simple_cmaker(gsims, ['PGA'])
         ctx = cmaker.new_ctx(4)
         ctx.mag = 6.
@@ -102,9 +103,9 @@ def test_coefficients_as_dictionary(self):
         self.assertAlmostEqual(output_coeffs["XYZ"][SA(3.0)]["XYZ"], 3.0)
 
     def get_mean_stds(self, **kw):
-        gmpe_name = 'AkkarEtAlRjb2014'
-        gmm1 = ModifiableGMPE(gmpe={gmpe_name: {}})
-        gmm2 = ModifiableGMPE(gmpe={gmpe_name: {}}, **kw)
+        gmm = valid.gsim('AkkarEtAlRjb2014')
+        gmm1 = valid.modified_gsim(gmm)
+        gmm2 = valid.modified_gsim(gmm, **kw)
         cmaker = simple_cmaker([gmm1, gmm2], ['PGA', 'SA(0.2)'])
         ctx = cmaker.new_ctx(4)
         ctx.mag = 6.
@@ -183,9 +184,8 @@ def test(self):
                           'BaumontEtAl2018High2210IAVGDC30n7',
                           'FaccioliCauzzi2006']:
 
-            gmm = ModifiableGMPE(gmpe={gmpe_name: {}},
-                                 apply_swiss_amplification={})
             gmpe = valid.gsim(gmpe_name)
+            gmm = valid.modified_gsim(gmpe, apply_swiss_amplification={})
             cmaker = simple_cmaker([gmm, gmpe], ['MMI'])
             ctx = self.get_ctx(cmaker)
             mea, _sig, _tau, phi = cmaker.get_mean_stds([ctx])
@@ -198,9 +198,8 @@ def test(self):
                           'EdwardsFah2013Foreland60Bars',
                           'ChiouYoungs2008SWISS01']:
 
-            gmm = ModifiableGMPE(gmpe={gmpe_name: {}},
-                                 apply_swiss_amplification_sa={})
             gmpe = valid.gsim(gmpe_name)
+            gmm = valid.modified_gsim(gmpe,apply_swiss_amplification_sa={})
             cmaker = simple_cmaker([gmm, gmpe], ['SA(0.3)'])
             ctx = self.get_ctx(cmaker)
             mea, _sig, _tau, phi = cmaker.get_mean_stds([ctx])