Merge pull request #10070 from gem/idl

ARISTOTLE: change the way a hypocenter is considered close to any mosaic models

openquake/calculators/event_based.py

@@ -600,7 +600,7 @@ def _read_scenario_ruptures(self):
                 elif oq.rupture_xml:
                     hypo = readinput.get_rupture(oq).hypocenter
                     lon, lat = [hypo.x, hypo.y]
-                mosaic_models = get_close_mosaic_models(lon, lat, 100)
+                mosaic_models = get_close_mosaic_models(lon, lat, 5)
                 # NOTE: using the first mosaic model
                 oq.mosaic_model = mosaic_models[0]
                 if len(mosaic_models) > 1:

openquake/engine/aristotle.py

@@ -24,6 +24,7 @@
 import logging
 from dataclasses import dataclass
 import numpy
+from shapely.geometry import Point
 from json.decoder import JSONDecodeError
 from urllib.error import HTTPError
 from openquake.baselib import config, hdf5, sap
@@ -54,32 +55,31 @@ class AristotleParam:
     ignore_shakemap: bool = False
 
 
-def get_close_mosaic_models(lon, lat, max_dist):
+def get_close_mosaic_models(lon, lat, buffer_radius):
     """
     :param lon: longitude
     :param lat: latitude
-    :param max_dist: max dist with respect to a model to consider it relevant
-    :returns: list of mosaic models with the max_dist distance
+    :param buffer_radius: radius of the buffer around the point.
+        This distance is in the same units as the point's
+        coordinates (i.e. degrees), and it defines how far from
+        the point the buffer should extend in all directions,
+        creating a circular buffer region around the point
+    :returns: list of mosaic models intersecting the circle
+        centered on the given coordinates having the specified radius
     """
     mosaic_df = readinput.read_mosaic_df(buffer=1)
-    hypocenter = geo.Point(lon, lat)
-    minx, miny, maxx, maxy = geo.utils.get_bounding_box([hypocenter], max_dist)
-    close_mosaic_models = set()
-    bbox_vertices = [(minx, maxy), (maxx, maxy), (maxx, miny), (minx, miny)]
-    for vertex in bbox_vertices:
-        lonlats = numpy.array([vertex[0], vertex[1]])
-        [mosaic_model] = geo.utils.geolocate([lonlats], mosaic_df)
-        close_mosaic_models.add(mosaic_model)
-    close_mosaic_models = sorted([model for model in close_mosaic_models
-                                  if model != '???'])
+    hypocenter = Point(lon, lat)
+    hypo_buffer = hypocenter.buffer(buffer_radius)
+    geoms = numpy.array([hypo_buffer])
+    [close_mosaic_models] = geo.utils.geolocate_geometries(geoms, mosaic_df)
     if not close_mosaic_models:
         raise ValueError(
-            f'({lon}, {lat}) is farther than {max_dist}km'
+            f'({lon}, {lat}) is farther than {buffer_radius} deg'
             f' from any mosaic model!')
     elif len(close_mosaic_models) > 1:
         logging.info(
-            '(%s, %s) is closer than %skm with respect to the following'
-            ' mosaic models: %s' % (lon, lat, max_dist, close_mosaic_models))
+            '(%s, %s) is closer than %s deg with respect to the following'
+            ' mosaic models: %s' % (lon, lat, buffer_radius, close_mosaic_models))
     return close_mosaic_models
 
 
@@ -166,7 +166,7 @@ def get_aristotle_params(arist):
         inputs['station_data'] = arist.station_data_file
     if not arist.mosaic_model:
         lon, lat = rupdic['lon'], rupdic['lat']
-        mosaic_models = get_close_mosaic_models(lon, lat, 100)
+        mosaic_models = get_close_mosaic_models(lon, lat, 5)
         # NOTE: using the first mosaic model
         arist.mosaic_model = mosaic_models[0]
         if len(mosaic_models) > 1:

openquake/engine/tests/aristotle4/aggrisk.txt

@@ -0,0 +1,4 @@
+| gsim             | weight | contents | nonstructural | structural | occupants | area      | number    | residents |
+|------------------+--------+----------+---------------+------------+-----------+-----------+-----------+-----------|
+| [McVerry2006Asc] | 1.0000 | 4.1815   | 0.0011        | 6.633E-04  | 6.700E-08 | 1.194E-06 | 7.000E-09 | 1.294E-07 |
+| Average          | 1.0000 | 4.1815   | 0.0011        | 6.633E-04  | 6.700E-08 | 1.194E-06 | 7.000E-09 | 1.294E-07 |

openquake/engine/tests/aristotle4/job.ini

@@ -0,0 +1,12 @@
+[general]
+description = Compute risk around a rupture close to the international date line
+calculation_mode = scenario_risk
+exposure_file = ../../../../exposure.hdf5
+
+[calculation]
+rupture_dict = {'lon': -179.8, 'lat': -39, 'dep': 30, 'mag': 8.6, 'rake': 63}
+truncation_level = 3.0
+maximum_distance = 600
+number_of_ground_motion_fields = 10
+asset_hazard_distance = 15
+tectonic_region_type = Active Shallow Crust

openquake/engine/tests/aristotle_test.py

@@ -1,18 +1,18 @@
 # -*- coding: utf-8 -*-
 # vim: tabstop=4 shiftwidth=4 softtabstop=4
-# 
+#
 # Copyright (C) 2024, GEM Foundation
-# 
+#
 # OpenQuake is free software: you can redistribute it and/or modify it
 # under the terms of the GNU Affero General Public License as published
 # by the Free Software Foundation, either version 3 of the License, or
 # (at your option) any later version.
-# 
+#
 # OpenQuake is distributed in the hope that it will be useful,
 # but WITHOUT ANY WARRANTY; without even the implied warranty of
 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 # GNU Affero General Public License for more details.
-# 
+#
 # You should have received a copy of the GNU Affero General Public License
 # along with OpenQuake.  If not, see <http://www.gnu.org/licenses/>.
 
@@ -44,7 +44,7 @@ def check_export_job(dstore):
                       'taxonomy_mapping.csv']
 
 
-@pytest.mark.parametrize('n', [1, 2, 3])
+@pytest.mark.parametrize('n', [1, 2, 3, 4])
 def test_aristotle(n):
     if not os.path.exists(cd.parent.parent.parent / 'exposure.hdf5'):
         raise unittest.SkipTest('Please download exposure.hdf5')

openquake/hazardlib/geo/utils.py

@@ -929,18 +929,38 @@ def geolocate(lonlats, geom_df, exclude=()):
     """
     :param lonlats: array of shape (N, 2) of (lon, lat)
     :param geom_df: DataFrame of geometries with a "code" field
+    :param exclude: List of codes to exclude from the results
     :returns: codes associated to the points
 
     NB: if the "code" field is not a primary key, i.e. there are
     different geometries with the same code, performs an "or", i.e.
     associates the code if at least one of the geometries matches
     """
     codes = numpy.array(['???'] * len(lonlats))
-    for code, df in geom_df.groupby('code'):
-        if code in exclude:
-            continue
+    filtered_geom_df = geom_df[~geom_df['code'].isin(exclude)]
+    for code, df in filtered_geom_df.groupby('code'):
         ok = numpy.zeros(len(lonlats), bool)
         for geom in df.geom:
             ok |= contains_xy(geom, lonlats)
         codes[ok] = code
     return codes
+
+
+def geolocate_geometries(geometries, geom_df, exclude=()):
+    """
+    :param geometries: NumPy array of Shapely geometries to check
+    :param geom_df: DataFrame of geometries with a "code" field
+    :param exclude: List of codes to exclude from the results
+    :returns: NumPy array where each element contains a list of codes
+        of geometries that intersect each input geometry
+    """
+    result_codes = numpy.empty(len(geometries), dtype=object)
+    filtered_geom_df = geom_df[~geom_df['code'].isin(exclude)]
+    for i, input_geom in enumerate(geometries):
+        intersecting_codes = set()  # to store intersecting codes for current geometry
+        for code, df in filtered_geom_df.groupby('code'):
+            target_geoms = df['geom'].values  # geometries associated with this code
+            if any(target_geom.intersects(input_geom) for target_geom in target_geoms):
+                intersecting_codes.add(code)
+        result_codes[i] = sorted(intersecting_codes)
+    return result_codes

openquake/server/tests/test_aristotle_mode.py

@@ -260,8 +260,17 @@ def test_get_rupture_data_from_shakemap_conversion_error(self):
         self.assertEqual(ret_dict['local_timestamp'],
                          '2023-02-06 04:17:34+03:00')
         self.assertEqual(ret_dict['time_event'], 'night')
-        self.assertEqual(ret_dict['mosaic_models'], ['MIE'])
+        self.assertEqual(ret_dict['mosaic_models'], ['ARB', 'MIE'])
         self.assertEqual(ret_dict['trts'], {
+            'ARB': ['TECTONIC_REGION_1',
+                    'TECTONIC_REGION_2',
+                    'TECTONIC_REGION_3',
+                    'TECTONIC_REGION_4',
+                    'Active Shallow Crust EMME',
+                    'Stable Shallow Crust EMME',
+                    'Subduction Interface EMME',
+                    'Subduction Inslab EMME',
+                    'Deep Seismicity EMME'],
             'MIE': ['Active Shallow Crust',
                     'Stable Shallow Crust',
                     'Subduction Interface',
@@ -347,7 +356,22 @@ def test_get_rupture_data_from_finite_fault(self):
         self.assertEqual(sorted(ret_dict.keys()), sorted(expected_keys))
         self.assertEqual(ret_dict['rupture_file'], None)
         self.assertEqual(ret_dict['usgs_id'], 'us6000jllz')
-        self.assertEqual(ret_dict['mosaic_models'], ['MIE'])
+        self.assertEqual(ret_dict['mosaic_models'], ['ARB', 'MIE'])
+        self.assertEqual(ret_dict['trts'], {
+            'ARB': ['TECTONIC_REGION_1',
+                    'TECTONIC_REGION_2',
+                    'TECTONIC_REGION_3',
+                    'TECTONIC_REGION_4',
+                    'Active Shallow Crust EMME',
+                    'Stable Shallow Crust EMME',
+                    'Subduction Interface EMME',
+                    'Subduction Inslab EMME',
+                    'Deep Seismicity EMME'],
+            'MIE': ['Active Shallow Crust',
+                    'Stable Shallow Crust',
+                    'Subduction Interface',
+                    'Subduction Inslab',
+                    'Deep Seismicity']})
 
     def test_run_by_usgs_id_then_remove_calc(self):
         data = dict(usgs_id='us6000jllz',

openquake/server/views.py

@@ -734,7 +734,7 @@ def aristotle_get_rupture_data(request):
             ' identifier "%s": %s' % (rupdic['usgs_id'], station_data_file))
         station_data_file = None
     try:
-        mosaic_models = get_close_mosaic_models(rupdic['lon'], rupdic['lat'], 100)
+        mosaic_models = get_close_mosaic_models(rupdic['lon'], rupdic['lat'], 5)
         for mosaic_model in mosaic_models:
             trts[mosaic_model] = get_trts_around(
                 mosaic_model,