sources.py: simplify_fault always absorbs a short interior plane into its right neighbour #114

Location: source_modelling/sources.py:1812-1830 (simplify_fault, interior-plane branch)

What is wrong

The comment says the code picks the neighbour to absorb into "by total deviation". The two numbers it compares are:

geo.point_to_segment_distance(plane.bounds[0],  left.bounds[0],  right.bounds[0])   # top-left of short plane vs top-left(left)->top-left(right)
geo.point_to_segment_distance(plane.bounds[-1], left.bounds[-1], right.bounds[-1])  # bottom-left of short plane vs bottom-left(left)->bottom-left(right)

Both measure the same thing: how far the short plane's left edge is from the line through the left edges of its neighbours. One is taken at the top of the fault and the other at the bottom. A Fault has constant dip and dip direction, so the bottom edge is the top edge translated, and the two distances are equal up to floating-point noise. The < test is therefore effectively always False, and the short plane is always merged into its right neighbour. The deviation of that choice is never measured, so the result ignores the fault's shape.

Reproducer (.ai/simp.py)

The fault trace runs north: P0=(0,0) → P1=(10000,0) → P2=(10500,0) → P3=(20500,1000), in NZTM offsets (m). The middle plane is 0.5 km long. P1 lies on the line P0–P2, so merging the short plane into the left plane changes nothing. P2 is a corner, so merging into the right plane moves the trace.

fault = sources.Fault([sources.Plane.from_nztm_trace(pts[i:i+2].copy(), 0, 10, dip, dip_dir_nztm=ddir) for i in range(3)])
s = sources.simplify_fault(fault, 1.0)

Output (the same for dip 90 and dip 60):

lengths before: [10.          0.5        10.04987562]
lengths after: [10.         10.54751155]
trace after (relative NZTM):
 [[    0.     0.]
 [10000.     0.]
 [10000.     0.]
 [20500.  1000.]]
kept corner P2 (correct: absorb into left plane)? False

The corner at P2 is removed, but the removable collinear vertex P1 is kept.

Suggested fix

Compare the deviation each merge would actually cause, i.e. the distance from the vertex that disappears to the new edge:

# absorb into left: vertex plane.bounds[0] disappears, new top edge is left.bounds[0] -> plane.bounds[1]
dev_left = geo.point_to_segment_distance(plane.bounds[0], left.bounds[0], plane.bounds[1])
# absorb into right: vertex plane.bounds[1] disappears, new top edge is plane.bounds[0] -> right.bounds[1]
dev_right = geo.point_to_segment_distance(plane.bounds[1], plane.bounds[0], right.bounds[1])
if dev_left < dev_right: ...absorb left... else: ...absorb right...

Add a deterministic test using the kinked trace above. test_simplify_fault only checks the resulting lengths and plane counts, so it cannot catch a wrong merge direction.

Filed by the weekly issue sweep. Kind: bug, severity: medium.