Report the real cause when a vertical plane is queried without depth

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.")
 

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]))