From 761b4ea856db09f406ac5f3865dbdfa449e8817b Mon Sep 17 00:00:00 2001 From: Marco Pagani Date: Mon, 10 Jun 2024 18:26:19 +0200 Subject: [PATCH 01/23] As in my laptop --- openquake/hazardlib/geo/surface/kite_fault.py | 72 +++++++++++++------ openquake/hazardlib/geo/utils.py | 8 ++- .../geo/surface/kite_fault_auxilliary_test.py | 5 +- 3 files changed, 59 insertions(+), 26 deletions(-) diff --git a/openquake/hazardlib/geo/surface/kite_fault.py b/openquake/hazardlib/geo/surface/kite_fault.py index af95498cc23f..ae1035e751b2 100644 --- a/openquake/hazardlib/geo/surface/kite_fault.py +++ b/openquake/hazardlib/geo/surface/kite_fault.py @@ -286,6 +286,8 @@ def _fix_right_hand(self): found = False irow = 0 icol = 0 + + # Search for a quadrilateral with finite vertexes while not found: if np.all(np.isfinite(self.mesh.lons[irow:irow + 2, icol:icol + 2])): @@ -297,6 +299,7 @@ def _fix_right_hand(self): icol = 1 if (irow + 1) >= self.mesh.lons.shape[0]: break + # Check strike if found: azi_strike = azimuth(self.mesh.lons[irow, icol], self.mesh.lats[irow, icol], @@ -307,7 +310,12 @@ def _fix_right_hand(self): self.mesh.lons[irow + 1, icol], self.mesh.lats[irow + 1, icol]) - if abs((azi_strike - 90) % 360 - azi_dip) < 40: + # Compare the dip direction from the strike against the one from + # the quadrilateral + # tmp = geo_utils._angles_diff(azi_strike, 90) + # if abs(geo_utils._angles_diff(tmp, azi_dip)) < 40: + tmp = (azi_strike + 90.0) % 360 + if abs(geo_utils._angles_diff(tmp, azi_dip)) > 40: tlo = np.fliplr(self.mesh.lons) tla = np.fliplr(self.mesh.lats) tde = np.fliplr(self.mesh.depths) @@ -831,32 +839,58 @@ def _dbg_plot_mesh(mesh): def _fix_right_hand(msh): # This function checks that the array describing the surface complies with - # the right hand rule. + # the right hand rule and, it required, flips the mesh to make it compliant # # :param msh: - # A :class:`numpy.ndarray` instance + # A :class:`numpy.ndarray` instance describing the mesh + # Check if the mesh complies with the right hand rule and flip it if + # required chk = _does_mesh_comply_with_right_hand_rule(msh) - - # Check if we need to flip the grid - # if not np.abs(_angles_diff(strike, top.azimuth)) < 60: if not chk: # Flip the grid to make it compliant with the right hand rule - nmsh = np.flip(msh, axis=1) - chk_flip = msh[0, 0, 0] == nmsh[0, -1, 0] + nmsh = np.empty_like(msh) + nmsh[:, :, :] = msh[:, ::-1, :] + chk_flip = ((msh[:, 0, 0] == nmsh[:, -1, 0]) & + (msh[:, 0, 2] == nmsh[:, -1, 2])) # Check again the average azimuth for the top edge of the surface msg = "The mesh still does not comply with the right hand rule" chk1 = _does_mesh_comply_with_right_hand_rule(nmsh) + chk1 = True + + # NOTE this is for debugging purposes + if True and not chk1: + _plot_mesh(nmsh) + _plot_mesh(msh) + assert chk1, msg return nmsh return msh +def _plot_mesh(nmsh): + scl = -0.01 + ax = plt.figure().add_subplot(projection='3d') + ax.set_aspect('equal') + for i in range(0, nmsh.shape[0]): + j = np.isfinite(nmsh[i, :, 0]) + plt.plot( + nmsh[i, j, 0], nmsh[i, j, 1], nmsh[i, j, 2] * scl, '-b', lw=.25) + for i in range(0, nmsh.shape[1]): + j = np.isfinite(nmsh[:, i, 0]) + plt.plot( + nmsh[j, i, 0], nmsh[j, i, 1], nmsh[j, i, 2] * scl, '-r', lw=.25) + assert np.sum(np.isfinite(nmsh[:, 0, 0])) > 0 + plt.plot(nmsh[:, 0, 0], nmsh[:, 0, 1], nmsh[:, 0, 2] * scl, '-g', lw=2.0) + plt.show() + + def _does_mesh_comply_with_right_hand_rule(msh, tolerance=45.): - # Given a mesh checks if it complies with the right hand rule + # Given a mesh, this function checks if it complies with the right hand + # rule # # :param msh: # A :class:`openquake.hazardlib.geo.mesh.Mesh` instance @@ -871,7 +905,7 @@ def _does_mesh_comply_with_right_hand_rule(msh, tolerance=45.): lr = np.isfinite(msh[1:, 1:, 0]) ia = np.nonzero(ur & ul & ll & lr) - idx = np.nonzero(msh[:, :, 0]) + # idx = np.nonzero(msh[:, :, 0]) # Coordinates of the cell with finite vertexes lons = [msh[ia[0][0], ia[1][0], 0], @@ -890,7 +924,6 @@ def _does_mesh_comply_with_right_hand_rule(msh, tolerance=45.): lats = np.array(lats) deps = np.array(deps) - # First we fit a plane through the four points and then we find the strike # of the plane _, vers = _get_plane(lons, lats, deps) @@ -898,11 +931,11 @@ def _does_mesh_comply_with_right_hand_rule(msh, tolerance=45.): # Next we find the azimuth of the top segment of the cell top = Line([Point(*msh[ia[0][0], ia[1][0], :]), - Point(*msh[ia[0][0], ia[1][0]+1, :])]) + Point(*msh[ia[0][0], ia[1][0] + 1, :])]) # Compute the difference between strike from the plane and cell top # direction - abs_angle_dff = np.abs(_angles_diff(top.azimuth, strike)) + abs_angle_dff = np.abs(geo_utils._angles_diff(top.azimuth, strike)) return abs_angle_dff < tolerance @@ -923,12 +956,7 @@ def _get_plane(lons, lats, deps): coo[:, 1] = tmp[:, 1] coo[:, 2] = deps coo[:, 2] *= -1 - return geo_utils.plane_fit(coo) - - -def _angles_diff(ang_a, ang_b): - dff = ang_a - ang_b - return (dff + 180) % 360 - 180 + return geo_utils.plane_fit(coo) def _align_profiles(prfr: list, prfl: list): @@ -1238,16 +1266,16 @@ def _dbg_plot(new_edges=None, profs=None, npr=None, ref_idx=None, if new_edges is not None: for key in new_edges: edg = np.array(new_edges[key]) - plt.plot(edg[:, 0], edg[:, 1], edg[:, 2], 'g-x', lw=1) + plt.plot(edg[:, 0], edg[:, 1], edg[:, 2], 'g-x', lw=.5) for ie, ee in enumerate(edg): ax.text(ee[0], ee[1], ee[2], s=f'{ie}') if profs is not None: for pro in profs: edg = np.array(pro) - plt.plot(edg[:, 0], edg[:, 1], edg[:, 2], 'r', lw=3) + plt.plot(edg[:, 0], edg[:, 1], edg[:, 2], 'r', lw=.5) edg = np.array(profs[ref_idx]) - plt.plot(edg[:, 0], edg[:, 1], edg[:, 2], 'y--', lw=3, zorder=100) + plt.plot(edg[:, 0], edg[:, 1], edg[:, 2], 'y--', lw=.5, zorder=100) if npr is not None: for pro in npr: diff --git a/openquake/hazardlib/geo/utils.py b/openquake/hazardlib/geo/utils.py index 378d6fdf54a2..f422b4ffcf05 100644 --- a/openquake/hazardlib/geo/utils.py +++ b/openquake/hazardlib/geo/utils.py @@ -832,7 +832,7 @@ def get_strike_from_plane_normal(nrml): nrml *= -1 # Get the strike - return numpy.rad2deg(numpy.arctan2(nrml[0], nrml[1])) - 90 + return _angles_diff(numpy.rad2deg(numpy.arctan2(nrml[0], nrml[1])), 90) def bbox2poly(bbox): @@ -935,3 +935,9 @@ def geolocate(lonlats, geom_df, exclude=()): continue codes[contains_xy(geom, lonlats)] = code return codes + +def _angles_diff(ang_a, ang_b): + # Computes the difference between the first and the second angle. Both are + # in decimal degrees. + dff = ang_a - ang_b + return (dff + 180.0) % 360.0 - 180.0 diff --git a/openquake/hazardlib/tests/geo/surface/kite_fault_auxilliary_test.py b/openquake/hazardlib/tests/geo/surface/kite_fault_auxilliary_test.py index 0b163b8abd81..17bf456e4316 100644 --- a/openquake/hazardlib/tests/geo/surface/kite_fault_auxilliary_test.py +++ b/openquake/hazardlib/tests/geo/surface/kite_fault_auxilliary_test.py @@ -17,8 +17,7 @@ import numpy as np import unittest from openquake.hazardlib.geo import utils as geo_utils -from openquake.hazardlib.geo.surface.kite_fault import ( - _angles_diff, _get_resampled_profs) +from openquake.hazardlib.geo.surface.kite_fault import _get_resampled_profs class AngleDifferenceTest(unittest.TestCase): @@ -28,7 +27,7 @@ def test_a(self): # Simple test angle_a = 30 angle_b = 60 - computed = _angles_diff(angle_a, angle_b) + computed = geo_utils._angles_diff(angle_a, angle_b) expected = -30 np.testing.assert_allclose(computed, expected) From 190d0cfc3fa0c954a20277c54ee197f3f3184849 Mon Sep 17 00:00:00 2001 From: Marco Pagani Date: Mon, 8 Jul 2024 18:29:18 +0200 Subject: [PATCH 02/23] More work --- openquake/hazardlib/geo/surface/kite_fault.py | 22 +++- openquake/hazardlib/geo/utils.py | 9 +- .../tests/geo/surface/kite_fault_test.py | 2 +- openquake/hazardlib/tests/geo/utils_test.py | 102 +++++++++++++++++- 4 files changed, 130 insertions(+), 5 deletions(-) diff --git a/openquake/hazardlib/geo/surface/kite_fault.py b/openquake/hazardlib/geo/surface/kite_fault.py index c8fb2fd21f3f..24a41a5a14f5 100644 --- a/openquake/hazardlib/geo/surface/kite_fault.py +++ b/openquake/hazardlib/geo/surface/kite_fault.py @@ -295,9 +295,10 @@ def _fix_right_hand(self): icol += 1 if (icol + 1) >= self.mesh.lons.shape[1]: irow += 1 - icol = 1 + icol = 0 if (irow + 1) >= self.mesh.lons.shape[0]: break + # Check strike if found: azi_strike = azimuth(self.mesh.lons[irow, icol], @@ -309,6 +310,25 @@ def _fix_right_hand(self): self.mesh.lons[irow + 1, icol], self.mesh.lats[irow + 1, icol]) + coo = np.empty((4, 3)) + coo[0, 0] = self.mesh.lons[irow, icol] + coo[0, 1] = self.mesh.lats[irow, icol] + coo[0, 1] = self.mesh.depths[irow, icol] + coo[1, 0] = self.mesh.lons[irow+1, icol] + coo[1, 1] = self.mesh.lats[irow+1, icol] + coo[1, 2] = self.mesh.depths[irow+1, icol] + coo[2, 0] = self.mesh.lons[irow+1, icol+1] + coo[2, 1] = self.mesh.lats[irow+1, icol+1] + coo[2, 2] = self.mesh.depths[irow+1, icol+1] + coo[3, 0] = self.mesh.lons[irow, icol+1] + coo[3, 1] = self.mesh.lats[irow, icol+1] + coo[3, 2] = self.mesh.depths[irow, icol+1] + + from openquake.hazardlib.geo.utils import ( + plane_fit, get_strike_from_plane_normal) + pnt_plane, nrml_plane = plane_fit(coo) + tmp_strike = get_strike_from_plane_normal(nrml_plane) + # Compare the dip direction from the strike against the one from # the quadrilateral # tmp = geo_utils._angles_diff(azi_strike, 90) diff --git a/openquake/hazardlib/geo/utils.py b/openquake/hazardlib/geo/utils.py index 89c64074f17b..561fd439f48f 100644 --- a/openquake/hazardlib/geo/utils.py +++ b/openquake/hazardlib/geo/utils.py @@ -840,7 +840,7 @@ def get_strike_from_plane_normal(nrml): plane. The positive z-direction is pointing upwards. :param nrml: - A vector with 3 elements + A vector with 3 elements defining the normal to the plane :returns: A float defining the strike direction """ @@ -850,7 +850,12 @@ def get_strike_from_plane_normal(nrml): nrml *= -1 # Get the strike - return _angles_diff(numpy.rad2deg(numpy.arctan2(nrml[0], nrml[1])), 90) + if nrml[1] >= 0: + tmp = (numpy.rad2deg(numpy.arctan2(nrml[0], nrml[1])) + 90.0) % 360.0 + else: + tmp = (numpy.rad2deg(numpy.arctan2(nrml[0], nrml[1])) + 270.0) % 360.0 + + return tmp def bbox2poly(bbox): diff --git a/openquake/hazardlib/tests/geo/surface/kite_fault_test.py b/openquake/hazardlib/tests/geo/surface/kite_fault_test.py index 159d56a6dd24..1771cf2e2dd2 100644 --- a/openquake/hazardlib/tests/geo/surface/kite_fault_test.py +++ b/openquake/hazardlib/tests/geo/surface/kite_fault_test.py @@ -34,7 +34,7 @@ NS = "{http://openquake.org/xmlns/nrml/0.5}" BASE_DATA_PATH = os.path.join(os.path.dirname(__file__), 'data') -PLOTTING = False +PLOTTING = True aae = np.testing.assert_almost_equal diff --git a/openquake/hazardlib/tests/geo/utils_test.py b/openquake/hazardlib/tests/geo/utils_test.py index c3a60b4d5040..870ca9a03543 100644 --- a/openquake/hazardlib/tests/geo/utils_test.py +++ b/openquake/hazardlib/tests/geo/utils_test.py @@ -21,6 +21,8 @@ from openquake.hazardlib import geo from openquake.hazardlib.geo import utils +from openquake.hazardlib.geo.utils import ( + plane_fit, get_strike_from_plane_normal) Point = collections.namedtuple("Point", 'lon lat') aac = numpy.testing.assert_allclose @@ -422,7 +424,105 @@ def test_plane_fit02(self): self.points[:, 2] += numpy.random.random(self.npts) * 0.01 pnt, par = utils.plane_fit(self.points) numpy.testing.assert_allclose(self.c[0:3], par, rtol=1e-3, atol=0) - self.assertAlmostEqual(self.c[-1], -sum(par*pnt), 2) + self.assertAlmostEqual(self.c[-1], -sum(par * pnt), 2) +class PlaneStrikeTest(unittest.TestCase): + + def test_get_strike01(self): + + coo = numpy.empty((4, 3)) + # Upper left + coo[0, :] = [10.0, 45.0, 0.0] + fnc = geo.geodetic.npoints_towards + # Lower right + tmp = fnc(coo[0, 0], coo[0, 1], coo[0, 2], 45, 20.0, 0.0, 2) + coo[1, :] = [tmp[0][1], tmp[1][1], tmp[2][1]] + # Lower left + tmp = fnc(coo[0, 0], coo[0, 1], coo[0, 2], 135, 20.0, 10.0, 2) + coo[2, :] = [tmp[0][1], tmp[1][1], tmp[2][1]] + # Lower right + tmp = fnc(coo[2, 0], coo[2, 1], coo[2, 2], 45, 20.0, 0.0, 2) + coo[3, :] = [tmp[0][1], tmp[1][1], tmp[2][1]] + + # Fit plane (without projecting the points) + ppnt, pnrm = plane_fit(coo) + + # Find strike + strike = get_strike_from_plane_normal(pnrm) + + # Test + expected = 45.0 + check = numpy.abs(strike - expected) < 15.0 + self.assertTrue(check) + + def test_get_strike02(self): + + expected = 135.0 + coo = _get_coo_poly(10., 45.0, expected) + + # Fit plane (without projecting the points) + ppnt, pnrm = plane_fit(coo) + + # Find strike + strike = get_strike_from_plane_normal(pnrm) + + # Test + check = numpy.abs(strike - expected) < 15.0 + self.assertTrue(check) + + def test_get_strike03(self): + + expected = 315.0 + coo = _get_coo_poly(10., 45.0, expected, True) + + # Fit plane (without projecting the points) + ppnt, pnrm = plane_fit(coo) + + # Find strike + strike = get_strike_from_plane_normal(pnrm) + + # Test + check = numpy.abs(strike - expected + 180.0) < 15.0 + self.assertTrue(check) + + _plot(coo) + + + +def _get_coo_poly(olon, olat, azim, opposite=False): + + coo = numpy.empty((4, 3)) + fnc = geo.geodetic.npoints_towards + + # Upper right + coo[0, :] = [olon, olat, 0.0] + + # Upper left + tmp = fnc(coo[0, 0], coo[0, 1], coo[0, 2], azim, 20.0, 0.0, 2) + coo[1, :] = [tmp[0][1], tmp[1][1], tmp[2][1]] + + # Lower right + factor = 1 if not opposite else -1 + tmp = (azim + 90.0 * factor) % 360. + tmp = fnc(coo[0, 0], coo[0, 1], coo[0, 2], tmp, 20.0, 10.0, 2) + coo[2, :] = [tmp[0][1], tmp[1][1], tmp[2][1]] + + # Lower left + tmp = fnc(coo[2, 0], coo[2, 1], coo[2, 2], azim, 20.0, 0.0, 2) + coo[3, :] = [tmp[0][1], tmp[1][1], tmp[2][1]] + + return coo + + +def _plot(coo): + + import matplotlib.pyplot as plt + ax = plt.figure().add_subplot(projection='3d') + idx = [0, 1] + plt.plot(coo[:, 0], coo[:, 1], coo[:, 2], 'o', color='red') + plt.plot(coo[idx, 0], coo[idx, 1], coo[idx, 2], '-', lw=2) + ax.invert_zaxis() + plt.show() + # NB: utils.assoc is tested in the engine From 86cd4faaeb9c1395c1da3a9ba8d7d4895ae6fd2f Mon Sep 17 00:00:00 2001 From: Kirsty Bayliss Date: Mon, 16 Sep 2024 10:47:23 +0200 Subject: [PATCH 03/23] trying to understand right hand rule problems --- openquake/hazardlib/geo/surface/kite_fault.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/openquake/hazardlib/geo/surface/kite_fault.py b/openquake/hazardlib/geo/surface/kite_fault.py index c8fb2fd21f3f..52e7db975be7 100644 --- a/openquake/hazardlib/geo/surface/kite_fault.py +++ b/openquake/hazardlib/geo/surface/kite_fault.py @@ -311,10 +311,10 @@ def _fix_right_hand(self): # Compare the dip direction from the strike against the one from # the quadrilateral - # tmp = geo_utils._angles_diff(azi_strike, 90) - # if abs(geo_utils._angles_diff(tmp, azi_dip)) < 40: - tmp = (azi_strike + 90.0) % 360 - if abs(geo_utils._angles_diff(tmp, azi_dip)) > 40: + tmp = geo_utils._angles_diff(azi_strike, 90) + if abs(geo_utils._angles_diff(tmp, azi_dip)) < 40: + #tmp = (azi_strike + 90.0) % 360 + #if abs(geo_utils._angles_diff(tmp, azi_dip)) > 40: tlo = np.fliplr(self.mesh.lons) tla = np.fliplr(self.mesh.lats) tde = np.fliplr(self.mesh.depths) From cae196f73662fabf71a44c50216f726246c51149 Mon Sep 17 00:00:00 2001 From: Kirsty Bayliss Date: Mon, 16 Sep 2024 11:18:41 +0200 Subject: [PATCH 04/23] removing unused chk_flip --- openquake/hazardlib/geo/surface/kite_fault.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/openquake/hazardlib/geo/surface/kite_fault.py b/openquake/hazardlib/geo/surface/kite_fault.py index 52e7db975be7..df8584b0ee26 100644 --- a/openquake/hazardlib/geo/surface/kite_fault.py +++ b/openquake/hazardlib/geo/surface/kite_fault.py @@ -38,7 +38,7 @@ from openquake.hazardlib.geo.geodetic import ( npoints_towards, distance, azimuth) from openquake.hazardlib.geo.surface import SimpleFaultSurface -from openquake.hazardlib.geo.surface.gridded import GriddedSurface +#from openquake.hazardlib.geo.surface.gridded import GriddedSurface TOL = 0.4 SMALL = 1e-5 @@ -835,7 +835,7 @@ def _dbg_plot_mesh(mesh): set_axes_equal(ax) plt.show() - + def _fix_right_hand(msh): # This function checks that the array describing the surface complies with # the right hand rule and, it required, flips the mesh to make it compliant @@ -851,8 +851,9 @@ def _fix_right_hand(msh): # Flip the grid to make it compliant with the right hand rule nmsh = np.empty_like(msh) nmsh[:, :, :] = msh[:, ::-1, :] - chk_flip = ((msh[:, 0, 0] == nmsh[:, -1, 0]) & - (msh[:, 0, 2] == nmsh[:, -1, 2])) + + #chk_flip = ((msh[:, 0, 0] == nmsh[:, -1, 0]) & + # (msh[:, 0, 2] == nmsh[:, -1, 2])) # Check again the average azimuth for the top edge of the surface msg = "The mesh still does not comply with the right hand rule" From c914a3da4687b8ad47092b93f2cf1a8d60500d2a Mon Sep 17 00:00:00 2001 From: Marco Pagani Date: Tue, 22 Oct 2024 08:04:46 +0200 Subject: [PATCH 05/23] fix comment --- openquake/hazardlib/geo/surface/kite_fault.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openquake/hazardlib/geo/surface/kite_fault.py b/openquake/hazardlib/geo/surface/kite_fault.py index 24a41a5a14f5..714c153642ec 100644 --- a/openquake/hazardlib/geo/surface/kite_fault.py +++ b/openquake/hazardlib/geo/surface/kite_fault.py @@ -858,7 +858,7 @@ def _dbg_plot_mesh(mesh): def _fix_right_hand(msh): # This function checks that the array describing the surface complies with - # the right hand rule and, it required, flips the mesh to make it compliant + # the right hand rule and, if required, flips the mesh to make it compliant # # :param msh: # A :class:`numpy.ndarray` instance describing the mesh From 4555eae34185214b19bcd82c486a0a5ca24f7221 Mon Sep 17 00:00:00 2001 From: Marco Pagani Date: Tue, 22 Oct 2024 12:43:08 +0200 Subject: [PATCH 06/23] More work --- openquake/hazardlib/geo/surface/kite_fault.py | 43 +++++++++++-------- .../tests/geo/surface/kite_fault_test.py | 4 ++ 2 files changed, 29 insertions(+), 18 deletions(-) diff --git a/openquake/hazardlib/geo/surface/kite_fault.py b/openquake/hazardlib/geo/surface/kite_fault.py index a70d0d1d634d..f5e7414c7671 100644 --- a/openquake/hazardlib/geo/surface/kite_fault.py +++ b/openquake/hazardlib/geo/surface/kite_fault.py @@ -301,6 +301,8 @@ def _fix_right_hand(self): # Check strike if found: + + # Get the azimuth direction for the strike and dip azi_strike = azimuth(self.mesh.lons[irow, icol], self.mesh.lats[irow, icol], self.mesh.lons[irow, icol + 1], @@ -314,27 +316,29 @@ def _fix_right_hand(self): coo[0, 0] = self.mesh.lons[irow, icol] coo[0, 1] = self.mesh.lats[irow, icol] coo[0, 1] = self.mesh.depths[irow, icol] - coo[1, 0] = self.mesh.lons[irow+1, icol] - coo[1, 1] = self.mesh.lats[irow+1, icol] - coo[1, 2] = self.mesh.depths[irow+1, icol] - coo[2, 0] = self.mesh.lons[irow+1, icol+1] - coo[2, 1] = self.mesh.lats[irow+1, icol+1] - coo[2, 2] = self.mesh.depths[irow+1, icol+1] - coo[3, 0] = self.mesh.lons[irow, icol+1] - coo[3, 1] = self.mesh.lats[irow, icol+1] - coo[3, 2] = self.mesh.depths[irow, icol+1] + coo[1, 0] = self.mesh.lons[irow + 1, icol] + coo[1, 1] = self.mesh.lats[irow + 1, icol] + coo[1, 2] = self.mesh.depths[irow + 1, icol] + coo[2, 0] = self.mesh.lons[irow + 1, icol + 1] + coo[2, 1] = self.mesh.lats[irow + 1, icol + 1] + coo[2, 2] = self.mesh.depths[irow + 1, icol + 1] + coo[3, 0] = self.mesh.lons[irow, icol + 1] + coo[3, 1] = self.mesh.lats[irow, icol + 1] + coo[3, 2] = self.mesh.depths[irow, icol + 1] from openquake.hazardlib.geo.utils import ( plane_fit, get_strike_from_plane_normal) - pnt_plane, nrml_plane = plane_fit(coo) + _, nrml_plane = plane_fit(coo) tmp_strike = get_strike_from_plane_normal(nrml_plane) + breakpoint() + # Compare the dip direction from the strike against the one from # the quadrilateral - tmp = geo_utils._angles_diff(azi_strike, 90) - if abs(geo_utils._angles_diff(tmp, azi_dip)) < 40: - #tmp = (azi_strike + 90.0) % 360 - #if abs(geo_utils._angles_diff(tmp, azi_dip)) > 40: + # tmp = geo_utils._angles_diff(azi_strike, 90) + #if abs(geo_utils._angles_diff(tmp, azi_dip)) < 40: + tmp = (azi_strike + 90) % 360 + if abs(geo_utils._angles_diff(tmp, azi_dip)) > 40: tlo = np.fliplr(self.mesh.lons) tla = np.fliplr(self.mesh.lats) tde = np.fliplr(self.mesh.depths) @@ -811,7 +815,7 @@ def _create_mesh(rprof, ref_idx, edge_sd, idl, align): msh = _fix_right_hand(msh) # INFO: this is just for debugging - # _dbg_plot_mesh(msh) + _dbg_plot_mesh(msh) return msh @@ -851,11 +855,14 @@ def _dbg_plot_mesh(mesh): ax.plot(mesh[:, :, 0].flatten(), mesh[:, :, 1].flatten(), mesh[:, :, 2].flatten() * scl, '.') - ax.invert_zaxis() + ax.plot(mesh[0, 0, 0].flatten(), + mesh[0, 0, 1].flatten(), + mesh[0, 0, 2].flatten() * scl, 'or') set_axes_equal(ax) + ax.invert_zaxis() plt.show() - + def _fix_right_hand(msh): # This function checks that the array describing the surface complies with # the right hand rule and, if required, flips the mesh to make it compliant @@ -871,7 +878,7 @@ def _fix_right_hand(msh): # Flip the grid to make it compliant with the right hand rule nmsh = np.empty_like(msh) nmsh[:, :, :] = msh[:, ::-1, :] - + #chk_flip = ((msh[:, 0, 0] == nmsh[:, -1, 0]) & # (msh[:, 0, 2] == nmsh[:, -1, 2])) diff --git a/openquake/hazardlib/tests/geo/surface/kite_fault_test.py b/openquake/hazardlib/tests/geo/surface/kite_fault_test.py index f52f6457b1eb..b475e50bd3fd 100644 --- a/openquake/hazardlib/tests/geo/surface/kite_fault_test.py +++ b/openquake/hazardlib/tests/geo/surface/kite_fault_test.py @@ -155,6 +155,8 @@ class KiteSurfaceFromMeshTest(unittest.TestCase): """ def setUp(self): + # The surface dips toward north. It doesn't comply with the right hand + # rule and must be flipped lons = [[np.nan, 0.05, 0.1, 0.15, 0.20], [0.00, 0.05, 0.1, 0.15, 0.20], [0.00, 0.05, 0.1, 0.15, 0.20], @@ -289,6 +291,8 @@ def test_get_tor(self): ax.plot(coo[:, 0], coo[:, 1], '-g', lw=4) plt.show() + breakpoint() + aae(elo, coo[:, 0]) aae(ela, coo[:, 1]) From cd0168ad352d9a77eb0e1dc0a072141b9750f465 Mon Sep 17 00:00:00 2001 From: Marco Pagani Date: Thu, 24 Oct 2024 08:05:01 +0200 Subject: [PATCH 07/23] Adding more tests --- openquake/hazardlib/geo/surface/kite_fault.py | 20 ++-- openquake/hazardlib/geo/utils.py | 17 +++ .../tests/geo/surface/kite_fault_test.py | 109 +++++++++++++++--- 3 files changed, 122 insertions(+), 24 deletions(-) diff --git a/openquake/hazardlib/geo/surface/kite_fault.py b/openquake/hazardlib/geo/surface/kite_fault.py index f5e7414c7671..fb10f9a4ca49 100644 --- a/openquake/hazardlib/geo/surface/kite_fault.py +++ b/openquake/hazardlib/geo/surface/kite_fault.py @@ -315,7 +315,7 @@ def _fix_right_hand(self): coo = np.empty((4, 3)) coo[0, 0] = self.mesh.lons[irow, icol] coo[0, 1] = self.mesh.lats[irow, icol] - coo[0, 1] = self.mesh.depths[irow, icol] + coo[0, 2] = self.mesh.depths[irow, icol] coo[1, 0] = self.mesh.lons[irow + 1, icol] coo[1, 1] = self.mesh.lats[irow + 1, icol] coo[1, 2] = self.mesh.depths[irow + 1, icol] @@ -329,16 +329,20 @@ def _fix_right_hand(self): from openquake.hazardlib.geo.utils import ( plane_fit, get_strike_from_plane_normal) _, nrml_plane = plane_fit(coo) + tmp_strike = get_strike_from_plane_normal(nrml_plane) + a_low = (tmp_strike + 10) % 360 + a_upp = (tmp_strike + 80) % 360 - breakpoint() + tmp = geo_utils.angles_within(a_low, a_upp, azi_dip) # Compare the dip direction from the strike against the one from # the quadrilateral # tmp = geo_utils._angles_diff(azi_strike, 90) - #if abs(geo_utils._angles_diff(tmp, azi_dip)) < 40: - tmp = (azi_strike + 90) % 360 - if abs(geo_utils._angles_diff(tmp, azi_dip)) > 40: + # if abs(geo_utils._angles_diff(tmp, azi_dip)) < 40: + # tmp = (azi_strike + 90) % 360 + #if abs(geo_utils._angles_diff(tmp, azi_dip)) > 40: + if not tmp: tlo = np.fliplr(self.mesh.lons) tla = np.fliplr(self.mesh.lats) tde = np.fliplr(self.mesh.depths) @@ -597,7 +601,9 @@ def geom_to_kite(geom): """ shape_y, shape_z = int(geom[1]), int(geom[2]) array = geom[3:].astype(np.float64).reshape(3, shape_y, shape_z) - return KiteSurface(RectangularMesh(*array)) + sfc = KiteSurface(RectangularMesh(*array)) + sfc._fix_right_hand() + return sfc def get_profiles_from_simple_fault_data( @@ -815,7 +821,7 @@ def _create_mesh(rprof, ref_idx, edge_sd, idl, align): msh = _fix_right_hand(msh) # INFO: this is just for debugging - _dbg_plot_mesh(msh) + # _dbg_plot_mesh(msh) return msh diff --git a/openquake/hazardlib/geo/utils.py b/openquake/hazardlib/geo/utils.py index 4cf880897cf2..fad2c4dbfe68 100644 --- a/openquake/hazardlib/geo/utils.py +++ b/openquake/hazardlib/geo/utils.py @@ -993,3 +993,20 @@ def geolocate_geometries(geometries, geom_df, exclude=()): intersecting_codes.add(code) result_codes[i] = sorted(intersecting_codes) return result_codes + + +def angles_within(a_first, a_second, angle): + """ + Check is 'angle' is within 'a_first' and 'a_second'. The interval + starts at 'a_first' and in clockwise direction goes to 'a_second' + + :param a_first: + :param a_second: + :param angle: + """ + out = numpy.zeros_like(angle) + if a_first < a_second: + out = (a_first <= angle) & (angle <= a_second) + else: + out = (a_first >= angle) & (angle >= a_second) + return out diff --git a/openquake/hazardlib/tests/geo/surface/kite_fault_test.py b/openquake/hazardlib/tests/geo/surface/kite_fault_test.py index b475e50bd3fd..baff04935149 100644 --- a/openquake/hazardlib/tests/geo/surface/kite_fault_test.py +++ b/openquake/hazardlib/tests/geo/surface/kite_fault_test.py @@ -68,7 +68,7 @@ def plot_mesh_3d(ax, smsh, zfa): smsh.mesh.depths[:, i] / zfa, '-r', lw=0.5) -def ppp(profiles: list, smsh: KiteSurface = None, title: str = '', +def ppp(profiles: list = None, smsh: KiteSurface = None, title: str = '', ax_equal=False, hold=False): """ Plots the 3D mesh @@ -86,17 +86,22 @@ def ppp(profiles: list, smsh: KiteSurface = None, title: str = '', ax = plt.figure().add_subplot(projection='3d') # Plotting original profiles - for i_pro, pro in enumerate(profiles): - coo = [[p.longitude, p.latitude, p.depth] for p in pro] - coo = np.array(coo) - ax.plot(coo[:, 0], coo[:, 1], coo[:, 2] * scl, '--g', lw=1) - ax.plot( - coo[:, 0], coo[:, 1], coo[:, 2] * scl, 'og', lw=1, markersize=3) - ax.text(coo[0, 0], coo[0, 1], coo[0, 2] * scl, s=f'{i_pro}') + if profiles is not None: + for i_pro, pro in enumerate(profiles): + coo = [[p.longitude, p.latitude, p.depth] for p in pro] + coo = np.array(coo) + ax.plot(coo[:, 0], coo[:, 1], coo[:, 2] * scl, '--g', lw=1) + ax.plot( + coo[:, 0], coo[:, 1], coo[:, 2] * scl, 'og', lw=1, ms=3) + ax.text(coo[0, 0], coo[0, 1], coo[0, 2] * scl, s=f'{i_pro}') # Plotting mesh if smsh is not None: + coo = smsh.tor.coo + ax.plot(coo[:, 0], coo[:, 1], coo[:, 2] * scl, '--', color='cyan', lw=2) + ax.plot(coo[0, 0], coo[0, 1], coo[0, 2] * scl, 'o', mfc='red') + # Plotting nodes idx = np.isfinite(smsh.mesh.lons) ax.plot(smsh.mesh.lons[idx].flatten(), @@ -137,6 +142,8 @@ def ppp(profiles: list, smsh: KiteSurface = None, title: str = '', smsh.mesh.depths[:, i] * scl, '-r', lw=0.5) """ + plt.xlabel('Longitude') + plt.ylabel('Latitude') plt.title(title) if ax_equal: @@ -151,7 +158,10 @@ def ppp(profiles: list, smsh: KiteSurface = None, title: str = '', class KiteSurfaceFromMeshTest(unittest.TestCase): """ - Tests the method that creates the external boundary of the rupture. + Test the method that creates the external boundary of the rupture. The + surface of the fault dips toward north. + + [x] test right hand / tor """ def setUp(self): @@ -199,8 +209,16 @@ def test_get_external_boundary(self): ax.invert_yaxis() plt.show() - def test_get_dip1(self): - self.ksfc.get_dip() + if PLOTTING: + title = 'Surface from Mesh' + ppp(None, self.ksfc, title=title, ax_equal=True) + + def test_get_strike(self): + strike = self.ksfc.get_strike() + msg = "The value of strike computed is wrong.\n" + msg += "computed: {strike:.3f} expected:" + self.assertTrue(abs(strike - 270) < 0.01, msg) + def test_get_cell_dimensions(self): ksfc = self.ksfc @@ -233,7 +251,9 @@ def test_get_tor(self): aae(coo[:, 1], [0.0, 0.0, 0.05]) def test_geom(self): + # Converts a KiteSurface to a numpy.ndarray instance geom = kite_to_geom(self.ksfc) + # Converts a numpy.ndarray instance back to a KiteSurface ksfc = geom_to_kite(geom) for par in ('lons', 'lats', 'depths'): orig = getattr(self.ksfc.mesh, par) @@ -246,6 +266,11 @@ class KiteSurfaceWithNaNs(unittest.TestCase): Test the creation of a surface which will contain NaNs. The :method:`openquake.hazardlib.geo.surface.kite_fault.Kite.KiteSurface._clean` removes rows and cols just containing NaNs. + + The profiles used to create this surface dip toward the south quadrant + hence the dip should also point in that direction. + + [x] test right hand """ NAME = 'KiteSurfaceWithNaNs' @@ -291,8 +316,6 @@ def test_get_tor(self): ax.plot(coo[:, 0], coo[:, 1], '-g', lw=4) plt.show() - breakpoint() - aae(elo, coo[:, 0]) aae(ela, coo[:, 1]) @@ -350,7 +373,7 @@ def test_rx_calculation(self): z = np.reshape(dst, self.mlons.shape) cs = plt.contour(self.mlons, self.mlats, z, 10, colors='k') _ = plt.clabel(cs) - coo = self.srfc.tor + coo = self.srfc.tor.coo ax.plot(coo[:, 0], coo[:, 1], '-g', lw=4) plt.title(f'{self.NAME} - Rx') plt.show() @@ -379,6 +402,10 @@ def test_get_dip2(self): class KiteSurfaceUCF1Tests(unittest.TestCase): + """ + Test the construction of a UCF3 surface + [ ] test right hand + """ def setUp(self): path = os.path.join(BASE_DATA_PATH, 'profiles10') @@ -402,6 +429,10 @@ def test_mesh_creationA(self): class KiteSurfaceUCF2Tests(unittest.TestCase): + """ + Test the construction of a UCF3 surface + [ ] test right hand + """ def setUp(self): path = os.path.join(BASE_DATA_PATH, 'profiles11') @@ -456,6 +487,7 @@ def set_axes_equal(ax): class KiteSurfaceSimpleTests(unittest.TestCase): """ Simple test for the creation of a KiteSurface + [ ] test right hand """ def setUp(self): @@ -520,7 +552,10 @@ def test_compute_joyner_boore_distance(self): class KinkedKiteSurfaceTestCase(unittest.TestCase): - """ Test the construction of a kinked kite fault surface. """ + """ + Test the construction of a kinked kite fault surface. + [ ] test right hand + """ def setUp(self): """ This creates a fault dipping to north """ @@ -565,7 +600,10 @@ def test_build_kinked_mesh_01(self): class KiteSurfaceTestCase(unittest.TestCase): - """ Test the construction of a Kite fault surface. """ + """ + Test the construction of a Kite fault surface + [ ] test right hand + """ def setUp(self): """ This creates a fault dipping to north """ @@ -669,13 +707,29 @@ class IdealisedSimpleMeshTest(unittest.TestCase): """ This is the simplest test implemented for the construction of the mesh. It uses just two parallel profiles gently dipping northward and it checks - that the size of the cells agrees with the input parameters + that the size of the cells agrees with the input parameters. The fault + surface dips to the north. + + [x] test right hand """ def setUp(self): path = os.path.join(BASE_DATA_PATH, 'profiles05') self.prf, _ = _read_profiles(path) + def test_right_hand(self): + # Create the mesh: two parallel profiles - no top alignment + hsmpl = 4 + vsmpl = 4 + idl = False + alg = False + srfc = KiteSurface.from_profiles(self.prf, hsmpl, vsmpl, idl, alg) + smsh = srfc.mesh + + expected = 270.002023 + strike = srfc.get_strike() + np.testing.assert_allclose([expected], [strike]) + def test_mesh_creationD(self): # Create the mesh: two parallel profiles - no top alignment hsmpl = 4 @@ -723,6 +777,8 @@ class IdealisedSimpleDisalignedMeshTest(unittest.TestCase): Similar to :class:`openquake.sub.tests.misc.mesh_test.IdealisedSimpleMeshTest` but with profiles at different depths + + [x] test right hand / tor """ def setUp(self): @@ -737,6 +793,11 @@ def setUp(self): self.smsh = KiteSurface.from_profiles(self.profiles, self.v_sampl, self.h_sampl, idl, alg) + def test_right_hand(self): + expected = 243.6459081 + strike = self.smsh.get_strike() + np.testing.assert_allclose([expected], [strike]) + def test_h_spacing(self): # Check h-spacing: two misaligned profiles - no top alignment @@ -786,6 +847,8 @@ class IdealisedAsimmetricMeshTest(unittest.TestCase): """ Tests the creation of a surface using profiles not 'aligned' at the top (i.e. they do not start at the same depth) and with different lenghts + + [x] test right hand / tor """ def setUp(self): @@ -847,6 +910,18 @@ def test_get_width(self): width = srfc.get_width() np.testing.assert_almost_equal(37.2501538, width) + def test_get_tor(self): + # test calculation of trace (i.e. surface projection of tor) + h_sampl = 2.5 + v_sampl = 2.5 + idl = False + alg = False + srfc = KiteSurface.from_profiles(self.profiles, v_sampl, h_sampl, + idl, alg) + if PLOTTING: + title = 'TOR' + ppp(self.profiles, srfc, title) + class IdealizedATest(unittest.TestCase): """ From 6157ca95566c9f4c73da7ea89f4d19aeee7511ca Mon Sep 17 00:00:00 2001 From: Marco Pagani Date: Thu, 24 Oct 2024 08:15:52 +0200 Subject: [PATCH 08/23] Removing unused variables --- openquake/hazardlib/geo/surface/kite_fault.py | 4 ---- openquake/hazardlib/tests/geo/surface/kite_fault_test.py | 1 - openquake/hazardlib/tests/geo/utils_test.py | 6 +++--- 3 files changed, 3 insertions(+), 8 deletions(-) diff --git a/openquake/hazardlib/geo/surface/kite_fault.py b/openquake/hazardlib/geo/surface/kite_fault.py index fb10f9a4ca49..797a907c0f60 100644 --- a/openquake/hazardlib/geo/surface/kite_fault.py +++ b/openquake/hazardlib/geo/surface/kite_fault.py @@ -303,10 +303,6 @@ def _fix_right_hand(self): if found: # Get the azimuth direction for the strike and dip - azi_strike = azimuth(self.mesh.lons[irow, icol], - self.mesh.lats[irow, icol], - self.mesh.lons[irow, icol + 1], - self.mesh.lats[irow, icol + 1]) azi_dip = azimuth(self.mesh.lons[irow, icol], self.mesh.lats[irow, icol], self.mesh.lons[irow + 1, icol], diff --git a/openquake/hazardlib/tests/geo/surface/kite_fault_test.py b/openquake/hazardlib/tests/geo/surface/kite_fault_test.py index baff04935149..6bb4b7124661 100644 --- a/openquake/hazardlib/tests/geo/surface/kite_fault_test.py +++ b/openquake/hazardlib/tests/geo/surface/kite_fault_test.py @@ -724,7 +724,6 @@ def test_right_hand(self): idl = False alg = False srfc = KiteSurface.from_profiles(self.prf, hsmpl, vsmpl, idl, alg) - smsh = srfc.mesh expected = 270.002023 strike = srfc.get_strike() diff --git a/openquake/hazardlib/tests/geo/utils_test.py b/openquake/hazardlib/tests/geo/utils_test.py index 870ca9a03543..f59ade474d1f 100644 --- a/openquake/hazardlib/tests/geo/utils_test.py +++ b/openquake/hazardlib/tests/geo/utils_test.py @@ -446,7 +446,7 @@ def test_get_strike01(self): coo[3, :] = [tmp[0][1], tmp[1][1], tmp[2][1]] # Fit plane (without projecting the points) - ppnt, pnrm = plane_fit(coo) + _, pnrm = plane_fit(coo) # Find strike strike = get_strike_from_plane_normal(pnrm) @@ -462,7 +462,7 @@ def test_get_strike02(self): coo = _get_coo_poly(10., 45.0, expected) # Fit plane (without projecting the points) - ppnt, pnrm = plane_fit(coo) + _, pnrm = plane_fit(coo) # Find strike strike = get_strike_from_plane_normal(pnrm) @@ -477,7 +477,7 @@ def test_get_strike03(self): coo = _get_coo_poly(10., 45.0, expected, True) # Fit plane (without projecting the points) - ppnt, pnrm = plane_fit(coo) + _, pnrm = plane_fit(coo) # Find strike strike = get_strike_from_plane_normal(pnrm) From f76b7c581ec384795400024cb251ef905b69b844 Mon Sep 17 00:00:00 2001 From: Marco Pagani Date: Thu, 24 Oct 2024 14:47:57 +0200 Subject: [PATCH 09/23] Fixing tests - In the QA 29 the polygon vtx go in the opposite dir --- openquake/calculators/tests/classical_test.py | 2 +- openquake/hazardlib/tests/source/multi_fault_test.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/openquake/calculators/tests/classical_test.py b/openquake/calculators/tests/classical_test.py index 1e306635da69..97fca7fd9eb1 100644 --- a/openquake/calculators/tests/classical_test.py +++ b/openquake/calculators/tests/classical_test.py @@ -283,7 +283,7 @@ def test_case_29(self): # non parametric source with 2 KiteSurfaces # check what QGIS will be seeing aw = extract(self.calc.datastore, 'rupture_info') poly = gzip.decompress(aw.boundaries).decode('ascii') - expected = '''POLYGON((0.17986 -0.00000, 0.13490 -0.00000, 0.08993 -0.00000, 0.04497 -0.00000, 0.00000 0.00000, 0.00000 0.04101, 0.00000 0.08202, 0.00000 0.12303, 0.00000 0.16404, 0.00000 0.20505, 0.00000 0.24606, 0.00000 0.28708, 0.04497 0.28708, 0.08993 0.28708, 0.13490 0.28708, 0.17987 0.28708, 0.17987 0.24606, 0.17987 0.20505, 0.17987 0.16404, 0.17986 0.12303, 0.17986 0.08202, 0.17986 0.04101, 0.17986 -0.00000, 0.17986 0.10000, 0.13490 0.10000, 0.08993 0.10000, 0.04497 0.10000, 0.00000 0.10000, 0.00000 0.14101, 0.00000 0.18202, 0.00000 0.22303, 0.00000 0.26404, 0.00000 0.30505, 0.00000 0.34606, 0.00000 0.38708, 0.04497 0.38708, 0.08993 0.38708, 0.13490 0.38708, 0.17987 0.38708, 0.17987 0.34606, 0.17987 0.30505, 0.17987 0.26404, 0.17987 0.22303, 0.17987 0.18202, 0.17986 0.14101, 0.17986 0.10000))''' + expected = '''POLYGON((0.00000 0.00000, 0.04497 -0.00000, 0.08993 -0.00000, 0.13490 -0.00000, 0.17986 -0.00000, 0.17986 0.04101, 0.17986 0.08202, 0.17986 0.12303, 0.17987 0.16404, 0.17987 0.20505, 0.17987 0.24606, 0.17987 0.28708, 0.13490 0.28708, 0.08993 0.28708, 0.04497 0.28708, 0.00000 0.28708, 0.00000 0.24606, 0.00000 0.20505, 0.00000 0.16404, 0.00000 0.12303, 0.00000 0.08202, 0.00000 0.04101, 0.00000 0.00000, 0.00000 0.10000, 0.04497 0.10000, 0.08993 0.10000, 0.13490 0.10000, 0.17986 0.10000, 0.17986 0.14101, 0.17987 0.18202, 0.17987 0.22303, 0.17987 0.26404, 0.17987 0.30505, 0.17987 0.34606, 0.17987 0.38708, 0.13490 0.38708, 0.08993 0.38708, 0.04497 0.38708, 0.00000 0.38708, 0.00000 0.34606, 0.00000 0.30505, 0.00000 0.26404, 0.00000 0.22303, 0.00000 0.18202, 0.00000 0.14101, 0.00000 0.10000))''' self.assertEqual(poly, expected) # This is for checking purposes. It creates a .txt file that can be diff --git a/openquake/hazardlib/tests/source/multi_fault_test.py b/openquake/hazardlib/tests/source/multi_fault_test.py index be0eb8aeb9b0..1eb7799e72ca 100644 --- a/openquake/hazardlib/tests/source/multi_fault_test.py +++ b/openquake/hazardlib/tests/source/multi_fault_test.py @@ -195,7 +195,7 @@ def test_slow(self): if col == 'probs_occur:2': continue print(col) - aac(df[col].to_numpy(), ctx[col], rtol=1E-5, equal_nan=1) + aac(df[col].to_numpy(), ctx[col], rtol=2E-5, equal_nan=1) def main100sites(): From 2326e02c69e236d227307cfa01d7bb909ae4e029 Mon Sep 17 00:00:00 2001 From: Marco Pagani Date: Fri, 25 Oct 2024 06:41:03 +0200 Subject: [PATCH 10/23] Fixing test --- .../qa_tests_data/classical/case_75/expected/context.org | 6 +++--- .../classical/case_75/expected/hcurve-mean.csv | 4 ++-- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/openquake/qa_tests_data/classical/case_75/expected/context.org b/openquake/qa_tests_data/classical/case_75/expected/context.org index 9ca04d670672..acc8687f80b0 100644 --- a/openquake/qa_tests_data/classical/case_75/expected/context.org +++ b/openquake/qa_tests_data/classical/case_75/expected/context.org @@ -1,3 +1,3 @@ -| rup_id | clat | clon | dip | grp_id | mag | occurrence_rate | probs_occur | rake | rjb | rrup | rx | ry0 | sids | src_id | weight | width | ztor | -|--------+--------+------+---------+--------+--------+-----------------+---------------+----------+--------+--------+--------+--------+------+--------+--------+---------+------| -| 0 | 0.0899 | 0.0 | 71.5521 | 0 | 6.0000 | NaN | 0.9600 0.0400 | 180.0000 | 6.6792 | 7.0397 | 2.2239 | 6.6792 | 0 | 0 | NaN | 10.0000 | 0.0 | \ No newline at end of file +| rup_id | clat | clon | dip | grp_id | mag | occurrence_rate | probs_occur | rake | rjb | rrup | rx | ry0 | sids | src_id | weight | width | ztor | +|--------+--------+------+---------+--------+--------+-----------------+---------------+----------+--------+--------+------------+--------+------+--------+--------+---------+------| +| 0 | 0.0899 | 0.0 | 71.5521 | 0 | 6.0000 | NaN | 0.9600 0.0400 | 180.0000 | 6.6792 | 7.0397 | -2.2238908 | 6.6792 | 0 | 0 | NaN | 10.0000 | 0.0 | \ No newline at end of file diff --git a/openquake/qa_tests_data/classical/case_75/expected/hcurve-mean.csv b/openquake/qa_tests_data/classical/case_75/expected/hcurve-mean.csv index 1eb2ad0f454a..0fcdf0b8a2ce 100644 --- a/openquake/qa_tests_data/classical/case_75/expected/hcurve-mean.csv +++ b/openquake/qa_tests_data/classical/case_75/expected/hcurve-mean.csv @@ -1,3 +1,3 @@ -#,,,,,,,,,,,,,,,,,,,,,,"generated_by='OpenQuake engine 3.21.0-git7ca6672329', start_date='2024-07-02T07:06:10', checksum=3814484434, kind='mean', investigation_time=1.0, imt='PGA'" +#,,,,,,,,,,,,,,,,,,,,,,"generated_by='OpenQuake engine 3.22.0-gitf76b7c581e', start_date='2024-10-24T16:46:01', checksum=3814484434, kind='mean', investigation_time=1.0, imt='PGA'" lon,lat,depth,poe-0.0050000,poe-0.0070015,poe-0.0098041,poe-0.0137286,poe-0.0192240,poe-0.0269192,poe-0.0376948,poe-0.0527836,poe-0.0739125,poe-0.1034991,poe-0.1449289,poe-0.2029427,poe-0.2841789,poe-0.3979333,poe-0.5572227,poe-0.7802743,poe-1.0926116,poe-1.5299748,poe-2.1424109,poe-3.0000000 -0.02000,0.15000,0.00000,1.677307E-01,1.677307E-01,1.677307E-01,1.677307E-01,1.677307E-01,1.676733E-01,1.672826E-01,1.657682E-01,1.615646E-01,1.523358E-01,1.359136E-01,1.119633E-01,8.321244E-02,5.474294E-02,3.144711E-02,1.561661E-02,6.615418E-03,2.333290E-03,6.587651E-04,1.411206E-04 +0.02000,0.15000,0.00000,1.677307E-01,1.677307E-01,1.677307E-01,1.677307E-01,1.677307E-01,1.676722E-01,1.672764E-01,1.657401E-01,1.614632E-01,1.520464E-01,1.352280E-01,1.106391E-01,8.111322E-02,5.203515E-02,2.862316E-02,1.325065E-02,5.035818E-03,1.502216E-03,3.221631E-04,3.778934E-05 From 48ccaebdd93e63675db115774a7a44bb58408ed3 Mon Sep 17 00:00:00 2001 From: Marco Pagani Date: Tue, 5 Nov 2024 22:45:42 +0900 Subject: [PATCH 11/23] Cleaning --- openquake/hazardlib/geo/surface/kite_fault.py | 11 ++--------- 1 file changed, 2 insertions(+), 9 deletions(-) diff --git a/openquake/hazardlib/geo/surface/kite_fault.py b/openquake/hazardlib/geo/surface/kite_fault.py index 797a907c0f60..f6b127ab67a7 100644 --- a/openquake/hazardlib/geo/surface/kite_fault.py +++ b/openquake/hazardlib/geo/surface/kite_fault.py @@ -330,14 +330,8 @@ def _fix_right_hand(self): a_low = (tmp_strike + 10) % 360 a_upp = (tmp_strike + 80) % 360 + # Check that the dip direction is within a range tmp = geo_utils.angles_within(a_low, a_upp, azi_dip) - - # Compare the dip direction from the strike against the one from - # the quadrilateral - # tmp = geo_utils._angles_diff(azi_strike, 90) - # if abs(geo_utils._angles_diff(tmp, azi_dip)) < 40: - # tmp = (azi_strike + 90) % 360 - #if abs(geo_utils._angles_diff(tmp, azi_dip)) > 40: if not tmp: tlo = np.fliplr(self.mesh.lons) tla = np.fliplr(self.mesh.lats) @@ -887,10 +881,9 @@ def _fix_right_hand(msh): # Check again the average azimuth for the top edge of the surface msg = "The mesh still does not comply with the right hand rule" chk1 = _does_mesh_comply_with_right_hand_rule(nmsh) - chk1 = True # NOTE this is for debugging purposes - if True and not chk1: + if False and not chk1: _plot_mesh(nmsh) _plot_mesh(msh) From 13e3055cd373f95f69a5730a32a31fdd2fff00cf Mon Sep 17 00:00:00 2001 From: Marco Pagani Date: Wed, 6 Nov 2024 20:18:33 +0900 Subject: [PATCH 12/23] Fixing tests --- openquake/calculators/tests/classical_test.py | 2 +- openquake/hazardlib/geo/surface/kite_fault.py | 63 +++++++++++-------- .../classical/case_75/expected/context.org | 6 +- .../case_75/expected/hcurve-mean.csv | 4 +- 4 files changed, 44 insertions(+), 31 deletions(-) diff --git a/openquake/calculators/tests/classical_test.py b/openquake/calculators/tests/classical_test.py index 97fca7fd9eb1..1e306635da69 100644 --- a/openquake/calculators/tests/classical_test.py +++ b/openquake/calculators/tests/classical_test.py @@ -283,7 +283,7 @@ def test_case_29(self): # non parametric source with 2 KiteSurfaces # check what QGIS will be seeing aw = extract(self.calc.datastore, 'rupture_info') poly = gzip.decompress(aw.boundaries).decode('ascii') - expected = '''POLYGON((0.00000 0.00000, 0.04497 -0.00000, 0.08993 -0.00000, 0.13490 -0.00000, 0.17986 -0.00000, 0.17986 0.04101, 0.17986 0.08202, 0.17986 0.12303, 0.17987 0.16404, 0.17987 0.20505, 0.17987 0.24606, 0.17987 0.28708, 0.13490 0.28708, 0.08993 0.28708, 0.04497 0.28708, 0.00000 0.28708, 0.00000 0.24606, 0.00000 0.20505, 0.00000 0.16404, 0.00000 0.12303, 0.00000 0.08202, 0.00000 0.04101, 0.00000 0.00000, 0.00000 0.10000, 0.04497 0.10000, 0.08993 0.10000, 0.13490 0.10000, 0.17986 0.10000, 0.17986 0.14101, 0.17987 0.18202, 0.17987 0.22303, 0.17987 0.26404, 0.17987 0.30505, 0.17987 0.34606, 0.17987 0.38708, 0.13490 0.38708, 0.08993 0.38708, 0.04497 0.38708, 0.00000 0.38708, 0.00000 0.34606, 0.00000 0.30505, 0.00000 0.26404, 0.00000 0.22303, 0.00000 0.18202, 0.00000 0.14101, 0.00000 0.10000))''' + expected = '''POLYGON((0.17986 -0.00000, 0.13490 -0.00000, 0.08993 -0.00000, 0.04497 -0.00000, 0.00000 0.00000, 0.00000 0.04101, 0.00000 0.08202, 0.00000 0.12303, 0.00000 0.16404, 0.00000 0.20505, 0.00000 0.24606, 0.00000 0.28708, 0.04497 0.28708, 0.08993 0.28708, 0.13490 0.28708, 0.17987 0.28708, 0.17987 0.24606, 0.17987 0.20505, 0.17987 0.16404, 0.17986 0.12303, 0.17986 0.08202, 0.17986 0.04101, 0.17986 -0.00000, 0.17986 0.10000, 0.13490 0.10000, 0.08993 0.10000, 0.04497 0.10000, 0.00000 0.10000, 0.00000 0.14101, 0.00000 0.18202, 0.00000 0.22303, 0.00000 0.26404, 0.00000 0.30505, 0.00000 0.34606, 0.00000 0.38708, 0.04497 0.38708, 0.08993 0.38708, 0.13490 0.38708, 0.17987 0.38708, 0.17987 0.34606, 0.17987 0.30505, 0.17987 0.26404, 0.17987 0.22303, 0.17987 0.18202, 0.17986 0.14101, 0.17986 0.10000))''' self.assertEqual(poly, expected) # This is for checking purposes. It creates a .txt file that can be diff --git a/openquake/hazardlib/geo/surface/kite_fault.py b/openquake/hazardlib/geo/surface/kite_fault.py index f6b127ab67a7..eb9de2b01076 100644 --- a/openquake/hazardlib/geo/surface/kite_fault.py +++ b/openquake/hazardlib/geo/surface/kite_fault.py @@ -866,26 +866,33 @@ def _fix_right_hand(msh): # :param msh: # A :class:`numpy.ndarray` instance describing the mesh + # Fit a plane through the four points of a non-null cell + lons, lats, deps = _get_non_null_cell(msh) + _, vers = _get_plane(lons, lats, deps) + + # Check if the plane is vertical + tmp = np.cross(vers, [0, 0, 1]) + if np.abs(tmp[2] < 0.01): + return msh + # Check if the mesh complies with the right hand rule and flip it if # required chk = _does_mesh_comply_with_right_hand_rule(msh) + if not chk: # Flip the grid to make it compliant with the right hand rule nmsh = np.empty_like(msh) nmsh[:, :, :] = msh[:, ::-1, :] - #chk_flip = ((msh[:, 0, 0] == nmsh[:, -1, 0]) & - # (msh[:, 0, 2] == nmsh[:, -1, 2])) - # Check again the average azimuth for the top edge of the surface msg = "The mesh still does not comply with the right hand rule" chk1 = _does_mesh_comply_with_right_hand_rule(nmsh) # NOTE this is for debugging purposes - if False and not chk1: - _plot_mesh(nmsh) - _plot_mesh(msh) + if True and not chk1: + _plot_mesh(nmsh, 'New Mesh') + _plot_mesh(msh, 'Old Mesh') assert chk1, msg return nmsh @@ -893,7 +900,7 @@ def _fix_right_hand(msh): return msh -def _plot_mesh(nmsh): +def _plot_mesh(nmsh, title=''): scl = -0.01 ax = plt.figure().add_subplot(projection='3d') ax.set_aspect('equal') @@ -907,10 +914,11 @@ def _plot_mesh(nmsh): nmsh[j, i, 0], nmsh[j, i, 1], nmsh[j, i, 2] * scl, '-r', lw=.25) assert np.sum(np.isfinite(nmsh[:, 0, 0])) > 0 plt.plot(nmsh[:, 0, 0], nmsh[:, 0, 1], nmsh[:, 0, 2] * scl, '-g', lw=2.0) + plt.title(title) plt.show() -def _does_mesh_comply_with_right_hand_rule(msh, tolerance=45.): +def _does_mesh_comply_with_right_hand_rule(msh, tolerance=45., vers=None): # Given a mesh, this function checks if it complies with the right hand # rule # @@ -921,14 +929,33 @@ def _does_mesh_comply_with_right_hand_rule(msh, tolerance=45.): # the top edge of the selected cell and the strike of the plane passingg # through the cell. + if vers is None: + lons, lats, deps = _get_non_null_cell(msh) + # Fit a plane through the four points + _, vers = _get_plane(lons, lats, deps) + + # Find the strike + strike = geo_utils.get_strike_from_plane_normal(vers) + + # Next we find the azimuth of the top segment of the cell + top = Line([Point(*msh[ia[0][0], ia[1][0], :]), + Point(*msh[ia[0][0], ia[1][0] + 1, :])]) + + # Compute the difference between strike from the plane and cell top + # direction + abs_angle_dff = np.abs(geo_utils._angles_diff(top.azimuth, strike)) + + return abs_angle_dff < tolerance + + +def _get_non_null_cell(msh): + ul = np.isfinite(msh[:-1, :-1, 0]) ur = np.isfinite(msh[:-1, 1:, 0]) ll = np.isfinite(msh[1:, :-1, 0]) lr = np.isfinite(msh[1:, 1:, 0]) ia = np.nonzero(ur & ul & ll & lr) - # idx = np.nonzero(msh[:, :, 0]) - # Coordinates of the cell with finite vertexes lons = [msh[ia[0][0], ia[1][0], 0], msh[ia[0][0], ia[1][0] + 1, 0], @@ -946,21 +973,7 @@ def _does_mesh_comply_with_right_hand_rule(msh, tolerance=45.): lats = np.array(lats) deps = np.array(deps) - # First we fit a plane through the four points and then we find the strike - # of the plane - _, vers = _get_plane(lons, lats, deps) - strike = geo_utils.get_strike_from_plane_normal(vers) - - # Next we find the azimuth of the top segment of the cell - top = Line([Point(*msh[ia[0][0], ia[1][0], :]), - Point(*msh[ia[0][0], ia[1][0] + 1, :])]) - - # Compute the difference between strike from the plane and cell top - # direction - abs_angle_dff = np.abs(geo_utils._angles_diff(top.azimuth, strike)) - - return abs_angle_dff < tolerance - + return lons, lats, deps def _get_plane(lons, lats, deps): diff --git a/openquake/qa_tests_data/classical/case_75/expected/context.org b/openquake/qa_tests_data/classical/case_75/expected/context.org index acc8687f80b0..9ca04d670672 100644 --- a/openquake/qa_tests_data/classical/case_75/expected/context.org +++ b/openquake/qa_tests_data/classical/case_75/expected/context.org @@ -1,3 +1,3 @@ -| rup_id | clat | clon | dip | grp_id | mag | occurrence_rate | probs_occur | rake | rjb | rrup | rx | ry0 | sids | src_id | weight | width | ztor | -|--------+--------+------+---------+--------+--------+-----------------+---------------+----------+--------+--------+------------+--------+------+--------+--------+---------+------| -| 0 | 0.0899 | 0.0 | 71.5521 | 0 | 6.0000 | NaN | 0.9600 0.0400 | 180.0000 | 6.6792 | 7.0397 | -2.2238908 | 6.6792 | 0 | 0 | NaN | 10.0000 | 0.0 | \ No newline at end of file +| rup_id | clat | clon | dip | grp_id | mag | occurrence_rate | probs_occur | rake | rjb | rrup | rx | ry0 | sids | src_id | weight | width | ztor | +|--------+--------+------+---------+--------+--------+-----------------+---------------+----------+--------+--------+--------+--------+------+--------+--------+---------+------| +| 0 | 0.0899 | 0.0 | 71.5521 | 0 | 6.0000 | NaN | 0.9600 0.0400 | 180.0000 | 6.6792 | 7.0397 | 2.2239 | 6.6792 | 0 | 0 | NaN | 10.0000 | 0.0 | \ No newline at end of file diff --git a/openquake/qa_tests_data/classical/case_75/expected/hcurve-mean.csv b/openquake/qa_tests_data/classical/case_75/expected/hcurve-mean.csv index 0fcdf0b8a2ce..9e0ebdc0bc3f 100644 --- a/openquake/qa_tests_data/classical/case_75/expected/hcurve-mean.csv +++ b/openquake/qa_tests_data/classical/case_75/expected/hcurve-mean.csv @@ -1,3 +1,3 @@ -#,,,,,,,,,,,,,,,,,,,,,,"generated_by='OpenQuake engine 3.22.0-gitf76b7c581e', start_date='2024-10-24T16:46:01', checksum=3814484434, kind='mean', investigation_time=1.0, imt='PGA'" +#,,,,,,,,,,,,,,,,,,,,,,"generated_by='OpenQuake engine 3.22.0-git48ccaebdd9', start_date='2024-11-06T20:14:54', checksum=3814484434, kind='mean', investigation_time=1.0, imt='PGA'" lon,lat,depth,poe-0.0050000,poe-0.0070015,poe-0.0098041,poe-0.0137286,poe-0.0192240,poe-0.0269192,poe-0.0376948,poe-0.0527836,poe-0.0739125,poe-0.1034991,poe-0.1449289,poe-0.2029427,poe-0.2841789,poe-0.3979333,poe-0.5572227,poe-0.7802743,poe-1.0926116,poe-1.5299748,poe-2.1424109,poe-3.0000000 -0.02000,0.15000,0.00000,1.677307E-01,1.677307E-01,1.677307E-01,1.677307E-01,1.677307E-01,1.676722E-01,1.672764E-01,1.657401E-01,1.614632E-01,1.520464E-01,1.352280E-01,1.106391E-01,8.111322E-02,5.203515E-02,2.862316E-02,1.325065E-02,5.035818E-03,1.502216E-03,3.221631E-04,3.778934E-05 +0.02000,0.15000,0.00000,1.677307E-01,1.677307E-01,1.677307E-01,1.677307E-01,1.677307E-01,1.676733E-01,1.672826E-01,1.657681E-01,1.615646E-01,1.523358E-01,1.359137E-01,1.119633E-01,8.321244E-02,5.474293E-02,3.144711E-02,1.561660E-02,6.615460E-03,2.333283E-03,6.587505E-04,1.410842E-04 From 61d542df25fe4c5581c09c15111f76e4ed3920bb Mon Sep 17 00:00:00 2001 From: Marco Pagani Date: Wed, 6 Nov 2024 21:42:53 +0900 Subject: [PATCH 13/23] Still fixing --- openquake/hazardlib/geo/surface/kite_fault.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/openquake/hazardlib/geo/surface/kite_fault.py b/openquake/hazardlib/geo/surface/kite_fault.py index eb9de2b01076..85506f60953c 100644 --- a/openquake/hazardlib/geo/surface/kite_fault.py +++ b/openquake/hazardlib/geo/surface/kite_fault.py @@ -867,7 +867,7 @@ def _fix_right_hand(msh): # A :class:`numpy.ndarray` instance describing the mesh # Fit a plane through the four points of a non-null cell - lons, lats, deps = _get_non_null_cell(msh) + lons, lats, deps, ia = _get_non_null_cell(msh) _, vers = _get_plane(lons, lats, deps) # Check if the plane is vertical @@ -930,7 +930,7 @@ def _does_mesh_comply_with_right_hand_rule(msh, tolerance=45., vers=None): # through the cell. if vers is None: - lons, lats, deps = _get_non_null_cell(msh) + lons, lats, deps, ia = _get_non_null_cell(msh) # Fit a plane through the four points _, vers = _get_plane(lons, lats, deps) @@ -973,7 +973,7 @@ def _get_non_null_cell(msh): lats = np.array(lats) deps = np.array(deps) - return lons, lats, deps + return lons, lats, deps, ia def _get_plane(lons, lats, deps): From 4638aeb7c58004c05ac58e8e96725a35c98bfa40 Mon Sep 17 00:00:00 2001 From: Marco Pagani Date: Wed, 6 Nov 2024 21:49:14 +0900 Subject: [PATCH 14/23] Still fixing --- openquake/hazardlib/geo/surface/kite_fault.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openquake/hazardlib/geo/surface/kite_fault.py b/openquake/hazardlib/geo/surface/kite_fault.py index 85506f60953c..4799a7f6200c 100644 --- a/openquake/hazardlib/geo/surface/kite_fault.py +++ b/openquake/hazardlib/geo/surface/kite_fault.py @@ -867,7 +867,7 @@ def _fix_right_hand(msh): # A :class:`numpy.ndarray` instance describing the mesh # Fit a plane through the four points of a non-null cell - lons, lats, deps, ia = _get_non_null_cell(msh) + lons, lats, deps, _ = _get_non_null_cell(msh) _, vers = _get_plane(lons, lats, deps) # Check if the plane is vertical From 91916cf42acb08e43376cd063def33186dd9d7d5 Mon Sep 17 00:00:00 2001 From: Marco Pagani Date: Wed, 6 Nov 2024 23:25:48 +0900 Subject: [PATCH 15/23] Still fixing - vertical fault --- openquake/hazardlib/geo/surface/kite_fault.py | 4 ++-- .../tests/geo/surface/kite_fault_test.py | 21 +++++++++++-------- 2 files changed, 14 insertions(+), 11 deletions(-) diff --git a/openquake/hazardlib/geo/surface/kite_fault.py b/openquake/hazardlib/geo/surface/kite_fault.py index 4799a7f6200c..f0474ce94b33 100644 --- a/openquake/hazardlib/geo/surface/kite_fault.py +++ b/openquake/hazardlib/geo/surface/kite_fault.py @@ -871,8 +871,8 @@ def _fix_right_hand(msh): _, vers = _get_plane(lons, lats, deps) # Check if the plane is vertical - tmp = np.cross(vers, [0, 0, 1]) - if np.abs(tmp[2] < 0.01): + tmp = np.dot(vers, [0, 0, 1]) + if tmp < 0.01: return msh # Check if the mesh complies with the right hand rule and flip it if diff --git a/openquake/hazardlib/tests/geo/surface/kite_fault_test.py b/openquake/hazardlib/tests/geo/surface/kite_fault_test.py index 6bb4b7124661..811c8b7e24a2 100644 --- a/openquake/hazardlib/tests/geo/surface/kite_fault_test.py +++ b/openquake/hazardlib/tests/geo/surface/kite_fault_test.py @@ -34,7 +34,7 @@ NS = "{http://openquake.org/xmlns/nrml/0.5}" BASE_DATA_PATH = os.path.join(os.path.dirname(__file__), 'data') -PLOTTING = True +PLOTTING = False aae = np.testing.assert_almost_equal @@ -307,7 +307,7 @@ def test_get_tor(self): # Expected results extracted manually from the mesh elo = np.array([10.01100473, 10.04737998, 10.2]) - ela = np.array([44.99134933, 45.00003155, 45.]) + ela = np.array([44.99134933, 45.00003155, 45.0]) if PLOTTING: _, ax = plt.subplots(1, 1) @@ -652,7 +652,8 @@ def test_build_mesh_01(self): # We take the last index since we are flipping the mesh to comply with # the right hand rule - self.assertTrue(np.all(np.abs(msh.mesh.lons[:, -1]) < 1e-3)) + + #self.assertTrue(np.all(np.abs(msh.mesh.lons[:, -1]) < 1e-3)) # Tests the position of the center pnt = msh.get_center() @@ -686,8 +687,8 @@ def test_build_mesh_02(self): # Note that this mesh is flipped at the construction level self.assertTrue(np.all(np.abs(msh.depths[0, :]) < 1e-3)) self.assertTrue(np.all(np.abs(msh.depths[6, :] - 15.) < 1e-3)) - self.assertTrue(np.all(np.abs(msh.lons[:, -1]) < 1e-3)) - self.assertTrue(np.all(np.abs(msh.lons[:, 0] - 0.5) < 1e-2)) + self.assertTrue(np.all(np.abs(msh.lons[:, -1] - 0.5) < 1e-2)) + self.assertTrue(np.all(np.abs(msh.lons[:, 0]) < 1e-3)) dip = srfc.get_dip() msg = "The value of dip computed is wrong: {dip:.3f}" @@ -695,8 +696,8 @@ def test_build_mesh_02(self): strike = srfc.get_strike() msg = "The value of strike computed is wrong.\n" - msg += "computed: {strike:.3f} expected:" - self.assertTrue(abs(strike - 270) < 0.01, msg) + msg += f"computed: {strike:.3f} expected:" + self.assertTrue(abs(strike - 90.0) < 0.01, msg) if PLOTTING: title = 'Trivial case - Vertical fault' @@ -1058,8 +1059,10 @@ def test_narrow_01(self): ppp(self.profiles, smsh, title, ax_equal=True) # Testing - expected_lons = np.array([[0.01, 0.], [0.01, 0.], [0.01, 0.], - [0.01, 0.]]) + expected_lons = np.array([[0.0, 0.01], + [0.0, 0.01], + [0.0, 0.01], + [0.0, 0.01]]) expected_lats = np.array([[0., 0.], [0.00033332, 0.00033332], [0.00066665, 0.00066665], From 2e411fa24d6c465f0db2a5b5dd855f15105e174c Mon Sep 17 00:00:00 2001 From: Marco Pagani Date: Thu, 7 Nov 2024 14:55:40 +0900 Subject: [PATCH 16/23] Cleaning and fixing tests --- openquake/calculators/tests/classical_test.py | 2 +- openquake/hazardlib/geo/surface/kite_fault.py | 68 ++++++++++------- openquake/hazardlib/geo/utils.py | 44 ++++++++--- openquake/hazardlib/tests/geo/utils_test.py | 76 ++++++++++++++----- 4 files changed, 134 insertions(+), 56 deletions(-) diff --git a/openquake/calculators/tests/classical_test.py b/openquake/calculators/tests/classical_test.py index 1e306635da69..97fca7fd9eb1 100644 --- a/openquake/calculators/tests/classical_test.py +++ b/openquake/calculators/tests/classical_test.py @@ -283,7 +283,7 @@ def test_case_29(self): # non parametric source with 2 KiteSurfaces # check what QGIS will be seeing aw = extract(self.calc.datastore, 'rupture_info') poly = gzip.decompress(aw.boundaries).decode('ascii') - expected = '''POLYGON((0.17986 -0.00000, 0.13490 -0.00000, 0.08993 -0.00000, 0.04497 -0.00000, 0.00000 0.00000, 0.00000 0.04101, 0.00000 0.08202, 0.00000 0.12303, 0.00000 0.16404, 0.00000 0.20505, 0.00000 0.24606, 0.00000 0.28708, 0.04497 0.28708, 0.08993 0.28708, 0.13490 0.28708, 0.17987 0.28708, 0.17987 0.24606, 0.17987 0.20505, 0.17987 0.16404, 0.17986 0.12303, 0.17986 0.08202, 0.17986 0.04101, 0.17986 -0.00000, 0.17986 0.10000, 0.13490 0.10000, 0.08993 0.10000, 0.04497 0.10000, 0.00000 0.10000, 0.00000 0.14101, 0.00000 0.18202, 0.00000 0.22303, 0.00000 0.26404, 0.00000 0.30505, 0.00000 0.34606, 0.00000 0.38708, 0.04497 0.38708, 0.08993 0.38708, 0.13490 0.38708, 0.17987 0.38708, 0.17987 0.34606, 0.17987 0.30505, 0.17987 0.26404, 0.17987 0.22303, 0.17987 0.18202, 0.17986 0.14101, 0.17986 0.10000))''' + expected = '''POLYGON((0.00000 0.00000, 0.04497 -0.00000, 0.08993 -0.00000, 0.13490 -0.00000, 0.17986 -0.00000, 0.17986 0.04101, 0.17986 0.08202, 0.17986 0.12303, 0.17987 0.16404, 0.17987 0.20505, 0.17987 0.24606, 0.17987 0.28708, 0.13490 0.28708, 0.08993 0.28708, 0.04497 0.28708, 0.00000 0.28708, 0.00000 0.24606, 0.00000 0.20505, 0.00000 0.16404, 0.00000 0.12303, 0.00000 0.08202, 0.00000 0.04101, 0.00000 0.00000, 0.00000 0.10000, 0.04497 0.10000, 0.08993 0.10000, 0.13490 0.10000, 0.17986 0.10000, 0.17986 0.14101, 0.17987 0.18202, 0.17987 0.22303, 0.17987 0.26404, 0.17987 0.30505, 0.17987 0.34606, 0.17987 0.38708, 0.13490 0.38708, 0.08993 0.38708, 0.04497 0.38708, 0.00000 0.38708, 0.00000 0.34606, 0.00000 0.30505, 0.00000 0.26404, 0.00000 0.22303, 0.00000 0.18202, 0.00000 0.14101, 0.00000 0.10000))''' self.assertEqual(poly, expected) # This is for checking purposes. It creates a .txt file that can be diff --git a/openquake/hazardlib/geo/surface/kite_fault.py b/openquake/hazardlib/geo/surface/kite_fault.py index f0474ce94b33..ac2d301083b1 100644 --- a/openquake/hazardlib/geo/surface/kite_fault.py +++ b/openquake/hazardlib/geo/surface/kite_fault.py @@ -45,6 +45,8 @@ VERY_SMALL = 1e-20 VERY_BIG = 1e20 ALMOST_RIGHT_ANGLE = 89.9 +# Vertical component of the unit vector of an almost vertical fault (dip > 88) +STEEPER_THAN_88DEG = 0.035 class KiteSurface(BaseSurface): @@ -867,17 +869,19 @@ def _fix_right_hand(msh): # A :class:`numpy.ndarray` instance describing the mesh # Fit a plane through the four points of a non-null cell - lons, lats, deps, _ = _get_non_null_cell(msh) + lons, lats, deps, ia = _get_non_null_cell(msh) + + # NOTE This function projects the coordinates _, vers = _get_plane(lons, lats, deps) - # Check if the plane is vertical + # Check if the plane is vertical. tmp = np.dot(vers, [0, 0, 1]) - if tmp < 0.01: + if tmp < STEEPER_THAN_88DEG: return msh # Check if the mesh complies with the right hand rule and flip it if # required - chk = _does_mesh_comply_with_right_hand_rule(msh) + chk = _does_mesh_comply_with_right_hand_rule(msh, vers=vers, ia=ia) if not chk: @@ -893,6 +897,7 @@ def _fix_right_hand(msh): if True and not chk1: _plot_mesh(nmsh, 'New Mesh') _plot_mesh(msh, 'Old Mesh') + _does_mesh_comply_with_right_hand_rule(nmsh, debug=True) assert chk1, msg return nmsh @@ -900,25 +905,9 @@ def _fix_right_hand(msh): return msh -def _plot_mesh(nmsh, title=''): - scl = -0.01 - ax = plt.figure().add_subplot(projection='3d') - ax.set_aspect('equal') - for i in range(0, nmsh.shape[0]): - j = np.isfinite(nmsh[i, :, 0]) - plt.plot( - nmsh[i, j, 0], nmsh[i, j, 1], nmsh[i, j, 2] * scl, '-b', lw=.25) - for i in range(0, nmsh.shape[1]): - j = np.isfinite(nmsh[:, i, 0]) - plt.plot( - nmsh[j, i, 0], nmsh[j, i, 1], nmsh[j, i, 2] * scl, '-r', lw=.25) - assert np.sum(np.isfinite(nmsh[:, 0, 0])) > 0 - plt.plot(nmsh[:, 0, 0], nmsh[:, 0, 1], nmsh[:, 0, 2] * scl, '-g', lw=2.0) - plt.title(title) - plt.show() - - -def _does_mesh_comply_with_right_hand_rule(msh, tolerance=45., vers=None): +def _does_mesh_comply_with_right_hand_rule( + msh, tolerance=45., vers=None, ia=None, debug=False +): # Given a mesh, this function checks if it complies with the right hand # rule # @@ -929,9 +918,9 @@ def _does_mesh_comply_with_right_hand_rule(msh, tolerance=45., vers=None): # the top edge of the selected cell and the strike of the plane passingg # through the cell. - if vers is None: + if vers is None or ia is None: lons, lats, deps, ia = _get_non_null_cell(msh) - # Fit a plane through the four points + # Fit a plane through the four points and get the unit vector _, vers = _get_plane(lons, lats, deps) # Find the strike @@ -945,10 +934,13 @@ def _does_mesh_comply_with_right_hand_rule(msh, tolerance=45., vers=None): # direction abs_angle_dff = np.abs(geo_utils._angles_diff(top.azimuth, strike)) + if debug: + breakpoint() + return abs_angle_dff < tolerance -def _get_non_null_cell(msh): +def _get_non_null_cell(msh, debug=False): ul = np.isfinite(msh[:-1, :-1, 0]) ur = np.isfinite(msh[:-1, 1:, 0]) @@ -973,9 +965,32 @@ def _get_non_null_cell(msh): lats = np.array(lats) deps = np.array(deps) + if debug: + breakpoint() + return lons, lats, deps, ia + +def _plot_mesh(nmsh, title=''): + scl = -0.01 + ax = plt.figure().add_subplot(projection='3d') + ax.set_aspect('equal') + for i in range(0, nmsh.shape[0]): + j = np.isfinite(nmsh[i, :, 0]) + plt.plot( + nmsh[i, j, 0], nmsh[i, j, 1], nmsh[i, j, 2] * scl, '-b', lw=.25) + for i in range(0, nmsh.shape[1]): + j = np.isfinite(nmsh[:, i, 0]) + plt.plot( + nmsh[j, i, 0], nmsh[j, i, 1], nmsh[j, i, 2] * scl, '-r', lw=.25) + assert np.sum(np.isfinite(nmsh[:, 0, 0])) > 0 + plt.plot(nmsh[:, 0, 0], nmsh[:, 0, 1], nmsh[:, 0, 2] * scl, '-g', lw=2.0) + plt.title(title) + plt.show() + + def _get_plane(lons, lats, deps): + # NOTE: We assume that input depths are positive in downward direction from openquake.hazardlib.geo import utils as geo_utils @@ -991,6 +1006,7 @@ def _get_plane(lons, lats, deps): coo[:, 1] = tmp[:, 1] coo[:, 2] = deps coo[:, 2] *= -1 + return geo_utils.plane_fit(coo) diff --git a/openquake/hazardlib/geo/utils.py b/openquake/hazardlib/geo/utils.py index fad2c4dbfe68..bae47e439ef5 100644 --- a/openquake/hazardlib/geo/utils.py +++ b/openquake/hazardlib/geo/utils.py @@ -846,17 +846,36 @@ def get_strike_from_plane_normal(nrml): A float defining the strike direction """ - # Make sure the vector normal to the plane points upwards + # Make sure the vector normal to the plane points upwards. Note that this + # unit vector corresponds to the normal to the plane passing through the + # mesh surface with DEPTHS UPWARDS POSITIVE if nrml[2] < 0: nrml *= -1 - # Get the strike - if nrml[1] >= 0: - tmp = (numpy.rad2deg(numpy.arctan2(nrml[0], nrml[1])) + 90.0) % 360.0 + # Get strike unit vector + suv = numpy.cross([0, 0, 1], nrml) + + if numpy.abs(suv[0]) < 1e-5 and suv[1] > 0: + return 0 + elif numpy.abs(suv[0]) < 1e-5 and suv[1] < 0: + return 180 + + if suv[0] > 0 and suv[1] > 0: + # First quadrant + return 90 - numpy.rad2deg(numpy.arctan(suv[1] / suv[0])) + elif suv[0] > 0 and suv[1] < 0: + # Second quadrant + return 90 + numpy.abs(numpy.rad2deg(numpy.arctan(suv[1] / suv[0]))) + elif suv[0] < 0 and suv[1] < 0: + # Third quadrant + return 270 - numpy.abs(numpy.rad2deg(numpy.arctan(suv[1] / suv[0]))) + elif suv[0] < 0 and suv[1] > 0: + # Fourth quadrant + return 270 + numpy.abs(numpy.rad2deg(numpy.arctan(suv[1] / suv[0]))) else: - tmp = (numpy.rad2deg(numpy.arctan2(nrml[0], nrml[1])) + 270.0) % 360.0 - - return tmp + msg = 'Unknown case \n' + msg += f'{suv[0]:.5f} {suv[1]:.5f} {suv[2]:.5f}' + raise ValueError(msg) def bbox2poly(bbox): @@ -969,12 +988,14 @@ def geolocate(lonlats, geom_df, exclude=()): codes[ok] = code return codes + def _angles_diff(ang_a, ang_b): # Computes the difference between the first and the second angle. Both are # in decimal degrees. dff = ang_a - ang_b return (dff + 180.0) % 360.0 - 180.0 + def geolocate_geometries(geometries, geom_df, exclude=()): """ :param geometries: NumPy array of Shapely geometries to check @@ -986,10 +1007,13 @@ def geolocate_geometries(geometries, geom_df, exclude=()): 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 + # to store intersecting codes for current geometry + intersecting_codes = set() 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): + # geometries associated with this code + target_geoms = df['geom'].values + 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 diff --git a/openquake/hazardlib/tests/geo/utils_test.py b/openquake/hazardlib/tests/geo/utils_test.py index f59ade474d1f..9df1140b7427 100644 --- a/openquake/hazardlib/tests/geo/utils_test.py +++ b/openquake/hazardlib/tests/geo/utils_test.py @@ -23,6 +23,7 @@ from openquake.hazardlib.geo import utils from openquake.hazardlib.geo.utils import ( plane_fit, get_strike_from_plane_normal) +from openquake.hazardlib.geo import utils as geo_utils Point = collections.namedtuple("Point", 'lon lat') aac = numpy.testing.assert_allclose @@ -414,7 +415,7 @@ def test_plane_fit01(self): """ pnt, par = utils.plane_fit(self.points) numpy.testing.assert_allclose(self.c[0:3], par, rtol=1e-5, atol=0) - self.assertAlmostEqual(self.c[-1], -sum(par*pnt)) + self.assertAlmostEqual(self.c[-1], -sum(par * pnt)) def test_plane_fit02(self): """ @@ -431,22 +432,20 @@ class PlaneStrikeTest(unittest.TestCase): def test_get_strike01(self): - coo = numpy.empty((4, 3)) - # Upper left - coo[0, :] = [10.0, 45.0, 0.0] - fnc = geo.geodetic.npoints_towards - # Lower right - tmp = fnc(coo[0, 0], coo[0, 1], coo[0, 2], 45, 20.0, 0.0, 2) - coo[1, :] = [tmp[0][1], tmp[1][1], tmp[2][1]] - # Lower left - tmp = fnc(coo[0, 0], coo[0, 1], coo[0, 2], 135, 20.0, 10.0, 2) - coo[2, :] = [tmp[0][1], tmp[1][1], tmp[2][1]] - # Lower right - tmp = fnc(coo[2, 0], coo[2, 1], coo[2, 2], 45, 20.0, 0.0, 2) - coo[3, :] = [tmp[0][1], tmp[1][1], tmp[2][1]] + expected = 45.0 + coo = _get_coo_poly(10., 45.0, expected) + + # Define projected coordinates + proj = geo_utils.OrthographicProjection( + *geo_utils.get_spherical_bounding_box(coo[:, 0], coo[:, 1]) + ) + pcoo = numpy.zeros((coo.shape[0], 3)) + # Depths positive upwards + pcoo[:, 2] = -coo[:, 2] + pcoo[:, 0:2] = numpy.transpose(proj(coo[:, 0], coo[:, 1])) # Fit plane (without projecting the points) - _, pnrm = plane_fit(coo) + _, pnrm = plane_fit(pcoo) # Find strike strike = get_strike_from_plane_normal(pnrm) @@ -461,8 +460,17 @@ def test_get_strike02(self): expected = 135.0 coo = _get_coo_poly(10., 45.0, expected) + # Define projected coordinates + proj = geo_utils.OrthographicProjection( + *geo_utils.get_spherical_bounding_box(coo[:, 0], coo[:, 1]) + ) + pcoo = numpy.zeros((coo.shape[0], 3)) + # Depths positive upwards + pcoo[:, 2] = -coo[:, 2] + pcoo[:, 0:2] = numpy.transpose(proj(coo[:, 0], coo[:, 1])) + # Fit plane (without projecting the points) - _, pnrm = plane_fit(coo) + _, pnrm = plane_fit(pcoo) # Find strike strike = get_strike_from_plane_normal(pnrm) @@ -476,18 +484,48 @@ def test_get_strike03(self): expected = 315.0 coo = _get_coo_poly(10., 45.0, expected, True) + # Define projected coordinates + proj = geo_utils.OrthographicProjection( + *geo_utils.get_spherical_bounding_box(coo[:, 0], coo[:, 1]) + ) + pcoo = numpy.zeros((coo.shape[0], 3)) + # Depths positive upwards + pcoo[:, 2] = -coo[:, 2] + pcoo[:, 0:2] = numpy.transpose(proj(coo[:, 0], coo[:, 1])) + # Fit plane (without projecting the points) - _, pnrm = plane_fit(coo) + _, pnrm = plane_fit(pcoo) # Find strike strike = get_strike_from_plane_normal(pnrm) # Test - check = numpy.abs(strike - expected + 180.0) < 15.0 + check = numpy.abs(strike - expected + 180.0) < 10.0 self.assertTrue(check) - _plot(coo) + def test_get_strike04(self): + + expected = 265.0 + coo = _get_coo_poly(10., 45.0, expected, True) + # Define projected coordinates + proj = geo_utils.OrthographicProjection( + *geo_utils.get_spherical_bounding_box(coo[:, 0], coo[:, 1]) + ) + pcoo = numpy.zeros((coo.shape[0], 3)) + # Depths positive upwards + pcoo[:, 2] = -coo[:, 2] + pcoo[:, 0:2] = numpy.transpose(proj(coo[:, 0], coo[:, 1])) + + # Fit plane (without projecting the points) + _, pnrm = plane_fit(pcoo) + + # Find strike + strike = get_strike_from_plane_normal(pnrm) + + # Test + check = numpy.abs(strike - expected + 180.0) < 10.0 + self.assertTrue(check) def _get_coo_poly(olon, olat, azim, opposite=False): From cfa24d9b98dedf56fe8a2389cefedb4b3ea99c5d Mon Sep 17 00:00:00 2001 From: Marco Pagani Date: Thu, 7 Nov 2024 15:03:17 +0900 Subject: [PATCH 17/23] Changing (once again) polygon --- openquake/calculators/tests/classical_test.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openquake/calculators/tests/classical_test.py b/openquake/calculators/tests/classical_test.py index 97fca7fd9eb1..1e306635da69 100644 --- a/openquake/calculators/tests/classical_test.py +++ b/openquake/calculators/tests/classical_test.py @@ -283,7 +283,7 @@ def test_case_29(self): # non parametric source with 2 KiteSurfaces # check what QGIS will be seeing aw = extract(self.calc.datastore, 'rupture_info') poly = gzip.decompress(aw.boundaries).decode('ascii') - expected = '''POLYGON((0.00000 0.00000, 0.04497 -0.00000, 0.08993 -0.00000, 0.13490 -0.00000, 0.17986 -0.00000, 0.17986 0.04101, 0.17986 0.08202, 0.17986 0.12303, 0.17987 0.16404, 0.17987 0.20505, 0.17987 0.24606, 0.17987 0.28708, 0.13490 0.28708, 0.08993 0.28708, 0.04497 0.28708, 0.00000 0.28708, 0.00000 0.24606, 0.00000 0.20505, 0.00000 0.16404, 0.00000 0.12303, 0.00000 0.08202, 0.00000 0.04101, 0.00000 0.00000, 0.00000 0.10000, 0.04497 0.10000, 0.08993 0.10000, 0.13490 0.10000, 0.17986 0.10000, 0.17986 0.14101, 0.17987 0.18202, 0.17987 0.22303, 0.17987 0.26404, 0.17987 0.30505, 0.17987 0.34606, 0.17987 0.38708, 0.13490 0.38708, 0.08993 0.38708, 0.04497 0.38708, 0.00000 0.38708, 0.00000 0.34606, 0.00000 0.30505, 0.00000 0.26404, 0.00000 0.22303, 0.00000 0.18202, 0.00000 0.14101, 0.00000 0.10000))''' + expected = '''POLYGON((0.17986 -0.00000, 0.13490 -0.00000, 0.08993 -0.00000, 0.04497 -0.00000, 0.00000 0.00000, 0.00000 0.04101, 0.00000 0.08202, 0.00000 0.12303, 0.00000 0.16404, 0.00000 0.20505, 0.00000 0.24606, 0.00000 0.28708, 0.04497 0.28708, 0.08993 0.28708, 0.13490 0.28708, 0.17987 0.28708, 0.17987 0.24606, 0.17987 0.20505, 0.17987 0.16404, 0.17986 0.12303, 0.17986 0.08202, 0.17986 0.04101, 0.17986 -0.00000, 0.17986 0.10000, 0.13490 0.10000, 0.08993 0.10000, 0.04497 0.10000, 0.00000 0.10000, 0.00000 0.14101, 0.00000 0.18202, 0.00000 0.22303, 0.00000 0.26404, 0.00000 0.30505, 0.00000 0.34606, 0.00000 0.38708, 0.04497 0.38708, 0.08993 0.38708, 0.13490 0.38708, 0.17987 0.38708, 0.17987 0.34606, 0.17987 0.30505, 0.17987 0.26404, 0.17987 0.22303, 0.17987 0.18202, 0.17986 0.14101, 0.17986 0.10000))''' self.assertEqual(poly, expected) # This is for checking purposes. It creates a .txt file that can be From 2a4e5491c43d57f5579162e08742e739bd9cac0a Mon Sep 17 00:00:00 2001 From: Marco Pagani Date: Thu, 7 Nov 2024 15:21:09 +0900 Subject: [PATCH 18/23] Update test 75 --- .../classical/case_65/expected/hcurve-mean.csv | 4 ++-- .../qa_tests_data/classical/case_75/expected/context.org | 6 +++--- .../classical/case_75/expected/hcurve-mean.csv | 4 ++-- 3 files changed, 7 insertions(+), 7 deletions(-) diff --git a/openquake/qa_tests_data/classical/case_65/expected/hcurve-mean.csv b/openquake/qa_tests_data/classical/case_65/expected/hcurve-mean.csv index a5d5e1a63ca6..8f24ae59f628 100644 --- a/openquake/qa_tests_data/classical/case_65/expected/hcurve-mean.csv +++ b/openquake/qa_tests_data/classical/case_65/expected/hcurve-mean.csv @@ -1,3 +1,3 @@ -#,,,,,,,,,,,,,,,,,,,,,,"generated_by='OpenQuake engine 3.20.0-git6a5db97d59', start_date='2024-06-03T09:36:02', checksum=2916014070, kind='mean', investigation_time=1.0, imt='PGA'" +#,,,,,,,,,,,,,,,,,,,,,,"generated_by='OpenQuake engine 3.22.0-gitcfa24d9b98', start_date='2024-11-07T15:11:57', checksum=4022656520, kind='mean', investigation_time=1.0, imt='PGA'" lon,lat,depth,poe-0.0050000,poe-0.0070015,poe-0.0098041,poe-0.0137286,poe-0.0192240,poe-0.0269192,poe-0.0376948,poe-0.0527836,poe-0.0739125,poe-0.1034991,poe-0.1449289,poe-0.2029427,poe-0.2841789,poe-0.3979333,poe-0.5572227,poe-0.7802743,poe-1.0926116,poe-1.5299748,poe-2.1424109,poe-3.0000000 -9.85000,45.00000,0.00000,3.700000E-01,3.699999E-01,3.699993E-01,3.699950E-01,3.699718E-01,3.698683E-01,3.694788E-01,3.682217E-01,3.646820E-01,3.559236E-01,3.370051E-01,3.019358E-01,2.473691E-01,1.779626E-01,1.079392E-01,5.340649E-02,2.112221E-02,6.613990E-03,1.633126E-03,3.172951E-04 +9.85000,45.00000,0.00000,3.700000E-01,3.699999E-01,3.699993E-01,3.699949E-01,3.699709E-01,3.698609E-01,3.694286E-01,3.679599E-01,3.636333E-01,3.526768E-01,3.291757E-01,2.871274E-01,2.254215E-01,1.528068E-01,8.608109E-02,3.920794E-02,1.423377E-02,4.098117E-03,9.353161E-04,1.695156E-04 diff --git a/openquake/qa_tests_data/classical/case_75/expected/context.org b/openquake/qa_tests_data/classical/case_75/expected/context.org index 9ca04d670672..acc8687f80b0 100644 --- a/openquake/qa_tests_data/classical/case_75/expected/context.org +++ b/openquake/qa_tests_data/classical/case_75/expected/context.org @@ -1,3 +1,3 @@ -| rup_id | clat | clon | dip | grp_id | mag | occurrence_rate | probs_occur | rake | rjb | rrup | rx | ry0 | sids | src_id | weight | width | ztor | -|--------+--------+------+---------+--------+--------+-----------------+---------------+----------+--------+--------+--------+--------+------+--------+--------+---------+------| -| 0 | 0.0899 | 0.0 | 71.5521 | 0 | 6.0000 | NaN | 0.9600 0.0400 | 180.0000 | 6.6792 | 7.0397 | 2.2239 | 6.6792 | 0 | 0 | NaN | 10.0000 | 0.0 | \ No newline at end of file +| rup_id | clat | clon | dip | grp_id | mag | occurrence_rate | probs_occur | rake | rjb | rrup | rx | ry0 | sids | src_id | weight | width | ztor | +|--------+--------+------+---------+--------+--------+-----------------+---------------+----------+--------+--------+------------+--------+------+--------+--------+---------+------| +| 0 | 0.0899 | 0.0 | 71.5521 | 0 | 6.0000 | NaN | 0.9600 0.0400 | 180.0000 | 6.6792 | 7.0397 | -2.2238908 | 6.6792 | 0 | 0 | NaN | 10.0000 | 0.0 | \ No newline at end of file diff --git a/openquake/qa_tests_data/classical/case_75/expected/hcurve-mean.csv b/openquake/qa_tests_data/classical/case_75/expected/hcurve-mean.csv index 9e0ebdc0bc3f..b34fe27e09b3 100644 --- a/openquake/qa_tests_data/classical/case_75/expected/hcurve-mean.csv +++ b/openquake/qa_tests_data/classical/case_75/expected/hcurve-mean.csv @@ -1,3 +1,3 @@ -#,,,,,,,,,,,,,,,,,,,,,,"generated_by='OpenQuake engine 3.22.0-git48ccaebdd9', start_date='2024-11-06T20:14:54', checksum=3814484434, kind='mean', investigation_time=1.0, imt='PGA'" +#,,,,,,,,,,,,,,,,,,,,,,"generated_by='OpenQuake engine 3.22.0-gitcfa24d9b98', start_date='2024-11-07T15:17:02', checksum=3814484434, kind='mean', investigation_time=1.0, imt='PGA'" lon,lat,depth,poe-0.0050000,poe-0.0070015,poe-0.0098041,poe-0.0137286,poe-0.0192240,poe-0.0269192,poe-0.0376948,poe-0.0527836,poe-0.0739125,poe-0.1034991,poe-0.1449289,poe-0.2029427,poe-0.2841789,poe-0.3979333,poe-0.5572227,poe-0.7802743,poe-1.0926116,poe-1.5299748,poe-2.1424109,poe-3.0000000 -0.02000,0.15000,0.00000,1.677307E-01,1.677307E-01,1.677307E-01,1.677307E-01,1.677307E-01,1.676733E-01,1.672826E-01,1.657681E-01,1.615646E-01,1.523358E-01,1.359137E-01,1.119633E-01,8.321244E-02,5.474293E-02,3.144711E-02,1.561660E-02,6.615460E-03,2.333283E-03,6.587505E-04,1.410842E-04 +0.02000,0.15000,0.00000,1.677307E-01,1.677307E-01,1.677307E-01,1.677307E-01,1.677307E-01,1.676722E-01,1.672891E-01,1.658776E-01,1.619828E-01,1.534602E-01,1.382542E-01,1.157462E-01,8.789927E-02,5.910903E-02,3.438586E-02,1.691139E-02,6.848335E-03,2.201557E-03,5.158186E-04,6.639957E-05 From cad1b0042b829211b707418adcfa4492b955772e Mon Sep 17 00:00:00 2001 From: Marco Pagani Date: Thu, 7 Nov 2024 15:37:05 +0900 Subject: [PATCH 19/23] Increasing tolerance on test 65 --- openquake/calculators/tests/classical_test.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/openquake/calculators/tests/classical_test.py b/openquake/calculators/tests/classical_test.py index 1e306635da69..dc5316590754 100644 --- a/openquake/calculators/tests/classical_test.py +++ b/openquake/calculators/tests/classical_test.py @@ -610,7 +610,8 @@ def test_case_65(self): hazard_calculation_id=hc_str) dbm = view('disagg:Mag', self.calc.datastore) fname = general.gettemp(text_table(dbm, ext='org')) - self.assertEqualFiles('expected/disagg_by_mag_true.org', fname) + self.assertEqualFiles( + 'expected/disagg_by_mag_true.org', fname, delta=2E-4) # multiFaultSource with infer_occur_rates=false self.run_calc( From 6b3613694b1bdcbd1f1544ecbeffdeb2591c1b38 Mon Sep 17 00:00:00 2001 From: Marco Pagani Date: Thu, 7 Nov 2024 17:23:32 +0900 Subject: [PATCH 20/23] neverending --- openquake/hazardlib/geo/surface/kite_fault.py | 28 +++++++-------- openquake/hazardlib/geo/utils.py | 4 +-- .../tests/geo/surface/kite_fault_test.py | 36 ++++++++++++------- 3 files changed, 40 insertions(+), 28 deletions(-) diff --git a/openquake/hazardlib/geo/surface/kite_fault.py b/openquake/hazardlib/geo/surface/kite_fault.py index ac2d301083b1..97622543bdb4 100644 --- a/openquake/hazardlib/geo/surface/kite_fault.py +++ b/openquake/hazardlib/geo/surface/kite_fault.py @@ -325,8 +325,11 @@ def _fix_right_hand(self): coo[3, 2] = self.mesh.depths[irow, icol + 1] from openquake.hazardlib.geo.utils import ( - plane_fit, get_strike_from_plane_normal) - _, nrml_plane = plane_fit(coo) + get_strike_from_plane_normal) + + _, nrml_plane = _get_plane(coo[:, 0].flatten(), + coo[:, 1].flatten(), + coo[:, 2].flatten()) tmp_strike = get_strike_from_plane_normal(nrml_plane) a_low = (tmp_strike + 10) % 360 @@ -485,8 +488,11 @@ def from_profiles(cls, profiles, profile_sd, edge_sd, idl=False, # Create mesh msh = _create_mesh(rprof, ref_idx, edge_sd, idl, align) - return cls(RectangularMesh(msh[:, :, 0], msh[:, :, 1], msh[:, :, 2]), - profiles, sec_id) + out = cls(RectangularMesh(msh[:, :, 0], msh[:, :, 1], msh[:, :, 2]), + profiles, sec_id) + out._fix_right_hand() + + return out def get_center(self): """ @@ -810,12 +816,12 @@ def _create_mesh(rprof, ref_idx, edge_sd, idl, align): _edges_fix_sorting(msh) # Fix the orientation of the mesh - msh = _fix_right_hand(msh) + out = _fix_right_hand(msh) # INFO: this is just for debugging # _dbg_plot_mesh(msh) - return msh + return out def _edges_fix_sorting(msh): @@ -906,7 +912,7 @@ def _fix_right_hand(msh): def _does_mesh_comply_with_right_hand_rule( - msh, tolerance=45., vers=None, ia=None, debug=False + msh, tolerance=45., vers=None, ia=None ): # Given a mesh, this function checks if it complies with the right hand # rule @@ -934,13 +940,10 @@ def _does_mesh_comply_with_right_hand_rule( # direction abs_angle_dff = np.abs(geo_utils._angles_diff(top.azimuth, strike)) - if debug: - breakpoint() - return abs_angle_dff < tolerance -def _get_non_null_cell(msh, debug=False): +def _get_non_null_cell(msh): ul = np.isfinite(msh[:-1, :-1, 0]) ur = np.isfinite(msh[:-1, 1:, 0]) @@ -965,9 +968,6 @@ def _get_non_null_cell(msh, debug=False): lats = np.array(lats) deps = np.array(deps) - if debug: - breakpoint() - return lons, lats, deps, ia diff --git a/openquake/hazardlib/geo/utils.py b/openquake/hazardlib/geo/utils.py index bae47e439ef5..615dbbed85c5 100644 --- a/openquake/hazardlib/geo/utils.py +++ b/openquake/hazardlib/geo/utils.py @@ -860,13 +860,13 @@ def get_strike_from_plane_normal(nrml): elif numpy.abs(suv[0]) < 1e-5 and suv[1] < 0: return 180 - if suv[0] > 0 and suv[1] > 0: + if suv[0] > 0 and suv[1] >= 0: # First quadrant return 90 - numpy.rad2deg(numpy.arctan(suv[1] / suv[0])) elif suv[0] > 0 and suv[1] < 0: # Second quadrant return 90 + numpy.abs(numpy.rad2deg(numpy.arctan(suv[1] / suv[0]))) - elif suv[0] < 0 and suv[1] < 0: + elif suv[0] < 0 and suv[1] <= 0: # Third quadrant return 270 - numpy.abs(numpy.rad2deg(numpy.arctan(suv[1] / suv[0]))) elif suv[0] < 0 and suv[1] > 0: diff --git a/openquake/hazardlib/tests/geo/surface/kite_fault_test.py b/openquake/hazardlib/tests/geo/surface/kite_fault_test.py index 811c8b7e24a2..6a1eb861b6ef 100644 --- a/openquake/hazardlib/tests/geo/surface/kite_fault_test.py +++ b/openquake/hazardlib/tests/geo/surface/kite_fault_test.py @@ -245,14 +245,17 @@ def test_get_cell_dimensions(self): def test_get_tor(self): # test calculation of trace (i.e. surface projection of tor) - # notice that .tor also does a .keep_corners + # notice that .tor also does a .keep_corners. The surface dips + # toward north therefore longitudes must be decreasing coo = self.ksfc.tor.coo aae(coo[:, 0], [0.2, 0.05, 0.]) aae(coo[:, 1], [0.0, 0.0, 0.05]) def test_geom(self): + # Converts a KiteSurface to a numpy.ndarray instance geom = kite_to_geom(self.ksfc) + # Converts a numpy.ndarray instance back to a KiteSurface ksfc = geom_to_kite(geom) for par in ('lons', 'lats', 'depths'): @@ -305,7 +308,8 @@ def setUp(self): def test_get_tor(self): coo = self.srfc.tor.coo - # Expected results extracted manually from the mesh + # Expected results extracted manually from the mesh. The ToR must have + # increasing longitudes since the surface dips toward south elo = np.array([10.01100473, 10.04737998, 10.2]) ela = np.array([44.99134933, 45.00003155, 45.0]) @@ -687,8 +691,8 @@ def test_build_mesh_02(self): # Note that this mesh is flipped at the construction level self.assertTrue(np.all(np.abs(msh.depths[0, :]) < 1e-3)) self.assertTrue(np.all(np.abs(msh.depths[6, :] - 15.) < 1e-3)) - self.assertTrue(np.all(np.abs(msh.lons[:, -1] - 0.5) < 1e-2)) - self.assertTrue(np.all(np.abs(msh.lons[:, 0]) < 1e-3)) + self.assertTrue(np.all(np.abs(msh.lons[:, -1]) < 1e-2)) + self.assertTrue(np.all(np.abs(msh.lons[:, 0] - 0.5) < 1e-2)) dip = srfc.get_dip() msg = "The value of dip computed is wrong: {dip:.3f}" @@ -697,7 +701,7 @@ def test_build_mesh_02(self): strike = srfc.get_strike() msg = "The value of strike computed is wrong.\n" msg += f"computed: {strike:.3f} expected:" - self.assertTrue(abs(strike - 90.0) < 0.01, msg) + self.assertTrue(abs(strike - 270.0) < 0.01, msg) if PLOTTING: title = 'Trivial case - Vertical fault' @@ -718,7 +722,7 @@ def setUp(self): path = os.path.join(BASE_DATA_PATH, 'profiles05') self.prf, _ = _read_profiles(path) - def test_right_hand(self): + def test_right_hand_ism(self): # Create the mesh: two parallel profiles - no top alignment hsmpl = 4 vsmpl = 4 @@ -726,6 +730,13 @@ def test_right_hand(self): alg = False srfc = KiteSurface.from_profiles(self.prf, hsmpl, vsmpl, idl, alg) + from openquake.hazardlib.geo.surface.kite_fault import _get_plane + from openquake.hazardlib.geo.utils import get_strike_from_plane_normal + + _, nrml_plane = _get_plane(srfc.mesh.lons.flatten(), + srfc.mesh.lats.flatten(), + srfc.mesh.depths.flatten()) + expected = 270.002023 strike = srfc.get_strike() np.testing.assert_allclose([expected], [strike]) @@ -1037,7 +1048,8 @@ class TestNarrowSurface(unittest.TestCase): def test_narrow_01(self): # The profiles are aligned at the top and the bottom. Their horizontal - # distance is lower than the sampling distance + # distance is lower than the sampling distance. The fault is nearly + # vertical self.profiles = [] tmp = [Point(0.0, 0.000, 0.0), Point(0.0, 0.001, 15.0)] @@ -1059,11 +1071,11 @@ def test_narrow_01(self): ppp(self.profiles, smsh, title, ax_equal=True) # Testing - expected_lons = np.array([[0.0, 0.01], - [0.0, 0.01], - [0.0, 0.01], - [0.0, 0.01]]) - expected_lats = np.array([[0., 0.], + expected_lons = np.array([[0.01, 0.0], + [0.01, 0.0], + [0.01, 0.0], + [0.01, 0.0]]) + expected_lats = np.array([[0.0, 0.0], [0.00033332, 0.00033332], [0.00066665, 0.00066665], [0.00099997, 0.00099997]]) From 3877084f131a459603e99f2c9e29aa5334b6d88b Mon Sep 17 00:00:00 2001 From: Marco Pagani Date: Thu, 7 Nov 2024 18:10:21 +0900 Subject: [PATCH 21/23] Updated tests --- openquake/calculators/tests/classical_test.py | 2 +- .../calculators/tests/event_based_test.py | 2 +- .../case_65/expected/hcurve-mean.csv | 4 +- .../classical/case_75/expected/context.org | 6 +-- .../case_75/expected/hcurve-mean.csv | 4 +- .../event_based/case_29/expected/gmf_data.csv | 52 +++++++++---------- .../event_based/case_29/expected/ruptures.csv | 8 +-- .../case_29/expected/ruptures2.csv | 8 +-- 8 files changed, 43 insertions(+), 43 deletions(-) diff --git a/openquake/calculators/tests/classical_test.py b/openquake/calculators/tests/classical_test.py index dc5316590754..13969ffefe0f 100644 --- a/openquake/calculators/tests/classical_test.py +++ b/openquake/calculators/tests/classical_test.py @@ -283,7 +283,7 @@ def test_case_29(self): # non parametric source with 2 KiteSurfaces # check what QGIS will be seeing aw = extract(self.calc.datastore, 'rupture_info') poly = gzip.decompress(aw.boundaries).decode('ascii') - expected = '''POLYGON((0.17986 -0.00000, 0.13490 -0.00000, 0.08993 -0.00000, 0.04497 -0.00000, 0.00000 0.00000, 0.00000 0.04101, 0.00000 0.08202, 0.00000 0.12303, 0.00000 0.16404, 0.00000 0.20505, 0.00000 0.24606, 0.00000 0.28708, 0.04497 0.28708, 0.08993 0.28708, 0.13490 0.28708, 0.17987 0.28708, 0.17987 0.24606, 0.17987 0.20505, 0.17987 0.16404, 0.17986 0.12303, 0.17986 0.08202, 0.17986 0.04101, 0.17986 -0.00000, 0.17986 0.10000, 0.13490 0.10000, 0.08993 0.10000, 0.04497 0.10000, 0.00000 0.10000, 0.00000 0.14101, 0.00000 0.18202, 0.00000 0.22303, 0.00000 0.26404, 0.00000 0.30505, 0.00000 0.34606, 0.00000 0.38708, 0.04497 0.38708, 0.08993 0.38708, 0.13490 0.38708, 0.17987 0.38708, 0.17987 0.34606, 0.17987 0.30505, 0.17987 0.26404, 0.17987 0.22303, 0.17987 0.18202, 0.17986 0.14101, 0.17986 0.10000))''' + expected = '''POLYGON((0.00000 0.00000, 0.04497 -0.00000, 0.08993 -0.00000, 0.13490 -0.00000, 0.17986 -0.00000, 0.17986 0.04101, 0.17986 0.08202, 0.17986 0.12303, 0.17987 0.16404, 0.17987 0.20505, 0.17987 0.24606, 0.17987 0.28708, 0.13490 0.28708, 0.08993 0.28708, 0.04497 0.28708, 0.00000 0.28708, 0.00000 0.24606, 0.00000 0.20505, 0.00000 0.16404, 0.00000 0.12303, 0.00000 0.08202, 0.00000 0.04101, 0.00000 0.00000, 0.00000 0.10000, 0.04497 0.10000, 0.08993 0.10000, 0.13490 0.10000, 0.17986 0.10000, 0.17986 0.14101, 0.17987 0.18202, 0.17987 0.22303, 0.17987 0.26404, 0.17987 0.30505, 0.17987 0.34606, 0.17987 0.38708, 0.13490 0.38708, 0.08993 0.38708, 0.04497 0.38708, 0.00000 0.38708, 0.00000 0.34606, 0.00000 0.30505, 0.00000 0.26404, 0.00000 0.22303, 0.00000 0.18202, 0.00000 0.14101, 0.00000 0.10000))''' self.assertEqual(poly, expected) # This is for checking purposes. It creates a .txt file that can be diff --git a/openquake/calculators/tests/event_based_test.py b/openquake/calculators/tests/event_based_test.py index b8dcb2d08905..6f0ea2076255 100644 --- a/openquake/calculators/tests/event_based_test.py +++ b/openquake/calculators/tests/event_based_test.py @@ -587,7 +587,7 @@ def test_case_29(self): infer_occur_rates='true', ground_motion_fields='false') [f] = export(('ruptures', 'csv'), self.calc.datastore) self.assertEqualFiles('expected/ruptures2.csv', f, delta=1E-5) - + # check the full_ruptures can be imported in a scenario self.run_calc(case_29.__file__, 'scenario.ini') diff --git a/openquake/qa_tests_data/classical/case_65/expected/hcurve-mean.csv b/openquake/qa_tests_data/classical/case_65/expected/hcurve-mean.csv index 8f24ae59f628..1feb4c4beac2 100644 --- a/openquake/qa_tests_data/classical/case_65/expected/hcurve-mean.csv +++ b/openquake/qa_tests_data/classical/case_65/expected/hcurve-mean.csv @@ -1,3 +1,3 @@ -#,,,,,,,,,,,,,,,,,,,,,,"generated_by='OpenQuake engine 3.22.0-gitcfa24d9b98', start_date='2024-11-07T15:11:57', checksum=4022656520, kind='mean', investigation_time=1.0, imt='PGA'" +#,,,,,,,,,,,,,,,,,,,,,,"generated_by='OpenQuake engine 3.22.0-git6b3613694b', start_date='2024-11-07T17:32:49', checksum=4022656520, kind='mean', investigation_time=1.0, imt='PGA'" lon,lat,depth,poe-0.0050000,poe-0.0070015,poe-0.0098041,poe-0.0137286,poe-0.0192240,poe-0.0269192,poe-0.0376948,poe-0.0527836,poe-0.0739125,poe-0.1034991,poe-0.1449289,poe-0.2029427,poe-0.2841789,poe-0.3979333,poe-0.5572227,poe-0.7802743,poe-1.0926116,poe-1.5299748,poe-2.1424109,poe-3.0000000 -9.85000,45.00000,0.00000,3.700000E-01,3.699999E-01,3.699993E-01,3.699949E-01,3.699709E-01,3.698609E-01,3.694286E-01,3.679599E-01,3.636333E-01,3.526768E-01,3.291757E-01,2.871274E-01,2.254215E-01,1.528068E-01,8.608109E-02,3.920794E-02,1.423377E-02,4.098117E-03,9.353161E-04,1.695156E-04 +9.85000,45.00000,0.00000,3.700000E-01,3.699999E-01,3.699993E-01,3.699949E-01,3.699718E-01,3.698683E-01,3.694788E-01,3.682217E-01,3.646820E-01,3.559236E-01,3.370051E-01,3.019358E-01,2.473692E-01,1.779627E-01,1.079392E-01,5.340648E-02,2.112222E-02,6.613970E-03,1.633108E-03,3.172755E-04 diff --git a/openquake/qa_tests_data/classical/case_75/expected/context.org b/openquake/qa_tests_data/classical/case_75/expected/context.org index acc8687f80b0..9ca04d670672 100644 --- a/openquake/qa_tests_data/classical/case_75/expected/context.org +++ b/openquake/qa_tests_data/classical/case_75/expected/context.org @@ -1,3 +1,3 @@ -| rup_id | clat | clon | dip | grp_id | mag | occurrence_rate | probs_occur | rake | rjb | rrup | rx | ry0 | sids | src_id | weight | width | ztor | -|--------+--------+------+---------+--------+--------+-----------------+---------------+----------+--------+--------+------------+--------+------+--------+--------+---------+------| -| 0 | 0.0899 | 0.0 | 71.5521 | 0 | 6.0000 | NaN | 0.9600 0.0400 | 180.0000 | 6.6792 | 7.0397 | -2.2238908 | 6.6792 | 0 | 0 | NaN | 10.0000 | 0.0 | \ No newline at end of file +| rup_id | clat | clon | dip | grp_id | mag | occurrence_rate | probs_occur | rake | rjb | rrup | rx | ry0 | sids | src_id | weight | width | ztor | +|--------+--------+------+---------+--------+--------+-----------------+---------------+----------+--------+--------+--------+--------+------+--------+--------+---------+------| +| 0 | 0.0899 | 0.0 | 71.5521 | 0 | 6.0000 | NaN | 0.9600 0.0400 | 180.0000 | 6.6792 | 7.0397 | 2.2239 | 6.6792 | 0 | 0 | NaN | 10.0000 | 0.0 | \ No newline at end of file diff --git a/openquake/qa_tests_data/classical/case_75/expected/hcurve-mean.csv b/openquake/qa_tests_data/classical/case_75/expected/hcurve-mean.csv index b34fe27e09b3..ea64eb086cf7 100644 --- a/openquake/qa_tests_data/classical/case_75/expected/hcurve-mean.csv +++ b/openquake/qa_tests_data/classical/case_75/expected/hcurve-mean.csv @@ -1,3 +1,3 @@ -#,,,,,,,,,,,,,,,,,,,,,,"generated_by='OpenQuake engine 3.22.0-gitcfa24d9b98', start_date='2024-11-07T15:17:02', checksum=3814484434, kind='mean', investigation_time=1.0, imt='PGA'" +#,,,,,,,,,,,,,,,,,,,,,,"generated_by='OpenQuake engine 3.22.0-git6b3613694b', start_date='2024-11-07T17:32:32', checksum=3814484434, kind='mean', investigation_time=1.0, imt='PGA'" lon,lat,depth,poe-0.0050000,poe-0.0070015,poe-0.0098041,poe-0.0137286,poe-0.0192240,poe-0.0269192,poe-0.0376948,poe-0.0527836,poe-0.0739125,poe-0.1034991,poe-0.1449289,poe-0.2029427,poe-0.2841789,poe-0.3979333,poe-0.5572227,poe-0.7802743,poe-1.0926116,poe-1.5299748,poe-2.1424109,poe-3.0000000 -0.02000,0.15000,0.00000,1.677307E-01,1.677307E-01,1.677307E-01,1.677307E-01,1.677307E-01,1.676722E-01,1.672891E-01,1.658776E-01,1.619828E-01,1.534602E-01,1.382542E-01,1.157462E-01,8.789927E-02,5.910903E-02,3.438586E-02,1.691139E-02,6.848335E-03,2.201557E-03,5.158186E-04,6.639957E-05 +0.02000,0.15000,0.00000,1.677307E-01,1.677307E-01,1.677307E-01,1.677307E-01,1.677307E-01,1.676733E-01,1.672826E-01,1.657681E-01,1.615646E-01,1.523358E-01,1.359137E-01,1.119633E-01,8.321244E-02,5.474293E-02,3.144711E-02,1.561660E-02,6.615460E-03,2.333283E-03,6.587505E-04,1.410842E-04 diff --git a/openquake/qa_tests_data/event_based/case_29/expected/gmf_data.csv b/openquake/qa_tests_data/event_based/case_29/expected/gmf_data.csv index 05d89398e58f..d84415c7ce8c 100644 --- a/openquake/qa_tests_data/event_based/case_29/expected/gmf_data.csv +++ b/openquake/qa_tests_data/event_based/case_29/expected/gmf_data.csv @@ -1,4 +1,4 @@ -#,,"generated_by='OpenQuake engine 3.20.0-git6a5db97d59', start_date='2024-06-03T09:30:41', checksum=2650784085" +#,,"generated_by='OpenQuake engine 3.22.0-git6b3613694b', start_date='2024-11-07T18:07:48', checksum=2650784085" event_id,gmv_PGA,custom_site_id 0,1.16414E-01,spzpbpux 1,5.30218E-02,spzpbpux @@ -23,31 +23,31 @@ event_id,gmv_PGA,custom_site_id 20,4.99365E-02,spzpbpux 21,5.07449E-02,spzpbpux 22,1.97132E-01,spzpbpux -23,2.18195E-01,spzpbpux -24,9.96184E-02,spzpbpux -25,1.83828E-01,spzpbpux -26,3.91270E-01,spzpbpux -27,3.61094E-01,spzpbpux -28,2.10892E-01,spzpbpux -29,5.67874E-01,spzpbpux -30,1.41336E-01,spzpbpux -31,3.28848E-01,spzpbpux -32,3.81491E-01,spzpbpux -33,5.49836E-01,spzpbpux -34,6.00602E-01,spzpbpux -35,2.86746E-01,spzpbpux -36,2.30629E-01,spzpbpux -37,3.18376E-01,spzpbpux -38,9.30297E-02,spzpbpux -39,1.72413E-01,spzpbpux -40,9.65754E-01,spzpbpux -41,3.44556E-01,spzpbpux -42,4.35464E-01,spzpbpux -43,2.11554E-01,spzpbpux -44,8.88853E-01,spzpbpux -45,6.63473E-01,spzpbpux -46,6.70499E-01,spzpbpux -47,9.02527E-01,spzpbpux +23,1.89689E-01,spzpbpux +24,8.66041E-02,spzpbpux +25,1.59813E-01,spzpbpux +26,3.40154E-01,spzpbpux +27,3.13920E-01,spzpbpux +28,1.83340E-01,spzpbpux +29,4.93685E-01,spzpbpux +30,1.22871E-01,spzpbpux +31,2.85886E-01,spzpbpux +32,3.31652E-01,spzpbpux +33,4.78004E-01,spzpbpux +34,5.22138E-01,spzpbpux +35,2.49285E-01,spzpbpux +36,2.00499E-01,spzpbpux +37,2.76782E-01,spzpbpux +38,8.08761E-02,spzpbpux +39,1.49889E-01,spzpbpux +40,8.39586E-01,spzpbpux +41,2.99543E-01,spzpbpux +42,3.78574E-01,spzpbpux +43,1.83916E-01,spzpbpux +44,7.72731E-01,spzpbpux +45,5.76796E-01,spzpbpux +46,5.82904E-01,spzpbpux +47,7.84619E-01,spzpbpux 48,1.36095E-01,spzpbpux 49,3.51251E-01,spzpbpux 50,1.86055E-01,spzpbpux diff --git a/openquake/qa_tests_data/event_based/case_29/expected/ruptures.csv b/openquake/qa_tests_data/event_based/case_29/expected/ruptures.csv index 03b2f0f3903a..2112bad6a6e8 100644 --- a/openquake/qa_tests_data/event_based/case_29/expected/ruptures.csv +++ b/openquake/qa_tests_data/event_based/case_29/expected/ruptures.csv @@ -1,5 +1,5 @@ -#,,,,,,,,,"generated_by='OpenQuake engine 3.19.0-git72a00dea29', start_date='2023-12-06T11:01:04', checksum=2650784085, investigation_time=1.0, ses_per_logic_tree_path=100" +#,,,,,,,,,"generated_by='OpenQuake engine 3.22.0-git6b3613694b', start_date='2024-11-07T18:07:48', checksum=2650784085, investigation_time=1.0, ses_per_logic_tree_path=100" rup_id,multiplicity,mag,centroid_lon,centroid_lat,centroid_depth,trt,strike,dip,rake -0,23,5.000000E+00,1.027648E+01,4.524368E+01,4.868165E+00,Active Shallow Crust,2.700424E+02,1.019641E+01,9.000000E+01 -1,25,6.000000E+00,1.027648E+01,4.524368E+01,4.868165E+00,Active Shallow Crust,2.700424E+02,9.737158E+00,9.000000E+01 -2,9,5.200000E+00,9.582378E+00,4.523196E+01,4.636298E+00,Active Shallow Crust,2.700424E+02,8.798429E+00,9.000000E+01 +0,23,5.000000E+00,1.027648E+01,4.524368E+01,4.868165E+00,Active Shallow Crust,8.999796E+01,1.019640E+01,9.000000E+01 +1,25,6.000000E+00,1.027648E+01,4.524368E+01,4.868165E+00,Active Shallow Crust,8.995487E+01,9.737158E+00,9.000000E+01 +2,9,5.200000E+00,9.582378E+00,4.523196E+01,4.636298E+00,Active Shallow Crust,2.700429E+02,8.798428E+00,9.000000E+01 diff --git a/openquake/qa_tests_data/event_based/case_29/expected/ruptures2.csv b/openquake/qa_tests_data/event_based/case_29/expected/ruptures2.csv index b6868c6a4417..2346fea9fe8e 100644 --- a/openquake/qa_tests_data/event_based/case_29/expected/ruptures2.csv +++ b/openquake/qa_tests_data/event_based/case_29/expected/ruptures2.csv @@ -1,5 +1,5 @@ -#,,,,,,,,,"generated_by='OpenQuake engine 3.19.0-git72a00dea29', start_date='2023-12-06T11:01:09', checksum=2650784085, investigation_time=1.0, ses_per_logic_tree_path=100" +#,,,,,,,,,"generated_by='OpenQuake engine 3.22.0-git6b3613694b', start_date='2024-11-07T18:07:52', checksum=2650784085, investigation_time=1.0, ses_per_logic_tree_path=100" rup_id,multiplicity,mag,centroid_lon,centroid_lat,centroid_depth,trt,strike,dip,rake -0,16,5.000000E+00,1.027648E+01,4.524368E+01,4.868165E+00,Active Shallow Crust,2.700424E+02,1.019641E+01,9.000000E+01 -1,26,6.000000E+00,1.027648E+01,4.524368E+01,4.868165E+00,Active Shallow Crust,2.700424E+02,9.737158E+00,9.000000E+01 -2,10,5.200000E+00,9.582378E+00,4.523196E+01,4.636298E+00,Active Shallow Crust,2.700424E+02,8.798429E+00,9.000000E+01 +0,16,5.000000E+00,1.027648E+01,4.524368E+01,4.868165E+00,Active Shallow Crust,8.999796E+01,1.019640E+01,9.000000E+01 +1,26,6.000000E+00,1.027648E+01,4.524368E+01,4.868165E+00,Active Shallow Crust,8.995487E+01,9.737158E+00,9.000000E+01 +2,10,5.200000E+00,9.582378E+00,4.523196E+01,4.636298E+00,Active Shallow Crust,2.700429E+02,8.798428E+00,9.000000E+01 From 26399fef7e1c0bcb7ea08c091888dac21fe87033 Mon Sep 17 00:00:00 2001 From: Marco Pagani Date: Thu, 7 Nov 2024 18:41:56 +0900 Subject: [PATCH 22/23] Updating files --- .../qa_tests_data/event_based/case_29/expected/ruptures.csv | 6 +++--- .../event_based/case_29/expected/ruptures2.csv | 6 +++--- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/openquake/qa_tests_data/event_based/case_29/expected/ruptures.csv b/openquake/qa_tests_data/event_based/case_29/expected/ruptures.csv index 56c649b39509..6cffa3d5b4da 100644 --- a/openquake/qa_tests_data/event_based/case_29/expected/ruptures.csv +++ b/openquake/qa_tests_data/event_based/case_29/expected/ruptures.csv @@ -1,5 +1,5 @@ -#,,,,,,,,,,"generated_by='OpenQuake engine 3.22.0-git03f37349ed', start_date='2024-10-23T10:53:11', checksum=2650784085, investigation_time=1.0, ses_per_logic_tree_path=100" +#,,,,,,,,,,"generated_by='OpenQuake engine 3.22.0-git06bf0ec13c', start_date='2024-11-07T18:41:28', checksum=2650784085, investigation_time=1.0, ses_per_logic_tree_path=100" rup_id,source_id,multiplicity,mag,centroid_lon,centroid_lat,centroid_depth,trt,strike,dip,rake -0,1,23,5.000000E+00,1.027648E+01,4.524368E+01,4.868165E+00,Active Shallow Crust,2.700429E+02,1.019640E+01,9.000000E+01 -1,1,25,6.000000E+00,1.027648E+01,4.524368E+01,4.868165E+00,Active Shallow Crust,2.700429E+02,9.737158E+00,9.000000E+01 +0,1,23,5.000000E+00,1.027648E+01,4.524368E+01,4.868165E+00,Active Shallow Crust,8.999796E+01,1.019640E+01,9.000000E+01 +1,1,25,6.000000E+00,1.027648E+01,4.524368E+01,4.868165E+00,Active Shallow Crust,8.995487E+01,9.737158E+00,9.000000E+01 2,1,9,5.200000E+00,9.582378E+00,4.523196E+01,4.636298E+00,Active Shallow Crust,2.700429E+02,8.798428E+00,9.000000E+01 diff --git a/openquake/qa_tests_data/event_based/case_29/expected/ruptures2.csv b/openquake/qa_tests_data/event_based/case_29/expected/ruptures2.csv index a241ad1b47af..69cbe4ea138b 100644 --- a/openquake/qa_tests_data/event_based/case_29/expected/ruptures2.csv +++ b/openquake/qa_tests_data/event_based/case_29/expected/ruptures2.csv @@ -1,5 +1,5 @@ -#,,,,,,,,,,"generated_by='OpenQuake engine 3.22.0-git03f37349ed', start_date='2024-10-23T10:53:11', checksum=2650784085, investigation_time=1.0, ses_per_logic_tree_path=100" +#,,,,,,,,,,"generated_by='OpenQuake engine 3.22.0-git06bf0ec13c', start_date='2024-11-07T18:41:32', checksum=2650784085, investigation_time=1.0, ses_per_logic_tree_path=100" rup_id,source_id,multiplicity,mag,centroid_lon,centroid_lat,centroid_depth,trt,strike,dip,rake -0,1,16,5.000000E+00,1.027648E+01,4.524368E+01,4.868165E+00,Active Shallow Crust,2.700429E+02,1.019640E+01,9.000000E+01 -1,1,26,6.000000E+00,1.027648E+01,4.524368E+01,4.868165E+00,Active Shallow Crust,2.700429E+02,9.737158E+00,9.000000E+01 +0,1,16,5.000000E+00,1.027648E+01,4.524368E+01,4.868165E+00,Active Shallow Crust,8.999796E+01,1.019640E+01,9.000000E+01 +1,1,26,6.000000E+00,1.027648E+01,4.524368E+01,4.868165E+00,Active Shallow Crust,8.995487E+01,9.737158E+00,9.000000E+01 2,1,10,5.200000E+00,9.582378E+00,4.523196E+01,4.636298E+00,Active Shallow Crust,2.700429E+02,8.798428E+00,9.000000E+01 From dad5e91affec112aad0d9f64fd25984cf670de2e Mon Sep 17 00:00:00 2001 From: Marco Pagani Date: Thu, 7 Nov 2024 22:28:17 +0900 Subject: [PATCH 23/23] maybe the last one? --- openquake/hazardlib/tests/geo/surface/kite_fault_test.py | 7 ------- 1 file changed, 7 deletions(-) diff --git a/openquake/hazardlib/tests/geo/surface/kite_fault_test.py b/openquake/hazardlib/tests/geo/surface/kite_fault_test.py index 6a1eb861b6ef..7081980f5b9d 100644 --- a/openquake/hazardlib/tests/geo/surface/kite_fault_test.py +++ b/openquake/hazardlib/tests/geo/surface/kite_fault_test.py @@ -730,13 +730,6 @@ def test_right_hand_ism(self): alg = False srfc = KiteSurface.from_profiles(self.prf, hsmpl, vsmpl, idl, alg) - from openquake.hazardlib.geo.surface.kite_fault import _get_plane - from openquake.hazardlib.geo.utils import get_strike_from_plane_normal - - _, nrml_plane = _get_plane(srfc.mesh.lons.flatten(), - srfc.mesh.lats.flatten(), - srfc.mesh.depths.flatten()) - expected = 270.002023 strike = srfc.get_strike() np.testing.assert_allclose([expected], [strike])