Correct moment and depth calculations for velocity model

pyproject.toml

@@ -22,7 +22,6 @@ dependencies = [
   "pandas[parquet, hdf5]",
   "pyyaml",
   "xarray[io]",
-
   # Numerics
   "numpy",
   "scipy",
@@ -35,7 +34,8 @@ dependencies = [
   "schema",    # For loading realisations
   "structlog", # Logging.
   "psutil",    # To get the CPU affinity for jobs
-
+  "parse>=1.21.0",
+  "rich>=14.3.2",
 ]
 
 [project.optional-dependencies]

workflow/scripts/gcmt_to_realisation.py

@@ -40,7 +40,7 @@
 
 from qcore import cli
 from qcore.uncertainties import distributions
-from source_modelling import community_fault_model, magnitude_scaling, sources
+from source_modelling import community_fault_model, magnitude_scaling, moment, sources
 from source_modelling.community_fault_model import NodalPlane
 from workflow import realisations
 from workflow.defaults import DefaultsVersion
@@ -89,9 +89,6 @@ class SourceType(StrEnum):
     """Use a point source approximation."""
 
 
-
-
-
 @cli.from_docstring(app)
 def gcmt_to_realisation(
     gcmt_event_id: Annotated[str, typer.Argument()],
@@ -166,17 +163,17 @@ def gcmt_to_realisation(
     if gcmt_event_id in gcmt_solutions.index:
         row = gcmt_solutions.loc[gcmt_event_id]
         latitude = float(row["Latitude"])  # type: ignore[invalid-argument-type]
-        longitude = float(gcmt_solutions.at[gcmt_event_id, "Longitude"]) # type: ignore[invalid-argument-type]
-        centroid_depth = float(gcmt_solutions.at[gcmt_event_id, "CD"]) # type: ignore[invalid-argument-type]
-        magnitude = float(gcmt_solutions.at[gcmt_event_id, "Mw"]) # type: ignore[invalid-argument-type]
-        
-        strike1 = float(gcmt_solutions.at[gcmt_event_id, "strike1"]) # type: ignore[invalid-argument-type]
-        dip1 = float(gcmt_solutions.at[gcmt_event_id, "dip1"]) # type: ignore[invalid-argument-type]
-        rake1 = float(gcmt_solutions.at[gcmt_event_id, "rake1"]) # type: ignore[invalid-argument-type]
-        
-        strike2 = float(gcmt_solutions.at[gcmt_event_id, "strike2"]) # type: ignore[invalid-argument-type]
-        dip2 = float(gcmt_solutions.at[gcmt_event_id, "dip2"]) # type: ignore[invalid-argument-type]
-        rake2 = float(gcmt_solutions.at[gcmt_event_id, "rake2"]) # type: ignore[invalid-argument-type]
+        longitude = float(gcmt_solutions.at[gcmt_event_id, "Longitude"])  # type: ignore[invalid-argument-type]
+        centroid_depth = float(gcmt_solutions.at[gcmt_event_id, "CD"])  # type: ignore[invalid-argument-type]
+        solution_moment = float(gcmt_solutions.at[gcmt_event_id, "Mo"])  # type: ignore[invalid-argument-type]
+
+        strike1 = float(gcmt_solutions.at[gcmt_event_id, "strike1"])  # type: ignore[invalid-argument-type]
+        dip1 = float(gcmt_solutions.at[gcmt_event_id, "dip1"])  # type: ignore[invalid-argument-type]
+        rake1 = float(gcmt_solutions.at[gcmt_event_id, "rake1"])  # type: ignore[invalid-argument-type]
+
+        strike2 = float(gcmt_solutions.at[gcmt_event_id, "strike2"])  # type: ignore[invalid-argument-type]
+        dip2 = float(gcmt_solutions.at[gcmt_event_id, "dip2"])  # type: ignore[invalid-argument-type]
+        rake2 = float(gcmt_solutions.at[gcmt_event_id, "rake2"])  # type: ignore[invalid-argument-type]
 
         nodal_plane_1 = NodalPlane(strike1, dip1, rake1)
         nodal_plane_2 = NodalPlane(strike2, dip2, rake2)
@@ -185,16 +182,18 @@ def gcmt_to_realisation(
         latitude = solution["location"]["latitude"]
         longitude = solution["location"]["longitude"]
         centroid_depth = solution["location"]["depth"]
-        magnitude = solution["magnitude"]
+        solution_moment = float(solution["moment"])
         nodal_plane_1 = NodalPlane(**solution["nodalPlanes"][0])
         nodal_plane_2 = NodalPlane(**solution["nodalPlanes"][1])
-        
+
     else:
         raise typer.BadParameter(
             f"GCMT event ID {gcmt_event_id} not found in either the published GCMT solutions or automated solutions.",
             param_hint="GCMT_EVENT_ID",
         )
 
+    magnitude = moment.moment_to_magnitude(solution_moment)
+
     model = community_fault_model.get_community_fault_model()
 
     match nodal_plane:

workflow/scripts/migrate.py

@@ -19,7 +19,7 @@
 from qcore import cli
 from workflow import realisations, utils
 from workflow.defaults import DefaultsVersion
-from workflow.realisations import RealisationMetadata, Seeds
+from workflow.realisations import Seeds
 
 app = typer.Typer()
 console = Console()

workflow/scripts/realisation_to_srf.py

@@ -652,8 +652,10 @@ def generate_point_source_srf(
     magnitude = params.magnitudes.magnitudes[name]
     moment_newton_metre = moment.magnitude_to_moment(magnitude)
 
-    velocity_model_df = params.velocity_model_1d.model
-    velocity_model_df["depth_km"] = velocity_model_df["thickness"].cumsum()
+    velocity_model_df = params.velocity_model_1d.model.copy()
+    velocity_model_df["depth_km"] = (
+        velocity_model_df["thickness"].cumsum() - velocity_model_df["thickness"]
+    )
 
     # Get the source depth
     # divide by 1000 to convert depth from meters to kilometers

workflow/utils.py

@@ -1,6 +1,5 @@
 """Miscellaneous workflow utilities that couldn't go anywhere else."""
 
-import inspect
 import hashlib
 import os
 import tempfile
@@ -15,7 +14,6 @@
 from shapely import Geometry, Polygon, geometry
 
 from qcore import coordinates
-from workflow import defaults
 
 NZ_COASTLINE_URL = "https://www.dropbox.com/scl/fi/zkohh794y0s2189t7b1hi/NZ.gmt?rlkey=02011f4morc4toutt9nzojrw1&st=vpz2ri8x&dl=1"