Initial implementation of logic trees for site models

doc/api-reference/openquake.hazardlib.rst

@@ -111,6 +111,14 @@ site
     :undoc-members:
     :show-inheritance:
 
+site_lt
+-------------------------------
+
+.. automodule:: openquake.hazardlib.site_lt
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
 sourceconverter
 ------------------------------------------
 

doc/user-guide/inputs/site-model-inputs.rst

@@ -127,5 +127,39 @@ We used this feature to split the ESHM20 model in two parts (Northern Europe and
 full hazard map was as trivial as joining the generated CSV files. Without the ``custom_site_id`` the site IDs would 
 overlap, thus making impossible to join the outputs.
 
-A geohash string (see https://en.wikipedia.org/wiki/Geohash) makes a good ``custom_site_id`` since it can enable the 
+A geohash string (see https://en.wikipedia.org/wiki/Geohash) makes a good ``custom_site_id`` since it can enable the
 unique identification of all potential sites across the globe.
+
+Site model logic tree
+---------------------
+
+The ``site_model_file`` can point to a NRML logic tree XML declaring alternative site models with weights, instead of
+a single site model file. This adds a third leg to the SSC × GMM logic tree; under full enumeration realizations become
+``R_SSC × R_GMM × R_SITE`` while under sampling all three legs are Monte-Carlo sampled ``num_samples`` times. Each
+realization path gains a third ``~``-separated branch ID (e.g. ``A~A~B``). All branches must reference the same sites
+(identical ``lon``/``lat`` and ``depth`` if present, in the same order, and identical field sets); only per-site parameter
+values (``vs30``, ``z1pt0``, ``z2pt5``, ...) may differ. Branch files may be CSV or NRML ``<siteModel>`` XML. *Site model
+logic trees are currently only supported in classical and disaggregation.*
+
+Example logic tree XML::
+
+    <?xml version="1.0" encoding="UTF-8"?>
+    <nrml xmlns="http://openquake.org/xmlns/nrml/0.5">
+      <logicTree logicTreeID="lt_site">
+        <logicTreeBranchSet uncertaintyType="siteModel" branchSetID="bs_site">
+          <logicTreeBranch branchID="rock">
+            <uncertaintyModel>site_model_rock.csv</uncertaintyModel>
+            <uncertaintyWeight>0.6</uncertaintyWeight>
+          </logicTreeBranch>
+          <logicTreeBranch branchID="soil">
+            <uncertaintyModel>site_model_soil.csv</uncertaintyModel>
+            <uncertaintyWeight>0.4</uncertaintyWeight>
+          </logicTreeBranch>
+        </logicTreeBranchSet>
+      </logicTree>
+    </nrml>
+
+Referenced from ``job.ini`` as::
+
+    [geometry]
+    site_model_file = site_model_lt.xml

openquake/calculators/base.py

@@ -157,8 +157,10 @@ def get_weights(oq, dstore):
     :returns: float32 array of realization weights
     """
     samples = oq.number_of_logic_tree_samples
-    if samples:
-        weights = numpy.ones(samples, dtype=F32)/samples
+    # Under a site-model LT weights may be non-uniform (late_weights),
+    # so always read them from the datastore in that case
+    if samples and 'full_lt/site_model_lt' not in dstore:
+        weights = numpy.ones(samples, dtype=F32) / samples
     else:
         weights = dstore['weights'][:]
     return weights
@@ -844,6 +846,24 @@ def R(self):
             return 1
         return len(get_weights(self.oqparam, self.datastore))
 
+    def _overlay_sitecol(self, arr):
+        """
+        Overlay ``arr``'s per-site params on ``self.sitecol`` and
+        ``dstore['sitecol/*']``; ``arr`` is a structured array with the
+        same rows as the sitecol (either a per-branch site model or a
+        copy used to restore it)
+        """
+        # Geometry fields are shared across branches - never overlay
+        skip = {'lon', 'lat', 'depth', 'sids'}
+        h5 = self.datastore.hdf5
+        for name in arr.dtype.names:
+            if name in skip:
+                continue
+            self.sitecol.array[name] = arr[name]
+            key = 'sitecol/' + name
+            if key in h5:
+                h5[key][:] = arr[name]
+
     def read_exposure(self, haz_sitecol):  # after load_crmodel
         """
         Read the exposure, the risk models and update the attributes

openquake/calculators/classical.py

@@ -36,7 +36,7 @@
 from openquake.hazardlib.calc import disagg
 from openquake.hazardlib.map_array import (
     RateMap, MapArray, rates_dt, check_hmaps, gen_chunks)
-from openquake.commonlib import calc
+from openquake.commonlib import calc, readinput
 from openquake.calculators import base, getters, preclassical, views
 
 get_weight = operator.attrgetter('weight')
@@ -303,6 +303,10 @@ def postclassical(pgetter, hstats, individual_rlzs, amplifier, monitor):
                 continue
             if R == 1 or individual_rlzs:
                 for r in range(R):
+                    # Under site LT each pgetter owns only its rlz_mask rlzs
+                    if (pgetter.rlz_mask is not None
+                            and not pgetter.rlz_mask[r]):
+                        continue
                     pmap_by_kind['hcurves-rlzs'][r].array[idx] = (
                         pc[:, r].reshape(M, L1))
             if hstats:
@@ -406,8 +410,10 @@ def agg_dicts(self, acc, dic):
         self.dparam_mb = max(dic.pop('dparam_mb'), self.dparam_mb)
         self.source_mb = max(dic.pop('source_mb'), self.source_mb)
 
-        # store rup_data if there are few sites
-        if self.few_sites and len(dic['rup_data']):
+        # Store rup_data if there are few sites; skip on site-LT re-runs
+        # since rupture data is source-only and would duplicate in rup/*
+        if (self.few_sites and len(dic['rup_data'])
+                and not getattr(self, '_skip_store_ctxs', False)):
             with self.monitor('saving rup_data'):
                 store_ctxs(self.datastore, dic['rup_data'], grp_id)
 
@@ -547,7 +553,10 @@ def execute(self):
                               oq.inputs)
         self.source_data = AccumDict(accum=[])
         sgs, ds = self._pre_execute()
-        self._execute(sgs, ds)
+        if getattr(self.full_lt, 'site_model_lt', None) is not None:
+            self._execute_epistemic_site(sgs, ds)
+        else:
+            self._execute(sgs, ds)
         if self.cfactor[0] == 0:
             if self.N == 1:
                 logging.error('The site is far from all seismic sources'
@@ -692,6 +701,48 @@ def check_mean_rates(self, mean_rates_by_src):
             ok = got[m] < 2.
             numpy.testing.assert_allclose(got[m, ok], exp[m, ok], atol=1E-5)
 
+    def _execute_epistemic_site(self, sgs, ds):
+        """
+        Run :meth:`_execute` once per site-model realization; each branch's
+        rates are stored in a separate ``_rates_site_i`` group
+        """
+        oq = self.oqparam
+        smep = readinput.get_site_models_epistemic(oq)
+        # Baseline site params from the canonical (first) branch
+        used = {r.site_rlz.ordinal
+                for r in self.full_lt.get_realizations()
+                if r.site_rlz is not None}
+        baseline = self.sitecol.array.copy()
+        for i, arr in enumerate(smep.arrays):
+            if i not in used:
+                continue
+            self._overlay_sitecol(arr)
+            self.rmap = {}
+            self.cfactor = numpy.zeros(2)
+            self.rel_ruptures = AccumDict(accum=0)
+            # Purge _rates so _execute starts clean; rup is source-only
+            # and is written once (see _skip_store_ctxs below)
+            h5 = self.datastore.hdf5
+            for key in ('_rates/slice_by_idx', '_rates/sid', '_rates/lid',
+                        '_rates/gid', '_rates/rate', '_rates', 'grp_keys'):
+                if key in h5:
+                    del h5[key]
+            self.datastore.create_df(
+                '_rates', [(n, rates_dt[n]) for n in rates_dt.names], GZIP)
+            self.datastore.create_dset(
+                '_rates/slice_by_idx', getters.slice_dt)
+            self._skip_store_ctxs = i > 0
+            self._execute(sgs, ds)
+            # Drop SWMR before renaming _rates
+            self.datastore.close()
+            self.datastore.open('a')
+            h5 = self.datastore.hdf5
+            target = '_rates_site_%d' % i
+            if target in h5:
+                del h5[target]
+            h5.move('_rates', target)
+        self._overlay_sitecol(baseline)
+
     def store_info(self):
         """
         Store full_lt, source_info and source_data
@@ -733,16 +784,22 @@ def collect_hazard(self, acc, pmap_by_kind):
         # this is practically instantaneous
         if pmap_by_kind is None:  # instead of a dict
             raise MemoryError('You ran out of memory!')
+        # Under site-model LT, sum across getters (each covers disjoint rlzs)
+        site_lt = getattr(self.full_lt, 'site_model_lt', None) is not None
         for kind in pmap_by_kind:  # hmaps-XXX, hcurves-XXX
             pmaps = pmap_by_kind[kind]
-            if kind in self.hazard:
-                array = self.hazard[kind]
-            else:
+            accum = site_lt and kind in self.hazard
+            if kind not in self.hazard:
                 dset = self.datastore.getitem(kind)
-                array = self.hazard[kind] = numpy.zeros(dset.shape, dset.dtype)
+                self.hazard[kind] = numpy.zeros(dset.shape, dset.dtype)
+            array = self.hazard[kind]
             for r, pmap in enumerate(pmaps):
                 for idx, sid in enumerate(pmap.sids):
-                    array[sid, r] = pmap.array[idx]  # shape (M, P)
+                    val = pmap.array[idx]  # shape (M, P) or (M, L1)
+                    if accum:
+                        array[sid, r] += val
+                    else:
+                        array[sid, r] = val
 
     def post_execute(self, dummy):
         """
@@ -806,7 +863,9 @@ def build_curves_maps(self):
         oq = self.oqparam
         hstats = oq.hazard_stats()
         N, S, M, P, L1 = self._create_hcurves_maps()
-        if '_rates' in set(self.datastore) or not self.datastore.parent:
+        has_rates = any(k == '_rates' or k.startswith('_rates_site_')
+                        for k in self.datastore)
+        if has_rates or not self.datastore.parent:
             dstore = self.datastore
         else:
             dstore = self.datastore.parent

openquake/calculators/disaggregation.py

@@ -31,7 +31,7 @@
 from openquake.hazardlib import stats, map_array, valid
 from openquake.hazardlib.calc import disagg, mean_rates
 from openquake.hazardlib.contexts import read_cmakers, read_ctx_by_grp
-from openquake.commonlib import util
+from openquake.commonlib import util, readinput
 from openquake.calculators import base, getters
 
 POE_TOO_BIG = '''\
@@ -47,7 +47,7 @@
 
 
 def compute_disagg(dstore, ctxt, sitecol, cmaker, bin_edges, src_mutex, rwdic,
-                   monitor):
+                   rlz_filter, monitor):
     """
     :param dstore:
         a DataStore instance
@@ -63,6 +63,9 @@ def compute_disagg(dstore, ctxt, sitecol, cmaker, bin_edges, src_mutex, rwdic,
         a dictionary src_id -> weight, usually empty
     :param rwdic:
         dictionary rlz -> weight, empty for individual realizations
+    :param rlz_filter:
+        set of rlz ordinals to keep, or ``None`` to disable filtering
+        (used by the site-model LT outer loop in :meth:`compute`)
     :param monitor:
         monitor of the currently running job
     :returns:
@@ -85,6 +88,10 @@ def compute_disagg(dstore, ctxt, sitecol, cmaker, bin_edges, src_mutex, rwdic,
 
             imtls = {imt: iml2[m] for m, imt in enumerate(cmaker.imts)}
             rlzs = dstore['best_rlzs'][dis.sid]
+        if rlz_filter is not None:
+            rlzs = numpy.array([r for r in rlzs if r in rlz_filter])
+            if len(rlzs) == 0:
+                continue
         res = dis.disagg_by_magi(imtls, rlzs, rwdic, src_mutex,
                                  mon0, mon1, mon2, mon3)
         out.extend(res)
@@ -112,12 +119,14 @@ def output_dict(shapedic, disagg_outputs, Z):
     return dic
 
 
-def submit(smap, dstore, ctxt, sitecol, cmaker, bin_edges, src_mutex, rwdic):
+def submit(smap, dstore, ctxt, sitecol, cmaker, bin_edges, src_mutex, rwdic,
+           rlz_filter=None):
     mags = list(numpy.unique(ctxt.mag))
     logging.debug('Sending %d/%d sites for grp_id=%d, mags=%s',
                   len(sitecol), len(sitecol.complete), ctxt.grp_id[0],
                   shortlist(mags))
-    smap.submit((dstore, ctxt, sitecol, cmaker, bin_edges, src_mutex, rwdic))
+    smap.submit((dstore, ctxt, sitecol, cmaker, bin_edges, src_mutex, rwdic,
+                 rlz_filter))
 
 
 def check_memory(N, Z, shape8D):
@@ -180,7 +189,8 @@ def full_disaggregation(self):
         try:
             full_lt = self.full_lt
         except AttributeError:
-            full_lt = self.datastore['full_lt'].init()
+            # Skip .init() here: map_getters will call it below
+            full_lt = self.datastore['full_lt']
         if oq.rlz_index is None and oq.num_rlzs_disagg == 0:
             oq.num_rlzs_disagg = self.R  # 0 means all rlzs
         self.oqparam.mags_by_trt = self.datastore['source_mags']
@@ -244,7 +254,38 @@ def full_disaggregation(self):
                      for z, r in enumerate(rlzs[s])}
         return self.compute()
 
-    def _submit_all(self, smap, cmakers, ctx_by_grp, src_mutex_by_grp):
+    def compute(self):
+        """
+        Submit disaggregation tasks and return the results
+        """
+        oq = self.oqparam
+        full_lt = getattr(self, 'full_lt', None) or self.datastore['full_lt']
+        site_lt = getattr(full_lt, 'site_model_lt', None)
+        if site_lt is None:
+            return self._compute_pass(rlz_filter=None)
+        # Site-model epistemic uncertainty: one pass per site branch
+        self.full_lt = full_lt
+        smep = readinput.get_site_models_epistemic(oq)
+        baseline = self.sitecol.array.copy()
+        all_rlzs = self.full_lt.get_realizations()
+        combined = None
+        for i, arr in enumerate(smep.arrays):
+            rlz_filter = {r.ordinal for r in all_rlzs
+                          if r.site_rlz.ordinal == i}
+            if not rlz_filter:
+                continue
+            self._overlay_sitecol(arr)
+            part = self._compute_pass(rlz_filter=rlz_filter)
+            if combined is None:
+                combined = part
+            else:
+                for k, v in part.items():
+                    combined[k] = combined.get(k, 0) + v
+        self._overlay_sitecol(baseline)
+        return combined
+
+    def _submit_all(self, smap, cmakers, ctx_by_grp, src_mutex_by_grp,
+                    rlz_filter):
         # compute the total weight of the contexts and the maxsize
         ct = self.oqparam.concurrent_tasks or 1
         totweight = sum(cmakers[grp_id].Z * len(ctx)
@@ -278,7 +319,8 @@ def _submit_all(self, smap, cmakers, ctx_by_grp, src_mutex_by_grp):
             if ntasks < 1 or len(src_mutex) or rup_mutex:
                 # do not split (test case_11)
                 submit(smap, self.datastore, ctxt, self.sitecol, cmaker,
-                       self.bin_edges, src_mutex, rwdic)
+                       self.bin_edges, src_mutex, rwdic,
+                       rlz_filter=rlz_filter)
                 continue
 
             # split by tiles
@@ -290,15 +332,19 @@ def _submit_all(self, smap, cmakers, ctx_by_grp, src_mutex_by_grp):
                     for c in disagg.split_by_magbin(
                             ctx, self.bin_edges[0]).values():
                         submit(smap, self.datastore, c, tile, cmaker,
-                               self.bin_edges, src_mutex, rwdic)
+                               self.bin_edges, src_mutex, rwdic,
+                               rlz_filter=rlz_filter)
                 elif len(ctx):
                     # see case_multi in the oq-risk-tests
                     submit(smap, self.datastore, ctx, tile, cmaker,
-                           self.bin_edges, src_mutex, rwdic)
+                           self.bin_edges, src_mutex, rwdic,
+                           rlz_filter=rlz_filter)
 
-    def compute(self):
+    def _compute_pass(self, rlz_filter):
         """
-        Submit disaggregation tasks and return the results
+        Single pass of the disaggregation compute loop; used both by
+        the regular path (``rlz_filter=None``) and by the site-model
+        epistemic outer loop (one pass per site rlz)
         """
         dstore = (self.datastore.parent if self.datastore.parent
                   else self.datastore)
@@ -326,7 +372,8 @@ def compute(self):
         # that would break the ordering of the indices causing an incredibly
         # worse performance, but visible only in extra-large calculations!
 
-        self._submit_all(smap, cmakers, ctx_by_grp, src_mutex_by_grp)
+        self._submit_all(smap, cmakers, ctx_by_grp, src_mutex_by_grp,
+                         rlz_filter)
         s = self.shapedic
         shape8D = (s['trt'], s['mag'], s['dist'], s['lon'], s['lat'], s['eps'],
                    s['M'], s['P'])

openquake/calculators/export/hazard.py

@@ -309,9 +309,12 @@ def get_metadata(rlzs, kind):
     """
     metadata = {}
     if kind.startswith('rlz-'):
-        smlt_path, gslt_path = rlzs[int(kind[4:])]['branch_path'].split('~')
-        metadata['smlt_path'] = smlt_path
-        metadata['gsimlt_path'] = gslt_path
+        # SSC~GMM or SSC~GMM~SITE when a site-model LT is active
+        parts = rlzs[int(kind[4:])]['branch_path'].split('~')
+        metadata['smlt_path'] = parts[0]
+        metadata['gsimlt_path'] = parts[1]
+        if len(parts) > 2:
+            metadata['site_lt_path'] = parts[2]
     elif kind.startswith('quantile-'):
         metadata['statistics'] = 'quantile'
         metadata['quantile_value'] = float(kind[9:])