Complex Fault Refactoring

openquake/calculators/tests/classical_test.py

@@ -276,7 +276,7 @@ def test_case_37(self):
                                "quantile_curve-0.16-PGA.csv",
                                "quantile_curve-0.5-PGA.csv",
                                "quantile_curve-0.84-PGA.csv"],
-                              case_37.__file__)
+                              case_37.__file__, delta=0.001)
 
     def test_case_38(self):
         # BC Hydro GMPEs with epistemic adjustments
@@ -309,7 +309,8 @@ def test_case_41(self):
     def test_case_42(self):
         # split/filter a long complex fault source with maxdist=1000 km
         self.assert_curves_ok(["hazard_curve-mean-PGA.csv",
-                               "hazard_map-mean-PGA.csv"], case_42.__file__)
+                               "hazard_map-mean-PGA.csv"],
+                              case_42.__file__, delta=0.0002)
 
         # check pandas readability of hmaps-stats
         df = self.calc.datastore.read_df('hmaps-stats', 'site_id',
@@ -506,13 +507,14 @@ def test_case_62(self):
         # multisurface with kite faults
         self.run_calc(case_62.__file__, 'job.ini')
         [f] = export(('hcurves/mean', 'csv'), self.calc.datastore)
-        self.assertEqualFiles('expected/hcurve-mean.csv', f, delta=1E-5)
+        self.assertEqualFiles('expected/hcurve-mean.csv', f, delta=2E-4)
 
     def test_case_63(self):
         # test soiltype
         self.run_calc(case_63.__file__, 'job.ini')
         [f] = export(('hcurves/mean', 'csv'), self.calc.datastore)
-        self.assertEqualFiles('expected/hazard_curve-mean-PGA.csv', f)
+        self.assertEqualFiles('expected/hazard_curve-mean-PGA.csv', f,
+                              delta=2E-3)
 
     def test_case_64(self):
         # LanzanoEtAl2016 with bas term

openquake/hazardlib/geo/line.py

@@ -25,16 +25,26 @@
 from openquake.hazardlib.geo import Point
 
 TOLERANCE = 0.1
+LARGE = 1e10
 
 
 def _update(rtra, rtra_prj, proj, pnt):
+    # Updates the arrays `rtra` and `rtra_prj` with the resampled trace (in
+    # projected and geographic coodinates). `proj` is the function for
+    # converting geographic and projected coordinates and `pnt` is the point to
+    # add
     xg, yg = proj(np.array([pnt[0]]), np.array([pnt[1]]), reverse=True)
     rtra.append(np.array([xg[0], yg[0], pnt[2]]))
     rtra_prj.append(pnt)
 
 
 def _resample(coo, sect_len, orig_extremes):
-    # returns array of resampled trace coordinates
+    # returns array of resampled trace coordinates where 'coo' is an array with
+    # N lines and 3 colums. `sect_len` is the length in km to be used for
+    # resampling and `orig_extremes` is a boolean controlling the inclusion or
+    # exclusion of the last original point
+
+    # Number of coordinates in the original polyline
     N = len(coo)
 
     # Project the coordinates of the trace
@@ -44,56 +54,87 @@ def _resample(coo, sect_len, orig_extremes):
     txy[:, 0], txy[:, 1] = proj(coo[:, 0], coo[:, 1])
 
     # Compute the total length of the original trace
-    tot_len = sum(utils.get_dist(txy[i], txy[i-1]) for i in range(1, N))
+    tot_len = sum(utils.get_dist(txy[i], txy[i - 1]) for i in range(1, N))
     inc_len = 0.
 
-    # Initialize the lists with the coordinates of the resampled trace
+    # Initialize the lists with the coordinates of the resampled trace in
+    # projected and geographic coordinates
     rtra_prj = [txy[0]]
     rtra = [coo[0]]
 
     # Resampling
     idx_vtx = -1
     while True:
 
-        # Computing distances from the reference point
-        dis = utils.get_dist(txy, rtra_prj[-1])
-        if idx_vtx > 0:
-            # Fixing distances for points before the index
-            dis[0:idx_vtx] = 100000
+        # Computing distances between the reference point and the points
+        # forming the original trace
+        tmp = [(np.sum((coo - rtra_prj[-1])**2)**0.5) for coo in txy]
+        dis = np.array(tmp).squeeze()
 
-        # Index of the point on the trace with a distance just below the
-        # sampling distance
-        idx = np.where(dis <= sect_len, dis, -np.inf).argmax()
+        # Fixing distances for points before the reference index
+        if idx_vtx > 0:
+            dis[0:idx_vtx] = LARGE
+
+        # Find the index of the point on the trace with a distance just below
+        # the sampling distance. This logic is required to cope with the test
+        # on the complex fault where the profiles have a sort of convex shape
+        max_dist = -1.0
+        idx = idx_vtx if idx_vtx >= 0 else 0
+        idx_pre = idx
+        while 1:
+            if dis[idx] <= sect_len and dis[idx] > max_dist:
+                idx_pre = copy.copy(idx)
+                max_dist = dis[idx]
+                if idx == len(dis) - 1:
+                    break
+                idx += 1
+            else:
+                idx = idx_pre
+                break
 
-        # If the pick a point that is not the last one on the trace we
-        # compute the new sample by interpolation
+        # If we pick a point that is not the last one on the trace we compute
+        # the new sample by interpolation
         if idx < len(dis) - 1:
+            # Find the point between two consecutive points on the original
+            # (projected) trace `txy` at a distance `sect_len` from the first
+            # one `txy[idx, :]`
             pnt = find_t(
                 txy[idx + 1, :], txy[idx, :], rtra_prj[-1], sect_len)
             if pnt is None:
                 raise ValueError('Did not find the intersection')
             _update(rtra, rtra_prj, proj, pnt)
             inc_len += sect_len
 
-            # Adding more points still on the same segment
+            # Adding more points still on the same segment. `delta` is the
+            # total increment along the x, y and z axes required to move from
+            # the last point on the resampled trace and the end of the current
+            # segment `txy[idx + 1]`
             delta = txy[idx + 1] - rtra_prj[-1]
-            chk_dst = utils.get_dist(txy[idx + 1], rtra_prj[-1])
+            chk_dst = (np.sum((txy[idx + 1] - rtra_prj[-1])**2))**0.5
             rat = delta / chk_dst
             while chk_dst > sect_len * 0.9999:
+                # Adding one more point
                 _update(rtra, rtra_prj, proj, rtra_prj[-1] + sect_len * rat)
                 inc_len += sect_len
                 # This is the distance between the last resampled point
                 # and the second vertex of the segment
-                chk_dst = utils.get_dist(txy[idx + 1], rtra_prj[-1])
+                chk_dst = (np.sum((txy[idx + 1] - rtra_prj[-1])**2))**0.5
         else:
-            # Adding one point
-            if tot_len - inc_len > 0.5 * sect_len and not orig_extremes:
-                # Adding more points still on the same segment
-                delta = txy[-1] - txy[-2]
-                chk_dst = utils.get_dist(txy[-1], txy[-2])
-                _update(rtra, rtra_prj, proj, rtra_prj[-1] +
-                       sect_len * delta / chk_dst)
-                inc_len += sect_len
+
+            # Adding more points still on the same segment, if needed
+            delta = txy[-1] - txy[-2]
+
+            chk_dst = (np.sum((txy[-1] - txy[-2])**2))**0.5
+
+            # New point
+            new_point = rtra_prj[-1] + sect_len * delta / chk_dst
+
+            # Distance between this point and the last one on the original
+            # section
+            dst = (np.sum((new_point - txy[-1])**2))**0.5
+
+            if dst < sect_len * 0.5:
+                _update(rtra, rtra_prj, proj, new_point)
             elif orig_extremes:
                 # Adding last point
                 rtra.append(coo[-1])
@@ -611,20 +652,22 @@ def find_t(pnt0, pnt1, ref_pnt, distance):
         A 1D :class:`numpy.ndarray` instance of length 3
     """
 
-    # First points on the line
+    # First point on the line
     x1 = pnt0[0]
     y1 = pnt0[1]
     z1 = pnt0[2]
 
-    #
+    # Second point on the line
     x2 = pnt1[0]
     y2 = pnt1[1]
     z2 = pnt1[2]
 
+    # Reference point
     x3 = ref_pnt[0]
     y3 = ref_pnt[1]
     z3 = ref_pnt[2]
 
+    # The sampling distance
     r = distance
 
     pa = (x2 - x1)**2 + (y2 - y1)**2 + (z2 - z1)**2

openquake/hazardlib/geo/surface/complex_fault.py

@@ -27,6 +27,7 @@
 from openquake.baselib.node import Node
 from openquake.hazardlib.geo.line import Line
 from openquake.hazardlib.geo.point import Point
+from openquake.hazardlib.geo.line import _resample
 from openquake.hazardlib.geo.surface.base import BaseSurface
 from openquake.hazardlib.geo.surface.planar import PlanarSurface
 from openquake.hazardlib.geo.mesh import Mesh, RectangularMesh
@@ -272,27 +273,47 @@ def from_fault_data(cls, edges, mesh_spacing):
         cls.check_fault_data(edges, mesh_spacing)
         surface_nodes = [complex_fault_node(edges)]
         mean_length = numpy.mean([edge.get_length() for edge in edges])
-        num_hor_points = int(round(mean_length / mesh_spacing)) + 1
+        num_hor_points = int(numpy.round(mean_length / mesh_spacing)) + 1
+
         if num_hor_points <= 1:
             raise ValueError(
                 'mesh spacing %.1f km is too big for mean length %.1f km' %
                 (mesh_spacing, mean_length)
             )
-        edges = [edge.resample_to_num_points(num_hor_points).points
-                 for i, edge in enumerate(edges)]
 
-        vert_edges = [Line(v_edge) for v_edge in zip(*edges)]
-        mean_width = numpy.mean([v_edge.get_length() for v_edge in vert_edges])
-        num_vert_points = int(round(mean_width / mesh_spacing)) + 1
+        lengths = [sum(edge.get_lengths()) for edge in edges]
+        redges = [_resample(edge.coo, tlen / (num_hor_points - 1) * 0.99,
+                            False)
+                  for edge, tlen in zip(edges, lengths)]
+        new_edges = [[Point(c[0], c[1], c[2]) for c in coo] for coo in redges]
+        vert_edges = [Line(v_edge) for v_edge in zip(*new_edges)]
+
+        vert_edges_lenghts = [v_edge.get_length() for v_edge in vert_edges]
+        mean_width = numpy.mean(vert_edges_lenghts)
+        num_vert_points = int(numpy.round(mean_width / mesh_spacing)) + 1
+
         if num_vert_points <= 1:
             raise ValueError(
                 'mesh spacing %.1f km is too big for mean width %.1f km' %
                 (mesh_spacing, mean_width)
             )
 
-        points = zip(*[v_edge.resample_to_num_points(num_vert_points).points
-                       for v_edge in vert_edges])
-        mesh = RectangularMesh.from_points_list(list(points))
+        vlengths = [lng / (num_vert_points - 1) for lng in vert_edges_lenghts]
+        fun = zip(vert_edges, vlengths)
+        vedges = []
+        for v_edge, lng in fun:
+            vedges.append(_resample(v_edge.coo, lng, False))
+
+        profiles = []
+        for i_row in range(len(vedges[0])):
+            tmp_profile = []
+            for edge in vedges:
+                tmp_profile.append(Point(edge[i_row][0], edge[i_row][1],
+                                         edge[i_row][2]))
+            profiles.append(tmp_profile)
+
+        mesh = RectangularMesh.from_points_list(profiles)
+
         assert 1 not in mesh.shape
         self = cls(mesh)
         self.surface_nodes = surface_nodes

openquake/hazardlib/tests/geo/line_test.py

@@ -18,21 +18,23 @@
 import matplotlib.pyplot as plt
 from openquake.hazardlib import geo
 
-PLOTTING = False
+PLOTTING = True
 
 
-def _plott(rtra_prj, txy):
+def _plott(rtra_prj, txy, scl=1.0):
     import matplotlib.pyplot as plt
-    fig, ax = plt.subplots(1, 1)
+    ax = plt.figure().add_subplot(projection='3d')
     tt = np.array(rtra_prj)
-    plt.plot(txy[:, 0], txy[:, 1], '-')
-    plt.plot(txy[:, 0], txy[:, 1], 'x', ms=2.0)
-    plt.plot(tt[:, 0], tt[:, 1], 'o')
+    ax.plot(txy[:, 0], txy[:, 1], txy[:, 2] * scl, '-', color='cyan')
+    ax.plot(tt[:, 0], tt[:, 1], tt[:, 2] * scl, 'or')
     for i, t in enumerate(tt):
-        plt.text(t[0], t[1], f'{i}')
-    ax.axis('equal')
+        ax.text(t[0], t[1], t[2] * scl, f'{i}')
+    ax.plot(txy[:, 0], txy[:, 1], txy[:, 2] * scl, 'xb', ms=2.0)
+    ax.set_aspect('equal')
+    ax.invert_zaxis()
     plt.show()
 
+
 class LineResampleTestCase(unittest.TestCase):
 
     def test_resample(self):
@@ -64,7 +66,7 @@ def test_resample_dense_trace(self):
                                     computed.coo[1:, 1], zro)
         np.testing.assert_allclose(dst, 2.0, atol=0.02)
         if PLOTTING:
-            _plott(computed.coo, line.coo)
+            _plott(computed.coo, line.coo, scl=0.01)
 
     def test_resample_2(self):
         """
@@ -90,6 +92,30 @@ def test_resample_3(self):
         with self.assertRaises(ValueError):
             geo.Line([p1, p2, p3]).resample(50.0)
 
+    def test_resample_complex_fault(self):
+        # This is a profile resampled while building a complex fault surface
+        coo = np.array([[0.00, 0.0198, 1.0],
+                        [0.02, 0.0099, 1.5],
+                        [0.00, 0.0198, 2.0]])
+
+        # Prepare the function for resampling
+        resampler = geo.line._resample
+        sampl = resampler(coo, 1.0, False)
+
+        # Expected array
+        expected = np.array([
+            [0.00000000e+00, 1.98000000e-02, 1.00000000e+00],
+            [7.90103510e-03, 1.58889880e-02, 1.19752587e+00],
+            [1.58020699e-02, 1.19779756e-02, 1.39505174e+00],
+            [8.39767622e-03, 1.56431506e-02, 1.79005810e+00],
+            [4.96641154e-04, 1.95541627e-02, 1.98758397e+00]])
+
+        if PLOTTING:
+            _plott(coo, sampl, scl=0.01)
+
+        # Test
+        np.testing.assert_almost_equal(sampl, expected)
+
 
 class LineCreationTestCase(unittest.TestCase):
 
@@ -579,6 +605,3 @@ def plot_pattern(lons, lats, z, plons, plats, label, num=5, show=True):
     if show:
         plt.show()
     return ax
-
-
-

openquake/hazardlib/tests/geo/surface/_utils.py

@@ -20,7 +20,7 @@
 from openquake.hazardlib.geo.mesh import Mesh
 
 
-def assert_mesh_is(testcase, surface, expected_mesh):
+def assert_mesh_is(testcase, surface, expected_mesh, delta):
     expected_mesh = list(itertools.chain(*expected_mesh))
     testcase.assertEqual(len(surface.mesh), len(expected_mesh))
     testcase.assertIsInstance(surface.mesh, Mesh)
@@ -29,11 +29,11 @@ def assert_mesh_is(testcase, surface, expected_mesh):
         expected_point = Point(*expected_mesh[i])
         distance = expected_point.distance(point) * 1e3
         testcase.assertAlmostEqual(
-            0, distance, delta=2,  # allow discrepancy of 2 meters
+            0, distance, delta=delta,  # allow discrepancy of `delta` meters
             msg="point %d is off: %s != %s (distance is %.3fm)"
                 % (i, point, expected_point, distance))
 
 
 class SurfaceTestCase(unittest.TestCase):
-    def assert_mesh_is(self, surface, expected_mesh):
-        return assert_mesh_is(self, surface, expected_mesh)
+    def assert_mesh_is(self, surface, expected_mesh, delta=2):
+        return assert_mesh_is(self, surface, expected_mesh, delta)