From 7b021491deeb5209ee95f441d45900581b4d00da Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Mon, 13 May 2024 10:50:53 +1200 Subject: [PATCH 01/15] add gsf module --- qcore/gsf.py | 60 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 60 insertions(+) create mode 100644 qcore/gsf.py diff --git a/qcore/gsf.py b/qcore/gsf.py new file mode 100644 index 00000000..26bd4dcf --- /dev/null +++ b/qcore/gsf.py @@ -0,0 +1,60 @@ +""" +Module providing utilities for reading GSF files. + +Functions: + read_gsf(gsf_filepath: str) -> pd.DataFrame: + Parse a GSF file into a pandas DataFrame. +""" + +import pandas as pd + + +def read_gsf(gsf_filepath: str) -> pd.DataFrame: + """Parse a GSF file into a pandas DataFrame. + + Parameters + ---------- + gsf_file_handle : TextIO + The file handle pointing to the GSF file to read. + + Returns + ------- + pd.DataFrame + A DataFrame containing all the points in the GSF file. The DataFrame's columns are + - lon (longitude) + - lat (latitude) + - depth (Kilometres below ground, i.e. depth = -10 indicates a point 10km underground). + - sub_dx (The subdivision size in the strike direction) + - sub_dy (The subdivision size in the dip direction) + - strike + - dip + - rake + - slip (nearly always -1) + - init_time (nearly always -1) + - seg_no (the fault segment this point belongs to) + """ + with open(gsf_filepath, mode="r", encoding="utf-8") as gsf_file_handle: + # we could use pd.read_csv with the skiprows argument, but it's not + # as versatile as simply skipping the first n rows with '#' + while gsf_file_handle.readline()[0] == "#": + pass + # NOTE: This skips the line after the last line beginning with #. + # This is ok as this line is always the number of points in the GSF + # file, which we do not need. + return pd.read_csv( + gsf_file_handle, + sep=r"\s+", + names=[ + "lon", + "lat", + "depth", + "sub_dx", + "sub_dy", + "strike", + "dip", + "rake", + "slip", + "init_time", + "seg_no", + ], + ) From e271c501d717a0cb01933585c932154d2ea073b8 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Mon, 13 May 2024 15:21:34 +1200 Subject: [PATCH 02/15] add documentation and gridpoint functions --- qcore/gsf.py | 158 +++++++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 154 insertions(+), 4 deletions(-) diff --git a/qcore/gsf.py b/qcore/gsf.py index 26bd4dcf..69d7f29f 100644 --- a/qcore/gsf.py +++ b/qcore/gsf.py @@ -1,13 +1,163 @@ -""" -Module providing utilities for reading GSF files. +"""This module provides functions for working with and generating GSF files. Functions: - read_gsf(gsf_filepath: str) -> pd.DataFrame: - Parse a GSF file into a pandas DataFrame. +- gridpoint_count_in_length: + Calculates the number of gridpoints that fit into a given length. +- coordinate_meshgrid: + Creates a meshgrid of points in a bounded plane region. +- write_fault_to_gsf_file: + Writes geometry data to a GSF file. +- read_gsf: + Parses a GSF file into a pandas DataFrame. + +The GSF file format is used to define the grid of the source model for a fault. +See https://wiki.canterbury.ac.nz/display/QuakeCore/File+Formats+Used+On+GM +for details on the GSF format. """ +from typing import TextIO + +import numpy as np import pandas as pd +from qcore import coordinates + + +def gridpoint_count_in_length(length: float, resolution: float) -> int: + """Calculate the number of gridpoints that fit into a given length. + + Computes the number of gridpoints that fit into a given length, if each + gridpoint is roughly resolution metres apart, and if the gridpoints + includes the endpoints. If length = 10, and resolution = 5, then the + function returns 3. + + 5m 5m + +-----+-----+ + + Parameters + ---------- + length : float + Length to distribute grid points along. + resolution : float + Resolution of the grid. + + Returns + ------- + int + The number of gridpoints that fit into length. + """ + return int(np.round(length / resolution + 2)) + + +def coordinate_meshgrid( + origin: np.ndarray, + x_upper: np.ndarray, + y_bottom: np.ndarray, + resolution: float, +) -> np.ndarray: + """Creates a meshgrid of points in a bounded plane region. + + Given the bounds of a rectangular planar region, create a meshgrid of + (lat, lon, depth) coordinates spaced at close to resolution metres apart + in the strike and dip directions. + + origin x_upper + ┌─────────────────────┐ + │. . . . . . . . . . .│ + │ │ + │. . . . . . . . . . .│ + │ │ + │. . . . . . . . . . .│ + │ │ + │. . . . . . . . . . .│ + │ │ + │. . . . . . . . . . .│ + │ ∧ ∧ │ + └────────────┼─┼──────┘ + y_bottom └─┘ + resolution + + Parameters + ---------- + origin : np.ndarray + Coordinates of the origin point (lat, lon, depth). + x_upper : np.ndarray + Coordinates of the upper x boundary (lat, lon, depth). + y_bottom : np.ndarray + Coordinates of the bottom y boundary (lat, lon, depth). + resolution : float + Resolution of the meshgrid. + + Returns + ------- + np.ndarray + The meshgrid of the rectangular planar region. + """ + origin = coordinates.wgs_depth_to_nztm(origin) + x_upper = coordinates.wgs_depth_to_nztm(x_upper) + y_bottom = coordinates.wgs_depth_to_nztm(y_bottom) + length_x = np.linalg.norm(x_upper - origin) + length_y = np.linalg.norm(y_bottom - origin) + nx = gridpoint_count_in_length(length_x, resolution) + ny = gridpoint_count_in_length(length_y, resolution) + x = np.linspace(0, length_x, nx) + y = np.linspace(0, length_y, ny) + xv, yv = np.meshgrid(x, y) + subdivision_coordinates = np.vstack([xv.ravel(), yv.ravel()]) + transformation_matrix = np.vstack( + [(x_upper - origin) / length_x, (y_bottom - origin) / length_y] + ).T + nztm_meshgrid = (transformation_matrix @ subdivision_coordinates).T + nztm_meshgrid += origin + return coordinates.nztm_to_wgs_depth(nztm_meshgrid) + + +def write_fault_to_gsf_file( + gsf_file_handle: TextIO, + meshgrids: list[np.ndarray], + lengths: np.ndarray, + widths: np.ndarray, + strikes: np.ndarray, + dips: np.ndarray, + rakes: np.ndarray, + resolution: int, +): + """Writes geometry data to a GSF file. + + Parameters + ---------- + gsf_file_handle : TextIO + The file handle pointing to the GSF file to write to. + meshgrids : list[np.ndarray] + List of meshgrid arrays representing fault segments (length: n). + lengths : np.ndarray + An (n x 1) array of segment lengths. + widths : np.ndarray + An (n x 1) array of segment widths. + strikes : np.ndarray + An (n x 1) array of segment strikes. + dips : np.ndarray + An (n x 1) array of segment dips. + rakes : np.ndarray + An (n x 1) array of segment rakes. + resolution : int + Resolution of the meshgrid. + """ + gsf_file_handle.write( + "# LON LAT DEP(km) SUB_DX SUB_DY LOC_STK LOC_DIP LOC_RAKE SLIP(cm) INIT_TIME SEG_NO\n" + ) + number_of_points = sum(meshgrid.shape[0] for meshgrid in meshgrids) + gsf_file_handle.write(f"{number_of_points}\n") + for i, (length, width, strike, dip, rake, meshgrid) in enumerate( + zip(lengths, widths, strikes, dips, rakes, meshgrids) + ): + strike_step = length / gridpoint_count_in_length(length * 1000, resolution) + dip_step = width / gridpoint_count_in_length(width * 1000, resolution) + for point in meshgrid: + gsf_file_handle.write( + f"{point[1]:11.5f} {point[0]:11.5f} {point[2] / 1000:11.5e} {strike_step:11.5e} {dip_step:11.5e} {strike:6.1f} {dip:6.1f} {rake:6.1f} {-1.0:8.2f} {-1.0:8.2f} {i:3d}\n" + ) + def read_gsf(gsf_filepath: str) -> pd.DataFrame: """Parse a GSF file into a pandas DataFrame. From 46fa924c45638f60ad749b25a5cfe643a08e2284 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Mon, 13 May 2024 15:21:51 +1200 Subject: [PATCH 03/15] include coordinate transformation module --- qcore/coordinates.py | 62 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 62 insertions(+) create mode 100644 qcore/coordinates.py diff --git a/qcore/coordinates.py b/qcore/coordinates.py new file mode 100644 index 00000000..f21b0a69 --- /dev/null +++ b/qcore/coordinates.py @@ -0,0 +1,62 @@ +""" +Module for coordinate conversions between WGS84 (latitude and longitude) and NZTM (New Zealand Transverse Mercator) coordinate systems. + +This module provides functions for converting coordinates between WGS84 and NZTM coordinate systems. +See linz.govt.nz for a description of the NZTM coordinate system. + +Functions: +- wgs_depth_to_nztm(wgs_depth_coordinates: np.ndarray) -> np.ndarray: + Converts WGS84 coordinates (latitude, longitude, depth) to NZTM coordinates. +- nztm_to_wgs_depth(nztm_coordinates: np.ndarray) -> np.ndarray: + Converts NZTM coordinates (x, y, depth) to WGS84 coordinates. +""" + +import numpy as np +import pyproj + +# Module level conversion constants for internal use only. +_WGS_CODE = 4326 +_NZTM_CODE = 2193 + +_WGS2NZTM = pyproj.Transformer.from_crs(_WGS_CODE, _NZTM_CODE) +_NZTM2WGS = pyproj.Transformer.from_crs(_NZTM_CODE, _WGS_CODE) + + +def wgs_depth_to_nztm(wgs_depth_coordinates: np.ndarray) -> np.ndarray: + """ + Convert WGS84 coordinates (latitude, longitude, depth) to NZTM coordinates. + + Parameters + ---------- + wgs_depth_coordinates : np.ndarray + An array of shape (N, 3) containing WGS84 coordinates (latitude, longitude, depth). + + Returns + ------- + np.ndarray + An array of shape (N, 3) containing NZTM coordinates (x, y, depth). + """ + nztm_coords = np.array( + _WGS2NZTM.transform(wgs_depth_coordinates[:, 0], wgs_depth_coordinates[:, 1]), + ).T + return np.append(nztm_coords, wgs_depth_coordinates[:, 2].reshape((-1, 1)), axis=-1) + + +def nztm_to_wgs_depth(nztm_coordinates: np.ndarray) -> np.ndarray: + """ + Convert NZTM coordinates (x, y, depth) to WGS84 coordinates. + + Parameters + ---------- + nztm_coordinates : np.ndarray + An array of shape (N, 3) containing NZTM coordinates (x, y, depth). + + Returns + ------- + np.ndarray + An array of shape (N, 3) containing WGS84 coordinates (latitude, longitude, depth). + """ + wgs_coords = np.array( + _NZTM2WGS.transform(nztm_coordinates[:, 0], nztm_coordinates[:, 1]), + ).T + return np.append(wgs_coords, nztm_coordinates[:, 2].reshape((-1, 1)), axis=-1) From 0bdd8f918df8fe1e79f7a4cc2a708f4b4fa926a7 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Tue, 14 May 2024 12:48:19 +1200 Subject: [PATCH 04/15] fix read_gsf docstring --- qcore/gsf.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/qcore/gsf.py b/qcore/gsf.py index 69d7f29f..eca579e2 100644 --- a/qcore/gsf.py +++ b/qcore/gsf.py @@ -164,7 +164,7 @@ def read_gsf(gsf_filepath: str) -> pd.DataFrame: Parameters ---------- - gsf_file_handle : TextIO + gsf_filepath : str The file handle pointing to the GSF file to read. Returns From fcede8d4910b31f28126770772906c2a7ea7054f Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Tue, 14 May 2024 12:52:13 +1200 Subject: [PATCH 05/15] use Path not str or TextIO --- qcore/gsf.py | 45 +++++++++++++++++++++++---------------------- 1 file changed, 23 insertions(+), 22 deletions(-) diff --git a/qcore/gsf.py b/qcore/gsf.py index eca579e2..d3575f6f 100644 --- a/qcore/gsf.py +++ b/qcore/gsf.py @@ -15,7 +15,7 @@ for details on the GSF format. """ -from typing import TextIO +from pathlib import Path import numpy as np import pandas as pd @@ -113,7 +113,7 @@ def coordinate_meshgrid( def write_fault_to_gsf_file( - gsf_file_handle: TextIO, + gsf_filepath: Path, meshgrids: list[np.ndarray], lengths: np.ndarray, widths: np.ndarray, @@ -126,8 +126,8 @@ def write_fault_to_gsf_file( Parameters ---------- - gsf_file_handle : TextIO - The file handle pointing to the GSF file to write to. + gsf_filepath : Path + The file path pointing to the GSF file to write to. meshgrids : list[np.ndarray] List of meshgrid arrays representing fault segments (length: n). lengths : np.ndarray @@ -143,28 +143,29 @@ def write_fault_to_gsf_file( resolution : int Resolution of the meshgrid. """ - gsf_file_handle.write( - "# LON LAT DEP(km) SUB_DX SUB_DY LOC_STK LOC_DIP LOC_RAKE SLIP(cm) INIT_TIME SEG_NO\n" - ) - number_of_points = sum(meshgrid.shape[0] for meshgrid in meshgrids) - gsf_file_handle.write(f"{number_of_points}\n") - for i, (length, width, strike, dip, rake, meshgrid) in enumerate( - zip(lengths, widths, strikes, dips, rakes, meshgrids) - ): - strike_step = length / gridpoint_count_in_length(length * 1000, resolution) - dip_step = width / gridpoint_count_in_length(width * 1000, resolution) - for point in meshgrid: - gsf_file_handle.write( - f"{point[1]:11.5f} {point[0]:11.5f} {point[2] / 1000:11.5e} {strike_step:11.5e} {dip_step:11.5e} {strike:6.1f} {dip:6.1f} {rake:6.1f} {-1.0:8.2f} {-1.0:8.2f} {i:3d}\n" - ) - - -def read_gsf(gsf_filepath: str) -> pd.DataFrame: + with open(gsf_filepath, "w", encoding="utf-8") as gsf_file_handle: + gsf_file_handle.write( + "# LON LAT DEP(km) SUB_DX SUB_DY LOC_STK LOC_DIP LOC_RAKE SLIP(cm) INIT_TIME SEG_NO\n" + ) + number_of_points = sum(meshgrid.shape[0] for meshgrid in meshgrids) + gsf_file_handle.write(f"{number_of_points}\n") + for i, (length, width, strike, dip, rake, meshgrid) in enumerate( + zip(lengths, widths, strikes, dips, rakes, meshgrids) + ): + strike_step = length / gridpoint_count_in_length(length * 1000, resolution) + dip_step = width / gridpoint_count_in_length(width * 1000, resolution) + for point in meshgrid: + gsf_file_handle.write( + f"{point[1]:11.5f} {point[0]:11.5f} {point[2] / 1000:11.5e} {strike_step:11.5e} {dip_step:11.5e} {strike:6.1f} {dip:6.1f} {rake:6.1f} {-1.0:8.2f} {-1.0:8.2f} {i:3d}\n" + ) + + +def read_gsf(gsf_filepath: Path) -> pd.DataFrame: """Parse a GSF file into a pandas DataFrame. Parameters ---------- - gsf_filepath : str + gsf_filepath : Path The file handle pointing to the GSF file to read. Returns From 7d6f13835ef50cbaab9dd907064a3164fb930e41 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Tue, 14 May 2024 14:48:17 +1200 Subject: [PATCH 06/15] simplify coordinate functions --- qcore/coordinates.py | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) diff --git a/qcore/coordinates.py b/qcore/coordinates.py index f21b0a69..0cf24fe4 100644 --- a/qcore/coordinates.py +++ b/qcore/coordinates.py @@ -36,10 +36,13 @@ def wgs_depth_to_nztm(wgs_depth_coordinates: np.ndarray) -> np.ndarray: np.ndarray An array of shape (N, 3) containing NZTM coordinates (x, y, depth). """ - nztm_coords = np.array( - _WGS2NZTM.transform(wgs_depth_coordinates[:, 0], wgs_depth_coordinates[:, 1]), + return np.array( + _WGS2NZTM.transform( + wgs_depth_coordinates[:, 0], + wgs_depth_coordinates[:, 1], + wgs_depth_coordinates[:, 2], + ) ).T - return np.append(nztm_coords, wgs_depth_coordinates[:, 2].reshape((-1, 1)), axis=-1) def nztm_to_wgs_depth(nztm_coordinates: np.ndarray) -> np.ndarray: @@ -56,7 +59,10 @@ def nztm_to_wgs_depth(nztm_coordinates: np.ndarray) -> np.ndarray: np.ndarray An array of shape (N, 3) containing WGS84 coordinates (latitude, longitude, depth). """ - wgs_coords = np.array( - _NZTM2WGS.transform(nztm_coordinates[:, 0], nztm_coordinates[:, 1]), + return np.array( + _WGS2NZTM.transform( + nztm_coordinates[:, 0], + nztm_coordinates[:, 1], + nztm_coordinates[:, 2], + ) ).T - return np.append(wgs_coords, nztm_coordinates[:, 2].reshape((-1, 1)), axis=-1) From 35f1a0ae6b05313f8923cfa40f6182a98203369c Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Wed, 15 May 2024 15:43:03 +1200 Subject: [PATCH 07/15] coordinate transform fixes --- qcore/coordinates.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/qcore/coordinates.py b/qcore/coordinates.py index 0cf24fe4..50e7bce8 100644 --- a/qcore/coordinates.py +++ b/qcore/coordinates.py @@ -36,6 +36,8 @@ def wgs_depth_to_nztm(wgs_depth_coordinates: np.ndarray) -> np.ndarray: np.ndarray An array of shape (N, 3) containing NZTM coordinates (x, y, depth). """ + if wgs_depth_coordinates.shape == (3,): + return np.array(_WGS2NZTM.transform(*wgs_depth_coordinates)) return np.array( _WGS2NZTM.transform( wgs_depth_coordinates[:, 0], @@ -59,8 +61,10 @@ def nztm_to_wgs_depth(nztm_coordinates: np.ndarray) -> np.ndarray: np.ndarray An array of shape (N, 3) containing WGS84 coordinates (latitude, longitude, depth). """ + if nztm_coordinates.shape == (3,): + return np.array(_NZTM2WGS.transform(*nztm_coordinates)) return np.array( - _WGS2NZTM.transform( + _NZTM2WGS.transform( nztm_coordinates[:, 0], nztm_coordinates[:, 1], nztm_coordinates[:, 2], From f97224c2b65646a6c8de78885c486178ec3b056c Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Tue, 21 May 2024 16:36:42 +1200 Subject: [PATCH 08/15] use a pandas dataframe for gsf writing --- qcore/gsf.py | 35 +++++++++++++---------------------- 1 file changed, 13 insertions(+), 22 deletions(-) diff --git a/qcore/gsf.py b/qcore/gsf.py index d3575f6f..6c803d59 100644 --- a/qcore/gsf.py +++ b/qcore/gsf.py @@ -114,12 +114,7 @@ def coordinate_meshgrid( def write_fault_to_gsf_file( gsf_filepath: Path, - meshgrids: list[np.ndarray], - lengths: np.ndarray, - widths: np.ndarray, - strikes: np.ndarray, - dips: np.ndarray, - rakes: np.ndarray, + gsf_df: pd.DataFrame, resolution: int, ): """Writes geometry data to a GSF file. @@ -128,18 +123,8 @@ def write_fault_to_gsf_file( ---------- gsf_filepath : Path The file path pointing to the GSF file to write to. - meshgrids : list[np.ndarray] - List of meshgrid arrays representing fault segments (length: n). - lengths : np.ndarray - An (n x 1) array of segment lengths. - widths : np.ndarray - An (n x 1) array of segment widths. - strikes : np.ndarray - An (n x 1) array of segment strikes. - dips : np.ndarray - An (n x 1) array of segment dips. - rakes : np.ndarray - An (n x 1) array of segment rakes. + gsf_df : pd.DataFrame + The GSF dataframe to write. resolution : int Resolution of the meshgrid. """ @@ -147,11 +132,17 @@ def write_fault_to_gsf_file( gsf_file_handle.write( "# LON LAT DEP(km) SUB_DX SUB_DY LOC_STK LOC_DIP LOC_RAKE SLIP(cm) INIT_TIME SEG_NO\n" ) - number_of_points = sum(meshgrid.shape[0] for meshgrid in meshgrids) + number_of_points = gsf_df.apply( + lambda row: row["meshgrid"].shape[0], axis=1 + ).sum() gsf_file_handle.write(f"{number_of_points}\n") - for i, (length, width, strike, dip, rake, meshgrid) in enumerate( - zip(lengths, widths, strikes, dips, rakes, meshgrids) - ): + for i, row in gsf_df.iterrows(): + length = row["length"] + width = row["width"] + strike = row["strike"] + dip = row["dip"] + rake = row["rake"] + meshgrid = row["meshgrid"] strike_step = length / gridpoint_count_in_length(length * 1000, resolution) dip_step = width / gridpoint_count_in_length(width * 1000, resolution) for point in meshgrid: From f5257de41be7b649601aa8d9961be975b03a7c0d Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Thu, 23 May 2024 16:19:00 +1200 Subject: [PATCH 09/15] Fix gsf meshgrid fault plane interleaving --- qcore/gsf.py | 48 +++++++++++++++++++++++++++++++----------------- 1 file changed, 31 insertions(+), 17 deletions(-) diff --git a/qcore/gsf.py b/qcore/gsf.py index 6c803d59..8472fbea 100644 --- a/qcore/gsf.py +++ b/qcore/gsf.py @@ -91,7 +91,9 @@ def coordinate_meshgrid( Returns ------- np.ndarray - The meshgrid of the rectangular planar region. + The meshgrid of the rectangular planar region. Has shape (ny, nx), where + ny is the number of points in the origin->y_bottom direction and nx the number of + points in the origin->x_upper direction. """ origin = coordinates.wgs_depth_to_nztm(origin) x_upper = coordinates.wgs_depth_to_nztm(x_upper) @@ -109,7 +111,7 @@ def coordinate_meshgrid( ).T nztm_meshgrid = (transformation_matrix @ subdivision_coordinates).T nztm_meshgrid += origin - return coordinates.nztm_to_wgs_depth(nztm_meshgrid) + return coordinates.nztm_to_wgs_depth(nztm_meshgrid).reshape((ny, nx, 3)) def write_fault_to_gsf_file( @@ -119,12 +121,17 @@ def write_fault_to_gsf_file( ): """Writes geometry data to a GSF file. + This code assumes that the dip is constant across all faults. + Parameters ---------- gsf_filepath : Path The file path pointing to the GSF file to write to. gsf_df : pd.DataFrame - The GSF dataframe to write. + The GSF dataframe to write. This dataframe must have the columns length, + width, strike, dip, rake, and meshgrid. Each row corresponds to one + fault plane, with the meshgrid column being the discretisation of the + fault planes. resolution : int Resolution of the meshgrid. """ @@ -133,22 +140,29 @@ def write_fault_to_gsf_file( "# LON LAT DEP(km) SUB_DX SUB_DY LOC_STK LOC_DIP LOC_RAKE SLIP(cm) INIT_TIME SEG_NO\n" ) number_of_points = gsf_df.apply( - lambda row: row["meshgrid"].shape[0], axis=1 + lambda row: np.prod(row["meshgrid"].shape[:2]), axis=1 ).sum() + + # Get the number of dip gridpoints by looking at the first dimension of + # the meshgrid of the first fault plane. See coordinate_meshgrid for an + # explanation of meshgrid dimensions. + n_dip = gsf_df.iloc[0]['meshgrid'].shape[0] + gsf_file_handle.write(f"{number_of_points}\n") - for i, row in gsf_df.iterrows(): - length = row["length"] - width = row["width"] - strike = row["strike"] - dip = row["dip"] - rake = row["rake"] - meshgrid = row["meshgrid"] - strike_step = length / gridpoint_count_in_length(length * 1000, resolution) - dip_step = width / gridpoint_count_in_length(width * 1000, resolution) - for point in meshgrid: - gsf_file_handle.write( - f"{point[1]:11.5f} {point[0]:11.5f} {point[2] / 1000:11.5e} {strike_step:11.5e} {dip_step:11.5e} {strike:6.1f} {dip:6.1f} {rake:6.1f} {-1.0:8.2f} {-1.0:8.2f} {i:3d}\n" - ) + for j in range(n_dip): + for i, row in gsf_df.iterrows(): + length = row["length"] + width = row["width"] + strike = row["strike"] + dip = row["dip"] + rake = row["rake"] + meshgrid = row["meshgrid"] + strike_step = length / gridpoint_count_in_length(length * 1000, resolution) + dip_step = width / gridpoint_count_in_length(width * 1000, resolution) + for point in meshgrid[j]: + gsf_file_handle.write( + f"{point[1]:11.5f} {point[0]:11.5f} {point[2] / 1000:11.5e} {strike_step:11.5e} {dip_step:11.5e} {strike:6.1f} {dip:6.1f} {rake:6.1f} {-1.0:8.2f} {-1.0:8.2f} {i:3d}\n" + ) def read_gsf(gsf_filepath: Path) -> pd.DataFrame: From 08990042f8ef1ef1ba38f4a84dec9d0eddb8aeb8 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Wed, 29 May 2024 16:27:00 +1200 Subject: [PATCH 10/15] improve docstrings --- qcore/coordinates.py | 21 +++++++++++++++++++-- 1 file changed, 19 insertions(+), 2 deletions(-) diff --git a/qcore/coordinates.py b/qcore/coordinates.py index 50e7bce8..f80c5a9e 100644 --- a/qcore/coordinates.py +++ b/qcore/coordinates.py @@ -29,12 +29,21 @@ def wgs_depth_to_nztm(wgs_depth_coordinates: np.ndarray) -> np.ndarray: Parameters ---------- wgs_depth_coordinates : np.ndarray - An array of shape (N, 3) containing WGS84 coordinates (latitude, longitude, depth). + An array of shape (N, 3) or (3,) containing WGS84 coordinates + (latitude, longitude, depth). Returns ------- np.ndarray An array of shape (N, 3) containing NZTM coordinates (x, y, depth). + + Examples + -------- + >>> wgs_depth_to_nztm(np.array([-43.5320, 172.6366, 0])) + array([5180040.61473068, 1570636.6812821 , 0. ]) + >>> wgs_depth_to_nztm(np.array([[-36.8509, 174.7645, 100], [-41.2924, 174.7787, 100]])) + array([[5.92021456e+06, 1.75731133e+06, 1.00000000e+02], + [5.42725716e+06, 1.74893148e+06, 0.00000000e+00]]) """ if wgs_depth_coordinates.shape == (3,): return np.array(_WGS2NZTM.transform(*wgs_depth_coordinates)) @@ -54,12 +63,20 @@ def nztm_to_wgs_depth(nztm_coordinates: np.ndarray) -> np.ndarray: Parameters ---------- nztm_coordinates : np.ndarray - An array of shape (N, 3) containing NZTM coordinates (x, y, depth). + An array of shape (N, 3) or (3,) containing NZTM coordinates (x, y, depth). Returns ------- np.ndarray An array of shape (N, 3) containing WGS84 coordinates (latitude, longitude, depth). + + Examples + -------- + >>> nztm_to_wgs_depth(np.array([5180040.61473068, 1570636.6812821, 0])) + array([-43.5320, 172.6366, 0]) + >>> wgs_depth_to_nztm(array([[5.92021456e+06, 1.75731133e+06, 100], + [5.42725716e+06, 1.74893148e+06, 0]])) + array([[-36.8509, 174.7645, 100], [-41.2924, 174.7787, 100]]) """ if nztm_coordinates.shape == (3,): return np.array(_NZTM2WGS.transform(*nztm_coordinates)) From e90f33d659f1666f3a7e1633455ad5f11e320269 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Wed, 29 May 2024 16:40:58 +1200 Subject: [PATCH 11/15] better document meshgrid generation --- qcore/gsf.py | 41 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 41 insertions(+) diff --git a/qcore/gsf.py b/qcore/gsf.py index 8472fbea..df72c2e8 100644 --- a/qcore/gsf.py +++ b/qcore/gsf.py @@ -95,22 +95,63 @@ def coordinate_meshgrid( ny is the number of points in the origin->y_bottom direction and nx the number of points in the origin->x_upper direction. """ + # These calculations are easier to do if the coordinates are in NZTM rather + # than (lat, lon, depth). origin = coordinates.wgs_depth_to_nztm(origin) x_upper = coordinates.wgs_depth_to_nztm(x_upper) y_bottom = coordinates.wgs_depth_to_nztm(y_bottom) + length_x = np.linalg.norm(x_upper - origin) length_y = np.linalg.norm(y_bottom - origin) + nx = gridpoint_count_in_length(length_x, resolution) ny = gridpoint_count_in_length(length_y, resolution) + + # We first create a meshgrid of coordinates across a flat rectangle like the following + # + # (0, 0) (length_x, 0) + # ┌─────────┐ + # │ │ + # │ │ + # │ │ + # │ │ + # │ │ + # │ │ + # │ │ + # └─────────┘ + # (0, length_y) + x = np.linspace(0, length_x, nx) y = np.linspace(0, length_y, ny) xv, yv = np.meshgrid(x, y) subdivision_coordinates = np.vstack([xv.ravel(), yv.ravel()]) + + # The subdivision coordinates lie on a rectangle that has the right size, + # but is not in the right orientation or position. The job of the + # transformation matrix is to rotate or shear the meshgrid to fit a plane + # with the same orientation as the desired plane. + # Diagramatically: + # + # ╱╲ + # ┌─────────┐ ╱ ╲ + # │ │ ╱ ╲ + # │ │ ╱ ╲ + # │ │ ╲ ╲ + # │ │ ─────> ╲ ╲ + # │ │ tr. matrix ╲ ╲ + # │ │ ╲ ╱ + # │ │ ╲ ╱ + # └─────────┘ ╲ ╱ + # ╲╱ transformation_matrix = np.vstack( [(x_upper - origin) / length_x, (y_bottom - origin) / length_y] ).T nztm_meshgrid = (transformation_matrix @ subdivision_coordinates).T + + # nztm_meshgrid is a grid of points along a plane with the same orientation + # as the desired plane, but it needs to be translated back to the origin. nztm_meshgrid += origin + return coordinates.nztm_to_wgs_depth(nztm_meshgrid).reshape((ny, nx, 3)) From d64e1ad4b927354977dfaeac0321b1bbe178aefd Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Wed, 29 May 2024 16:41:13 +1200 Subject: [PATCH 12/15] better variable names for write_fault_to_gsf --- qcore/gsf.py | 26 ++++++++++++++------------ 1 file changed, 14 insertions(+), 12 deletions(-) diff --git a/qcore/gsf.py b/qcore/gsf.py index df72c2e8..af2d0853 100644 --- a/qcore/gsf.py +++ b/qcore/gsf.py @@ -187,22 +187,24 @@ def write_fault_to_gsf_file( # Get the number of dip gridpoints by looking at the first dimension of # the meshgrid of the first fault plane. See coordinate_meshgrid for an # explanation of meshgrid dimensions. - n_dip = gsf_df.iloc[0]['meshgrid'].shape[0] + n_dip = gsf_df.iloc[0]["meshgrid"].shape[0] gsf_file_handle.write(f"{number_of_points}\n") - for j in range(n_dip): - for i, row in gsf_df.iterrows(): - length = row["length"] - width = row["width"] - strike = row["strike"] - dip = row["dip"] - rake = row["rake"] - meshgrid = row["meshgrid"] - strike_step = length / gridpoint_count_in_length(length * 1000, resolution) + for dip_index in range(n_dip): + for plane_index, plane in gsf_df.iterrows(): + length = plane["length"] + width = plane["width"] + strike = plane["strike"] + dip = plane["dip"] + rake = plane["rake"] + meshgrid = plane["meshgrid"] + strike_step = length / gridpoint_count_in_length( + length * 1000, resolution + ) dip_step = width / gridpoint_count_in_length(width * 1000, resolution) - for point in meshgrid[j]: + for point in meshgrid[dip_index]: gsf_file_handle.write( - f"{point[1]:11.5f} {point[0]:11.5f} {point[2] / 1000:11.5e} {strike_step:11.5e} {dip_step:11.5e} {strike:6.1f} {dip:6.1f} {rake:6.1f} {-1.0:8.2f} {-1.0:8.2f} {i:3d}\n" + f"{point[1]:11.5f} {point[0]:11.5f} {point[2] / 1000:11.5e} {strike_step:11.5e} {dip_step:11.5e} {strike:6.1f} {dip:6.1f} {rake:6.1f} {-1.0:8.2f} {-1.0:8.2f} {plane_index:3d}\n" ) From 2af37cdb89dfc085eec0c27cc8afa899dd0c144b Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Wed, 29 May 2024 16:42:46 +1200 Subject: [PATCH 13/15] improve docstrings --- qcore/gsf.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/qcore/gsf.py b/qcore/gsf.py index af2d0853..df88c704 100644 --- a/qcore/gsf.py +++ b/qcore/gsf.py @@ -29,7 +29,7 @@ def gridpoint_count_in_length(length: float, resolution: float) -> int: Computes the number of gridpoints that fit into a given length, if each gridpoint is roughly resolution metres apart, and if the gridpoints includes the endpoints. If length = 10, and resolution = 5, then the - function returns 3. + function returns 3 grid points spaced as follows: 5m 5m +-----+-----+ From 59acb57676a9f36c67c8e58e3403158792331893 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Wed, 29 May 2024 16:45:39 +1200 Subject: [PATCH 14/15] provide better formatting for module-level docstrings --- qcore/coordinates.py | 13 +++++++++---- qcore/gsf.py | 5 ++++- 2 files changed, 13 insertions(+), 5 deletions(-) diff --git a/qcore/coordinates.py b/qcore/coordinates.py index f80c5a9e..a1db925e 100644 --- a/qcore/coordinates.py +++ b/qcore/coordinates.py @@ -1,14 +1,19 @@ """ Module for coordinate conversions between WGS84 (latitude and longitude) and NZTM (New Zealand Transverse Mercator) coordinate systems. -This module provides functions for converting coordinates between WGS84 and NZTM coordinate systems. -See linz.govt.nz for a description of the NZTM coordinate system. - -Functions: +Functions +---------- - wgs_depth_to_nztm(wgs_depth_coordinates: np.ndarray) -> np.ndarray: Converts WGS84 coordinates (latitude, longitude, depth) to NZTM coordinates. - nztm_to_wgs_depth(nztm_coordinates: np.ndarray) -> np.ndarray: Converts NZTM coordinates (x, y, depth) to WGS84 coordinates. + +References +---------- +This module provides functions for converting coordinates between WGS84 and NZTM coordinate systems. +See LINZ[0] for a description of the NZTM coordinate system. + +[0]: https://www.linz.govt.nz/guidance/geodetic-system/coordinate-systems-used-new-zealand/projections/new-zealand-transverse-mercator-2000-nztm2000 """ import numpy as np diff --git a/qcore/gsf.py b/qcore/gsf.py index df88c704..d82cfc40 100644 --- a/qcore/gsf.py +++ b/qcore/gsf.py @@ -1,6 +1,7 @@ """This module provides functions for working with and generating GSF files. -Functions: +Functions +--------- - gridpoint_count_in_length: Calculates the number of gridpoints that fit into a given length. - coordinate_meshgrid: @@ -10,6 +11,8 @@ - read_gsf: Parses a GSF file into a pandas DataFrame. +References +---------- The GSF file format is used to define the grid of the source model for a fault. See https://wiki.canterbury.ac.nz/display/QuakeCore/File+Formats+Used+On+GM for details on the GSF format. From 9177c4fe37c0a56a3ed15884badd01d4c5391ac5 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Wed, 29 May 2024 16:45:39 +1200 Subject: [PATCH 15/15] provide better formatting for module-level docstrings --- qcore/coordinates.py | 16 +++++++++++----- qcore/gsf.py | 5 ++++- 2 files changed, 15 insertions(+), 6 deletions(-) diff --git a/qcore/coordinates.py b/qcore/coordinates.py index f80c5a9e..4b9db200 100644 --- a/qcore/coordinates.py +++ b/qcore/coordinates.py @@ -1,14 +1,20 @@ """ -Module for coordinate conversions between WGS84 (latitude and longitude) and NZTM (New Zealand Transverse Mercator) coordinate systems. +Module for coordinate conversions between WGS84 (latitude and longitude) and +NZTM (New Zealand Transverse Mercator) coordinate systems. -This module provides functions for converting coordinates between WGS84 and NZTM coordinate systems. -See linz.govt.nz for a description of the NZTM coordinate system. - -Functions: +Functions +---------- - wgs_depth_to_nztm(wgs_depth_coordinates: np.ndarray) -> np.ndarray: Converts WGS84 coordinates (latitude, longitude, depth) to NZTM coordinates. - nztm_to_wgs_depth(nztm_coordinates: np.ndarray) -> np.ndarray: Converts NZTM coordinates (x, y, depth) to WGS84 coordinates. + +References +---------- +This module provides functions for converting coordinates between WGS84 and NZTM coordinate systems. +See LINZ[0] for a description of the NZTM coordinate system. + +[0]: https://www.linz.govt.nz/guidance/geodetic-system/coordinate-systems-used-new-zealand/projections/new-zealand-transverse-mercator-2000-nztm2000 """ import numpy as np diff --git a/qcore/gsf.py b/qcore/gsf.py index df88c704..d82cfc40 100644 --- a/qcore/gsf.py +++ b/qcore/gsf.py @@ -1,6 +1,7 @@ """This module provides functions for working with and generating GSF files. -Functions: +Functions +--------- - gridpoint_count_in_length: Calculates the number of gridpoints that fit into a given length. - coordinate_meshgrid: @@ -10,6 +11,8 @@ - read_gsf: Parses a GSF file into a pandas DataFrame. +References +---------- The GSF file format is used to define the grid of the source model for a fault. See https://wiki.canterbury.ac.nz/display/QuakeCore/File+Formats+Used+On+GM for details on the GSF format.