energyml_section/
erosion.rs

1//! Erosions suspected from the surfaces' ages (ADR-0015):
2//! detect erosions from the ages of the surfaces, propose local corrections or
3//! suggestions of corrections when something is wrong, and colour the portions
4//! identified as erosions).
5//!
6//! # The rule
7//!
8//! On one section, take two **horizons** that both state an age
9//! ([`TraceItem::age`], an interval `[old, young]` in Ma). `A` is older than `B`
10//! only when the two intervals are **disjoint** - `A.young > B.old`; an
11//! overlap, or two equal ages, orders nothing and is not judged. An older
12//! horizon lies **above** a younger one nowhere, so wherever, at the same `s`,
13//! the younger one is below the older one (`z` up: `z_young < z_old`), the
14//! order of the surfaces contradicts the order of the ages. Two readings fit:
15//! the younger surface truncates the older one (an erosion: the older one should
16//! stop against it), or one of the two ages is wrong. Which one is the
17//! geologist's to say, so the module reports the **portion** `[s0, s1]` and
18//! proposes the first reading in words. **A proposal is never applied.**
19//!
20//! What is not judged: a horizon without an age (the rule - most files date
21//! nothing, and nothing is inferred from a name or a depth), a fault (it does
22//! not stack), a span where only one of the two lines exists (a surface that
23//! stops is not a surface that crosses), and a grid's envelope.
24//!
25//! Pure: [`detect`] takes traces, [`of_section`] a section, [`notes`] a
26//! catalogue; nothing is read.
27
28use crate::model::{Catalogue, Item, Note, ObjectId, Section, SectionKey, TraceItem, TraceRole};
29
30/// The note's code ([`Note::code`]).
31pub const CODE: &str = "erosion";
32
33/// How far below, in the section's units, counts as below: float noise is not
34/// an erosion.
35const Z_TOLERANCE: f64 = 1e-6;
36/// A portion shorter than this, along `s`, is a crossing point, not a portion.
37const MIN_LENGTH: f64 = 1e-6;
38
39/// A stretch of a section where a younger horizon passes below an older one.
40#[derive(Clone, Debug, PartialEq)]
41pub struct Portion {
42    /// On which section.
43    pub section: SectionKey,
44    /// The horizon the ages say is older - the one that should stop.
45    pub older: ObjectId,
46    /// Its age, `[old, young]` in Ma.
47    pub older_age: [f64; 2],
48    /// The horizon the ages say is younger - the one that truncates.
49    pub younger: ObjectId,
50    /// Its age, `[old, young]` in Ma.
51    pub younger_age: [f64; 2],
52    /// Where it starts along the section.
53    pub s0: f64,
54    /// Where it ends.
55    pub s1: f64,
56}
57
58impl Portion {
59    /// The correction proposed, in words - and only in words.
60    pub fn proposal(&self) -> String {
61        format!(
62            "`{older}` should stop against `{younger}` between s = {s0} and s = {s1} (an erosion: the younger surface \
63             truncates the older one), or one of the two ages is inverted; nothing was changed",
64            older = self.older.title,
65            younger = self.younger.title,
66            s0 = round1(self.s0),
67            s1 = round1(self.s1),
68        )
69    }
70
71    /// The sentence of the note.
72    pub fn said(&self) -> String {
73        format!(
74            "erosion suspected on {section}: `{younger}` ({ya}) passes below `{older}` ({oa}) between s = {s0} and s = {s1}, \
75             though it is the younger - proposal: {proposal}",
76            section = section_name(self.section),
77            younger = self.younger.title,
78            ya = age_text(self.younger_age),
79            older = self.older.title,
80            oa = age_text(self.older_age),
81            s0 = round1(self.s0),
82            s1 = round1(self.s1),
83            proposal = self.proposal(),
84        )
85    }
86}
87
88/// The portions of two or more dated traces, by pair of horizons: for each
89/// ordered pair (the older first in the item order), each stretch where the
90/// younger passes below the older, left to right. Traces with no age, faults,
91/// and pairs the ages do not order give nothing.
92pub fn detect(section: SectionKey, traces: &[&TraceItem]) -> Vec<Portion> {
93    let dated: Vec<(&TraceItem, [f64; 2])> = traces
94        .iter()
95        .filter(|t| t.role == TraceRole::Horizon)
96        .filter_map(|t| t.age.map(|age| (*t, age)))
97        .collect();
98    let mut out = Vec::new();
99    for (i, (a, a_age)) in dated.iter().enumerate() {
100        for (b, b_age) in dated.iter().skip(i + 1) {
101            // The older first: disjoint intervals only.
102            let (older, o_age, younger, y_age) = if a_age[1] > b_age[0] {
103                (*a, *a_age, *b, *b_age)
104            } else if b_age[1] > a_age[0] {
105                (*b, *b_age, *a, *a_age)
106            } else {
107                continue;
108            };
109            for (s0, s1) in below(&older.spans, &younger.spans) {
110                out.push(Portion {
111                    section,
112                    older: older.object.clone(),
113                    older_age: o_age,
114                    younger: younger.object.clone(),
115                    younger_age: y_age,
116                    s0,
117                    s1,
118                });
119            }
120        }
121    }
122    out
123}
124
125/// [`detect`] over what a section draws.
126pub fn of_section(section: &Section) -> Vec<Portion> {
127    let traces: Vec<&TraceItem> = section
128        .items
129        .iter()
130        .filter_map(|item| match item {
131            Item::Trace(t) => Some(t),
132            _ => None,
133        })
134        .collect();
135    detect(section.key, &traces)
136}
137
138/// One note per portion, over the sections given (the long one first).
139pub fn notes_of<'a>(sections: impl IntoIterator<Item = &'a Section>) -> Vec<Note> {
140    sections
141        .into_iter()
142        .flat_map(of_section)
143        .map(|portion| Note {
144            code: CODE.to_string(),
145            sections: vec![portion.section],
146            objects: vec![portion.older.clone(), portion.younger.clone()],
147            said: portion.said(),
148        })
149        .collect()
150}
151
152/// [`notes_of`] over a finished catalogue.
153pub fn notes(catalogue: &Catalogue) -> Vec<Note> {
154    notes_of(std::iter::once(&catalogue.long).chain(&catalogue.crosses))
155}
156
157/// The `z` of a polyline set at `s`, when one of its spans covers `s`.
158fn z_at(spans: &[Vec<[f64; 2]>], s: f64) -> Option<f64> {
159    spans.iter().find_map(|span| {
160        let (first, last) = (span.first()?, span.last()?);
161        if s < first[0] || s > last[0] {
162            return None;
163        }
164        span.windows(2).find(|w| s >= w[0][0] && s <= w[1][0]).map(|w| {
165            let run = w[1][0] - w[0][0];
166            if run <= 0.0 {
167                w[0][1]
168            } else {
169                w[0][1] + (w[1][1] - w[0][1]) * (s - w[0][0]) / run
170            }
171        })
172    })
173}
174
175/// The stretches of `s` where `lower` is below `upper` - both exist there, and
176/// `z_lower < z_upper` by more than the tolerance. Between two breakpoints (a
177/// vertex of either) the difference is linear, so its sign change is exact.
178fn below(upper: &[Vec<[f64; 2]>], lower: &[Vec<[f64; 2]>]) -> Vec<(f64, f64)> {
179    let mut at: Vec<f64> = upper
180        .iter()
181        .chain(lower)
182        .flatten()
183        .map(|p| p[0])
184        .filter(|s| s.is_finite())
185        .collect();
186    at.sort_by(f64::total_cmp);
187    at.dedup();
188    let mut found: Vec<(f64, f64)> = Vec::new();
189    for pair in at.windows(2) {
190        let (a, b) = (pair[0], pair[1]);
191        let mid = 0.5 * (a + b);
192        // Both lines must run across the whole segment: probe the middle too,
193        // so a gap between two spans is not bridged.
194        let (Some(ua), Some(ub), Some(la), Some(lb), Some(_), Some(_)) = (
195            z_at(upper, a),
196            z_at(upper, b),
197            z_at(lower, a),
198            z_at(lower, b),
199            z_at(upper, mid),
200            z_at(lower, mid),
201        ) else {
202            continue;
203        };
204        let (fa, fb) = (la - ua, lb - ub);
205        let neg = |f: f64| f < -Z_TOLERANCE;
206        let part = match (neg(fa), neg(fb)) {
207            (true, true) => Some((a, b)),
208            (true, false) => Some((a, a + (b - a) * fa / (fa - fb))),
209            (false, true) => Some((a + (b - a) * fa / (fa - fb), b)),
210            (false, false) => None,
211        };
212        if let Some((s0, s1)) = part.filter(|(s0, s1)| s1 - s0 > MIN_LENGTH) {
213            match found.last_mut() {
214                Some(last) if (s0 - last.1).abs() <= MIN_LENGTH => last.1 = s1,
215                _ => found.push((s0, s1)),
216            }
217        }
218    }
219    found
220}
221
222fn section_name(key: SectionKey) -> String {
223    match key {
224        SectionKey::Long => "the long section".to_string(),
225        SectionKey::Cross { index } => format!("cross section {index}"),
226    }
227}
228
229/// `10 Ma`, or `10-5 Ma` for a span.
230fn age_text(age: [f64; 2]) -> String {
231    if (age[0] - age[1]).abs() < 1e-9 {
232        format!("{} Ma", trim(age[0]))
233    } else {
234        format!("{}-{} Ma", trim(age[0]), trim(age[1]))
235    }
236}
237
238fn trim(v: f64) -> String {
239    let text = format!("{v:.3}");
240    text.trim_end_matches('0').trim_end_matches('.').to_string()
241}
242
243fn round1(v: f64) -> String {
244    format!("{v:.1}")
245}
energyml_section/erosion.rs

energyml_section/
erosion.rs
