sources.py: vertical planes (dip == 90) reject 2D coordinates as "not on fault"

[lispandfound]

Location: source_modelling/sources.py:791

What happens: Plane.wgs_depth_coordinates_to_fault_coordinates computes

coordinate_length = (
    3 if global_coordinates.shape[-1] == 3 or self.dip == 90 else 2
)

When dip == 90 and the caller passes a 2D (lat, lon) point, coordinate_length becomes 3 while the input only has 2 components. The offset expression at sources.py:799-802 then subtracts a 3-vector from a 2-vector and raises ValueError: operands could not be broadcast together with shapes (2,) (3,). Fault.wgs_depth_coordinates_to_fault_coordinates catches that ValueError at sources.py:1352 and converts it into raise ValueError("Given coordinates are not on fault.") (sources.py:1354) — so a point that lies exactly on the fault is silently reported as off it.

Why that is wrong: The method docstring explicitly supports 2D input — sources.py:786-789: "While not passing depth information is supported, depth information greatly improves the accuracy of the estimation." The or self.dip == 90 clause is presumably meant to force the 3D path when depth is available (a vertical plane is degenerate in plan view), but it fires regardless of the input dimensionality. The generic except ValueError at sources.py:1352 then masks the shape error as a geometric one.

How to reproduce:

import numpy as np
from qcore import coordinates
from source_modelling.sources import Plane, Fault

o = coordinates.wgs_depth_to_nztm(np.array([-43.5, 172.6, 0.0]))
along = np.array([10000.0, 10000.0, 0.0])   # along strike, 45 deg
down  = np.array([0.0, 0.0, 10000.0])       # straight down => dip 90
p = Plane(np.array([o, o + along, o + along + down, o + down]))
assert p.dip == 90.0

p.wgs_depth_coordinates_to_fault_coordinates(p.centroid)       # [0.5 0.5]  ok
p.wgs_depth_coordinates_to_fault_coordinates(p.centroid[:2])   # ValueError: shapes (2,) (3,)
Fault([p]).wgs_depth_coordinates_to_fault_coordinates(p.centroid[:2])
# ValueError: Given coordinates are not on fault.   <- the centroid IS on the fault

Control: rebuilding the same plane with down = [700.0, -700.0, 10000.0] gives dip == 84.35 and the identical 2D query returns [0.5 0.5]. Only the exactly-vertical case fails, and dip == 90 is common in NZ fault models.

Suggested direction: Decide coordinate_length from the input shape alone, and handle the vertical-plane degeneracy separately (e.g. project onto the trace when no depth is supplied). Narrowing the except ValueError at sources.py:1352 so it only swallows genuine off-fault errors would stop shape bugs from masquerading as geometry.

Confidence: high

[lispandfound]

Intended behaviour as inversion is ill-conditioned for vertical fault with 2D input.