Merge branch 'master' of https://github.com/gem/oq-engine into us_202…

…3_basin_terms

.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 }}

debian/changelog

@@ -1,4 +1,6 @@
   [Michele Simionato]
+  * Internal: changed the ordering in the composite risk model from
+    (loss_type, riskid) -> (riskid, loss_type)
   * Added an exporter for trt_gsim
   * Internal: added utility function readinput.read_source_models
 

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
 --------------------------------
 
@@ -399,10 +455,8 @@ Running the engine tests
 ------------------------
 
 If you are a hazard scientist contributing a bug fix to a GMPE (or any other kind of bug fix) you may need to run the
-engine tests and possibly change the expected files if there is a change in the numbers. The way to do it is to start
-the dbserver and then run the tests from the repository root::
+engine tests and possibly change the expected files if there is a change in the numbers. The way to do it is to give the following command from the repository root::
 
-	$ oq dbserver start
 	$ pytest -vx openquake/calculators
 
 If you get an error like the following::

docker/scripts/oq-start.sh

@@ -19,7 +19,7 @@
 # This is required to load a custom  local_settings.py when 'oq webui' is used.
 export PYTHONPATH=$HOME
 
-oq engine --upgrade-db &
+oq dbserver start &
 
 # Wait the DbServer to come up; may be replaced with a "oq dbserver wait"
 while :

openquake/baselib/general.py

@@ -45,7 +45,7 @@
 import numpy
 import pandas
 from decorator import decorator
-from openquake.baselib import __version__
+from openquake.baselib import __version__, config
 from openquake.baselib.python3compat import decode
 
 U8 = numpy.uint8
@@ -436,7 +436,8 @@ def gettemp(content=None, dir=None, prefix="tmp", suffix="tmp", remove=True):
     if dir is not None:
         if not os.path.exists(dir):
             os.makedirs(dir)
-    fh, path = tempfile.mkstemp(dir=dir, prefix=prefix, suffix=suffix)
+    fh, path = tempfile.mkstemp(dir=dir or config.directory.custom_tmp,
+                                prefix=prefix, suffix=suffix)
     if remove:
         _tmp_paths.append(path)
     with os.fdopen(fh, "wb") as fh:

openquake/calculators/multi_risk.py

@@ -58,7 +58,7 @@ def get_dmg_csq(crm, assets_by_site, gmf, time_event):
                         for k, w in zip(df.risk_id, df.weight)]
                 # NB: risk logic trees are not yet supported in multi_risk
                 [peril] = rm.imt_by_lt.values()
-                fracs[li] = rm.scenario_damage(loss_type, assets,
+                fracs[li] = rm.scenario_damage('earthquake', loss_type, assets,
                                                pandas.DataFrame({peril: [gmv]}))
             csq = crm.compute_csq(assets, fracs, df, crm.oqparam)
             number = assets['value-number']

openquake/commands/compare.py

@@ -469,8 +469,8 @@ def compare_oqparam(calc_ids: int):
     ds1 = datastore.read(calc_ids[1])
     dic0 = vars(ds0['oqparam'])
     dic1 = vars(ds1['oqparam'])
-    common = sorted(set(dic0) & set(dic1))
-    for key in common:
+    common = set(dic0) & set(dic1) - {'hdf5path'}
+    for key in sorted(common):
         if dic0[key] != dic1[key]:
             print('%s: %s != %s' % (key, dic0[key], dic1[key]))
 

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/engine/tests/aristotle1/aggrisk.txt

@@ -1,6 +1,6 @@
 | gsim              | weight | contents    | nonstructural | structural | occupants | area   | number   | residents |
 |-------------------+--------+-------------+---------------+------------+-----------+--------+----------+-----------|
-| [ChiouYoungs2014] | 0.3500 | 47_928_608  | 86_326_184    | 49_106_732 | 1.3142    | 31_790 | 107.8218 | 1_063     |
-| [BooreEtAl2014]   | 0.3500 | 71_780_840  | 76_948_504    | 55_768_280 | 0.4272    | 43_355 | 250.2    | 1_185     |
-| [ZhaoEtAl2016Asc] | 0.3000 | 123_750_848 | 88_760_720    | 40_181_644 | 0.2262    | 61_196 | 349.9    | 1_304     |
-| Average           | 1.0000 | 79_023_561  | 83_774_357    | 48_760_747 | 0.6773    | 44_660 | 230.3    | 1_178     |
+| [ChiouYoungs2014] | 0.3500 | 68_505_824  | 77_741_872    | 48_760_420 | 1.3142    | 31_790 | 107.8218 | 1_063     |
+| [BooreEtAl2014]   | 0.3500 | 81_108_936  | 104_466_288   | 45_423_496 | 0.4272    | 43_355 | 250.2    | 1_185     |
+| [ZhaoEtAl2016Asc] | 0.3000 | 105_540_112 | 68_223_000    | 39_220_012 | 0.2262    | 61_196 | 349.9    | 1_304     |
+| Average           | 1.0000 | 84_027_200  | 84_239_756    | 44_730_374 | 0.6773    | 44_660 | 230.3    | 1_178     |

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,