ARISTOTLE: add to level 2 interface the possibility to choose among different workflows

bin/arist

@@ -3,4 +3,6 @@ import sys
 from openquake.engine.aristotle import main_cmd
 if __name__ == '__main__':
     # us7000k9f0
-    main_cmd(sys.argv[1], maximum_distance='200', loglevel='info', userlevel=1)
+    main_cmd(sys.argv[1], maximum_distance='200', msr='',
+             approach='use_shakemap_from_usgs', loglevel='info',
+             userlevel=1)

openquake/engine/aristotle.py

@@ -64,6 +64,7 @@ def main_cmd(usgs_id, rupture_file=None,
              ses_seed='42', local_timestamp='',
              exposure_hdf5=None, station_data_file=None,
              maximum_distance_stations='',
+             msr='', approach='use_shakemap_from_usgs',
              loglevel='warn',
              userlevel=2):
     """
@@ -80,6 +81,8 @@ def main_cmd(usgs_id, rupture_file=None,
         'rupture_dict', 'rupture_file', 'station_data_file'}
     post = {f: loc.get(f) for f in fields}
     post['usgs_id'] = usgs_id
+    post['msr'] = msr
+    post['approach'] = approach
     monitor = performance.Monitor()
     _rup, rupdic, oqparams, err = aristotle_validate(
         post, User(level=userlevel), rupture_file, station_data_file, monitor)
@@ -99,6 +102,7 @@ def main_cmd(usgs_id, rupture_file=None,
     else:
         callback(job.calc_id, oqparams, exc=None)
 
+
 main_cmd.usgs_id = 'ShakeMap ID'  # i.e. us6000m0xl
 main_cmd.rupture_file = 'XML file with the rupture model (optional)'
 main_cmd.rupture_dict = 'Used by the command `oq mosaic aristotle`'
@@ -116,6 +120,8 @@ def main_cmd(usgs_id, rupture_file=None,
                           'and vulnerability functions')
 main_cmd.station_data_file = 'CSV file with the station data'
 main_cmd.maximum_distance_stations = 'Maximum distance from stations in km'
+main_cmd.msr = 'Magnitude scaling relationship'
+main_cmd.approach = 'For instance use_shakemap_from_usgs'
 main_cmd.loglevel = 'Log level'
 main_cmd.userlevel = 'User level'
 

openquake/hazardlib/scalerel/__init__.py

@@ -36,7 +36,7 @@ def _get_available_class(base_class):
     scalerel submodule with classes that derives from `base_class`,
     keyed by class name.
     '''
-    gsims = {}
+    classes = {}  # class_name -> class
     for fname in os.listdir(os.path.dirname(__file__)):
         if fname.endswith('.py'):
             modname, _ext = os.path.splitext(fname)
@@ -46,16 +46,23 @@ def _get_available_class(base_class):
                 if inspect.isclass(cls) and issubclass(cls, base_class) \
                         and cls != base_class \
                         and not inspect.isabstract(cls):
-                    gsims[cls.__name__] = cls
-    return dict((k, gsims[k]) for k in sorted(gsims))
+                    classes[cls.__name__] = cls
+    return dict((k, classes[k]) for k in sorted(classes))
 
 
 def get_available_magnitude_scalerel():
     '''
-    Return an ordered dictionary with the available Magnitude ScaleRel
-    classes, keyed by class name.
+    Return a list of the available Magnitude ScaleRel instances.
     '''
-    return _get_available_class(BaseMSR)
+    dic = _get_available_class(BaseMSR)
+    msrs = []
+    for name, cls in dic.items():
+        if inspect.signature(cls).parameters:
+            # ignore CScaling
+            pass
+        else:
+            msrs.append(cls())
+    return msrs
 
 
 def get_available_sigma_magnitude_scalerel():
@@ -88,6 +95,6 @@ def get_available_scalerel():
     keyed by class name.
     '''
     ret = get_available_area_scalerel()
-    ret.update(get_available_magnitude_scalerel())
+    ret.update(_get_available_class(BaseMSR))
 
     return ret

openquake/hazardlib/shakemap/parsers.py

@@ -43,6 +43,7 @@
 from openquake.baselib.node import node_from_xml
 from openquake.hazardlib import nrml, sourceconverter
 from openquake.hazardlib.source.rupture import get_multiplanar, is_matrix
+from openquake.hazardlib.scalerel import get_available_magnitude_scalerel
 
 NOT_FOUND = 'No file with extension \'.%s\' file found'
 US_GOV = 'https://earthquake.usgs.gov'
@@ -71,6 +72,7 @@
 }
 REQUIRED_IMTS = {'PGA', 'PSA03', 'PSA10'}
 
+
 @dataclass
 class User:
     level: int = 0
@@ -257,6 +259,8 @@ def utc_to_local_time(utc_timestamp, lon, lat):
     utc_time = utc_zone.localize(utc_time)
     local_zone = pytz.timezone(timezone_str)
     local_timestamp = utc_time.astimezone(local_zone)
+    # NOTE: the validated timestamp format has no microseconds
+    local_timestamp = local_timestamp.replace(microsecond=0)
     return local_timestamp
 
 
@@ -402,18 +406,19 @@ def usgs_to_ecd_format(stations, exclude_imts=()):
     return df_seismic_non_null
 
 
-def _get_preferred_shakemap(shakemaps):
-    preferred_weights = [shakemap['preferredWeight'] for shakemap in shakemaps]
+def _get_preferred_item(items):
+    # items can be for instance shakemaps, moment tensors or finite faults
+    preferred_weights = [item['preferredWeight'] for item in items]
     preferred_idxs = [idx for idx, val in enumerate(preferred_weights)
                       if val == max(preferred_weights)]
-    preferred_shakemaps = [shakemaps[idx] for idx in preferred_idxs]
-    if len(preferred_shakemaps) > 1:
-        update_times = [shakemap['updateTime'] for shakemap in preferred_shakemaps]
+    preferred_items = [items[idx] for idx in preferred_idxs]
+    if len(preferred_items) > 1:
+        update_times = [item['updateTime'] for item in preferred_items]
         latest_idx = update_times.index(max(update_times))
-        shakemap = preferred_shakemaps[latest_idx]
+        item = preferred_items[latest_idx]
     else:
-        shakemap = preferred_shakemaps[0]
-    return shakemap
+        item = preferred_items[0]
+    return item
 
 
 def download_station_data_file(usgs_id, contents, user):
@@ -472,17 +477,34 @@ def load_rupdic_from_finite_fault(usgs_id, mag, products):
     """
     Extract the finite fault properties from products.
     NB: if the finite-fault list contains multiple elements we take the
-    first one.
+    preferred one.
     """
     logging.info('Getting finite-fault properties')
-    try:
-        ff = products['finite-fault']
-    except KeyError:
+    if 'finite-fault' not in products:
+        # e.g. us6000phrk
         # FIXME: not tested
-        raise MissingLink('There is no finite-fault info for %s' % usgs_id)
-    if isinstance(ff, list):
-        ff = ff[0]
+        err_msg = f'There is no finite-fault info for {usgs_id}'
+        err = {"status": "failed", "error_msg": err_msg}
+        return None, err
+    ffs = products['finite-fault']
+    ff = _get_preferred_item(ffs)
     p = ff['properties']
+    # TODO: we probably need to get the rupture coordinates from shakemap_polygon.txt
+    # if 'shakemap_polygon.txt' in ff['contents']:
+    #     # with open(f'/tmp/{usgs_id}-shakemap_polygon.txt', 'wb') as f:
+    #     #       f.write(urlopen(url).read())
+    #     if user.testdir:  # in parsers_test
+    #         fname = os.path.join(user.testdir, f'{usgs_id}-shakemap_polygon.txt')
+    #         text = open(fname).read()
+    #     else:
+    #         url = ff['contents']['shakemap_polygon.txt']['url']
+    #         logging.info('Downloading shakemap_polygon.txt')
+    #         text = urlopen(url).read()
+    #         lines = text.decode('utf8').split("\n")
+    #         numerical_data = [line for line in lines
+    #                           if line and not line.startswith("#")]
+    #         coords = [tuple(map(float, line.split()))
+    #                   for line in numerical_data if line]
     lon = float(p['longitude'])
     lat = float(p['latitude'])
     utc_time = p['eventtime']
@@ -494,7 +516,41 @@ def load_rupdic_from_finite_fault(usgs_id, mag, products):
               'require_dip_strike': True,
               'pga_map_png': None, 'mmi_map_png': None,
               'usgs_id': usgs_id, 'rupture_file': None}
-    return rupdic
+    return rupdic, None
+
+
+def load_rupdic_from_origin(usgs_id, products):
+    """
+    Extract the origin properties from products.
+    NB: if the origin list contains multiple elements we take the
+    preferred one.
+    """
+    # TODO: we may try to unify this with the very similar
+    # load_rupdic_from_finite_fault
+    logging.info('Getting origin properties')
+    if 'origin' not in products:
+        # FIXME: not tested
+        err_msg = f'There is no origin info for {usgs_id}'
+        err = {"status": "failed", "error_msg": err_msg}
+        return None, err
+    origins = products['origin']
+    origin = _get_preferred_item(origins)
+    p = origin['properties']
+    mag = float(p['magnitude'])
+    lon = float(p['longitude'])
+    lat = float(p['latitude'])
+    dep = float(p['depth'])
+    rake = 0.
+    utc_time = p['eventtime']
+    local_time = utc_to_local_time(utc_time, lon, lat)
+    time_event = local_time_to_time_event(local_time)
+    rupdic = {'lon': lon, 'lat': lat, 'dep': dep,
+              'mag': mag, 'rake': rake,
+              'local_timestamp': str(local_time), 'time_event': time_event,
+              'require_dip_strike': True,
+              'pga_map_png': None, 'mmi_map_png': None,
+              'usgs_id': usgs_id, 'rupture_file': None}
+    return rupdic, None
 
 
 # NB: not used right now
@@ -638,7 +694,7 @@ def _contents_properties_shakemap(usgs_id, user, use_shakemap, monitor):
     properties = js['properties']
 
     # NB: currently we cannot find a case with missing shakemap
-    shakemap = _get_preferred_shakemap(properties['products']['shakemap'])
+    shakemap = _get_preferred_item(properties['products']['shakemap'])
     contents = shakemap['contents']
 
     if (user.level == 1 or use_shakemap) and 'download/grid.xml' in contents:
@@ -657,15 +713,54 @@ def _contents_properties_shakemap(usgs_id, user, use_shakemap, monitor):
     return contents, properties, shakemap_array, err
 
 
-def get_rup_dic(usgs_id, user=User(), use_shakemap=False,
-                rupture_file=None, station_data_file=None,
-                monitor=performance.Monitor()):
+def _get_nodal_planes(properties):
+    # in parsers_test
+    err = None
+    if 'moment-tensor' not in properties['products']:
+        err = 'Unable to retrieve information about the nodal options'
+        return None, err
+    moment_tensor = _get_preferred_item(properties['products']['moment-tensor'])
+    props = moment_tensor['properties']
+    nodal_planes = {}
+    for key, value in props.items():
+        if key.startswith('nodal-plane-'):
+            parts = key.split('-')
+            plane = f'NP{parts[2]}'
+            attr = parts[3]  # Get the attribute (i.e. 'dip', 'rake' or 'strike')
+            if plane not in nodal_planes:
+                nodal_planes[plane] = {}
+            nodal_planes[plane][attr] = float(value)
+    return nodal_planes, err
+
+
+def _get_rup_dic_from_xml(usgs_id, user, rupture_file, station_data_file):
+    [rup_node] = nrml.read(os.path.join(user.testdir, rupture_file)
+                           if user.testdir else rupture_file)
+    rup = sourceconverter.RuptureConverter(
+        rupture_mesh_spacing=5.).convert_node(rup_node)
+    rup.tectonic_region_type = '*'
+    hp = rup.hypocenter
+    rupdic = dict(lon=hp.x, lat=hp.y, dep=hp.z,
+                  mag=rup.mag, rake=rup.rake,
+                  strike=rup.surface.get_strike(),
+                  dip=rup.surface.get_dip(),
+                  usgs_id=usgs_id,
+                  rupture_file=rupture_file,
+                  station_data_file=station_data_file)
+    return rup, rupdic
+
+
+def get_rup_dic(dic, user=User(), approach='use_shakemap_from_usgs',
+                use_shakemap=False, rupture_file=None,
+                station_data_file=None, monitor=performance.Monitor()):
     """
     If the rupture_file is None, download a rupture from the USGS site given
     the ShakeMap ID, else build the rupture locally with the given usgs_id.
 
-    :param usgs_id: ShakeMap ID
+    :param dic: dictionary with ShakeMap ID and other parameters
     :param user: User instance
+    :param approach: the workflow selected by the user
+        (default: 'use_shakemap_from_usgs')
     :param use_shakemap: download the ShakeMap only if True
     :param rupture_file: None
     :param station_data_file: None
@@ -674,20 +769,18 @@ def get_rup_dic(usgs_id, user=User(), use_shakemap=False,
     rupdic = {}
     rup_data = {}
     err = {}
+    usgs_id = dic['usgs_id']
+    rup = None
+    if approach == 'provide_rup_params':
+        rupdic = dic.copy()
+        rupdic.update(rupture_file=rupture_file,
+                      station_data_file=station_data_file,
+                      require_dip_strike=True)
+        return rup, rupdic, err
+
     if rupture_file and rupture_file.endswith('.xml'):
-        [rup_node] = nrml.read(os.path.join(user.testdir, rupture_file)
-                               if user.testdir else rupture_file)
-        rup = sourceconverter.RuptureConverter(
-            rupture_mesh_spacing=5.).convert_node(rup_node)
-        rup.tectonic_region_type = '*'
-        hp = rup.hypocenter
-        rupdic = dict(lon=hp.x, lat=hp.y, dep=hp.z,
-                      mag=rup.mag, rake=rup.rake,
-                      strike=rup.surface.get_strike(),
-                      dip=rup.surface.get_dip(),
-                      usgs_id=usgs_id,
-                      rupture_file=rupture_file,
-                      station_data_file=station_data_file)
+        rup, rupdic = _get_rup_dic_from_xml(
+            usgs_id, user, rupture_file, station_data_file)
         if usgs_id == 'FromFile':
             return rup, rupdic, err
     elif rupture_file and rupture_file.endswith('.json'):
@@ -704,18 +797,31 @@ def get_rup_dic(usgs_id, user=User(), use_shakemap=False,
         usgs_id, user, use_shakemap, monitor)
     if err:
         return None, None, err
-
-    if 'download/rupture.json' not in contents:
+    if approach in ['use_pnt_rup_from_usgs', 'build_rup_from_usgs']:
+        rupdic, err = load_rupdic_from_origin(usgs_id, properties['products'])
+        if err:
+            return None, None, err
+    elif ('download/rupture.json' not in contents
+          or approach == 'use_finite_rup_from_usgs'):
         # happens for us6000f65h in parsers_test
-        rupdic = load_rupdic_from_finite_fault(
+        rupdic, err = load_rupdic_from_finite_fault(
             usgs_id, properties['mag'], properties['products'])
+        if err:
+            return None, None, err
+
+    if approach == 'build_rup_from_usgs':
+        rupdic['nodal_planes'], err = _get_nodal_planes(properties)
+        rupdic['msrs'] = [msr.__class__.__name__
+                          for msr in get_available_magnitude_scalerel()]
+        if err:
+            return None, rupdic, err
+
     if not rupdic:
         if not rup_data:
             with monitor('Downloading rupture json'):
                 rup_data, rupture_file = download_rupture_data(
                     usgs_id, contents, user)
         rupdic = convert_rup_data(rup_data, usgs_id, rupture_file, shakemap)
-
     if user.level == 2 and not station_data_file:
         with monitor('Downloading stations'):
             rupdic['station_data_file'], rupdic['station_data_issue'] = (