Merge pull request #10115 from gem/doc

Documented the GmfComputer

.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/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/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])

openquake/risklib/scientific.py

@@ -45,6 +45,7 @@
                       'losses', 'collapsed',
                       'injured', 'fatalities', 'homeless', 'non_operational']
 
+PERILS = 'earthquake', 'liquefaction', 'landslide'
 LOSSTYPE = numpy.array('''\
 business_interruption contents nonstructural structural
 occupants occupants_day occupants_night occupants_transit