From ccd6ca8b2f36bd181040d862de7ff5ea39dca188 Mon Sep 17 00:00:00 2001 From: Jake Faulkner Date: Tue, 8 Sep 2026 23:22:01 +1200 Subject: [PATCH] Report the real cause when a vertical plane is queried without depth coordinate_length was forced to 3 whenever dip == 90, regardless of the input's dimensionality, so a 2D (lat, lon) query against a vertical plane produced a broadcast ValueError. Fault's per-plane loop caught that and re-raised "Given coordinates are not on fault" -- for a point that is on the fault. A plane dipping 84 degrees handled the same query fine, so only the exactly-vertical case failed, and dip == 90 is common in NZ fault models. Sizes the slices from the input alone, and raises a specific error for the vertical-without-depth case: a vertical plane projects onto a line in plan view, so a (lat, lon) pair maps to every depth on the plane and the dip coordinate is genuinely undetermined. Fabricating one would trade a confusing error for a silently wrong answer. Adds CoordinatesNotOnPlaneError (a ValueError subclass, so external callers catching ValueError are unaffected) for the genuine off-plane miss, so Fault's loop only swallows real misses and the undetermined-dip error surfaces to the caller instead of being reported as a geometric one. Fixes #85 Co-Authored-By: Claude Opus 5 --- source_modelling/sources.py | 38 +++++++++++++++++++++----- tests/test_sources.py | 54 +++++++++++++++++++++++++++++++++++++ 2 files changed, 85 insertions(+), 7 deletions(-) diff --git a/source_modelling/sources.py b/source_modelling/sources.py index 162971c..ad02ffe 100644 --- a/source_modelling/sources.py +++ b/source_modelling/sources.py @@ -28,6 +28,16 @@ _KM_TO_M = 1000 +class CoordinatesNotOnPlaneError(ValueError): + """Raised when global coordinates do not lie within a plane. + + Subclasses :class:`ValueError` for backwards compatibility. It exists so + that "this point is not on this plane" can be told apart from "this point + cannot be located on this plane", which callers iterating over planes must + not treat as a miss. + """ + + @dataclasses.dataclass class Point: """A representation of a point source. @@ -778,19 +788,31 @@ def wgs_depth_coordinates_to_fault_coordinates( Raises ------ - ValueError + CoordinatesNotOnPlaneError If the given coordinates do not lie in the fault plane. + ValueError + If the plane is vertical (``dip == 90``) and no depth is given, + because the dip coordinate is then undetermined. Notes ----- While not passing depth information is supported, depth information *greatly* improves the accuracy of the estimation. No guarantees are made about the accuracy of the inversion if you do not pass - depth information. + depth information. Vertical planes are the exception: they project + onto a line in plan view, so depth is required rather than merely + recommended. """ - coordinate_length = ( - 3 if global_coordinates.shape[-1] == 3 or self.dip == 90 else 2 - ) + coordinate_length = 3 if global_coordinates.shape[-1] == 3 else 2 + if coordinate_length == 2 and self.dip == 90: + # A vertical plane projects onto a line in plan view, so a + # (lat, lon) pair maps to every depth on the plane and the dip + # coordinate is genuinely undetermined. Fail loudly rather than + # invent one. + raise ValueError( + "Depth is required to locate coordinates on a vertical plane " + "(dip == 90); the dip coordinate is undetermined without it." + ) strike_direction = ( self.bounds[1, :coordinate_length] - self.bounds[0, :coordinate_length] ) @@ -815,7 +837,9 @@ def wgs_depth_coordinates_to_fault_coordinates( | np.isclose(fault_local_coordinates, 1, atol=tolerance) ) ): - raise ValueError("Specified coordinates do not lie in plane") + raise CoordinatesNotOnPlaneError( + "Specified coordinates do not lie in plane" + ) return np.clip(fault_local_coordinates, 0, 1) def rrup_distance(self, points: np.ndarray) -> np.ndarray | float: @@ -1349,7 +1373,7 @@ def wgs_depth_coordinates_to_fault_coordinates( return np.array([left_edges[i], 0]) + plane_coordinates * np.array( [left_edges[i + 1] - left_edges[i], 1] ) - except ValueError: + except CoordinatesNotOnPlaneError: continue raise ValueError("Given coordinates are not on fault.") diff --git a/tests/test_sources.py b/tests/test_sources.py index 72b304a..17567f5 100644 --- a/tests/test_sources.py +++ b/tests/test_sources.py @@ -1490,3 +1490,57 @@ def test_multi_fault_rx_ry(): assert rx == pytest.approx(0.0) assert ry == pytest.approx(0.0) + + +def _vertical_plane() -> Plane: + """Build a strictly vertical plane (dip == 90) from NZTM corners.""" + origin = coordinates.wgs_depth_to_nztm(np.array([-43.5, 172.6, 0.0])) + along_strike = np.array([10000.0, 10000.0, 0.0]) + down_dip = np.array([0.0, 0.0, 10000.0]) + return Plane( + np.array( + [ + origin, + origin + along_strike, + origin + along_strike + down_dip, + origin + down_dip, + ] + ) + ) + + +def test_vertical_plane_with_depth_is_unaffected(): + """A vertical plane still inverts exactly when depth is supplied.""" + plane = _vertical_plane() + assert plane.dip == 90.0 + assert np.allclose( + plane.wgs_depth_coordinates_to_fault_coordinates(plane.centroid), [0.5, 0.5] + ) + + +def test_vertical_plane_without_depth_reports_the_real_cause(): + """A vertical plane queried without depth must say depth is required. + + Regression test: ``coordinate_length`` was forced to 3 whenever + ``dip == 90``, regardless of the input's dimensionality, so a 2D query + raised a broadcast error ("operands could not be broadcast together with + shapes (2,) (3,)") which ``Fault`` then swallowed and reported as "not on + fault" -- for a point that is on the fault. + """ + plane = _vertical_plane() + centroid_2d = plane.centroid[:2] + + with pytest.raises(ValueError, match="Depth is required"): + plane.wgs_depth_coordinates_to_fault_coordinates(centroid_2d) + + # and the cause must survive Fault's per-plane loop rather than being + # reported as a geometric miss + with pytest.raises(ValueError, match="Depth is required"): + Fault([plane]).wgs_depth_coordinates_to_fault_coordinates(centroid_2d) + + +def test_fault_still_reports_genuine_misses_as_not_on_fault(): + """Points genuinely off the fault must still raise "not on fault".""" + fault = Fault([_vertical_plane()]) + with pytest.raises(ValueError, match="not on fault"): + fault.wgs_depth_coordinates_to_fault_coordinates(np.array([-41.0, 174.0, 0.0]))