Merge pull request #9034 from gem/line_resampling_new

Fixed resampling bug in simple fault sources

debian/changelog

@@ -1,3 +1,7 @@
+  [Marco Pagani]
+  * Fixed a bug in the code that resamples a line, with a small effect
+    on the hazard of models containing simple fault sources
+
   [Michele Simionato]
   * Optimized the calculation of quantile hazard maps (7x)
   * Internal: stored `_rates` instead of `_poes` in classical calculations
@@ -36,14 +40,14 @@
   * Adjusted the Swiss-specific implementations of the GMPEs used
     in the Swiss national seismic hazard (SUIhaz15) and risk (ERM-CH23) models.
     The new implementation returns all of total, intra- and inter-event sigmas,
-    rather than just the total one.
+    rather than just the total one
   * Extended the ModifiableGMPE class by adding an
-    apply_swiss_amplification_sa method that is used in ERM-CH23 to
+    `apply_swiss_amplification_sa` method that is used in ERM-CH23 to
     apply site specific adjustments for site effects and intra-event
-    uncertainty.
+    uncertainty
+
   * Added ch_ampl03, ch_ampl06, ch_phis2s03, ch_phis2s06,
-    ch_phisss03, ch_phis2s06 site parameters to the site.py file.
-    These parameters are required by the apply_swiss_amplification_sa method.
+    ch_phisss03, ch_phis2s06 site parameters to the site.py file
 
   [Paolo Tormene]
   * Upgraded requirements: Shapely to version 2.0.1 and Pandas to version 2.0.3

openquake/calculators/export/hazard.py

@@ -280,7 +280,7 @@ def export_uhs_xml(ekey, dstore):
 
 class Location(object):
     def __init__(self, xyz):
-        self.x, self.y = tuple(xyz)[:2]
+        self.x, self.y = xyz['lon'], xyz['lat']
         self.wkt = 'POINT(%s %s)' % (self.x, self.y)
 
 

openquake/engine/tests/aelo_test.py

@@ -32,13 +32,14 @@
 MOSAIC_DIR = os.path.dirname(mosaic.__file__)
 aac = numpy.testing.assert_allclose
 
+
 SITES = ['far -90.071 16.60'.split(), 'close -85.071 10.606'.split()]
-EXPECTED = [[0.318181, 0.661788, 0.758431], [0.763345, 1.84953, 1.28969]]
-ASCE7 = ['0.78388', '0.50000', '0.50000', '1.84953', '0.95911', '1.50000',
-         '1.50000', 'Very High', '1.28969', '0.95669', '0.60000', '0.60000',
+EXPECTED = [[0.320349, 0.667217, 0.761118], [0.76334, 1.84954, 1.28972]]
+ASCE7 = ['0.78388', '0.50000', '0.50000', '1.84954', '0.95911', '1.50000',
+         '1.50000', 'Very High', '1.28972', '0.95670', '0.60000', '0.60000',
          'Very High']
-ASCE41 = [1.5, 1.45971, 1.45971, 0.83824, 0.83824, 1., 0.6,
-          1.00811, 0.6 , 0.4, 0.57329, 0.4]
+ASCE41 = [1.5, 1.45971, 1.45971, 0.83825, 0.83825, 1., 0.6,
+          1.00814, 0.6 , 0.4, 0.57332, 0.4]
 
 
 def test_CCA():

openquake/hazardlib/geo/line.py

@@ -145,7 +145,7 @@ def average_azimuth(self):
                                                lons[1:], lats[1:])
         return get_average_azimuth(azimuths, distances)
 
-    def resample(self, section_length):
+    def resample_old(self, section_length):
         """
         Resample this line into sections.  The first point in the resampled
         line corresponds to the first point in the original line.  Starting
@@ -160,7 +160,7 @@ def resample(self, section_length):
         general smaller or greater (depending on the rounding) than the length
         of the original line.  For a straight line, the difference between the
         resulting length and the original length is at maximum half of the
-        ``section_length``.  For a curved line, the difference my be larger,
+        ``section_length``.  For a curved line, the difference may be larger,
         because of corners getting cut.
 
         :param section_length:
@@ -172,8 +172,10 @@ def resample(self, section_length):
         :rtype:
             An instance of :class:`Line`
         """
+
         if len(self.points) < 2:
             return Line(self.points)
+
         resampled_points = []
         # 1. Resample the first section. 2. Loop over the remaining points
         # in the line and resample the remaining sections.
@@ -193,6 +195,157 @@ def resample(self, section_length):
 
         return Line(resampled_points)
 
+    def resample(self, sect_len: float, orig_extremes=False):
+        """
+        Resample this line into sections.  The first point in the resampled
+        line corresponds to the first point in the original line.  Starting
+        from the first point in the original line, a line segment is defined as
+        the line connecting the last point in the resampled line and the next
+        point in the original line.
+
+
+        :param float sect_len:
+            The length of the section, in km.
+        :param bool original_extremes:
+            A boolean controlling the way in which the last point is added.
+            When true the first and last point match the original extremes.
+            When false the last point is at a `sect_len` distance from the
+            previous one, before or after the last point.
+        :returns:
+            A new line resampled into sections based on the given length.
+        """
+        if len(self.points) < 2:
+            raise ValueError('The line contains less than two points')
+
+        # Project the coordinates
+        sbb = utils.get_spherical_bounding_box
+        west, east, north, south = sbb(self.coo[:, 0], self.coo[:, 1])
+        proj = utils.OrthographicProjection(west, east, north, south)
+
+        tmp = self.points
+        coo = np.array([[p.longitude, p.latitude, p.depth] for p in tmp])
+
+        # Project the coordinates of the trace
+        txy = copy.copy(coo)
+        txy[:, 0], txy[:, 1] = proj(coo[:, 0], coo[:, 1])
+
+        # Initialise the list where we store the coordinates of the resampled
+        # trace
+        rtra_prj = [[txy[0, 0], txy[0, 1], txy[0, 2]]]
+        rtra = [self.points[0]]
+
+        # Compute the total length of the original trace
+        tot_len = np.sum(((txy[:-1, 0] - txy[1:, 0])**2 +
+                          (txy[:-1, 1] - txy[1:, 1])**2 +
+                          (txy[:-1, 2] - txy[1:, 2])**2)**0.5)
+        inc_len = 0.0
+
+        # Resampling
+        idx_vtx = -1
+        while 1:
+
+            # Computing distances from the reference point
+            dis = ((txy[:, 0] - rtra_prj[-1][0])**2 +
+                   (txy[:, 1] - rtra_prj[-1][1])**2 +
+                   (txy[:, 2] - rtra_prj[-1][2])**2)**0.5
+
+            # Fixing distances for points before the index
+            if idx_vtx > 0:
+                dis[0:idx_vtx] = 100000
+
+            # Index of the point on the trace with a distance just below the
+            # sampling distance
+            tmp = dis - sect_len
+            idx = np.where(tmp <= 0, dis, -np.inf).argmax()
+
+            # If the pick a point that is not the last one on the trace we
+            # compute the new sample by interpolation
+            if idx < len(dis) - 1:
+
+                tmp = [rtra_prj[-1][0], rtra_prj[-1][1], rtra_prj[-1][2]]
+                pnt = find_t(txy[idx + 1, :], txy[idx, :], tmp, sect_len)
+
+                if pnt is None:
+                    raise ValueError('Did not find the intersection')
+
+                # Update the list of coordinates
+                xg, yg = proj(np.array([pnt[0]]), np.array([pnt[1]]),
+                              reverse=True)
+                rtra.append(Point(xg, yg, pnt[2]))
+                rtra_prj.append(list(pnt))
+
+                # Updating incremental length
+                inc_len += sect_len
+
+                # Adding more points still on the same segment
+                Δx = (txy[idx + 1, 0] - rtra_prj[-1][0])
+                Δy = (txy[idx + 1, 1] - rtra_prj[-1][1])
+                Δz = (txy[idx + 1, 2] - rtra_prj[-1][2])
+                chk_dst = (Δx**2 + Δy**2 + Δz**2)**0.5
+                Δx_ratio = Δx / chk_dst
+                Δy_ratio = Δy / chk_dst
+                Δz_ratio = Δz / chk_dst
+
+                while chk_dst > sect_len * 0.9999:
+
+                    x_pnt = rtra_prj[-1][0] + Δx_ratio * sect_len
+                    y_pnt = rtra_prj[-1][1] + Δy_ratio * sect_len
+                    z_pnt = rtra_prj[-1][2] + Δz_ratio * sect_len
+
+                    # New point
+                    xg, yg = proj(np.array([x_pnt]), np.array([y_pnt]),
+                                  reverse=True)
+                    rtra.append(Point(xg[0], yg[0], z_pnt))
+                    rtra_prj.append([x_pnt, y_pnt, z_pnt])
+
+                    # Updating incremental length
+                    inc_len += sect_len
+
+                    # This is the distance between the last resampled point
+                    # and the second vertex of the segment
+                    chk_dst = ((txy[idx + 1, 0] - rtra_prj[-1][0])**2 +
+                               (txy[idx + 1, 1] - rtra_prj[-1][1])**2 +
+                               (txy[idx + 1, 2] - rtra_prj[-1][2])**2)**0.5
+
+            else:
+
+                # Adding one point
+                if tot_len - inc_len > 0.5 * sect_len and not orig_extremes:
+
+                    # Adding more points still on the same segment
+                    dx = (txy[-1, 0] - txy[-2, 0])
+                    dy = (txy[-1, 1] - txy[-2, 1])
+                    dz = (txy[-1, 2] - txy[-2, 2])
+                    chk_dst = (dx**2 + dy**2 + dz**2)**0.5
+                    dx_ratio = dx / chk_dst
+                    dy_ratio = dy / chk_dst
+                    dz_ratio = dz / chk_dst
+
+                    x_pnt = rtra_prj[-1][0] + dx_ratio * sect_len
+                    y_pnt = rtra_prj[-1][1] + dy_ratio * sect_len
+                    z_pnt = rtra_prj[-1][2] + dz_ratio * sect_len
+
+                    # New point
+                    xg, yg = proj(np.array([x_pnt]), np.array([y_pnt]),
+                                  reverse=True)
+                    rtra.append(Point(xg[0], yg[0], z_pnt))
+                    rtra_prj.append([x_pnt, y_pnt, z_pnt])
+
+                    # Updating incremental length
+                    inc_len += sect_len
+
+                elif orig_extremes:
+
+                    # Adding last point
+                    rtra.append(Point(coo[-1, 0], coo[-1, 1], coo[-1, 2]))
+
+                break
+
+            # Updating index
+            idx_vtx = idx + 1
+
+        return Line(rtra)
+
     def get_lengths(self) -> np.ndarray:
         """
         Calculate a numpy.ndarray instance with the length
@@ -230,9 +383,9 @@ def keep_corners(self, delta):
         # Compute the azimuth of all the segments
         azim = geodetic.azimuth(coo[:-1, 0], coo[:-1, 1],
                                 coo[1:, 0], coo[1:, 1])
-        pidx = set([0, coo.shape[0]-1])
+        pidx = set([0, coo.shape[0] - 1])
         idx = np.nonzero(np.abs(np.diff(azim)) > delta)[0]
-        pidx = sorted(list(pidx.union(set(idx+1))))
+        pidx = sorted(list(pidx.union(set(idx + 1))))
         self.points = [Point(coo[c, 0], coo[c, 1]) for c in pidx]
         self.coo = coo[pidx, :]
 
@@ -275,7 +428,8 @@ def resample_to_num_points(self, num_points):
 
     def get_tu(self, mesh):
         """
-        Computes the U and T coordinates of the GC2 method for a mesh of points.
+        Computes the U and T coordinates of the GC2 method for a mesh of
+        points.
 
         :param mesh:
             An instance of :class:`openquake.hazardlib.geo.mesh.Mesh`
@@ -340,8 +494,8 @@ def get_ui_ti(self, mesh, uhat, that):
         txy[:, 0], txy[:, 1] = proj(self.coo[:, 0], self.coo[:, 1])
 
         # Initializing ti and ui coordinates
-        ui = np.zeros((txy.shape[0]-1, sxy.shape[0]))
-        ti = np.zeros((txy.shape[0]-1, sxy.shape[0]))
+        ui = np.zeros((txy.shape[0] - 1, sxy.shape[0]))
+        ti = np.zeros((txy.shape[0] - 1, sxy.shape[0]))
 
         # For each section
         for i in range(ui.shape[0]):
@@ -375,10 +529,10 @@ def get_tu_hat(self):
         # Projected coordinates
         sx, sy = self.proj(self.coo[:, 0], self.coo[:, 1])
 
-        slen = ((sx[1:]-sx[:-1])**2 + (sy[1:]-sy[:-1])**2)**0.5
-        sg = np.zeros((len(sx)-1, 3))
-        sg[:, 0] = sx[1:]-sx[:-1]
-        sg[:, 1] = sy[1:]-sy[:-1]
+        slen = ((sx[1:] - sx[:-1])**2 + (sy[1:] - sy[:-1])**2)**0.5
+        sg = np.zeros((len(sx) - 1, 3))
+        sg[:, 0] = sx[1:] - sx[:-1]
+        sg[:, 1] = sy[1:] - sy[:-1]
         uhat = get_versor(sg)
         that = get_versor(np.cross(sg, np.array([0, 0, 1])))
         return slen, uhat, that
@@ -469,13 +623,14 @@ def get_ti_weights(ui, ti, segments_len):
                               ui[i, :] > (segments_len[i] + TOLERANCE))
         iii = np.logical_and(cond1, cond2)
         if len(iii):
-            weights[i, iii] = 1./(ui[i, iii] - segments_len[i]) - 1./ui[i, iii]
+            weights[i, iii] = (1. / (ui[i, iii] - segments_len[i])
+                               - 1. / ui[i, iii])
 
         # Case for sites on one segment
         cond3 = np.logical_and(ui[i, :] >= (0. - TOLERANCE),
                                ui[i, :] <= (segments_len[i] + TOLERANCE))
         jjj = np.logical_and(cond1, cond3)
-        weights[i, jjj] = 1/(-0.01-segments_len[i])+1/0.01
+        weights[i, jjj] = 1 / (-0.01 - segments_len[i]) + 1 / 0.01
         idx_on_trace[jjj] = 1.0
 
     return weights, idx_on_trace
@@ -486,3 +641,69 @@ def get_versor(arr):
     Returns the versor (i.e. a unit vector) of a vector
     """
     return np.divide(arr.T, np.linalg.norm(arr, axis=1)).T
+
+
+def find_t(pnt0, pnt1, ref_pnt, distance):
+    """
+    Find the point on the segment within `pnt0` and `pnt1` at `distance` from
+    `ref_pnt`. See https://tinyurl.com/meyt4ft3
+
+    :param pnt0:
+        A 1D :class:`numpy.ndarray` instance of length 3
+    :param pnt1:
+        A 1D :class:`numpy.ndarray` instance of length 3
+    :param ref_pnt:
+        A 1D :class:`numpy.ndarray` instance of length 3
+    :param distance:
+        A float with the distance in km from `ref_pnt` to the point on the
+        segment.
+    :returns:
+        A 1D :class:`numpy.ndarray` instance of length 3
+    """
+
+    # First points on the line
+    x1 = pnt0[0]
+    y1 = pnt0[1]
+    z1 = pnt0[2]
+
+    #
+    x2 = pnt1[0]
+    y2 = pnt1[1]
+    z2 = pnt1[2]
+
+    x3 = ref_pnt[0]
+    y3 = ref_pnt[1]
+    z3 = ref_pnt[2]
+
+    r = distance
+
+    pa = (x2 - x1)**2 + (y2 - y1)**2 + (z2 - z1)**2
+    pb = 2 * ((x2 - x1) * (x1 - x3) + (y2 - y1) *
+              (y1 - y3) + (z2 - z1) * (z1 - z3))
+    pc = (x3**2 + y3**2 + z3**2 + x1**2 + y1**2 + z1**2 -
+          2 * (x3 * x1 + y3 * y1 + z3 * z1) - r**2)
+
+    chk = pb * pb - 4 * pa * pc
+
+    # In this case the line is not intersecting the sphere
+    if chk < 0:
+        return None
+
+    # Computing the points of intersection
+    pu = (-pb + (pb**2 - 4 * pa * pc)**0.5) / (2 * pa)
+    x = x1 + pu * (x2 - x1)
+    y = y1 + pu * (y2 - y1)
+    z = z1 + pu * (z2 - z1)
+
+    if (x >= np.min([x1, x2]) and x <= np.max([x1, x2]) and
+            y >= np.min([y1, y2]) and y <= np.max([y1, y2]) and
+            z >= np.min([z1, z2]) and z <= np.max([z1, z2])):
+        out = [x, y, z]
+    else:
+        pu = (-pb - (pb**2 - 4 * pa * pc)**0.5) / (2 * pa)
+        x = x1 + pu * (x2 - x1)
+        y = y1 + pu * (y2 - y1)
+        z = z1 + pu * (z2 - z1)
+        out = [x, y, z]
+
+    return np.array(out)