parseSEG 0.9.0: survey joins match by stations; receiver-domain (O) relations

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

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 = '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. |
 | `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. |

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 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")                     # O 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

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.

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"

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 <file> 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()