From e1d231114403038ee63734231bed7892d0b99332 Mon Sep 17 00:00:00 2001 From: xrfez Date: Wed, 11 Nov 2026 18:29:42 +0000 Subject: [PATCH 1/3] Survey join: match traces by their stations; O (receiver-domain) relations Field record numbers and channel numbers in the data often differ from the relation file's (a recorder numbers all its live channels; renumbered or reshot records), while the header stations are right. The join now matches by stations and never replaces a header station: - both stations in the header: coordinates by station only; - source station only: the receiver from the X records of that source station (the one with the trace's record number among repeats, else they must agree); - receiver station only: the source from the O gathers of that receiver, or the X relation of record and channel if its receiver is the header's; - no station: record and channel, the fallback. `overwrite` now replaces coordinates only. Auxiliary traces are left alone. `joinStations` exposes the station step; `JoinMatch` and `JoinStats` (from new `applySurvey` overloads) say how traces matched. The CLI prints them. Relations of both domains (`RelationDomain`): O records (record letter O, or X records in an `.ops` file) are the X layout with source and receiver swapped. New lookups: `records`, `relations`, `expectedChannels`, `recordsAt`, `channelOf`, `sameStation`. `receiverRelations` builds the O records of an X survey, one receiver line in memory at a time. Also fixes writeGeometry/copyWithGeometry reporting a value that rounds to the bytes already in the file (shotpoint 99.5 into an espnum holding 100) as a write the journal dropped (a fatal error): such a trace is no change. Docs (conversion guide, api.md, SPS spec notes) and tests follow; tests that relied on the old (record, channel)-first join now use headers without the station they fill. Co-Authored-By: Claude Opus 5.5 --- docs/api.md | 13 +- docs/guide/conversion.md | 76 +++-- docs/spec/sps.md | 17 + src/parseSEG/cli.nim | 15 +- src/parseSEG/survey/join.nim | 512 +++++++++++++++++++++++++----- src/parseSEG/survey/sps.nim | 29 +- src/parseSEG/survey/writeback.nim | 4 +- tests/tcli.nim | 50 +-- tests/tgeomwrite.nim | 35 +- tests/tsurvey.nim | 238 +++++++++++++- 10 files changed, 828 insertions(+), 161 deletions(-) diff --git a/docs/api.md b/docs/api.md index eb986f6..5f211c2 100644 --- a/docs/api.md +++ b/docs/api.md @@ -100,18 +100,21 @@ without keys, SEG-D channel types other than seismic or EM and SEG-2 `TEST_DATA` traces (SEG-Y traces then count as seismic). ## Survey files and external geometry (Phase 10) -SPS (Rev 1.0, 2.1) and SEG P1 (1983) give the coordinates that SEG-D 2.x headers lack, and SPS -X records name the stations each channel was recorded from. Spec notes: `docs/spec/sps.md`, +SPS (Rev 1.0, 2.1) and SEG P1 (1983) give the coordinates that SEG-D 2.x headers lack, SPS +X records name the stations each channel was recorded from, and O records (`.ops`) do the same +in the receiver domain. Traces are joined by their stations; record numbers are a fallback. Spec notes: `docs/spec/sps.md`, `docs/spec/p1_p2_p3.md`. | Proc | What | |---|---| -| `readSps(path or stream, SpsOptions)`, `parseSps(text)` → `SpsFile` | H, R/S (`points`, `kind` `pkReceiver`/`pkSource`), X (`relations`) and C records; revision from H00 or the columns. Every record keeps its raw text. Blank integers are `MissingInt`, blank reals NaN. | +| `readSps(path or stream, SpsOptions)`, `parseSps(text)` → `SpsFile` | H, R/S (`points`, `kind` `pkReceiver`/`pkSource`), X and O (`relations`, `recordType`) and C records; revision from H00 or the columns. Every record keeps its raw text. Blank integers are `MissingInt`, blank reals NaN. | | `writeSps(f, path or stream, SpsWriteOptions(revision, lineEnding, lenient))`, `spsText(f)` (a `SurveyText`: `text`, `log`), `encodeSps(record, rev)` | Column-exact; unchanged records byte for byte; 1.0 ↔ 2.1. A value that does not fit is an error and nothing is written (unless `lenient`). | | `readP1(path or stream, P1Options(gridDecimals, depthDecimals, angleUnit))`, `parseP1` → `P1File` | Logical files (`files`), each with its header text and records; P1 decoded, P2-P4 kept raw. | | `writeP1(f, path or stream, P1WriteOptions)`, `p1Text(f)` (a `SurveyText`), `encodeP1(record, opts)` | Header blocks padded to 20 records; unchanged records byte for byte. | -| `SurveyGeometry` (`initSurveyGeometry()`): `addSps(g, sps or path)`, `addP1(g, p1 or path[, P1Options], elevation)`, `addSource`/`addReceiver`, `relation(g, record, channel)` (an `Option[StationRelation]`), `resolveLine(g, name, numeric)` (an `Option`) | Station coordinates (`table`, `sourceCount`, `receiverCount`) and relation runs (`relationCount`). Text line names go through `mapLine(g, name, line)` (`lineNames`), then `trailingDigits`. `multiComponent`, `tape` refine X records; problems are in `log`. | -| `applySurvey(f, keys or index, g, overwrite = false)` | Stations from relations (on `TraceKeys.record` and `.channel`), then coordinates, offset, midpoint; in memory only, the file does not change. | +| `SurveyGeometry` (`initSurveyGeometry()`): `addSps(g, sps, xDomain = rdSource)`, `addSps(g, path, opts, xDomain = relationDomainOf(path))`, `addP1(g, p1 or path[, P1Options], elevation)`, `addSource`/`addReceiver`, `relation(g, record, channel, domain = rdSource)` (an `Option[StationRelation]`), `resolveLine(g, name, numeric)` (an `Option`) | Station coordinates (`table`, `sourceCount`, `receiverCount`) and relation runs (`relationCount(g, domain)`) of both domains (`RelationDomain`: `rdSource` X, `rdReceiver` O; X records of an `.ops` file are O by default). Text line names go through `mapLine(g, name, line)` (`lineNames`), then `trailingDigits`. `multiComponent`, `tape` refine relations; problems are in `log`. | +| `records(g, domain)`, `relations(g, record, domain)` (iterators: ascending records; `(channel, relation)` per channel, a channel two runs cover once), `expectedChannels(g, record, domain)`, `recordsAt(g, station, domain)`, `channelOf(g, record, station, domain)` (an `Option`), `sameStation(a, b)` | Enumerate and look up the relations: the records of a source station (X) or gathers of a receiver station (O), the channel of a station in a record. Index 0 matches any. | +| `receiverRelations(g, firstGather = 1, recordType = 'X')` → `seq[SpsRelation]` | The O records of the X relations: one gather per receiver station (receiver line, point order), its sources by line and point as runs, repeats after; write with `writeSps`. One receiver line in memory at a time. | +| `joinStations(keys, g)` → `JoinMatch`; `applySurvey(f, keys or index, g[, stats], overwrite = false)` | Fill the stations a header lacks, matched **by station**: both in the header (`jmStations`); source only → the X records of that source (`jmSource`); receiver only → the O gathers of that receiver (`jmReceiver`); no station → record and channel (`jmRecord`); also `jmAmbiguous`, `jmConflict`, `jmUnmatched`, `jmAuxiliary` (left alone). Header stations are never replaced; a relation must agree with the stations the header has. `applySurvey` then fills coordinates by station, offset, midpoint (`overwrite`: replace known coordinates); in memory only. `JoinStats` (`traces`, `changed`, `matches`, `noSourceXy`, `noReceiverXy`; `$`). | | `geometryChanges(f, fields[, log])` | What writing the joined keys of `f` (opened with keys, `applySurvey` done) would change, as `GeometryChanges` (traces per kind, `unwritable` values; `$` gives a one-line summary), with the reports in `log`; nothing is written. | | `writeGeometry(editor, fields = GeometryFields())`, `writeGeometry(editor, f.index, fields)` | Writes the joined geometry in place through the journaled `HeaderEditor`: from the editor's own index (an editor made with `openEditor(f)` from the joined file; `ValueError` without keys), or from the index of the same file opened separately. Same-size fields; undo with `rollbackJournal`. SEG-Y: through the layout (`sht_x`, `rec_x`, `cdp_x` with a coordinate scalar chosen to keep the precision in 4 bytes; `selev`/`relev`; `offset`; `espnum`/`sp`; Ext1 copies in use; `sht_line`, `rec_point`, ... where a Layout stanza defines them). SEG-D Rev 1/2: THE#1 receiver line/point/index. Rev 3: THE#1 and the Position / Source Description blocks already in the trace header extensions. Values with no slot or too large for theirs are reported in `editor.log`. `GeometryFields` selects source/receiver X/Y, elevations, midpoints, offsets, stations. SEG-2 is refused for now. | | `copyWithGeometry(src, dst, g, fields, overwrite)`, `convertWithGeometry(src, dst, g, opts, fields, overwrite)` | A copy in the same standard (multiplexed SEG-D demultiplexed), or a `convertFile` SEG-Y output, with `g` joined and written. `copyWithGeometry` returns a `GeometryCopy` (`changes`, `log`), `convertWithGeometry` a `DiagLog`; both logs carry an info summary. | diff --git a/docs/guide/conversion.md b/docs/guide/conversion.md index ba6e201..a8d1ede 100644 --- a/docs/guide/conversion.md +++ b/docs/guide/conversion.md @@ -89,21 +89,38 @@ tape devices included. SEG-D 2.x headers usually carry station numbers but no coordinates, and some files not even station numbers per channel. SPS files (Rev 1.0, 2.1) give the coordinates of every source (S records) and receiver (R records), and X records say which stations each channel of each -field record was recorded from. SEG P1 files give shotpoint coordinates. +field record was recorded from. O records (often in `.ops` files) are the same relations in +the receiver domain: the X layout with source and receiver swapped, one gather per receiver. +SEG P1 files give shotpoint coordinates. A `SurveyGeometry` collects them; `applySurvey` joins them to a file's keys **in memory** -(the file does not change): stations from the relations (matched on `TraceKeys.record` and -`.channel`), then coordinates from the stations, then offset and midpoint. +(the file does not change). **Traces are matched by their stations** (line, point, index), +not by record number: field record and channel numbers often differ from the relation +file's. Stations a header has are never replaced; the relations only fill the ones it lacks: + +- source and receiver station in the header: coordinates by station, nothing else; +- source station only: the X records **of that source station** give the receiver of the + trace's channel (the one with the trace's record number when there are several repeats, + else they must agree); +- receiver station only: the O gathers of that receiver give the source likewise (without + O records, the X relation of the trace's record and channel, if its receiver is the + header's); +- no station: the trace's record and channel, the fallback. + +Then coordinates come from the stations, then offset and midpoint. `overwrite` replaces +coordinates already known, never stations. Auxiliary traces are left alone. A `JoinStats` +says how the traces matched. ```nim import std/math import parseSEG -# A small SPS 2.1 file: one receiver, one source, and the relation of record 63, channel 1. +# A small SPS 2.1 file: record.sgd's receiver and source, and the relation of +# record 63, channel 1. writeFile("line.sps", """H00 SPS format version number SPS 2.1; -R 10.00 101.00 1G1 500100.0 6000000.0 12.5 18213250 -S 20.00 205.00 1V6 500000.0 6000400.0 10.0 19001150 -XT1 6311 20.00 205.001 1 11 10.00 101.00 101.001 +R 10.00 104.00 1G1 500100.0 6000000.0 12.5 18213250 +S 20.00 99.50 1V6 500000.0 6000400.0 10.0 19001150 +XT1 6311 20.00 99.501 1 11 10.00 104.00 104.001 """) var g = initSurveyGeometry() @@ -114,13 +131,29 @@ echo g.sourceCount, " source(s), ", g.receiverCount, " receiver(s), ", var f = openSeismic("record.sgd", keys = true) doAssert f.traceKeys(0).rcvX.isNaN # the SEG-D 2.0 headers have no coordinates -let changed = f.applySurvey(g, overwrite = true) +var stats: JoinStats +let changed = f.applySurvey(g, stats) +echo stats # 1 traces, 1 changed; by stations 1 let k = f.traceKeys(0) echo changed, " trace(s) joined: receiver ", k.rcvLine, "/", k.rcvPoint, " at (", k.rcvX, ", ", k.rcvY, "), offset ", round(k.offset, 1) -doAssert changed == 1 and k.rcvX == 500100 and k.srcPoint == 205 +doAssert changed == 1 and k.rcvX == 500100 and k.srcPoint == 99.5 +doAssert stats.matches[jmStations] == 1 # the header's own stations were the key + +# The receiver domain: O records built from the X relations, written as an +# .ops file (X records, as such files hold them) and read back. +let ops = SpsFile(revision: sps21, headers: readSps("line.sps").headers, + relations: g.receiverRelations()) +doAssert not writeSps(ops, "line.ops").hasErrors +var go = initSurveyGeometry() +go.addSps("line.ops") # `.ops`: its X records are receiver gathers +doAssert go.relationCount(rdReceiver) == 1 +doAssert go.recordsAt(stationKey(10, 104, 1), rdReceiver) == @[1'i64] ``` +- `records(g)`, `relations(g, record)`, `expectedChannels`, `recordsAt(g, station)` (the + records of a source station, or the gathers of a receiver) and `channelOf` enumerate and + look up the relations, in either domain (`rdSource`, `rdReceiver`). - Line names that are not numbers (SPS 1.0, P1) are mapped with `g.mapLine(name, number)`, or by their trailing digits with `g.trailingDigits = true`. - `g.multiComponent` makes X records with a channel increment cover consecutive component @@ -142,37 +175,34 @@ import std/os import parseSEG writeFile("line.sps", """H00 SPS format version number SPS 2.1; -R 10.00 101.00 1G1 500100.0 6000000.0 12.5 18213250 -S 20.00 205.00 1V6 500000.0 6000400.0 10.0 19001150 -XT1 6311 20.00 205.001 1 11 10.00 101.00 101.001 +R 10.00 104.00 1G1 500100.0 6000000.0 12.5 18213250 +S 20.00 99.50 1V6 500000.0 6000400.0 10.0 19001150 +XT1 6311 20.00 99.501 1 11 10.00 104.00 104.001 """) var g = initSurveyGeometry() g.addSps("line.sps") -# 1. In place, journaled and undoable: same-size header fields only. +# 1. In place, journaled and undoable: same-size header fields only. record.sgd's +# SEG-D 2.0 headers already hold its stations, and have no coordinate fields. copyFile("record.sgd", "field.sgd") var f = openSeismic("field.sgd", keys = true) -discard f.applySurvey(g, overwrite = true) +discard f.applySurvey(g) var report: DiagLog let plan = f.geometryChanges(GeometryFields(), report) # a dry run: nothing is written echo "would change ", plan # traces per kind, values with no slot for d in report: echo " ", d.msg +doAssert plan.traces == 0 and plan.unwritable > 0 var ed = openEditor(move f) # takes over the open file and its joined keys let done = ed.writeGeometry() ed.close() -doAssert done.receiverStations == 1 -var check = openSeismic("field.sgd", keys = true) -doAssert check.traceKeys(0).rcvPoint == 101 -doAssert rollbackJournal("field.sgd", "field.sgd.edits.nif").len == 0 # undo, if wanted +doAssert done == plan # rollbackJournal undoes a write, if wanted # 2. Into a copy in the same standard. -let copied = copyWithGeometry("record.sgd", "record_geo.sgd", g, - overwrite = true) -doAssert not copied.log.hasErrors and copied.changes.receiverStations == 1 +let copied = copyWithGeometry("record.sgd", "record_geo.sgd", g) +doAssert not copied.log.hasErrors and copied.changes.traces == 0 # 3. Converted to SEG-Y, which has coordinate fields for everything. -let convLog = convertWithGeometry("record.sgd", "record_geo.sgy", g, - overwrite = true) +let convLog = convertWithGeometry("record.sgd", "record_geo.sgy", g) doAssert not convLog.isFatal var y = openSeismic("record_geo.sgy", keys = true) doAssert y.traceKeys(0).rcvX == 500100 diff --git a/docs/spec/sps.md b/docs/spec/sps.md index 202f0bc..01131c1 100644 --- a/docs/spec/sps.md +++ b/docs/spec/sps.md @@ -543,6 +543,23 @@ Not verifiable from the examples: the exact columns of 1.0 static/depth/datum/up 5. **Units and CRS** (H12-H259, H20/H201) have no slot in `GeometryTable`. Keep them on the SPS header model. 6. **Elevation vs depth.** SPS gives the surface `elevation`, `pointDepth` (below surface) and `waterDepth`. `Coordinates.z` gets the surface elevation. A source depth is `pointDepth`, which the caller applies. +### 10.1b Join order (as implemented) + +Field record numbers and channel numbers in the data often differ from the X file's: a +recorder numbers its traces its own way (all live channels, not the X file's patch), and a +reshoot or renumbered post-plot moves records. Seen on a 3D land survey: a SEG-D file whose +FFID was another shot's X record, and whose trace numbers matched the X channels for 23 of +1127 traces, while its header stations matched the R/S/X stations exactly. So +`survey/join` matches **by stations**: header stations are never replaced; a missing +receiver comes from the X records of the header's source station, a missing source from the +O gathers of the header's receiver station, and only a trace with no station uses +(record, channel). See the module documentation. + +**O records** (receiver domain): not in either SPS revision, but in use (`.ops` files): the X +layout with the source and receiver fields swapped, one gather per receiver station, its +traces numbered across the sources that recorded it. Read with record letter `O`, or `X` in +an `.ops` file; `receiverRelations` builds them from X. + ### 10.2 Reader/writer guidance - Parse each line with a per-revision column table (`[(first, last, fmt)]`). Pad short lines to 80. Lines longer than 80: keep the extra in a `tail` field and warn. diff --git a/src/parseSEG/cli.nim b/src/parseSEG/cli.nim index 9e14100..321f544 100644 --- a/src/parseSEG/cli.nim +++ b/src/parseSEG/cli.nim @@ -44,7 +44,8 @@ Survey files (SPS 1.0/2.1, SEG P1 1983; docs/spec/sps.md, p1_p2_p3.md) --line-map=NAME=N[,...] text line names -> line numbers --trailing-digits unmapped names use their trailing digits --multi-component SPS channel increment k: k components/station - --overwrite replace stations and coordinates already known + --overwrite replace coordinates already known (header + stations are never replaced) survey-apply SURVEY [-o out | --in-place [--journal=J]] [--values=LIST] write the joined geometry into the headers: into a copy (-o; multiplexed SEG-D is written @@ -264,9 +265,11 @@ proc joinSurvey(a: Args, f: var SeismicFile): DiagLog = ## Joins the SURVEY options' geometry to `f` (opened with keys) and prints ## what was joined. let g = loadSurvey(a, result) - let n = f.index.applySurvey(g, a.has("overwrite")) + var stats: JoinStats + discard f.index.applySurvey(g, stats, a.has("overwrite")) echo "survey: ", g.table.sources.len, " sources, ", g.table.receivers.len, - " receivers, ", g.relationCount, " relation runs; ", n, " traces updated" + " receivers, ", g.relationCount, " X + ", g.relationCount(rdReceiver), " O relation runs" + echo "joined: ", stats proc info(a: Args): int = a.onlyOptions(@["keys"] & @SurveyOptions) @@ -922,9 +925,11 @@ proc surveyApply(a: Args): int = var log: DiagLog let g = loadSurvey(a, log) let overwrite = a.has("overwrite") - let n = f.index.applySurvey(g, overwrite) + var stats: JoinStats + discard f.index.applySurvey(g, stats, overwrite) echo "survey: ", g.sourceCount, " sources, ", g.receiverCount, " receivers, ", - g.relationCount, " relation runs; ", n, " traces joined" + g.relationCount, " X + ", g.relationCount(rdReceiver), " O relation runs" + echo "joined: ", stats var c: GeometryChanges if output.len > 0: f.close() diff --git a/src/parseSEG/survey/join.nim b/src/parseSEG/survey/join.nim index 0864775..6ca8f1f 100644 --- a/src/parseSEG/survey/join.nim +++ b/src/parseSEG/survey/join.nim @@ -1,25 +1,46 @@ ## Joins survey geometry (SPS, SEG P1) to seismic traces. ## -## A `SurveyGeometry` collects station coordinates (a `GeometryTable`) and, -## from SPS X records, which stations every (field record, channel) was -## recorded from. `applySurvey` then fills each trace's `TraceKeys`: +## A `SurveyGeometry` collects station coordinates (a `GeometryTable`) and +## the relations between field records and stations: +## - **X records** (source domain, `rdSource`): channels of a source's field +## record → receivers; +## - **O records** (receiver domain, `rdReceiver`, `.ops` files): the X layout +## with the roles swapped, traces of a receiver's gather → sources. See +## `receiverRelations` to build them from X. ## -## 1. When the trace's field record and channel match an X record, its -## source and receiver stations (line, point, index) are set from the -## relation, where unknown (or always, with `overwrite`). This is how -## SEG-D 2.x traces without station numbers in their headers get them. -## 2. The stations' coordinates come from the table (`keys.applyGeometry`), -## then the offset and midpoint. +## **Traces are matched by their stations** (line, point, index); the record +## number is a fallback for traces whose headers carry no station. +## `joinStations` fills the stations a header lacks, and never replaces one it +## has (header stations are the truth: field record numbers and channel +## numbers often differ from the relation file's). Per trace: ## -## Records are matched on `TraceKeys.record` (SEG-D file number, SEG-Y field -## record) and channels on `TraceKeys.channel`. Relations are stored as runs -## (one per X record), never expanded per channel. +## 1. source and receiver station in the header: nothing to fill; +## 2. source station only: the X records **of that source station** give the +## receiver of the trace's channel (the record equal to the trace's own +## record number when there are several; else they must agree); +## 3. receiver station only: the O gathers of that receiver station give the +## source likewise; without O records, the X relation of the trace's +## record and channel, only if its receiver is the header's; +## 4. no station: the relation of the trace's record and channel (X, then O). +## +## A relation is used only when it agrees with every station the header has. +## A station is known when both its line and its point are; a header with +## only one of them (e.g. a SEG-Y shotpoint without a line) has it filled +## whole. +## `applySurvey` then takes coordinates from the table by station +## (`keys.applyGeometry`), then the offset and midpoint. Auxiliary traces are +## left alone. The channel → station step (2-4) assumes the recorder numbers +## channels as the relation file does; when headers carry both stations, +## channel numbers are never used. +## +## Relations are stored as runs (one per X/O record), never expanded per +## channel. ## ## Line names that are not numbers (SPS 1.0, P1) go through the geometry's ## `lineNames` (a caller-supplied name → number map, see `mapLine`), then ## (with `trailingDigits`) the trailing digits of the name (`CLT4960` → 4960). -import std/[options, tables, strutils, math, sets] +import std/[options, tables, strutils, math, sets, algorithm, os] import ../common/diagnostics import ../model/[keys, traceindex] import ./[textrec, cards, sps, p1] @@ -28,18 +49,29 @@ import ../private/checks export keys type + RelationDomain* = enum + ## The domain of a relation record. + rdSource = "source" ## X: a source's field record; channels → receivers. + rdReceiver = "receiver" ## O: a receiver's gather; traces → sources. + RelationRun = object fieldTape: string fromChannel, toChannel: int64 step: int64 - src: StationKey - rcvLine, fromRcv, toRcv: float64 - rcvIndex: int16 + fixed: StationKey # X: the source; O: the receiver + line, first, last: float64 # X: receiver line and points; O: source line and points + index: int16 # point index of `line`'s stations + + StationRecord = object + index: int16 + record: int64 SurveyGeometry* = object ## Station coordinates and channel relations from survey files. table*: GeometryTable - runs: Table[int64, seq[RelationRun]] + runs: array[RelationDomain, Table[int64, seq[RelationRun]]] + byStation: array[RelationDomain, Table[StationKey, seq[StationRecord]]] + # (line, point, index 0) of a run's fixed station → its records multiComponent*: bool ## X records with a channel increment k > 1 cover k consecutive ## channels per station (components 1..k, SPS 1.0 §"X" note); else @@ -60,17 +92,39 @@ type stations* = true ## Source and receiver line, point and point index. StationRelation* = object - ## The stations one channel of a field record was recorded from. + ## The stations one channel of a field record (X) or one trace of a + ## receiver gather (O) was recorded from. source*, receiver*: StationKey component*: int ## 1-based within the station (multi-component), else 0. + JoinMatch* = enum + ## How `joinStations` matched a trace (see the module documentation). + jmStations = "by stations" ## Both stations were in the header. + jmSource = "via the source" ## Receiver from the X records of the source station. + jmReceiver = "via the receiver" ## Source from the O gathers of the receiver station. + jmRecord = "by record number" ## Record and channel (no station, or checked by one). + jmAmbiguous = "ambiguous" ## The station's records disagree on the channel. + jmConflict = "conflicting" ## The record's relation contradicts a header station. + jmUnmatched = "unmatched" ## No relation found. + jmAuxiliary = "auxiliary" ## Not seismic: left alone. + + JoinStats* = object + ## What `applySurvey` did to a file's traces. + traces*: int ## Traces seen. + changed*: int ## Traces whose keys changed. + matches*: array[JoinMatch, int] ## Traces by how they matched. + noSourceXy*, noReceiverXy*: int + ## Seismic traces left without source / receiver coordinates (a + ## station missing from the survey, or no station at all). + func log*(g: SurveyGeometry): lent DiagLog {.inline.} = ## Problems met while adding survey files. g.log -func relationCount*(g: SurveyGeometry): int = - ## The number of X runs loaded. - for runs in g.runs.values: result += runs.len + +func relationCount*(g: SurveyGeometry, domain = rdSource): int = + ## The number of relation runs (X records, or O records) loaded. + for runs in g.runs[domain].values: result += runs.len func initSurveyGeometry*(): SurveyGeometry = ## An empty geometry: no stations, no relations. @@ -130,23 +184,71 @@ proc lineOf(g: var SurveyGeometry, unresolved: var HashSet[string], file, name: g.note(sevWarning, file, "line '" & name & "' is not a number and has no mapping; " & "its records are skipped (use a line name map)", rec) +func `==`(a, b: StationRecord): bool = a.index == b.index and a.record == b.record + +func anyIndex(k: StationKey): StationKey {.inline.} = + StationKey(line: k.line, point: k.point, index: 0) + +func sameStation*(a, b: StationKey): bool = + ## Same line and point (to 1/100), and the same point index unless either + ## is 0 (unknown: SEG-D headers often leave it 0 where survey files write 1). + let (x, y) = (stationKey(a.line, a.point, a.index), stationKey(b.line, b.point, b.index)) + x.line == y.line and x.point == y.point and + (x.index == 0 or y.index == 0 or x.index == y.index) + +func lastChannel(run: RelationRun, multi: bool): int64 {.inline.} = + if multi: run.toChannel + run.step - 1 else: run.toChannel + +func covers(run: RelationRun, channel: int64, multi: bool): bool {.inline.} = + channel >= run.fromChannel and channel <= run.lastChannel(multi) and + (multi or (channel - run.fromChannel) mod run.step == 0) + +func at(run: RelationRun, channel: int64, multi: bool, domain: RelationDomain): StationRelation = + ## The relation of a channel `run` covers. + let d = channel - run.fromChannel + let i = d div run.step + let n = (run.toChannel - run.fromChannel) div run.step + let p = if n == 0: run.first + else: run.first + float64(i) * (run.last - run.first) / float64(n) + let moving = stationKey(run.line, p, run.index) + result.component = if multi and run.step > 1: int(d mod run.step) + 1 else: 0 + if domain == rdSource: (result.source, result.receiver) = (run.fixed, moving) + else: (result.source, result.receiver) = (moving, run.fixed) + +func sameRelation(a, b: StationRelation): bool = + a.source == b.source and a.receiver == b.receiver and a.component == b.component + func overlaps(runs: seq[RelationRun], run: RelationRun): bool = for o in runs: if run.fromChannel <= o.toChannel and o.fromChannel <= run.toChannel: return true -proc addRelation(g: var SurveyGeometry, file: string, i: int, r: SpsRelation, sl, rl: float64) = +proc addRelation(g: var SurveyGeometry, domain: RelationDomain, file: string, i: int, + r: SpsRelation, fixedLine, line: float64) = + # In both domains the record's source columns hold the fixed station and + # its receiver columns the run of moving stations. let run = RelationRun(fieldTape: r.fieldTape, fromChannel: r.fromChannel, toChannel: r.toChannel, step: r.channelIncrement, - src: stationKey(sl, r.srcPoint, idx(r.srcIndex)), rcvLine: rl, - fromRcv: r.fromReceiver, toRcv: r.toReceiver, rcvIndex: idx(r.rcvIndex)) - if r.fieldRecord in g.runs and g.runs[r.fieldRecord].overlaps(run): - g.note(sevWarning, file, "field record " & $r.fieldRecord & ": channels " & - $run.fromChannel & ".." & $run.toChannel & " overlap an earlier X record " & - "(different tapes or line sequences?); the earlier one wins", i) - g.runs.mgetOrPut(r.fieldRecord, @[]).add run - -proc addSps*(g: var SurveyGeometry, f: SpsFile) = - ## Adds the stations (R and S records) and relations (X records) of `f`. + fixed: stationKey(fixedLine, r.srcPoint, idx(r.srcIndex)), line: line, + first: r.fromReceiver, last: r.toReceiver, index: idx(r.rcvIndex)) + if r.fieldRecord in g.runs[domain] and g.runs[domain][r.fieldRecord].overlaps(run): + g.note(sevWarning, file, (if domain == rdSource: "field record " else: "gather ") & + $r.fieldRecord & ": channels " & $run.fromChannel & ".." & $run.toChannel & + " overlap an earlier record (different tapes or line sequences?); the earlier one wins", i) + g.runs[domain].mgetOrPut(r.fieldRecord, @[]).add run + let entry = StationRecord(index: run.fixed.index, record: r.fieldRecord) + let key = run.fixed.anyIndex + if key notin g.byStation[domain] or entry notin g.byStation[domain][key]: + g.byStation[domain].mgetOrPut(key, @[]).add entry + +func relationDomainOf*(path: string): RelationDomain = + ## The domain of the X-type records of a file, by its name: `rdReceiver` + ## for an `.ops` file (any case), else `rdSource`. + if path.splitFile.ext.toLowerAscii == ".ops": rdReceiver else: rdSource + +proc addSps*(g: var SurveyGeometry, f: SpsFile, xDomain = rdSource) = + ## Adds the stations (R and S records) and relations of `f`. `O` records + ## are receiver-domain relations; `X` records are of `xDomain` (pass + ## `rdReceiver` for an O file written with `X` records, e.g. `.ops`). ## Grid coordinates are taken as they are (units of H20); z is the surface ## elevation. Records whose line or point cannot be made numeric are ## skipped, with one warning per line name. @@ -164,19 +266,21 @@ proc addSps*(g: var SurveyGeometry, f: SpsFile) = else: g.table.addReceiver(l.get, p.point, idx(p.index), p.easting, p.northing, p.elevation) for i, r in f.relations: + let domain = if r.recordType == 'O': rdReceiver else: xDomain if g.tape.len > 0 and r.fieldTape != g.tape: discard elif r.fieldRecord == MissingInt or r.fromChannel == MissingInt or r.toChannel == MissingInt or r.channelIncrement < 1 or r.toChannel < r.fromChannel: - g.note(sevWarning, f.name, "X record without a usable record or channel range; skipped", i) + g.note(sevWarning, f.name, r.recordType & " record without a usable record or channel " & + "range; skipped", i) else: let sl = g.lineOf(unresolved, f.name, r.srcLineName, r.srcLine, i) let rl = if sl.isSome: g.lineOf(unresolved, f.name, r.rcvLineName, r.rcvLine, i) else: none(float64) if rl.isNone: discard elif r.srcPoint.isNaN or r.fromReceiver.isNaN or r.toReceiver.isNaN: - g.note(sevWarning, f.name, "X record with a non-numeric point; skipped", i) + g.note(sevWarning, f.name, r.recordType & " record with a non-numeric point; skipped", i) else: - g.addRelation(f.name, i, r, sl.get, rl.get) + g.addRelation(domain, f.name, i, r, sl.get, rl.get) proc addP1*(g: var SurveyGeometry, f: P1File, elevation = false) = ## Adds the shotpoints of the P1 logical files in `f` as sources. The @@ -207,12 +311,14 @@ proc addLog(g: var SurveyGeometry, log: DiagLog) = if log.isFatal: raise newException(IOError, log[log.len - 1].msg) g.log.add log -proc addSps*(g: var SurveyGeometry, path: string, opts = SpsOptions()) = +proc addSps*(g: var SurveyGeometry, path: string, opts = SpsOptions(), + xDomain = relationDomainOf(path)) = ## Reads the SPS file (Rev 1.0 or 2.1) at `path` and adds it; its - ## diagnostics go to `g.log`. Raises `IOError` if it cannot be read. + ## diagnostics go to `g.log`. `X` records of an `.ops` file are receiver + ## domain (O) by default. Raises `IOError` if it cannot be read. let f = readSps(path, opts) g.addLog(f.log) - g.addSps(f) + g.addSps(f, xDomain) proc addP1*(g: var SurveyGeometry, path: string, opts = P1Options(), elevation = false) = ## Reads the SEG P1 file at `path` (`readP1` with `opts`: the implicit @@ -226,43 +332,150 @@ proc addP1*(g: var SurveyGeometry, path: string, opts = P1Options(), elevation = g.addLog(f.log) g.addP1(f, elevation) -proc relation*(g: SurveyGeometry, record, channel: int64): Option[StationRelation] = - ## The stations channel `channel` of field record `record` was recorded - ## from, per the X records; none if no run covers it. - if record notin g.runs: return none(StationRelation) - for run in g.runs[record]: - let last = if g.multiComponent: run.toChannel + run.step - 1 else: run.toChannel - if channel >= run.fromChannel and channel <= last: - let d = channel - run.fromChannel - if g.multiComponent or d mod run.step == 0: - let i = d div run.step - let n = (run.toChannel - run.fromChannel) div run.step - let rp = if n == 0: run.fromRcv - else: run.fromRcv + float64(i) * (run.toRcv - run.fromRcv) / float64(n) - return some(StationRelation(source: run.src, - receiver: stationKey(run.rcvLine, rp, run.rcvIndex), - component: if g.multiComponent and run.step > 1: int(d mod run.step) + 1 else: 0)) +proc relation*(g: SurveyGeometry, record, channel: int64, + domain = rdSource): Option[StationRelation] = + ## The stations channel `channel` of field record `record` (X, `rdSource`) + ## or trace `channel` of gather `record` (O, `rdReceiver`) was recorded + ## from; none if no run covers it. + if record notin g.runs[domain]: return none(StationRelation) + for run in g.runs[domain][record]: + if run.covers(channel, g.multiComponent): + return some(run.at(channel, g.multiComponent, domain)) none(StationRelation) -proc applySurvey*(k: var TraceKeys, g: SurveyGeometry, overwrite = false): bool = - ## Fills `k` from the survey: stations from the X relations, then - ## coordinates, offset and midpoint. True if anything changed. - result = false - let found = if k.record.isNaN or k.channel.isNaN: none(StationRelation) - else: g.relation(int64(k.record), int64(k.channel)) - if found.isSome: - let rel = found.get - if overwrite or k.srcLine.isNaN or k.srcPoint.isNaN: - (k.srcLine, k.srcPoint, k.srcIndex) = (rel.source.line, rel.source.point, rel.source.index) - result = true - if overwrite or k.rcvLine.isNaN or k.rcvPoint.isNaN: - (k.rcvLine, k.rcvPoint, k.rcvIndex) = - (rel.receiver.line, rel.receiver.point, rel.receiver.index) - result = true - if rel.component > 0 and (overwrite or k.component == 0): - k.component = uint8(min(rel.component, 255)) - result = true - if k.applyGeometry(g.table, overwrite): result = true +iterator records*(g: SurveyGeometry, domain = rdSource): int64 = + ## The field records (X) or gathers (O) that have relations, ascending. + var all = newSeqOfCap[int64](g.runs[domain].len) + for r in g.runs[domain].keys: all.add r + all.sort + for r in all: yield r + +iterator relations*(g: SurveyGeometry, record: int64, + domain = rdSource): tuple[channel: int64, relation: StationRelation] = + ## Every channel of `record` with its stations, run by run in file order; + ## a channel two runs cover is yielded once, for the earlier run (as + ## `relation` does). + if record in g.runs[domain]: + let runs = g.runs[domain][record] + let multi = g.multiComponent + for j, run in runs: + var c = run.fromChannel + while c <= run.lastChannel(multi): + if run.covers(c, multi): + var earlier = false + for o in 0 ..< j: + if runs[o].covers(c, multi): + earlier = true + break + if not earlier: yield (c, run.at(c, multi, domain)) + inc c + +proc expectedChannels*(g: SurveyGeometry, record: int64, domain = rdSource): int = + ## The number of channels (X) or traces (O) `record` has relations for. + for _ in g.relations(record, domain): inc result + +proc recordsAt*(g: SurveyGeometry, station: StationKey, domain = rdSource): seq[int64] = + ## The field records of a source station (X, `rdSource`) or the gathers of + ## a receiver station (O, `rdReceiver`), ascending: several for repeated or + ## reshot points. Index 0 on either side matches any index. + let key = stationKey(station.line, station.point, 0) + if key in g.byStation[domain]: + for e in g.byStation[domain][key]: + if station.index == 0 or e.index == 0 or e.index == station.index: + if e.record notin result: result.add e.record + result.sort + +proc channelOf*(g: SurveyGeometry, record: int64, station: StationKey, + domain = rdSource): Option[int64] = + ## The channel of field record `record` recorded at receiver `station` + ## (X), or the trace of gather `record` from source `station` (O); the + ## first component's channel for multi-component runs. None if the record + ## does not list the station. + if record notin g.runs[domain]: return none(int64) + let st = stationKey(station.line, station.point, station.index) + for run in g.runs[domain][record]: + if run.line == st.line and (st.index == 0 or run.index == 0 or st.index == run.index): + let n = (run.toChannel - run.fromChannel) div run.step + let i = if n == 0: 0'i64 + elif run.last == run.first: -1'i64 + else: int64(round((st.point - run.first) / (run.last - run.first) * float64(n))) + if i in 0'i64 .. n: + let p = if n == 0: run.first + else: run.first + float64(i) * (run.last - run.first) / float64(n) + if stationKey(run.line, p, run.index).point == st.point: + let c = run.fromChannel + i * run.step + let won = g.relation(record, c, domain) # an earlier run may cover `c` + if won.isSome and sameRelation(won.get, run.at(c, g.multiComponent, domain)): + return some(c) + none(int64) + +type ViaStation = enum vsNone, vsFound, vsAmbiguous + +proc viaStation(g: SurveyGeometry, domain: RelationDomain, station: StationKey, + record: float64, channel: int64, rel: var StationRelation): ViaStation = + ## The relation of `channel` in the records of `station`: the record equal + ## to the trace's record number when it is one of them, else the one + ## relation all of them agree on. + let recs = g.recordsAt(station, domain) + if recs.len == 0: return vsNone + if not record.isNaN and int64(record) in recs: + let r = g.relation(int64(record), channel, domain) + if r.isNone: return vsNone + rel = r.get + return vsFound + result = vsNone + for rec in recs: + let r = g.relation(rec, channel, domain) + if r.isSome: + if result == vsNone: + rel = r.get + result = vsFound + elif not sameRelation(rel, r.get): return vsAmbiguous + +proc joinStations*(k: var TraceKeys, g: SurveyGeometry): JoinMatch = + ## Fills the source and receiver stations (line, point, index) and the + ## component a trace's header lacks, from the relations, by the rules in + ## the module documentation. Stations the header has are never changed. + ## Auxiliary traces are left alone. + if k.status == tsAuxiliary: return jmAuxiliary + let hasSrc = not (k.srcLine.isNaN or k.srcPoint.isNaN) + let hasRcv = not (k.rcvLine.isNaN or k.rcvPoint.isNaN) + if hasSrc and hasRcv: return jmStations + if k.channel.isNaN: return jmUnmatched + let channel = int64(k.channel) + var rel: StationRelation + if hasSrc: + case g.viaStation(rdSource, stationKey(k.srcLine, k.srcPoint, k.srcIndex), k.record, + channel, rel) + of vsNone: return jmUnmatched + of vsAmbiguous: return jmAmbiguous + of vsFound: result = jmSource + elif hasRcv: + let st = stationKey(k.rcvLine, k.rcvPoint, k.rcvIndex) + case g.viaStation(rdReceiver, st, k.record, channel, rel) + of vsAmbiguous: return jmAmbiguous + of vsFound: result = jmReceiver + of vsNone: + let r = if k.record.isNaN: none(StationRelation) + else: g.relation(int64(k.record), channel) + if r.isNone: return jmUnmatched + if not sameStation(r.get.receiver, st): return jmConflict + rel = r.get + result = jmRecord + else: + if k.record.isNaN: return jmUnmatched + var r = g.relation(int64(k.record), channel) + if r.isNone: r = g.relation(int64(k.record), channel, rdReceiver) + if r.isNone: return jmUnmatched + rel = r.get + result = jmRecord + if not hasSrc: + (k.srcLine, k.srcPoint, k.srcIndex) = (rel.source.line, rel.source.point, rel.source.index) + if not hasRcv: + (k.rcvLine, k.rcvPoint, k.rcvIndex) = + (rel.receiver.line, rel.receiver.point, rel.receiver.index) + if rel.component > 0 and k.component == 0: + k.component = uint8(min(rel.component, 255)) func coordinatesDiffer(a, b: TraceKeys): bool = func same(x, y: float64): bool = x == y or (x.isNaN and y.isNaN) @@ -276,27 +489,71 @@ func stationsDiffer(a, b: TraceKeys): bool = same(a.rcvLine, b.rcvLine) and same(a.rcvPoint, b.rcvPoint) and a.srcIndex == b.srcIndex and a.rcvIndex == b.rcvIndex and a.component == b.component) -proc applySurvey*(ix: var TraceIndex, g: SurveyGeometry, overwrite = false): int = - ## `applySurvey` on every trace; returns the traces changed. Raises - ## `ValueError` for an index built without `keys = true`. +proc join(k: var TraceKeys, g: SurveyGeometry, overwrite: bool, match: var JoinMatch): bool = + let before = k + match = k.joinStations(g) + result = stationsDiffer(before, k) + if match != jmAuxiliary and k.applyGeometry(g.table, overwrite): result = true + +proc applySurvey*(k: var TraceKeys, g: SurveyGeometry, overwrite = false): bool = + ## Fills `k` from the survey: the stations its header lacks + ## (`joinStations`), then coordinates by station, offset and midpoint. + ## `overwrite` replaces known coordinates, offset and midpoint; header + ## stations are never replaced. True if anything changed. + var match: JoinMatch + k.join(g, overwrite, match) + +proc applySurvey*(ix: var TraceIndex, g: SurveyGeometry, stats: var JoinStats, + overwrite = false): int = + ## `applySurvey` on every trace; returns the traces changed, and how they + ## matched in `stats`. Raises `ValueError` for an index built without + ## `keys = true`. ix.requireKeys() ix.joined = true + stats = JoinStats(traces: ix.keys.len) for i, k in ix.keys.mpairs: let before = k - if k.applySurvey(g, overwrite): + var match: JoinMatch + if k.join(g, overwrite, match): inc result var groups: set[JoinGroup] if coordinatesDiffer(before, k): groups.incl jgCoordinates if stationsDiffer(before, k): groups.incl jgStations ix.markJoined(i, groups) + inc stats.matches[match] + if match != jmAuxiliary: + if not k.hasSourceXY: inc stats.noSourceXy + if not k.hasReceiverXY: inc stats.noReceiverXy + stats.changed = result if result > 0: inc ix.keysRevision +proc applySurvey*(ix: var TraceIndex, g: SurveyGeometry, overwrite = false): int = + ## `applySurvey` on every trace; returns the traces changed. Raises + ## `ValueError` for an index built without `keys = true`. + var stats: JoinStats + ix.applySurvey(g, stats, overwrite) + +proc applySurvey*(f: var SeismicFile, g: SurveyGeometry, stats: var JoinStats, + overwrite = false): int = + ## Fills the trace keys from the survey, in memory only: the file does not + ## change. Returns the traces changed, and how they matched in `stats`. + ## Raises `ValueError` for a file opened without `keys = true`. + f.index.applySurvey(g, stats, overwrite) + proc applySurvey*(f: var SeismicFile, g: SurveyGeometry, overwrite = false): int = ## Fills the trace keys from the survey, in memory only: the file does not ## change. Returns the traces changed. Raises `ValueError` for a file opened ## without `keys = true`. f.index.applySurvey(g, overwrite) +func `$`*(s: JoinStats): string = + ## One line: traces, changed, the matches that occurred, missing coordinates. + result = $s.traces & " traces, " & $s.changed & " changed" + for m in JoinMatch: + if s.matches[m] > 0: result.add "; " & $m & " " & $s.matches[m] + if s.noSourceXy > 0: result.add "; " & $s.noSourceXy & " without source X/Y" + if s.noReceiverXy > 0: result.add "; " & $s.noReceiverXy & " without receiver X/Y" + proc applySurvey*(k: var TraceKeys, g: SurveyGeometry, fields: GeometryFields, overwrite = false): bool = ## `applySurvey`, keeping only the joined values `fields` selects (the @@ -322,3 +579,104 @@ proc applySurvey*(k: var TraceKeys, g: SurveyGeometry, fields: GeometryFields, same(k.srcPoint, before.srcPoint) and same(k.rcvLine, before.rcvLine) and same(k.rcvPoint, before.rcvPoint) and k.srcIndex == before.srcIndex and k.rcvIndex == before.rcvIndex and k.component == before.component) + +# ------------------------------------------------------- receiver domain --- + +type SourceAt = object + line, point: float64 + index: int16 + +func `<`(a, b: SourceAt): bool = + if a.line != b.line: a.line < b.line + elif a.point != b.point: a.point < b.point + else: a.index < b.index + +func stationLess(a, b: StationKey): int = + if a.line != b.line: cmp(a.line, b.line) + elif a.point != b.point: cmp(a.point, b.point) + else: cmp(a.index, b.index) + +proc receiverRelations*(g: SurveyGeometry, firstGather = 1'i64, + recordType = 'X'): seq[SpsRelation] = + ## The receiver-domain (O) records of the X relations: one gather per + ## receiver station, numbered from `firstGather` in receiver line, point + ## order; its traces 1..n are the sources that recorded it, by source + ## line and point, as runs along source lines; a repeated source point + ## (one X record per repeat) is one more trace, after the first ones, so + ## each X pair appears once. Components are not + ## represented (one trace per source and receiver station). Write them with + ## `writeSps` (records start with `recordType`, `X` as in `.ops` files, or + ## `O`); `addSps(..., xDomain = rdReceiver)` reads them back. Works one + ## receiver line at a time, so memory follows the largest line, not the + ## survey. Raises `ValueError` for a `recordType` other than `X` or `O`. + if recordType notin {'X', 'O'}: + raise newException(ValueError, "record type must be 'X' or 'O'") + let multi = g.multiComponent + # X runs by receiver line: (record, run). + var byLine: Table[float64, seq[(int64, int)]] + for rec, runs in g.runs[rdSource]: + for j, run in runs: byLine.mgetOrPut(run.line, @[]).add (rec, j) + var lines = newSeqOfCap[float64](byLine.len) + for l in byLine.keys: lines.add l + lines.sort + var gather = firstGather + for line in lines: + var at: Table[StationKey, seq[SourceAt]] + for (rec, j) in byLine[line]: + let runs = g.runs[rdSource][rec] + let run = runs[j] + var c = run.fromChannel + while c <= run.lastChannel(multi): + if run.covers(c, multi): + var earlier = false + for o in 0 ..< j: + if runs[o].covers(c, multi): + earlier = true + break + let rel = run.at(c, multi, rdSource) + if not earlier and rel.component <= 1: + at.mgetOrPut(rel.receiver, @[]).add SourceAt(line: rel.source.line, + point: rel.source.point, index: rel.source.index) + inc c + var stations = newSeqOfCap[StationKey](at.len) + for k in at.keys: stations.add k + stations.sort(stationLess) + for st in stations: + var srcs = at[st] + srcs.sort(proc (a, b: SourceAt): int = (if a < b: -1 elif b < a: 1 else: 0)) + # Each source's first occurrence, then its repeats (pass by pass), so a + # repeated shot does not cut its line's run in two. + var passes: seq[SourceAt] + var rest = srcs + while rest.len > 0: + var later: seq[SourceAt] + for j, a in rest: + if j > 0 and not (rest[j - 1] < a): later.add a + else: passes.add a + rest = later + srcs = passes + var trace = 1'i64 + var i = 0 + while i < srcs.len: + var j = i + if i + 1 < srcs.len and srcs[i + 1].line == srcs[i].line and + srcs[i + 1].index == srcs[i].index and srcs[i + 1].point > srcs[i].point: + let step = srcs[i + 1].point - srcs[i].point + j = i + 1 + while j + 1 < srcs.len and srcs[j + 1].line == srcs[i].line and + srcs[j + 1].index == srcs[i].index and + abs(srcs[j + 1].point - srcs[j].point - step) < 0.005: + inc j + var r = initSpsRelation() + r.recordType = recordType + r.fieldRecord = gather + (r.srcLine, r.srcPoint, r.srcIndex) = (st.line, st.point, int(st.index)) + r.fromChannel = trace + r.toChannel = trace + int64(j - i) + r.channelIncrement = 1 + (r.rcvLine, r.fromReceiver, r.toReceiver, r.rcvIndex) = + (srcs[i].line, srcs[i].point, srcs[j].point, int(srcs[i].index)) + result.add r + trace += int64(j - i + 1) + i = j + 1 + inc gather diff --git a/src/parseSEG/survey/sps.nim b/src/parseSEG/survey/sps.nim index 98d78d9..d5818a1 100644 --- a/src/parseSEG/survey/sps.nim +++ b/src/parseSEG/survey/sps.nim @@ -3,8 +3,8 @@ ## ## An SPS survey is a set of files that share one block of header (H) ## records: receiver points (R), source points (S), relations (X: field -## record and channel → source and receiver stations) and, optionally, -## comments (C). `readSps` reads any of them (or a file mixing them) into one +## record and channel → source and receiver stations; O: the same layout in +## the receiver domain, see `SpsRelation`) and, optionally, comments (C). `readSps` reads any of them (or a file mixing them) into one ## `SpsFile`; `writeSps` writes the records back, column-exact. ## ## Reading never stops on a bad record: each problem is a diagnostic in @@ -70,6 +70,14 @@ type ## `channelIncrement`) of field record `fieldRecord` were recorded from ## receivers `fromReceiver..toReceiver` on `rcvLine`, shot at ## (`srcLine`, `srcPoint`, `srcIndex`). + ## + ## An O record (receiver domain, `recordType` `'O'`) has the same layout + ## with the roles swapped: traces `fromChannel..toChannel` of gather + ## `fieldRecord` were recorded at the receiver in the source fields, from + ## the sources `fromReceiver..toReceiver` on `rcvLine` (see + ## `survey/join`). Files named `.ops` often hold O relations as `X` + ## records. + recordType*: char = 'X' ## `X`, or `O` for a receiver-domain record. fieldTape*: string fieldRecord*: int64 recordIncrement*: int @@ -130,7 +138,7 @@ func initSpsPoint*(kind = pkSource): SpsPoint = func initSpsRelation*(): SpsRelation = ## A relation with every value missing (increments and indices 1). - SpsRelation(fieldRecord: MissingInt, recordIncrement: 1, instrumentCode: '1', srcLine: NaN, + SpsRelation(recordType: 'X', fieldRecord: MissingInt, recordIncrement: 1, instrumentCode: '1', srcLine: NaN, srcPoint: NaN, srcIndex: 1, fromChannel: MissingInt, toChannel: MissingInt, channelIncrement: 1, rcvLine: NaN, fromReceiver: NaN, toReceiver: NaN, rcvIndex: 1) @@ -251,6 +259,7 @@ proc parsePoint(x: Ctx, log: var DiagLog, line: string, rev: SpsRevision): SpsPo proc parseRelation(x: Ctx, log: var DiagLog, line: string, rev: SpsRevision): SpsRelation = let cols = relCols[rev] result = initSpsRelation() + result.recordType = line[0] result.fieldTape = line.c(cols.tape).strip result.fieldRecord = x.intField(log, line.c(cols.rec), "X.fieldRecord", MissingInt) result.recordIncrement = x.smallField(log, line.c(cols.recInc), "X.recordIncrement", 1) @@ -316,7 +325,7 @@ func layoutScore(line: string): int = not line.col(12, 21).isBlank: inc result else: dec result - of 'X': + of 'X', 'O': if parseFortranFloat(line.col(18, 27)).isSome and parseFortranFloat(line.col(28, 37)).isSome: inc result else: dec result @@ -341,7 +350,7 @@ proc readSps*(s: Stream, name = "", opts = SpsOptions()): SpsFile = rev = rv result.revisionDeclared = true break - elif r.text[0] in {'R', 'S', 'X'}: + elif r.text[0] in {'R', 'S', 'X', 'O'}: score += layoutScore(r.text) inc seen if seen >= 50: break @@ -374,7 +383,7 @@ proc readSps*(s: Stream, name = "", opts = SpsOptions()): SpsFile = of 'R', 'S': result.order.add SpsEntry(kind: sePoint, at: result.points.len) result.points.add x.parsePoint(result.log, line, rev) - of 'X': + of 'X', 'O': result.order.add SpsEntry(kind: seRelation, at: result.relations.len) result.relations.add x.parseRelation(result.log, line, rev) of 'C': @@ -408,7 +417,7 @@ proc looksLikeSps*(s: Stream, maxRecords = 200): bool = let (e, n) = (parseFortranFloat(t.col(47, 55), 1), parseFortranFloat(t.col(56, 65), 1)) if e.isSome and n.isSome and not (e.get.isNaN or n.get.isNaN): inc sps - of 'X': + of 'X', 'O': if t.col(39, 48).strip.len > 0 and t.col(39, 48).strip.allCharsInSet(Digits + {' '}): inc sps of ' ': @@ -516,8 +525,10 @@ proc encodeSps*(r: SpsRelation, rev: SpsRevision): string = ## when a value does not fit its field, `ValueError` for `spsAuto`. checkTarget(rev) let cols = relCols[rev] + if r.recordType notin {'X', 'O'}: + raiseFormat("X.recordType: '" & $r.recordType & "' is not X or O") result = repeat(' ', 80) - result[0] = 'X' + result[0] = r.recordType result.putCol(cols.tape.a, alignLeft(fit(r.fieldTape & (if r.fieldTape.len == 0: " " else: ""), cols.tape.len, "X.fieldTape"), cols.tape.len)) result.putCol(cols.rec.a, fit(formatInt(r.fieldRecord, cols.rec.len), cols.rec.len, @@ -564,7 +575,7 @@ proc pointLine(p: SpsPoint, rev, fileRev: SpsRevision, strict: bool): string = encodeSps(p, rev) proc relationLine(r: SpsRelation, rev, fileRev: SpsRevision, strict: bool): string = - if r.raw.len > 0 and rev == fileRev and r.raw[0] == 'X': + if r.raw.len > 0 and rev == fileRev and r.raw[0] in {'X', 'O'}: var scratch: DiagLog if sameValues(quietCtx(strict).parseRelation(scratch, r.raw, rev), r): return r.raw encodeSps(r, rev) diff --git a/src/parseSEG/survey/writeback.nim b/src/parseSEG/survey/writeback.nim index 7cd84d6..6f64312 100644 --- a/src/parseSEG/survey/writeback.nim +++ b/src/parseSEG/survey/writeback.nim @@ -657,7 +657,9 @@ template eachChange(f: SeismicFile, keys: openArray[TraceKeys], fields: Geometry if e.offset >= 0 and e.offset + n <= f.mappedSize and n > 0: buf.setLen n copyMem(addr buf[0], addr mem[e.offset], n) - let changed = p.planTrace(f.index, i, keys[i], buf) + var changed = p.planTrace(f.index, i, keys[i], buf) + if changed.card > 0 and equalMem(addr buf[0], addr mem[e.offset], n): + changed = {} # the values round to the bytes already there (e.g. 99.5 → 100) p.count(changed) if changed.card > 0: body diff --git a/tests/tcli.nim b/tests/tcli.nim index 222ea65..1c6f1f3 100644 --- a/tests/tcli.nim +++ b/tests/tcli.nim @@ -262,12 +262,13 @@ sampleBlock vendormaps: # ------------------------------------------- Phase 7 review regressions -- sampleBlock survey: + # 63.sgd's header stations: source 13/99.5, receiver 13/104. let sps = tmp / "t.sps" writeFile(sps, [ "H00 SPS format version number SPS 2.1;", - "R 10.00 101.00 1G1 500100.0 6000000.0 12.5 18213250", - "S 20.00 205.00 1V6 500000.0 6000400.0 10.0 19001150", - "XT1 6311 20.00 205.001 1 11 10.00 101.00 101.001", + "R 13.00 104.00 1G1 500100.0 6000000.0 12.5 18213250", + "S 13.00 99.50 1V6 500000.0 6000400.0 10.0 19001150", + "XT1 6311 13.00 99.501 1 11 13.00 104.00 104.001", ].join("\r\n") & "\r\n") let s = ok("survey", sps) doAssert "SPS 2.1 (H00)" in s and "CR LF" in s and "1 H, 1 S, 1 R, 1 X, 0 C" in s @@ -277,11 +278,13 @@ sampleBlock survey: doAssert readFile(tmp / "t21.sps") == readFile(sps), "1.0 -> 2.1 restores the file" discard fails(2, "survey", sps, "-o", sps) discard fails(2, "survey", sps, "--rev=3") - let k0 = ok("keys", s63, "--sps=" & sps) - doAssert "1 relation runs; 0 traces updated" in k0, "the headers' stations are kept" - let k = ok("keys", s63, "--sps=" & sps, "--overwrite") - doAssert "1 traces updated" in k and "\t63\t1\t20\t205\t500000\t6000400\t10\t101\t500100" in k - doAssert "412.31" in k + for extra in ["", "--overwrite"]: + let k = if extra.len == 0: ok("keys", s63, "--sps=" & sps) + else: ok("keys", s63, "--sps=" & sps, extra) + doAssert "1 X + 0 O relation runs" in k, k + doAssert "joined: 1 traces, 1 changed; by stations 1" in k, k + doAssert "\t63\t1\t13\t99.50\t500000\t6000400\t13\t104\t500100" in k, k + doAssert "412.31" in k doAssert "relation runs" in ok("info", s63, "--sps", sps) discard fails(2, "keys", s63, "--p1=" & sps, "--decimals=9,1") discard fails(2, "keys", s63, "--line-map=oops") @@ -291,24 +294,31 @@ sampleBlock survey: discard fails(2, "keys", s63, "--overwrite") # needs --sps or --p1 sampleBlock survey_apply: - let sps = tmp / "t.sps" # from the survey block: record 63 -> stations 20/205, 10/101 - let path = fresh(s63) + # 63.sgd converted to SEG-Y: its keys have no stations (the layout's fields + # are not read into them), so the join fills them by record and channel, + # and SEG-Y has fields for the coordinates. + let sps = tmp / "t.sps" # from the survey block: record 63 -> stations 13/99.5, 13/104 + let path = tmp / "63cli.sgy" + discard ok("convert", s63, path) let original = readFile(path) - let dry = ok("survey-apply", path, "--sps=" & sps, "--overwrite") - doAssert "would change" in dry and "rcv station: 1 traces" in dry, dry - doAssert "receiver X/Y: SEG-D Rev 1/2 trace headers have no coordinate fields" in dry + let dry = ok("survey-apply", path, "--sps=" & sps) + doAssert "would change" in dry and "by record number 1" in dry, dry + doAssert "receiver X/Y:1 traces" in dry and "rcv station: 0 traces" in dry, dry # 13/104 already doAssert readFile(path) == original, "a dry run writes nothing" - let inPlace = ok("survey-apply", path, "--sps=" & sps, "--overwrite", "--in-place") + let inPlace = ok("survey-apply", path, "--sps=" & sps, "--in-place") doAssert "changed in place" in inPlace and "journal:" in inPlace - doAssert "\t10\t101\t" in ok("keys", path) + doAssert "\t500000\t6000400\t-\t-\t500100\t6000000\t412" in ok("keys", path) # stations: layout discard ok("journal", "rollback", path, path & ".edits.nif") doAssert readFile(path) == original - let copy = tmp / "63geo.sgd" - doAssert "changed in " & copy in ok("survey-apply", s63, "--sps=" & sps, "--overwrite", - "-o", copy, "--values=stations") - doAssert "\t10\t101\t" in ok("keys", copy) + let copy = tmp / "63geo.sgy" + let c = ok("survey-apply", path, "--sps=" & sps, "-o", copy, "--values=stations") + doAssert "changed in " & copy in c and "src station: 1 traces" in c, c + doAssert "500100" notin ok("keys", copy), "stations only" + # SEG-D 2.0 headers already carry the stations and have no coordinate fields. + let segd = ok("survey-apply", s63, "--sps=" & sps) + doAssert "rcv station: 0 traces" in segd, segd discard fails(2, "survey-apply", s63) # needs --sps - discard fails(2, "survey-apply", s63, "--sps=" & sps, "-o", copy, "--in-place") + discard fails(2, "survey-apply", s63, "--sps=" & sps, "-o", tmp / "x.sgd", "--in-place") discard fails(2, "survey-apply", s63, "--sps=" & sps, "--journal=j") discard fails(2, "survey-apply", s63, "--sps=" & sps, "--values=nope") discard fails(2, "survey-apply", s63, "--sps=" & sps, "-o", s63) diff --git a/tests/tgeomwrite.nim b/tests/tgeomwrite.nim index 1fee24f..bf442cc 100644 --- a/tests/tgeomwrite.nim +++ b/tests/tgeomwrite.nim @@ -46,7 +46,7 @@ proc copySample(name: string): string = result = tmp / name copyFile(exampleDir / name, result) -proc segyFile(path: string, rev = 1, traces = 4) = +proc segyFile(path: string, rev = 1, traces = 4, espnum = 0'i32) = ## A small SEG-Y file: field record 10, channels 1..n, no coordinates. var f = SegyFile(text: initTextHeader(teEbcdic), byteOrder: boBig, format: some(sfFloat32), revision: SegyRevision(major: uint8(rev)), extended: initExtendedTextHeader(), @@ -58,7 +58,7 @@ proc segyFile(path: string, rev = 1, traces = 4) = for i in 0.. 0 @@ -240,6 +246,25 @@ block convert_with_geometry_from_segy: k.close() doAssert not fileExists(dst & ".tmp") and not fileExists(dst & ".tmp.geometry.nif") +block values_that_round_to_the_bytes_there_are_no_change: + # Shotpoint 99.5 into a 4-byte espnum that already holds 100: the planned + # bytes equal the file's, so nothing changes (it used to count as a write + # the journal had dropped, a fatal error). + let path = tmp / "r.sgy" + segyFile(path, rev = 0, traces = 1, espnum = 100) # Rev 0: no scaled shotpoint field + let g = sps("r.sps", r21('S', "20", "99.5", '1', "500000.0", "6000400.0", "10.0"), + r21('R', "10", "101", '1', "500100.0", "6000000.0", "12.5"), + x21(10, "20", "99.5", 1, 1, "10", "101", "101")) + let r = copyWithGeometry(path, tmp / "r2.sgy", g, GeometryFields(sourceXy: false, + receiverXy: false, elevations: false, midpoints: false, offsets: false)) + doAssert not r.log.isFatal and not r.log.hasErrors, $r.log.len + doAssert r.changes.traces == 0 and r.changes.sourceStations == 0, $r.changes + var f = openSeismic(path, keys = true) + doAssert f.applySurvey(g) == 1 + doAssert f.geometryChanges(GeometryFields(sourceXy: false, receiverXy: false, + elevations: false, midpoints: false, offsets: false)).traces == 0 + f.close() + block refusals: let path = tmp / "g.sgy" segyFile(path) diff --git a/tests/tsurvey.nim b/tests/tsurvey.nim index 0aef37a..0778699 100644 --- a/tests/tsurvey.nim +++ b/tests/tsurvey.nim @@ -3,11 +3,11 @@ ## (`docs/spec/sps.md` §5.3, `docs/spec/p1_p2_p3.md` §3.4); the extracted ## spec text collapses blanks, so its lines cannot be used verbatim. -import std/[assertions, options, streams, strutils, math, os, tables] +import std/[assertions, options, streams, strutils, math, os, tables, algorithm] import parseSEG import parseSEG/common/ebcdic import parseSEG/survey/cards # column helpers (internal) -import ./utils/testutils +import ../tools/demofiles proc rec(first: char, fields: openArray[(int, string)]): string = ## An 80-column record with `fields` placed at their 1-based columns. @@ -271,26 +271,232 @@ block join_line_names: doAssert stationKey(4960, 12340, 3) in g4.table.sources doAssert g4.table.sources[stationKey(4960, 12340, 3)].z == 857 -sampleBlock join_to_a_real_file: - # 63.sgd (SEG-D 2.0): file 63, one trace (channel 1), stations in THE#1. +block join_to_a_real_file: + # The demo record.sgd (SEG-D 2.0): file 63, one trace (channel 1), stations + # in THE#1: source 20/99.5, receiver 10/104. The X record of file 63 names + # other stations (as when field record numbers differ from the relation + # file's): it is never used, since the header has both stations. let dir = getTempDir() / "tsurvey2" createDir(dir) - let x = x21("T1", 63, "20.00", "205.00", '1', 1, 1, '1', "10.00", "101.00", "101.00", '1') - let r = r21('R', "10.00", "101.00", '1', "G1", "500100.0", "6000000.0", "12.5") - let s = r21('S', "20.00", "205.00", '1', "V6", "500000.0", "6000400.0", "10.0", " 19001150") - writeFile(dir / "t.sps", [H00, r, s, x].join("\n")) - var f = openSeismic("SEG Spec/exampleFiles/63.sgd", keys = true) - doAssert f.index.keys[0].record == 63 and f.index.keys[0].channel == 1 + writeRecordSegd(dir / "record.sgd") + let x = x21("T1", 63, "30.00", "205.00", '1', 1, 1, '1', "40.00", "101.00", "101.00", '1') + let r = r21('R', "10.00", "104.00", '1', "G1", "500100.0", "6000000.0", "12.5") + let r2 = r21('R', "40.00", "101.00", '1', "G1", "700100.0", "6000000.0", "12.5") + let s = r21('S', "20.00", "99.50", '1', "V6", "500000.0", "6000400.0", "10.0", " 19001150") + let s2 = r21('S', "30.00", "205.00", '1', "V6", "700000.0", "6000400.0", "10.0", " 19001150") + writeFile(dir / "t.sps", [H00, r, r2, s, s2, x].join("\n")) var g: SurveyGeometry g.addSps(readSps(dir / "t.sps")) - doAssert f.index.applySurvey(g) == 0, "the headers' own stations are kept" - doAssert f.index.applySurvey(g, overwrite = true) == 1 - let k = f.index.keys[0] - doAssert k.srcPoint == 205 and k.rcvPoint == 101 and k.rcvX == 500100 - doAssert abs(k.offset - hypot(100.0, 400.0)) < 1e-9 - f.close() + for overwrite in [false, true]: + var f = openSeismic(dir / "record.sgd", keys = true) + doAssert f.index.keys[0].record == 63 and f.index.keys[0].channel == 1 + var stats: JoinStats + doAssert f.applySurvey(g, stats, overwrite) == 1 + let k = f.index.keys[0] + doAssert k.srcLine == 20 and k.srcPoint == 99.5 and k.rcvLine == 10 and k.rcvPoint == 104, + "header stations are never replaced" + doAssert k.rcvX == 500100 and k.srcX == 500000 + doAssert abs(k.offset - hypot(100.0, 400.0)) < 1e-9 + doAssert stats.traces == 1 and stats.changed == 1 and stats.matches[jmStations] == 1 + doAssert stats.noSourceXy == 0 and stats.noReceiverXy == 0 + doAssert $stats == "1 traces, 1 changed; by stations 1", $stats + f.close() removeDir(dir) +# ------------------------------------------------------ station-first join --- + +# Two shots on one receiver line 100 (101..104). Shot 129/113 is record 77 +# (channels 1-4 → 101..104) and repeated as record 78; shot 511/117 is +# record 151 (channels 1-4 → 104..101). A field file numbered 151 that holds +# shot 129/113 must still join to record 77's stations. +let joinSurvey = [H00, + r21('S', "129.00", "113.00", '1', "E1", "1000.0", "5000.0", "10.0", " 19001150"), + r21('S', "511.00", "117.00", '1', "V1", "2000.0", "5000.0", "10.0", " 19001150"), + r21('R', "100.00", "101.00", '1', "G1", "1100.0", "6000.0", "1.0"), + r21('R', "100.00", "102.00", '1', "G1", "1200.0", "6000.0", "2.0"), + r21('R', "100.00", "103.00", '1', "G1", "1300.0", "6000.0", "3.0"), + r21('R', "100.00", "104.00", '1', "G1", "1400.0", "6000.0", "4.0"), + x21("T1", 77, "129.00", "113.00", '1', 1, 4, '1', "100.00", "101.00", "104.00", '1'), + x21("T1", 78, "129.00", "113.00", '1', 1, 4, '1', "100.00", "101.00", "104.00", '1'), + x21("T1", 151, "511.00", "117.00", '1', 1, 4, '1', "100.00", "104.00", "101.00", '1')] + +proc traceAt(record, channel: float64, src, rcv: (float64, float64), srcIndex = 0'i16): + TraceKeys = + result = initTraceKeys() + result.status = tsLive + (result.record, result.channel) = (record, channel) + (result.srcLine, result.srcPoint, result.srcIndex) = (src[0], src[1], srcIndex) + (result.rcvLine, result.rcvPoint) = rcv + +block join_station_first: + var g: SurveyGeometry + g.addSps(parseSps(joinSurvey.join("\n"))) + const none = (NaN, NaN) + # Both stations in the header: the record number (151 = another shot) is + # not used, with or without overwrite; coordinates come by station. + for overwrite in [false, true]: + var k = traceAt(151, 2, (129.0, 113.0), (100.0, 102.0)) + doAssert k.applySurvey(g, overwrite) + doAssert k.srcPoint == 113 and k.rcvPoint == 102 and k.srcX == 1000 and k.rcvX == 1200 + doAssert k.srcIndex == 0, "the header's index stays; 0 matches the survey's 1" + var k = traceAt(151, 2, (129.0, 113.0), (100.0, 102.0)) + doAssert k.joinStations(g) == jmStations + # Source only: the receiver comes from the source's records (77 and 78 + # agree), not from record 151 (which would give 103). + k = traceAt(151, 2, (129.0, 113.0), none) + doAssert k.joinStations(g) == jmSource and k.rcvLine == 100 and k.rcvPoint == 102 + k = traceAt(151, 2, (129.0, 113.0), none) + doAssert k.applySurvey(g) and k.rcvX == 1200 + # Receiver only: no O records, so record and channel, checked by the receiver. + k = traceAt(151, 2, none, (100.0, 103.0)) + doAssert k.joinStations(g) == jmRecord and k.srcLine == 511 and k.srcPoint == 117 + k = traceAt(151, 2, none, (100.0, 102.0)) + doAssert k.joinStations(g) == jmConflict and k.srcLine.isNaN, "record 151 ch 2 is 100/103" + # No station: the record number is the only key. + k = traceAt(151, 2, none, none) + doAssert k.joinStations(g) == jmRecord and k.srcPoint == 117 and k.rcvPoint == 103 + k = traceAt(152, 2, none, none) + doAssert k.joinStations(g) == jmUnmatched and k.srcLine.isNaN + k = traceAt(77, NaN, none, none) + doAssert k.joinStations(g) == jmUnmatched + # A source not in the relations, or a channel its records lack. + k = traceAt(77, 2, (130.0, 113.0), none) + doAssert k.joinStations(g) == jmUnmatched and k.rcvLine.isNaN + k = traceAt(77, 9, (129.0, 113.0), none) + doAssert k.joinStations(g) == jmUnmatched + # Auxiliary traces are left alone. + k = traceAt(151, 2, none, none) + k.status = tsAuxiliary + doAssert k.joinStations(g) == jmAuxiliary and not k.applySurvey(g) and k.srcLine.isNaN + # Repeats that disagree: ambiguous, unless the record number picks one. + var g2: SurveyGeometry + g2.addSps(parseSps((joinSurvey[0 .. ^4] & @[ + x21("T1", 77, "129.00", "113.00", '1', 1, 4, '1', "100.00", "101.00", "104.00", '1'), + x21("T1", 78, "129.00", "113.00", '1', 1, 4, '1', "100.00", "104.00", "101.00", '1')]).join("\n"))) + k = traceAt(151, 2, (129.0, 113.0), none) + doAssert k.joinStations(g2) == jmAmbiguous and k.rcvLine.isNaN + k = traceAt(78, 2, (129.0, 113.0), none) + doAssert k.joinStations(g2) == jmSource and k.rcvPoint == 103 + # A point index that does not match is another station. + k = traceAt(151, 2, (129.0, 113.0), none, srcIndex = 2) + doAssert k.joinStations(g) == jmUnmatched + doAssert sameStation(stationKey(1, 2, 0), stationKey(1, 2.000001, 1)) + doAssert not sameStation(stationKey(1, 2, 2), stationKey(1, 2, 1)) + +block join_enumeration_and_lookups: + var g: SurveyGeometry + g.addSps(parseSps(joinSurvey.join("\n"))) + var recs: seq[int64] + for r in g.records: recs.add r + doAssert recs == @[77'i64, 78, 151] + var n = 0 + for (c, rel) in g.relations(151): + inc n + doAssert g.relation(151, c).get == rel + doAssert rel.receiver.point == 105 - float64(c) and rel.source.point == 117 + doAssert n == 4 and g.expectedChannels(151) == 4 and g.expectedChannels(5) == 0 + doAssert g.recordsAt(stationKey(129, 113, 0)) == @[77'i64, 78] + doAssert g.recordsAt(stationKey(129, 113, 1)) == @[77'i64, 78] + doAssert g.recordsAt(stationKey(129, 113, 2)).len == 0 + doAssert g.recordsAt(stationKey(100, 101, 1), rdReceiver).len == 0 + doAssert g.channelOf(77, stationKey(100, 103, 0)) == some(3'i64) + doAssert g.channelOf(151, stationKey(100, 103, 1)) == some(2'i64) + doAssert g.channelOf(151, stationKey(100, 105, 1)).isNone + doAssert g.channelOf(151, stationKey(101, 103, 1)).isNone + # Overlapping runs: a channel is yielded once, for the earlier run. + var g3: SurveyGeometry + g3.addSps(parseSps([H00, + x21("T1", 5, "1.00", "1.00", '1', 1, 4, '1', "100.00", "101.00", "104.00", '1'), + x21("T1", 5, "1.00", "1.00", '1', 3, 6, '1', "200.00", "101.00", "104.00", '1')].join("\n"))) + var lines: seq[float64] + for (c, rel) in g3.relations(5): lines.add rel.receiver.line + doAssert lines == @[100.0, 100, 100, 100, 200, 200] and g3.log.count(sevWarning) == 1 + doAssert g3.channelOf(5, stationKey(200, 101, 1)).isNone, "channel 3 belongs to line 100" + # Multi-component: k channels per station. + var g4 = SurveyGeometry(multiComponent: true) + g4.addSps(parseSps([H00, x21("T1", 7, "1.00", "1.00", '1', 1, 7, '3', "100.00", "101.00", + "103.00", '1')].join("\n"))) + var comps: seq[int] + for (c, rel) in g4.relations(7): comps.add rel.component + doAssert comps == @[1, 2, 3, 1, 2, 3, 1, 2, 3] + doAssert g4.channelOf(7, stationKey(100, 102, 1)) == some(4'i64) + +# --------------------------------------------------------- receiver domain --- + +proc pairs(g: SurveyGeometry, domain: RelationDomain): seq[(float64, float64, float64, float64)] = + for r in g.records(domain): + for (_, rel) in g.relations(r, domain): + result.add (rel.source.line, rel.source.point, rel.receiver.line, rel.receiver.point) + result.sort + +block receiver_domain_o_records: + # Read: an O record is a receiver's gather; the X layout with the roles swapped. + let o = rec('O', [(2, "T1"), (8, align("12", 8)), (16, "1"), (17, "1"), + (18, "100.00".align(10)), (28, "102.00".align(10)), (38, "1"), (39, align("1", 5)), + (44, align("3", 5)), (49, "1"), (50, "129.00".align(10)), (60, "113.00".align(10)), + (70, "115.00".align(10)), (80, "1")]) + let f = parseSps([H00, o].join("\n") & "\n") + doAssert f.log.len == 0 and f.relations.len == 1 and f.relations[0].recordType == 'O' + doAssert spsText(f).text == [H00, o].join("\n") & "\n" + var edited = f + edited.relations[0].raw = "" + doAssert spsText(edited).text == [H00, o].join("\n") & "\n", "encoded with its letter" + var g: SurveyGeometry + g.addSps(f) + doAssert g.relationCount(rdReceiver) == 1 and g.relationCount == 0 + let rel = g.relation(12, 2, rdReceiver).get + doAssert rel.receiver == stationKey(100, 102, 1) and rel.source == stationKey(129, 114, 1) + doAssert g.recordsAt(stationKey(100, 102, 0), rdReceiver) == @[12'i64] + # A receiver gather trace with only its receiver station joins via the O gathers, + # whatever its record number. + var k = traceAt(NaN, 3, (NaN, NaN), (100.0, 102.0)) + doAssert k.joinStations(g) == jmReceiver and k.srcLine == 129 and k.srcPoint == 115 + k = traceAt(NaN, 4, (NaN, NaN), (100.0, 102.0)) + doAssert k.joinStations(g) == jmUnmatched + # `.ops` files may hold O relations as X records. + doAssert relationDomainOf("a/b.OPS") == rdReceiver and relationDomainOf("b.xps") == rdSource + var g2: SurveyGeometry + g2.addSps(parseSps([H00, "X" & o[1 .. ^1]].join("\n")), xDomain = rdReceiver) + doAssert g2.relation(12, 2, rdReceiver) == some(rel) + var bad = f.relations[0] + bad.recordType = 'Y' + doAssertRaises(SegFormatError): discard encodeSps(bad, sps21) + +block receiver_domain_from_x: + # Sources 20/201..203 record receivers 10/101..105; 21/201 records 10/103..107 + # (descending); a repeat of 20/202 (record 4) records 10/101..102; line 11 too. + let x = [H00, + x21("T1", 1, "20.00", "201.00", '1', 1, 5, '1', "10.00", "101.00", "105.00", '1'), + x21("T1", 2, "20.00", "202.00", '1', 1, 5, '1', "10.00", "101.00", "105.00", '1'), + x21("T1", 3, "20.00", "203.00", '1', 1, 5, '1', "10.00", "101.00", "105.00", '1'), + x21("T1", 3, "20.00", "203.00", '1', 6, 7, '1', "11.00", "101.00", "102.00", '1'), + x21("T1", 4, "20.00", "202.00", '1', 1, 2, '1', "10.00", "101.00", "102.00", '1'), + x21("T1", 5, "21.00", "201.00", '1', 1, 5, '1', "10.00", "107.00", "103.00", '1')] + var g: SurveyGeometry + g.addSps(parseSps(x.join("\n"))) + let rels = g.receiverRelations() + doAssert rels.len > 0 and rels[0].fieldRecord == 1 and rels[0].recordType == 'X' + # Gather 1 is receiver 10/101: sources 20/201..203 as one run, then the + # repeat of 20/202 as a run of its own. + doAssert rels[0].srcLine == 10 and rels[0].srcPoint == 101 + doAssert rels[0].rcvLine == 20 and rels[0].fromReceiver == 201 and rels[0].toReceiver == 203 + doAssert rels[0].fromChannel == 1 and rels[0].toChannel == 3 + let o = SpsFile(revision: sps21, headers: parseSps(H00).headers, relations: rels) + let text = spsText(o) + doAssert not text.log.hasErrors + let back = parseSps(text.text) + doAssert back.log.len == 0 and spsText(back).text == text.text, "reads back byte-identical" + var go: SurveyGeometry + go.addSps(back, xDomain = rdReceiver) + doAssert go.pairs(rdReceiver) == g.pairs(rdSource), "X → O → pairs equals X's pairs" + var gathers: seq[int64] + for r in go.records(rdReceiver): gathers.add r + doAssert gathers == @[1'i64, 2, 3, 4, 5, 6, 7, 8, 9], $gathers # 10/101..107, 11/101..102 + doAssert go.recordsAt(stationKey(11, 102, 1), rdReceiver) == @[9'i64] + let oo = g.receiverRelations(firstGather = 100, recordType = 'O') + doAssert oo[0].recordType == 'O' and oo[0].fieldRecord == 100 + doAssertRaises(ValueError): discard g.receiverRelations(recordType = 'Z') + # --------------------------------------------------------- review (10.4) --- block review_p1_non_p1_lines_stay_raw: From d2c7efeea775a7916f00d021a7227c68ce81c72d Mon Sep 17 00:00:00 2001 From: xrfez Date: Wed, 11 Nov 2026 19:11:33 +0000 Subject: [PATCH 2/3] receiverRelations: write O records by default (X on request) Co-Authored-By: Claude Opus 5.5 --- docs/api.md | 2 +- docs/guide/conversion.md | 4 ++-- src/parseSEG/survey/join.nim | 7 ++++--- tests/tsurvey.nim | 12 ++++++++---- 4 files changed, 15 insertions(+), 10 deletions(-) diff --git a/docs/api.md b/docs/api.md index 5f211c2..7cd5df8 100644 --- a/docs/api.md +++ b/docs/api.md @@ -113,7 +113,7 @@ in the receiver domain. Traces are joined by their stations; record numbers are | `writeP1(f, path or stream, P1WriteOptions)`, `p1Text(f)` (a `SurveyText`), `encodeP1(record, opts)` | Header blocks padded to 20 records; unchanged records byte for byte. | | `SurveyGeometry` (`initSurveyGeometry()`): `addSps(g, sps, xDomain = rdSource)`, `addSps(g, path, opts, xDomain = relationDomainOf(path))`, `addP1(g, p1 or path[, P1Options], elevation)`, `addSource`/`addReceiver`, `relation(g, record, channel, domain = rdSource)` (an `Option[StationRelation]`), `resolveLine(g, name, numeric)` (an `Option`) | Station coordinates (`table`, `sourceCount`, `receiverCount`) and relation runs (`relationCount(g, domain)`) of both domains (`RelationDomain`: `rdSource` X, `rdReceiver` O; X records of an `.ops` file are O by default). Text line names go through `mapLine(g, name, line)` (`lineNames`), then `trailingDigits`. `multiComponent`, `tape` refine relations; problems are in `log`. | | `records(g, domain)`, `relations(g, record, domain)` (iterators: ascending records; `(channel, relation)` per channel, a channel two runs cover once), `expectedChannels(g, record, domain)`, `recordsAt(g, station, domain)`, `channelOf(g, record, station, domain)` (an `Option`), `sameStation(a, b)` | Enumerate and look up the relations: the records of a source station (X) or gathers of a receiver station (O), the channel of a station in a record. Index 0 matches any. | -| `receiverRelations(g, firstGather = 1, recordType = 'X')` → `seq[SpsRelation]` | The O records of the X relations: one gather per receiver station (receiver line, point order), its sources by line and point as runs, repeats after; write with `writeSps`. One receiver line in memory at a time. | +| `receiverRelations(g, firstGather = 1, recordType = 'O')` → `seq[SpsRelation]` | The O records of the X relations: one gather per receiver station (receiver line, point order), its sources by line and point as runs, repeats after; write with `writeSps`. One receiver line in memory at a time. | | `joinStations(keys, g)` → `JoinMatch`; `applySurvey(f, keys or index, g[, stats], overwrite = false)` | Fill the stations a header lacks, matched **by station**: both in the header (`jmStations`); source only → the X records of that source (`jmSource`); receiver only → the O gathers of that receiver (`jmReceiver`); no station → record and channel (`jmRecord`); also `jmAmbiguous`, `jmConflict`, `jmUnmatched`, `jmAuxiliary` (left alone). Header stations are never replaced; a relation must agree with the stations the header has. `applySurvey` then fills coordinates by station, offset, midpoint (`overwrite`: replace known coordinates); in memory only. `JoinStats` (`traces`, `changed`, `matches`, `noSourceXy`, `noReceiverXy`; `$`). | | `geometryChanges(f, fields[, log])` | What writing the joined keys of `f` (opened with keys, `applySurvey` done) would change, as `GeometryChanges` (traces per kind, `unwritable` values; `$` gives a one-line summary), with the reports in `log`; nothing is written. | | `writeGeometry(editor, fields = GeometryFields())`, `writeGeometry(editor, f.index, fields)` | Writes the joined geometry in place through the journaled `HeaderEditor`: from the editor's own index (an editor made with `openEditor(f)` from the joined file; `ValueError` without keys), or from the index of the same file opened separately. Same-size fields; undo with `rollbackJournal`. SEG-Y: through the layout (`sht_x`, `rec_x`, `cdp_x` with a coordinate scalar chosen to keep the precision in 4 bytes; `selev`/`relev`; `offset`; `espnum`/`sp`; Ext1 copies in use; `sht_line`, `rec_point`, ... where a Layout stanza defines them). SEG-D Rev 1/2: THE#1 receiver line/point/index. Rev 3: THE#1 and the Position / Source Description blocks already in the trace header extensions. Values with no slot or too large for theirs are reported in `editor.log`. `GeometryFields` selects source/receiver X/Y, elevations, midpoints, offsets, stations. SEG-2 is refused for now. | diff --git a/docs/guide/conversion.md b/docs/guide/conversion.md index a8d1ede..3a66d78 100644 --- a/docs/guide/conversion.md +++ b/docs/guide/conversion.md @@ -141,12 +141,12 @@ doAssert changed == 1 and k.rcvX == 500100 and k.srcPoint == 99.5 doAssert stats.matches[jmStations] == 1 # the header's own stations were the key # The receiver domain: O records built from the X relations, written as an -# .ops file (X records, as such files hold them) and read back. +# .ops file and read back. let ops = SpsFile(revision: sps21, headers: readSps("line.sps").headers, relations: g.receiverRelations()) doAssert not writeSps(ops, "line.ops").hasErrors var go = initSurveyGeometry() -go.addSps("line.ops") # `.ops`: its X records are receiver gathers +go.addSps("line.ops") # O records are receiver gathers doAssert go.relationCount(rdReceiver) == 1 doAssert go.recordsAt(stationKey(10, 104, 1), rdReceiver) == @[1'i64] ``` diff --git a/src/parseSEG/survey/join.nim b/src/parseSEG/survey/join.nim index 6ca8f1f..84dd534 100644 --- a/src/parseSEG/survey/join.nim +++ b/src/parseSEG/survey/join.nim @@ -597,7 +597,7 @@ func stationLess(a, b: StationKey): int = else: cmp(a.index, b.index) proc receiverRelations*(g: SurveyGeometry, firstGather = 1'i64, - recordType = 'X'): seq[SpsRelation] = + recordType = 'O'): seq[SpsRelation] = ## The receiver-domain (O) records of the X relations: one gather per ## receiver station, numbered from `firstGather` in receiver line, point ## order; its traces 1..n are the sources that recorded it, by source @@ -605,8 +605,9 @@ proc receiverRelations*(g: SurveyGeometry, firstGather = 1'i64, ## (one X record per repeat) is one more trace, after the first ones, so ## each X pair appears once. Components are not ## represented (one trace per source and receiver station). Write them with - ## `writeSps` (records start with `recordType`, `X` as in `.ops` files, or - ## `O`); `addSps(..., xDomain = rdReceiver)` reads them back. Works one + ## `writeSps` (records start with `recordType`: `O`, or `X` for software + ## that expects X records in `.ops` files; `addSps(..., xDomain = + ## rdReceiver)` reads those back). Works one ## receiver line at a time, so memory follows the largest line, not the ## survey. Raises `ValueError` for a `recordType` other than `X` or `O`. if recordType notin {'X', 'O'}: diff --git a/tests/tsurvey.nim b/tests/tsurvey.nim index 0778699..45cfbdd 100644 --- a/tests/tsurvey.nim +++ b/tests/tsurvey.nim @@ -475,7 +475,7 @@ block receiver_domain_from_x: var g: SurveyGeometry g.addSps(parseSps(x.join("\n"))) let rels = g.receiverRelations() - doAssert rels.len > 0 and rels[0].fieldRecord == 1 and rels[0].recordType == 'X' + doAssert rels.len > 0 and rels[0].fieldRecord == 1 and rels[0].recordType == 'O' # Gather 1 is receiver 10/101: sources 20/201..203 as one run, then the # repeat of 20/202 as a run of its own. doAssert rels[0].srcLine == 10 and rels[0].srcPoint == 101 @@ -487,14 +487,18 @@ block receiver_domain_from_x: let back = parseSps(text.text) doAssert back.log.len == 0 and spsText(back).text == text.text, "reads back byte-identical" var go: SurveyGeometry - go.addSps(back, xDomain = rdReceiver) + go.addSps(back) # O records: receiver domain doAssert go.pairs(rdReceiver) == g.pairs(rdSource), "X → O → pairs equals X's pairs" var gathers: seq[int64] for r in go.records(rdReceiver): gathers.add r doAssert gathers == @[1'i64, 2, 3, 4, 5, 6, 7, 8, 9], $gathers # 10/101..107, 11/101..102 doAssert go.recordsAt(stationKey(11, 102, 1), rdReceiver) == @[9'i64] - let oo = g.receiverRelations(firstGather = 100, recordType = 'O') - doAssert oo[0].recordType == 'O' and oo[0].fieldRecord == 100 + let oo = g.receiverRelations(firstGather = 100, recordType = 'X') + doAssert oo[0].recordType == 'X' and oo[0].fieldRecord == 100 + var gx: SurveyGeometry + gx.addSps(SpsFile(revision: sps21, relations: oo), xDomain = rdReceiver) + doAssert gx.pairs(rdReceiver) == g.pairs(rdSource), "X-lettered O records, read as O" + doAssertRaises(ValueError): discard g.receiverRelations(recordType = 'Z') # --------------------------------------------------------- review (10.4) --- From 601494334576301021e1bd8886584bdc2588a06d Mon Sep 17 00:00:00 2001 From: xrfez Date: Wed, 11 Nov 2026 19:11:33 +0000 Subject: [PATCH 3/3] Release 0.9.0 New since 0.8.0: survey joins match traces by their stations (header stations are never replaced; record and channel only for traces without stations; overwrite refreshes coordinates only); JoinMatch / JoinStats; receiver-domain O relations (.ops), receiverRelations to build them from X; relation lookups (records, relations, expectedChannels, recordsAt, channelOf). Fix: a geometry value that rounds to the bytes already in the file is no change (it was reported as a dropped journal write). Co-Authored-By: Claude Opus 5.5 --- README.md | 4 ++-- parseSEG.nimble | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/README.md b/README.md index 182e375..3ba0d05 100644 --- a/README.md +++ b/README.md @@ -26,10 +26,10 @@ atlas use https://github.com/xrfez/parseSEG With Nimble, in your `.nimble` file: ```nim -requires "https://github.com/xrfez/parseSEG >= 0.8.0" +requires "https://github.com/xrfez/parseSEG >= 0.9.0" ``` -Releases are git tags (`v0.8.0`). Then: +Releases are git tags (`v0.9.0`). Then: ```nim import parseSEG # the whole public API diff --git a/parseSEG.nimble b/parseSEG.nimble index b65da4f..f8484b8 100644 --- a/parseSEG.nimble +++ b/parseSEG.nimble @@ -1,6 +1,6 @@ # Package -version = "0.8.0" +version = "0.9.0" author = "Joshua Fenner" description = "Read, write, index and edit SEG-Y (Rev 0-2.1) and SEG-D (Rev 2.0-3.0) seismic data on disk and tape" license = "MIT"