docs: def _find_turning_points(mesh, tol=1.0):
      """
      Identifies the turning points in a rectangular mesh based on the
      deviation in the azimuth between successive points on the upper edge.
      A turning point is flagged if the change in azimuth change is greater than
      the specified tolerance (in degrees)
  
      :param mesh:
          Mesh for downsampling as instance of :class:
          openquake.hazardlib.geo.mesh.RectangularMesh
  
      :param float tol:
          Maximum difference in azimuth (decimal degrees) between successive
          points to identify a turning point
  
      :returns:
          Column indices of turning points (as numpy array)
      """
      assert isinstance(mesh, RectangularMesh)
      azimuths = geodetic.azimuth(mesh.lons[0, :-1], mesh.lats[0, :-1],
                                  mesh.lons[0, 1:], mesh.lats[0, 1:])
      naz = len(azimuths)
      azim = azimuths[0]
      # Retain initial point
      idx = [0]
      for i in range(1, naz):
          if numpy.fabs(azimuths[i] - azim) > tol:
              idx.append(i)
              azim = azimuths[i]
      # Add on last point - if not already in the set
      if idx[-1] != mesh.lons.shape[1] - 1:
          idx.append(mesh.lons.shape[1] - 1)
      return numpy.array(idx)