rsqsim_api.io.write_utils ========================= .. py:module:: rsqsim_api.io.write_utils .. autoapi-nested-parse:: Utilities for writing RSQSim catalogues and generating structured quad meshes from fault surfaces. Functions --------- .. autoapisummary:: rsqsim_api.io.write_utils.write_catalogue_dataframe_and_arrays rsqsim_api.io.write_utils.create_quad_mesh_from_fault rsqsim_api.io.write_utils.fit_plane_to_points rsqsim_api.io.write_utils.get_fault_rotation_matrix rsqsim_api.io.write_utils.fault_global_to_local rsqsim_api.io.write_utils.fault_local_to_global rsqsim_api.io.write_utils.get_quad_mesh_edges rsqsim_api.io.write_utils.get_edge rsqsim_api.io.write_utils.create_local_grid rsqsim_api.io.write_utils.create_cells_from_dims rsqsim_api.io.write_utils.tri_interpolate_zcoords rsqsim_api.io.write_utils.project_2d_coords rsqsim_api.io.write_utils.find_projected_coords rsqsim_api.io.write_utils.get_mesh_boundary rsqsim_api.io.write_utils.array_to_mesh Module Contents --------------- .. py:function:: write_catalogue_dataframe_and_arrays(prefix, catalogue, directory = None, write_index = True) Write an RSQSim catalogue to a CSV file and associated NumPy arrays. Creates five files with the given prefix: ``_catalogue.csv``, ``_events.npy``, ``_patches.npy``, ``_slip.npy``, and ``_slip_time.npy``. :param prefix: File name prefix (with or without trailing underscore). :param catalogue: RSQSim catalogue object exposing ``catalogue_df``, ``event_list``, ``patch_list``, ``patch_slip``, and ``patch_time_list`` attributes. :param directory: Directory in which to write the files. Defaults to the current working directory. :param write_index: If ``True`` (default), write the DataFrame index to the CSV. :raises AssertionError: If ``directory`` is specified but does not exist, or if ``prefix`` is not a non-empty string. .. py:function:: create_quad_mesh_from_fault(points, edges, triangles, resolution = 5000.0, num_search_tris = 10, fit_plane_epsilon = 1e-05, cutoff_rotation_vecmag = 0.98, is_plane_epsilon = 1.0) Create a structured quad mesh from a triangulated fault surface. Fits a best-fit plane through the input points, rotates to local fault coordinates, identifies the four boundary edges, generates a regular grid in local coordinates, and transforms back to global coordinates. :param points: Array of shape (n_points, 3) with (x, y, z) coordinates of the triangulated fault surface vertices in NZTM (metres). :param edges: Array of shape (n_edge_points, 3) with coordinates of the ordered boundary edge points. :param triangles: Integer array of shape (n_triangles, 3) with vertex indices into ``points``. :param resolution: Target grid spacing in metres. Defaults to 5000.0 m. :param num_search_tris: Number of nearest triangles to search when interpolating z coordinates for non-planar faults. Defaults to 10. :param fit_plane_epsilon: Threshold below which plane-normal components are zeroed. Defaults to 1e-5. :param cutoff_rotation_vecmag: Dot-product threshold above which the reference direction for the rotation matrix is switched from vertical to north. Defaults to 0.98. :param is_plane_epsilon: Tolerance in metres for deciding whether the fault is planar based on mean out-of-plane displacement. Defaults to 1.0. :returns: **fault_mesh_info** -- Dictionary with keys: - ``plane_normal`` : numpy.ndarray of shape (3,) - ``plane_origin`` : numpy.ndarray of shape (3,) - ``rotation_matrix`` : numpy.ndarray of shape (3, 3) - ``points_local`` : numpy.ndarray — input points in local coords - ``edges_local`` : numpy.ndarray — boundary edges in local coords - ``fault_is_plane`` : bool - ``quad_edges`` : dict with keys ``left_edge``, ``right_edge``, ``bottom_edge``, ``top_edge`` - ``mesh_points_local`` : numpy.ndarray of shape (n_pts, 3) - ``num_horiz_points`` : int - ``num_vert_points`` : int - ``mesh_points_global`` : numpy.ndarray of shape (n_pts, 3) :rtype: dict .. py:function:: fit_plane_to_points(points, eps = 1e-05) Fit a best-fit plane through a set of 3-D points using SVD. Computes the centroid of ``points``, subtracts it, and uses SVD on the moment matrix to find the normal direction. The plane is constrained to pass through the centroid. :param points: Input 3-D point coordinates (x, y, z) in NZTM (metres). :type points: numpy.ndarray of shape (n_points, 3) :param eps: Components of the plane normal smaller than this value are zeroed before normalisation. Defaults to 1e-5. :returns: * **plane_normal** (*numpy.ndarray of shape (3,)*) -- Unit normal vector to the best-fit plane. Always oriented with a non-negative z-component. * **plane_origin** (*numpy.ndarray of shape (3,)*) -- Centroid of the input points, used as the plane origin. .. py:function:: get_fault_rotation_matrix(plane_normal, cutoff_vecmag = 0.98) Compute a rotation matrix that aligns the fault plane normal with the z-axis. Builds an orthonormal frame from the plane normal and a reference direction. If the normal is nearly vertical (dot product with ``[0,0,1] > cutoff_vecmag``), uses north ``[0,1,0]`` as the reference direction instead. :param plane_normal: Unit normal vector to the fault plane. :type plane_normal: numpy.ndarray of shape (3,) :param cutoff_vecmag: Dot-product threshold above which the north reference direction is used. Defaults to 0.98. :returns: **rotation_matrix** -- Rotation matrix whose columns are ``(tan_dir1, tan_dir2, plane_normal)``. :rtype: numpy.ndarray of shape (3, 3) .. py:function:: fault_global_to_local(points, edges, rotation_matrix, plane_origin, plane_epsilon = 1.0) Transform global fault coordinates to the local fault-plane frame. Translates by ``plane_origin`` and rotates using ``rotation_matrix``. If the mean out-of-plane (z-local) displacement is below ``plane_epsilon``, the fault is considered planar. :param points: Global (x, y, z) coordinates of fault surface vertices. :type points: numpy.ndarray of shape (n_points, 3) :param edges: Global (x, y, z) coordinates of boundary edge points. :type edges: numpy.ndarray of shape (n_edge_points, 3) :param rotation_matrix: Rotation matrix from :func:`get_fault_rotation_matrix`. :type rotation_matrix: numpy.ndarray of shape (3, 3) :param plane_origin: Origin point subtracted before rotation. :type plane_origin: numpy.ndarray of shape (3,) :param plane_epsilon: Planarity tolerance in metres. Defaults to 1.0. :returns: * **points_local** (*numpy.ndarray of shape (n_points, 3)*) -- Fault-local coordinates of the input points. * **edges_local** (*numpy.ndarray of shape (n_edge_points, 3)*) -- Fault-local coordinates of the boundary edges. * **fault_is_plane** (*bool*) -- ``True`` if the mean absolute out-of-plane displacement is below ``plane_epsilon``. .. py:function:: fault_local_to_global(points, edges, rotation_matrix, plane_origin) Transform local fault-plane coordinates back to global coordinates. :param points: Fault-local (x, y, z) point coordinates. :type points: numpy.ndarray of shape (n_points, 3) :param edges: Fault-local boundary edge coordinates. :type edges: numpy.ndarray of shape (n_edge_points, 3) :param rotation_matrix: Rotation matrix from :func:`get_fault_rotation_matrix`. :type rotation_matrix: numpy.ndarray of shape (3, 3) :param plane_origin: Origin point added back after rotation. :type plane_origin: numpy.ndarray of shape (3,) :returns: * **points_global** (*numpy.ndarray of shape (n_points, 3)*) -- Global (x, y, z) point coordinates. * **edges_global** (*numpy.ndarray of shape (n_edge_points, 3)*) -- Global (x, y, z) boundary edge coordinates. .. py:function:: get_quad_mesh_edges(edges, corner_separation = 3000.0) Identify the four boundary edges of a semi-quadrilateral fault surface. Uses dot products of successive boundary segments to locate the four corners (smallest dot products = sharpest turns), then partitions the boundary into left, right, top, and bottom edges. :param edges: Ordered boundary edge coordinates in fault-local coordinates. :type edges: numpy.ndarray of shape (n_points, 3) :param corner_separation: Reserved for future use to handle chopped-off corners. Currently unused. :returns: **edge_sides** -- Dictionary with keys ``"left_edge"``, ``"right_edge"``, ``"bottom_edge"``, and ``"top_edge"``, each an array of shape (n_edge_pts, 3) in fault-local coordinates sorted along the primary axis for that edge. :rtype: dict .. admonition:: Notes All coordinates are assumed to be in fault-local coordinates. .. py:function:: get_edge(edges, ind1, ind2, ind_min, ind_max, sort_ind) Extract a contiguous edge segment from the boundary array. Handles the wrap-around case where the segment crosses the start/end of the ``edges`` array. :param edges: Ordered boundary edge coordinates. :type edges: numpy.ndarray of shape (n_points, 3) :param ind1: Index of the first corner endpoint. :param ind2: Index of the second corner endpoint. :param ind_min: Minimum index among all four corners (used for wrap detection). :param ind_max: Maximum index among all four corners (used for wrap detection). :param sort_ind: Column index (0 for x, 1 for y) used to sort the extracted edge. :returns: **edge** -- Edge coordinates sorted along ``sort_ind``. :rtype: numpy.ndarray of shape (n_edge_pts, 3) .. py:function:: create_local_grid(points, edge_sides, triangles, fault_is_plane, resolution = 5000.0, num_search_tris = 10) Generate a structured grid of points in fault-local coordinates. Interpolates grid point positions from the four boundary edges and, for non-planar faults, interpolates z-coordinates from the triangular mesh. :param points: Fault-local vertex coordinates. :type points: numpy.ndarray of shape (n_points, 3) :param edge_sides: Dictionary of the four boundary edges returned by :func:`get_quad_mesh_edges`. :type edge_sides: dict :param triangles: Triangle connectivity (vertex indices into ``points``). :type triangles: numpy.ndarray of shape (n_triangles, 3) :param fault_is_plane: If ``True``, z-coordinates are set to zero (planar fault). :param resolution: Target grid spacing in metres. Defaults to 5000.0. :param num_search_tris: Number of nearest triangle centroids to search when interpolating z. Defaults to 10. :returns: * **mesh_points** (*numpy.ndarray of shape (n_horiz * n_vert, 3)*) -- Grid point coordinates in fault-local frame. * **num_horiz_points** (*int*) -- Number of grid points in the along-strike direction. * **num_vert_points** (*int*) -- Number of grid points in the down-dip direction. .. py:function:: create_cells_from_dims(num_verts_x, num_verts_y) Build a quad cell connectivity array from grid dimensions. :param num_verts_x: Number of vertices in the x (along-strike) direction. :param num_verts_y: Number of vertices in the y (down-dip) direction. :returns: **cell_array** -- Each row gives the four vertex indices (in row-major order) composing one quad cell. :rtype: numpy.ndarray of shape (n_cells, 4), dtype int .. py:function:: tri_interpolate_zcoords(points, triangles, mesh_points, is_mesh_edge, num_search_tris = 10) Interpolate z-coordinates onto a 2-D point set using a 3-D triangular mesh. For each interior grid point (where ``is_mesh_edge`` is ``False``), finds the enclosing triangle using a KD-tree on triangle centroids and uses barycentric coordinates to interpolate the z-value. :param points: 3-D vertex coordinates of the triangular mesh. :type points: numpy.ndarray of shape (n_points, 3) :param triangles: Triangle vertex indices into ``points``. :type triangles: numpy.ndarray of shape (n_triangles, 3) :param mesh_points: 2-D (x, y) coordinates of the query points. :type mesh_points: numpy.ndarray of shape (n_grid, 2) :param is_mesh_edge: Mask that is ``True`` for boundary grid points (z is not interpolated for these). :type is_mesh_edge: numpy.ndarray of shape (n_grid,), dtype bool :param num_search_tris: Number of nearest triangle centroids to try when searching for the enclosing triangle. Defaults to 10. :returns: **z** -- Interpolated z-values. Boundary points retain their zero initialisation value (caller is responsible for overwriting them). :rtype: numpy.ndarray of shape (n_grid,) :raises ValueError: If no enclosing triangle is found for any interior point (via :func:`project_2d_coords`). .. py:function:: project_2d_coords(tri_coords, coord, tree, num_search_tris = 10) Interpolate the z-coordinate of a 2-D point from a set of triangle vertices. Searches the ``num_search_tris`` nearest triangle centroids and returns the barycentric interpolation of z from the first triangle that contains the query point. :param tri_coords: 3-D vertex coordinates for each triangle. :type tri_coords: numpy.ndarray of shape (n_triangles, 3, 3) :param coord: 2-D (x, y) query point. :type coord: numpy.ndarray of shape (2,) :param tree: KD-tree built from triangle centroid coordinates. :type tree: scipy.spatial.cKDTree :param num_search_tris: Number of nearest centroids to search. Defaults to 10. :returns: Interpolated z-value at ``coord``. :rtype: float :raises ValueError: If none of the ``num_search_tris`` nearest triangles contains the query point. .. py:function:: find_projected_coords(tri_coord, point) Determine whether a 2-D point lies within a triangle and compute its 3-D projection. Uses Shapely for the containment test and barycentric coordinates for the interpolation. :param tri_coord: 3-D vertex coordinates of the triangle. :type tri_coord: numpy.ndarray of shape (3, 3) :param point: 2-D (x, y) query point. :type point: array-like of shape (2,) :returns: * **in_tri** (*bool*) -- ``True`` if ``point`` lies within (or on the boundary of) the triangle projected to 2-D. * **projected_coords** (*numpy.ndarray of shape (3,) or None*) -- Barycentric interpolation of ``tri_coord`` at ``point`` giving the (x, y, z) position. ``None`` if ``in_tri`` is ``False``. .. py:function:: get_mesh_boundary(triangles) Find the ordered outer boundary vertices of a triangulated mesh. Identifies edges that appear in only one triangle (boundary edges) and traverses them to produce an ordered sequence of vertex indices. :param triangles: Triangle vertex index array. :type triangles: numpy.ndarray of shape (n_triangles, 3) :returns: **vert_inds** -- Ordered vertex indices forming the closed outer boundary. :rtype: numpy.ndarray of shape (n_boundary_verts,), dtype int :raises ValueError: If the traversal cannot find a next connected edge, indicating that the mesh has holes or is otherwise non-manifold. .. admonition:: Notes Assumes the mesh has no holes. Boundary is determined purely from connectivity (edges shared by only one triangle). .. py:function:: array_to_mesh(triangle_array) Convert a flat triangle vertex array to a meshio Mesh object. :param triangle_array: Each row contains the (x, y, z) coordinates of the three corner vertices: ``[x1,y1,z1, x2,y2,z2, x3,y3,z3]``. :type triangle_array: numpy.ndarray of shape (n_triangles, 9) :returns: Mesh with deduplicated ``points`` and triangular cells. :rtype: meshio.Mesh :raises AssertionError: If ``triangle_array`` does not have exactly 9 columns.