rsqsim_api.io.write_utils

Utilities for writing RSQSim catalogues and generating structured quad meshes from fault surfaces.

Functions

write_catalogue_dataframe_and_arrays(prefix, catalogue)

Write an RSQSim catalogue to a CSV file and associated NumPy arrays.

create_quad_mesh_from_fault(points, edges, triangles)

Create a structured quad mesh from a triangulated fault surface.

fit_plane_to_points(points[, eps])

Fit a best-fit plane through a set of 3-D points using SVD.

get_fault_rotation_matrix(plane_normal[, cutoff_vecmag])

Compute a rotation matrix that aligns the fault plane normal with the z-axis.

fault_global_to_local(points, edges, rotation_matrix, ...)

Transform global fault coordinates to the local fault-plane frame.

fault_local_to_global(points, edges, rotation_matrix, ...)

Transform local fault-plane coordinates back to global coordinates.

get_quad_mesh_edges(edges[, corner_separation])

Identify the four boundary edges of a semi-quadrilateral fault surface.

get_edge(edges, ind1, ind2, ind_min, ind_max, sort_ind)

Extract a contiguous edge segment from the boundary array.

create_local_grid(points, edge_sides, triangles, ...)

Generate a structured grid of points in fault-local coordinates.

create_cells_from_dims(num_verts_x, num_verts_y)

Build a quad cell connectivity array from grid dimensions.

tri_interpolate_zcoords(points, triangles, ...[, ...])

Interpolate z-coordinates onto a 2-D point set using a 3-D triangular mesh.

project_2d_coords(tri_coords, coord, tree[, ...])

Interpolate the z-coordinate of a 2-D point from a set of triangle vertices.

find_projected_coords(tri_coord, point)

Determine whether a 2-D point lies within a triangle and compute its 3-D projection.

get_mesh_boundary(triangles)

Find the ordered outer boundary vertices of a triangulated mesh.

array_to_mesh(triangle_array)

Convert a flat triangle vertex array to a meshio Mesh object.

Module Contents

rsqsim_api.io.write_utils.write_catalogue_dataframe_and_arrays(prefix, catalogue, directory=None, write_index=True)[source]

Write an RSQSim catalogue to a CSV file and associated NumPy arrays.

Creates five files with the given prefix: <prefix>_catalogue.csv, <prefix>_events.npy, <prefix>_patches.npy, <prefix>_slip.npy, and <prefix>_slip_time.npy.

Parameters:
  • prefix (str) – File name prefix (with or without trailing underscore).

  • catalogue – RSQSim catalogue object exposing catalogue_df, event_list, patch_list, patch_slip, and patch_time_list attributes.

  • directory (str) – Directory in which to write the files. Defaults to the current working directory.

  • write_index (bool) – 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.

rsqsim_api.io.write_utils.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)[source]

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.

Parameters:
  • points (numpy.ndarray) – Array of shape (n_points, 3) with (x, y, z) coordinates of the triangulated fault surface vertices in NZTM (metres).

  • edges (numpy.ndarray) – Array of shape (n_edge_points, 3) with coordinates of the ordered boundary edge points.

  • triangles (numpy.ndarray) – Integer array of shape (n_triangles, 3) with vertex indices into points.

  • resolution (float) – Target grid spacing in metres. Defaults to 5000.0 m.

  • num_search_tris (int) – Number of nearest triangles to search when interpolating z coordinates for non-planar faults. Defaults to 10.

  • fit_plane_epsilon (float) – Threshold below which plane-normal components are zeroed. Defaults to 1e-5.

  • cutoff_rotation_vecmag (float) – Dot-product threshold above which the reference direction for the rotation matrix is switched from vertical to north. Defaults to 0.98.

  • is_plane_epsilon (float) – 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)

Return type:

dict

rsqsim_api.io.write_utils.fit_plane_to_points(points, eps=1e-05)[source]

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.

Parameters:
  • points (numpy.ndarray of shape (n_points, 3)) – Input 3-D point coordinates (x, y, z) in NZTM (metres).

  • eps (float) – 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.

rsqsim_api.io.write_utils.get_fault_rotation_matrix(plane_normal, cutoff_vecmag=0.98)[source]

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.

Parameters:
  • plane_normal (numpy.ndarray of shape (3,)) – Unit normal vector to the fault plane.

  • cutoff_vecmag (float) – 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).

Return type:

numpy.ndarray of shape (3, 3)

rsqsim_api.io.write_utils.fault_global_to_local(points, edges, rotation_matrix, plane_origin, plane_epsilon=1.0)[source]

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.

Parameters:
  • points (numpy.ndarray of shape (n_points, 3)) – Global (x, y, z) coordinates of fault surface vertices.

  • edges (numpy.ndarray of shape (n_edge_points, 3)) – Global (x, y, z) coordinates of boundary edge points.

  • rotation_matrix (numpy.ndarray of shape (3, 3)) – Rotation matrix from get_fault_rotation_matrix().

  • plane_origin (numpy.ndarray of shape (3,)) – Origin point subtracted before rotation.

  • plane_epsilon (float) – 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.

rsqsim_api.io.write_utils.fault_local_to_global(points, edges, rotation_matrix, plane_origin)[source]

Transform local fault-plane coordinates back to global coordinates.

Parameters:
  • points (numpy.ndarray of shape (n_points, 3)) – Fault-local (x, y, z) point coordinates.

  • edges (numpy.ndarray of shape (n_edge_points, 3)) – Fault-local boundary edge coordinates.

  • rotation_matrix (numpy.ndarray of shape (3, 3)) – Rotation matrix from get_fault_rotation_matrix().

  • plane_origin (numpy.ndarray of shape (3,)) – Origin point added back after rotation.

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.

rsqsim_api.io.write_utils.get_quad_mesh_edges(edges, corner_separation=3000.0)[source]

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.

Parameters:
  • edges (numpy.ndarray of shape (n_points, 3)) – Ordered boundary edge coordinates in fault-local coordinates.

  • corner_separation (float) – 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.

Return type:

dict

Notes

All coordinates are assumed to be in fault-local coordinates.

rsqsim_api.io.write_utils.get_edge(edges, ind1, ind2, ind_min, ind_max, sort_ind)[source]

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.

Parameters:
  • edges (numpy.ndarray of shape (n_points, 3)) – Ordered boundary edge coordinates.

  • ind1 (int) – Index of the first corner endpoint.

  • ind2 (int) – Index of the second corner endpoint.

  • ind_min (int) – Minimum index among all four corners (used for wrap detection).

  • ind_max (int) – Maximum index among all four corners (used for wrap detection).

  • sort_ind (int) – Column index (0 for x, 1 for y) used to sort the extracted edge.

Returns:

edge – Edge coordinates sorted along sort_ind.

Return type:

numpy.ndarray of shape (n_edge_pts, 3)

rsqsim_api.io.write_utils.create_local_grid(points, edge_sides, triangles, fault_is_plane, resolution=5000.0, num_search_tris=10)[source]

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.

Parameters:
  • points (numpy.ndarray of shape (n_points, 3)) – Fault-local vertex coordinates.

  • edge_sides (dict) – Dictionary of the four boundary edges returned by get_quad_mesh_edges().

  • triangles (numpy.ndarray of shape (n_triangles, 3)) – Triangle connectivity (vertex indices into points).

  • fault_is_plane (bool) – If True, z-coordinates are set to zero (planar fault).

  • resolution (float) – Target grid spacing in metres. Defaults to 5000.0.

  • num_search_tris (int) – 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.

rsqsim_api.io.write_utils.create_cells_from_dims(num_verts_x, num_verts_y)[source]

Build a quad cell connectivity array from grid dimensions.

Parameters:
  • num_verts_x (int) – Number of vertices in the x (along-strike) direction.

  • num_verts_y (int) – 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.

Return type:

numpy.ndarray of shape (n_cells, 4), dtype int

rsqsim_api.io.write_utils.tri_interpolate_zcoords(points, triangles, mesh_points, is_mesh_edge, num_search_tris=10)[source]

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.

Parameters:
  • points (numpy.ndarray of shape (n_points, 3)) – 3-D vertex coordinates of the triangular mesh.

  • triangles (numpy.ndarray of shape (n_triangles, 3)) – Triangle vertex indices into points.

  • mesh_points (numpy.ndarray of shape (n_grid, 2)) – 2-D (x, y) coordinates of the query points.

  • is_mesh_edge (numpy.ndarray of shape (n_grid,), dtype bool) – Mask that is True for boundary grid points (z is not interpolated for these).

  • num_search_tris (int) – 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).

Return type:

numpy.ndarray of shape (n_grid,)

Raises:

ValueError – If no enclosing triangle is found for any interior point (via project_2d_coords()).

rsqsim_api.io.write_utils.project_2d_coords(tri_coords, coord, tree, num_search_tris=10)[source]

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.

Parameters:
  • tri_coords (numpy.ndarray of shape (n_triangles, 3, 3)) – 3-D vertex coordinates for each triangle.

  • coord (numpy.ndarray of shape (2,)) – 2-D (x, y) query point.

  • tree (scipy.spatial.cKDTree) – KD-tree built from triangle centroid coordinates.

  • num_search_tris (int) – Number of nearest centroids to search. Defaults to 10.

Returns:

Interpolated z-value at coord.

Return type:

float

Raises:

ValueError – If none of the num_search_tris nearest triangles contains the query point.

rsqsim_api.io.write_utils.find_projected_coords(tri_coord, point)[source]

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.

Parameters:
  • tri_coord (numpy.ndarray of shape (3, 3)) – 3-D vertex coordinates of the triangle.

  • point (array-like of shape (2,)) – 2-D (x, y) query point.

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.

rsqsim_api.io.write_utils.get_mesh_boundary(triangles)[source]

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.

Parameters:

triangles (numpy.ndarray of shape (n_triangles, 3)) – Triangle vertex index array.

Returns:

vert_inds – Ordered vertex indices forming the closed outer boundary.

Return type:

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.

Notes

Assumes the mesh has no holes. Boundary is determined purely from connectivity (edges shared by only one triangle).

rsqsim_api.io.write_utils.array_to_mesh(triangle_array)[source]

Convert a flat triangle vertex array to a meshio Mesh object.

Parameters:

triangle_array (numpy.ndarray of shape (n_triangles, 9)) – Each row contains the (x, y, z) coordinates of the three corner vertices: [x1,y1,z1, x2,y2,z2, x3,y3,z3].

Returns:

Mesh with deduplicated points and triangular cells.

Return type:

meshio.Mesh

Raises:

AssertionError – If triangle_array does not have exactly 9 columns.