Use extension 1 in rotation calculations

Author: achaikou

lib/src/segy.c

@@ -3257,7 +3257,7 @@ static int scaled_standard_header_cdp(
     const segy_entry_definition* map,
     int cdp_offset,
     int scalar_offset,
-    float* cdp
+    double* cdp
 ) {
 
     segy_field_data fd = {0};
@@ -3294,39 +3294,86 @@ static int scaled_standard_header_cdp(
     return SEGY_OK;
 }
 
-static int scaled_cdp(
+static int scaled_ext1_header_cdp(
+    const char* header,
+    const segy_entry_definition* map,
+    int cdp_offset,
+    double* cdp
+) {
+
+    segy_field_data fd = { 0 };
+    int err = segy_get_tracefield( header, map, cdp_offset, &fd );
+    if( err != SEGY_OK ) return err;
+
+    switch( fd.entry_type ) {
+        case SEGY_ENTRY_TYPE_IEEE64:
+            *cdp = fd.value.f64;
+            return SEGY_OK;
+            break;
+        default:
+            return SEGY_INVALID_FIELD_DATATYPE;
+    }
+}
+
+static int scaled_cdp_in_dimension(
     segy_datasource* ds,
     int traceno,
-    float* cdpx,
-    float* cdpy
+    int dimension_standard_name,
+    int dimension_ext1_name,
+    double* cdp
 ) {
 
     char trheader[SEGY_TRACE_HEADER_SIZE];
 
-    int err = segy_read_standard_traceheader( ds, traceno, trheader );
-    if( err != 0 ) return err;
+    const segy_entry_definition* ext1_map =
+        ds->traceheader_mapping_extension1.offset_to_entry_definition;
+    const int cdp_ext1_offset =
+        ds->traceheader_mapping_extension1.name_to_offset[dimension_ext1_name];
+
+    if( ext1_map[cdp_ext1_offset - 1].entry_type != SEGY_ENTRY_TYPE_UNDEFINED ) {
+        int err = segy_read_traceheader( ds, traceno, 1, ext1_map, trheader );
+        if( err != SEGY_OK ) return err;
+
+        err = scaled_ext1_header_cdp(
+            trheader, ext1_map, cdp_ext1_offset, cdp
+        );
+        if( err != SEGY_OK ) return err;
+        if( *cdp != 0 || !ext1_map[cdp_ext1_offset - 1].requires_nonzero_value ) {
+            return SEGY_OK;
+        }
+    }
 
     const segy_entry_definition* standard_map =
         ds->traceheader_mapping_standard.offset_to_entry_definition;
 
-    const int cdp_x_offset =
-        ds->traceheader_mapping_standard.name_to_offset[SEGY_TR_CDP_X];
-    const int cdp_y_offset =
-        ds->traceheader_mapping_standard.name_to_offset[SEGY_TR_CDP_Y];
+    const int cdp_standard_offset =
+        ds->traceheader_mapping_standard.name_to_offset[dimension_standard_name];
     const int scalar_offset =
         ds->traceheader_mapping_standard.name_to_offset[SEGY_TR_SOURCE_GROUP_SCALAR];
 
-    err = scaled_standard_header_cdp(
-        trheader, standard_map, cdp_x_offset, scalar_offset, cdpx
-    );
+    int err = segy_read_standard_traceheader( ds, traceno, trheader );
     if( err != SEGY_OK ) return err;
 
     err = scaled_standard_header_cdp(
-        trheader, standard_map, cdp_y_offset, scalar_offset, cdpy
+        trheader, standard_map, cdp_standard_offset, scalar_offset, cdp
     );
     return err;
 }
 
+static int scaled_cdp(
+    segy_datasource* ds,
+    int traceno,
+    double* cdpx,
+    double* cdpy
+) {
+    int err;
+    err = scaled_cdp_in_dimension( ds, traceno, SEGY_TR_CDP_X, SEGY_EXT1_CDP_X, cdpx );
+    if( err != SEGY_OK ) return err;
+
+    err = scaled_cdp_in_dimension( ds, traceno, SEGY_TR_CDP_Y, SEGY_EXT1_CDP_Y, cdpy );
+    return err;
+}
+
 int segy_rotation_cw( segy_datasource* ds,
                       int line_length,
                       int stride,
@@ -3335,7 +3382,7 @@ int segy_rotation_cw( segy_datasource* ds,
                       int linenos_sz,
                       float* rotation) {
 
-    struct coord { float x, y; } nw, sw;
+    struct coord { double x, y; } nw, sw;
 
     int err;
     int traceno;
@@ -3355,8 +3402,8 @@ int segy_rotation_cw( segy_datasource* ds,
     err = scaled_cdp( ds, traceno, &nw.x, &nw.y );
     if( err != 0 ) return err;
 
-    float x = nw.x - sw.x;
-    float y = nw.y - sw.y;
+    double x = nw.x - sw.x;
+    double y = nw.y - sw.y;
     double radians = x || y ? atan2( x, y ) : 0;
     if( radians < 0 ) radians += 2 * acos(-1);
 

python/examples/make-rotated-rev2.py

@@ -0,0 +1,60 @@
+import sys
+import itertools as itr
+import segyio
+import segyio.su as su
+
+
+def product(f):
+    return itr.product(range(len(f.ilines)), range(len(f.xlines)))
+
+# this program turns the rev1 file into rev2 file with CDP values set to
+# different set of cdp-x and cdp-y coordinates in standard and ext1 headers.
+# The formulas are taken from make-rotated-copies.py script
+
+def main():
+    if len(sys.argv) != 3:
+        sys.exit(
+            "Usage: {} [source-file] [destination-file]".format(sys.argv[0]))
+
+    srcfile = sys.argv[1]
+    dstfile = sys.argv[2] if len(sys.argv) > 2 else srcfile
+
+    with segyio.open(srcfile) as src:
+        spec = segyio.spec()
+        spec.format = int(src.format)
+        spec.sorting = int(src.sorting)
+        spec.samples = src.samples
+        spec.ilines = src.ilines
+        spec.xlines = src.xlines
+        spec.traceheader_count = 2
+        spec.tracecount = src.tracecount
+
+        with segyio.create(dstfile, spec) as dst:
+            for i in range(1 + src.ext_headers):
+                dst.text[i] = src.text[i]
+            dst.bin = src.bin
+            dst.trace = src.trace
+            dst.traceheader = src.traceheader
+
+            dst.bin[segyio.BinField.SEGYRevision] = 2
+            dst.bin[segyio.BinField.MaxAdditionalTraceHeaders] = 1
+
+            for i in range(src.tracecount):
+                dst.traceheader[i][1].header_name = bytes("SEG00001", "ascii")
+
+    with segyio.open(dstfile, 'r+') as dst:
+        for i, (x, y) in enumerate(product(src)):
+            trh = dst.traceheader[i][0]
+            # "right" rotation formula
+            trh.cdp_x = y
+            trh.cdp_y = 100 - x
+            trh.co_scal = 1
+
+            trh = dst.traceheader[i][1]
+            # "left" rotation formula
+            trh.cdp_x = (100 - y) * 21
+            trh.cdp_y = x * 21
+
+
+if __name__ == '__main__':
+    main()

python/test/tools.py

@@ -174,6 +174,20 @@ def test_rotation_lsb():
             assert rotation(msb, line = 'fast') == rotation(lsb, line = 'fast')
             assert rotation(msb, line = 'slow') == rotation(lsb, line = 'slow')
 
+def test_rotation_ext1():
+    # note that in the test file ext1 values are taken from "left" and standard
+    # header values are taken from "right", yet the result does not correspond
+    # to "left" results. This happens because some of the values in the ext1
+    # header are 0 and 'use-if-non-zero' flag is set for ext1 header. So these
+    # values fall down to 'right' standard header instead.
+    #
+    # While in reality it is likely that standard header would be 0-ed, for us
+    # it is good opportunity to test 'use-if-non-zero' flag. It was checked that
+    # test result values correspond with the mixed left-right picked data.
+    with segyio.open(testdata / 'rotated-small-rev2.sgy') as f:
+        assert 4.712 == approx(segyio.tools.rotation(f, line = 'fast')[0], abs = 1e-3)
+        assert 3.142 == approx(segyio.tools.rotation(f, line = 'slow')[0], abs = 1e-3)
+
 def test_metadata():
     spec = segyio.spec()
     spec.ilines = [1, 2, 3, 4, 5]

test-data/README.md

@@ -57,6 +57,12 @@ To recreate:
 | left-small.sgy      | 270°                                                                  |
 | inv-acute-small.sgy | 315°                                                                  |
 
+File below is used to test revision 2 features and is created with:
+`python python/examples/make-rotated-rev2.py test-data/small.sgy test-data/rotated-small-rev2.sgy`
+
+| File                   | Purpose                                                                    |
+|------------------------|--------------------------------------------------------------------------- |
+| rotated-small-rev2.sgy | Revision 2 file with different values on standard and extension 1 headers. |
 
 ## Dimensions sorting test files