diff --git a/openquake/calculators/tests/classical_test.py b/openquake/calculators/tests/classical_test.py index cd23d56fc719..38ded6e483b0 100644 --- a/openquake/calculators/tests/classical_test.py +++ b/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 diff --git a/openquake/hazardlib/geo/line.py b/openquake/hazardlib/geo/line.py index acf7a899d999..45d7cfc7188b 100644 --- a/openquake/hazardlib/geo/line.py +++ b/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,10 +54,11 @@ 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]] @@ -55,19 +66,38 @@ def _resample(coo, sect_len, orig_extremes): 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: @@ -75,25 +105,36 @@ def _resample(coo, sect_len, orig_extremes): _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 diff --git a/openquake/hazardlib/geo/surface/complex_fault.py b/openquake/hazardlib/geo/surface/complex_fault.py index 0234bcce30de..f2d7dc52839f 100644 --- a/openquake/hazardlib/geo/surface/complex_fault.py +++ b/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 diff --git a/openquake/hazardlib/tests/geo/line_test.py b/openquake/hazardlib/tests/geo/line_test.py index 5221781b358d..24d9c1753870 100644 --- a/openquake/hazardlib/tests/geo/line_test.py +++ b/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 - - - diff --git a/openquake/hazardlib/tests/geo/surface/_utils.py b/openquake/hazardlib/tests/geo/surface/_utils.py index b0d814552928..1e004f778fe9 100644 --- a/openquake/hazardlib/tests/geo/surface/_utils.py +++ b/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) diff --git a/openquake/hazardlib/tests/geo/surface/complex_fault_test.py b/openquake/hazardlib/tests/geo/surface/complex_fault_test.py index 768d3ef087a3..5338f57aa806 100644 --- a/openquake/hazardlib/tests/geo/surface/complex_fault_test.py +++ b/openquake/hazardlib/tests/geo/surface/complex_fault_test.py @@ -142,7 +142,25 @@ def test_invalid_surface_polygon_topo(self): ) +def plot(edges, mesh=None): + import matplotlib.pyplot as plt + ax = plt.figure().add_subplot(projection='3d') + for edge in edges: + coo = edge.coo + plt.plot(coo[:, 0], coo[:, 1], coo[:, 2], '-r') + + if mesh is not None: + coo = mesh.array + plt.plot(coo[0].flatten(), coo[1].flatten(), coo[2].flatten(), '.b') + + ax.set_box_aspect([1, 1, 1]) + ax.set_xlabel('lon') + ax.set_ylabel('lat') + plt.show() + + class ComplexFaultFromFaultDataTestCase(utils.SurfaceTestCase): + def test_1(self): edge1 = Line([Point(0, 0), Point(0.03, 0)]) edge2 = Line([Point(0, 0, 2.224), Point(0.03, 0, 2.224)]) @@ -155,7 +173,8 @@ def test_1(self): (0.02, 0, 1.112), (0.03, 0, 1.112)], [(0, 0, 2.224), (0.01, 0, 2.224), (0.02, 0, 2.224), (0.03, 0, 2.224)], - ]) + ], delta=40) + # plot([edge1, edge2], mesh=surface.mesh) def test_2(self): # this is a regression test. Reference values have been obtained @@ -163,33 +182,28 @@ def test_2(self): edge1 = Line([Point(0, 0, 1), Point(0, 0.02, 1)]) edge2 = Line([Point(0.02, 0, 1.5), Point(0.02, 0.01, 1.5)]) edge3 = Line([Point(0, 0, 2), Point(0, 0.02, 2)]) - surface = ComplexFaultSurface.from_fault_data([edge1, edge2, edge3], - mesh_spacing=1) - self.assert_mesh_is(surface=surface, expected_mesh=[ - [(0.0, 0.0, 1.0), - (0.0, 0.01, 1.0), - (0.0, 0.02, 1.0)], - - [(0.008, 0.0, 1.2), - (0.008, 0.008, 1.2), - (0.008, 0.016, 1.2)], - - [(0.016, 0.0, 1.4), - (0.016, 0.006, 1.4), - (0.016, 0.012, 1.4)], + # plot([edge1, edge2, edge3]) - [(0.016, 0.0, 1.6), - (0.016, 0.006, 1.6), - (0.016, 0.012, 1.6)], - - [(0.008, 0, 1.8), - (0.008, 0.008, 1.8), - (0.008, 0.016, 1.8)], - - [(0.0, 0.0, 2.0), - (0.0, 0.01, 2.0), - (0.0, 0.02, 2.0)], - ]) + surface = ComplexFaultSurface.from_fault_data([edge1, edge2, edge3], + mesh_spacing=1.0) + expected = numpy.array([ + [[0.00000000e+00, 0.00000000e+00, 0.00000000e+00], + [8.00000005e-03, 8.00000011e-03, 8.00000028e-03], + [1.60000001e-02, 1.60000001e-02, 1.60000002e-02], + [8.57028427e-03, 8.53924247e-03, 8.46354157e-03], + [5.70284230e-04, 5.39242373e-04, 4.63541314e-04]], + [[0.00000000e+00, 9.90000005e-03, 1.98000001e-02], + [-2.68595095e-26, 7.92000013e-03, 1.58400003e-02], + [-1.79063396e-26, 5.94000007e-03, 1.18800001e-02], + [-2.74067477e-26, 7.78653765e-03, 1.56105473e-02], + [-3.10016364e-27, 9.76653758e-03, 1.95705472e-02]], + [[1.00000000e+00, 1.00000000e+00, 1.00000000e+00], + [1.20000000e+00, 1.20000000e+00, 1.20000000e+00], + [1.40000000e+00, 1.40000000e+00, 1.40000000e+00], + [1.78574289e+00, 1.78651894e+00, 1.78841147e+00], + [1.98574289e+00, 1.98651894e+00, 1.98841147e+00]]]) + + numpy.testing.assert_almost_equal(surface.mesh.array, expected) def test_mesh_spacing_more_than_two_lengths(self): edge1 = Line([Point(0, 0, 0), Point(0, 0.1, 0)]) @@ -221,7 +235,7 @@ def test_surface_crossing_international_date_line(self): self.assert_mesh_is(surface=surface, expected_mesh=[ [(179.95, 0., 0.), (-179.95, 0., 0.)], [(179.95, 0., 10.), (-179.95, 0., 10.)] - ]) + ], delta=150) class ComplexFaultSurfaceProjectionTestCase(unittest.TestCase):