Merge pull request #10093 from gem/read_sections

Added functionality to read multi fault sources

Author: micheles

debian/changelog

@@ -1,3 +1,6 @@
+  [Michele Simionato]
+  * Internal: added utility function readinput.read_source_models
+
   [Manuela Villani]
   * Improved the "Governing MCE" plot
 

openquake/calculators/preclassical.py

@@ -26,7 +26,6 @@
 from openquake.hazardlib import pmf, geo, source_reader
 from openquake.baselib.general import AccumDict, groupby, block_splitter
 from openquake.hazardlib.contexts import read_cmakers
-from openquake.hazardlib.geo.surface.multi import build_secparams
 from openquake.hazardlib.source.point import grid_point_sources
 from openquake.hazardlib.source.base import get_code2cls
 from openquake.hazardlib.sourceconverter import SourceGroup
@@ -278,9 +277,8 @@ def populate_csm(self):
             else:
                 normal_sources.extend(sg)
         if multifaults:
-            # this is ultra-fast
-            sections = multifaults[0].get_sections()
-            secparams = build_secparams(sections)
+            with hdf5.File(multifaults[0].hdf5path, 'r') as h5:
+                secparams = h5['secparams'][:]
             logging.warning(
                 'There are %d multiFaultSources (secparams=%s)',
                 len(multifaults), general.humansize(secparams.nbytes))

openquake/commonlib/readinput.py

@@ -1638,3 +1638,24 @@ def read_countries_df(buffer=0.1):
     fname = os.path.join(os.path.dirname(global_risk.__file__),
                          'geoBoundariesCGAZ_ADM0.shp')
     return read_geometries(fname, 'shapeGroup', buffer)
+
+
+def read_source_models(fnames, hdf5path='', **converterparams):
+    """
+    :param fnames: a list of source model files
+    :param hdf5path: auxiliary .hdf5 file used to store the multifault sources
+    :param converterparams: a dictionary of parameters like rupture_mesh_spacing
+    :returns: a list of SourceModel instances
+    """
+    converter = sourceconverter.SourceConverter()
+    vars(converter).update(converterparams)
+    smodels = list(nrml.read_source_models(fnames, converter))
+    smdict = dict(zip(fnames, smodels))
+    src_groups = [sg for sm in smdict.values() for sg in sm.src_groups]
+    secparams = source_reader.fix_geometry_sections(smdict, src_groups, hdf5path)
+    for smodel in smodels:
+        for sg in smodel.src_groups:
+            for src in sg:
+                if src.code == b'F':  # multifault
+                    src.set_msparams(secparams)
+    return smodels

openquake/commonlib/tests/readinput_test.py

@@ -31,7 +31,7 @@
 from openquake.risklib import asset
 from openquake.commonlib import readinput, datastore
 from openquake.qa_tests_data.logictree import case_02, case_15, case_21
-from openquake.qa_tests_data.classical import case_34
+from openquake.qa_tests_data.classical import case_34, case_65
 from openquake.qa_tests_data.event_based import case_16
 from openquake.qa_tests_data.event_based_risk import case_2, case_caracas
 from openquake.qa_tests_data import mosaic
@@ -564,3 +564,18 @@ def test_read_station_data(self):
             readinput.get_station_data(oq, sitecol)
         self.assertIn("Stations_NIED.csv: has duplicate sites ['GIF001', 'GIF013']",
                       str(ctx.exception))
+
+
+class ReadSourceModelsTestCase(unittest.TestCase):
+    def test(self):
+        base = os.path.dirname(case_65.__file__)
+        hdf5path = general.gettemp(suffix='.hdf5')
+        fnames = [os.path.join(base, 'ssm.xml'), os.path.join(base, 'sections.xml')]
+        smodels = readinput.read_source_models(fnames, hdf5path, investigation_time=1.)
+        nrups = 0
+        for smodel in smodels:
+            for sg in smodel.src_groups:
+                for src in sg:
+                    for rup in src.iter_ruptures():
+                        nrups += 1
+        self.assertEqual(nrups, 3)

openquake/hazardlib/contexts.py

@@ -198,6 +198,7 @@ class Oq(object):
     """
     mea_tau_phi = False
     split_sources = True
+    use_rates = False
 
     def __init__(self, **hparams):
         vars(self).update(hparams)

openquake/hazardlib/source/multi_fault.py

@@ -31,7 +31,8 @@
     NonParametricSeismicSource as NP)
 from openquake.hazardlib.geo.surface.kite_fault import (
     geom_to_kite, kite_to_geom)
-from openquake.hazardlib.geo.surface.multi import MultiSurface, build_msparams
+from openquake.hazardlib.geo.surface.multi import (
+    MultiSurface, build_msparams, build_secparams)
 from openquake.hazardlib.geo.utils import (
     angular_distance, KM_TO_DEGREES, get_spherical_bounding_box)
 from openquake.hazardlib.source.base import BaseSeismicSource
@@ -335,6 +336,7 @@ def save_and_split(mfsources, sectiondict, hdf5path, site1=None,
                 split_dic[src.source_id].append(split)
         h5.save_vlen('multi_fault_sections',
                      [kite_to_geom(sec) for sec in sectiondict.values()])
+        h5['secparams'] = build_secparams(src.get_sections())
 
     return split_dic
 
@@ -346,7 +348,7 @@ def load(hdf5path):
     srcs = []
     with hdf5.File(hdf5path, 'r') as h5:
         for key in list(h5):
-            if key == 'multi_fault_sections':
+            if key in ('multi_fault_sections', 'secparams'):
                 continue
             data = h5[key]
             name = data.attrs['name']

openquake/hazardlib/source_reader.py

@@ -310,7 +310,13 @@ def get_csm(oq, full_lt, dstore=None):
                            fillvalue=None)
 
     # must be called *after* _fix_dupl_ids
-    fix_geometry_sections(smdict, csm, dstore)
+    if oq.sites and len(oq.sites) == 1 and oq.use_rates:
+        lon, lat, _dep = oq.sites[0]
+        site1 = site.SiteCollection.from_points([lon], [lat])
+    else:
+        site1 = None
+    hdf5path = dstore.tempname if dstore else ''
+    fix_geometry_sections(smdict, csm.src_groups, hdf5path, site1)
     return csm
 
 
@@ -377,7 +383,7 @@ def replace(lst, splitdic, key):
     lst[:] = new
 
 
-def fix_geometry_sections(smdict, csm, dstore):
+def fix_geometry_sections(smdict, src_groups, hdf5path='', site1=None):
     """
     If there are MultiFaultSources, fix the sections according to the
     GeometryModels (if any).
@@ -399,22 +405,20 @@ def fix_geometry_sections(smdict, csm, dstore):
 
     if sections:
         # save in the temporary file
-        assert dstore, ('You forgot to pass the dstore to '
+        assert hdf5path, ('You forgot to pass the dstore to '
                         'get_composite_source_model')
-        oq = dstore['oqparam']
+
         mfsources = []
-        for sg in csm.src_groups:
+        for sg in src_groups:
             for src in sg:
                 if src.code == b'F':
                     mfsources.append(src)
-        if oq.sites and len(oq.sites) == 1 and oq.use_rates:
-            lon, lat, _dep = oq.sites[0]
-            site1 = site.SiteCollection.from_points([lon], [lat])
-        else:
-            site1 = None
-        split_dic = save_and_split(mfsources, sections, dstore.tempname, site1)
-        for sg in csm.src_groups:
+        split_dic = save_and_split(mfsources, sections, hdf5path, site1)
+        for sg in src_groups:
             replace(sg.sources, split_dic, 'source_id')
+        with hdf5.File(hdf5path, 'r') as h5:
+            return h5['secparams'][:]
+    return ()
 
 
 def _groups_ids(smlt_dir, smdict, fnames):