Work in progress on extraction to numpy data structures.

examples/extract_trace_header_field.py

@@ -0,0 +1,113 @@
+#!/usr/bin/env python3
+
+"""Extract a single field from the trace headers of  3D seismic volume to a 2D Numpy array.
+
+This utility assumes the inline and crossline numbers are evenly spaced.
+Each inline of the source data will be represented as a single row, and
+each crossline as a single column in the resulting 2D array.
+
+Usage: extract_trace_header_field.py [-h] [--null NULL]
+                    segy-file npy-file field-name
+
+Positional arguments:
+  segy-file      Path to an existing SEG Y file of 3D seismic data
+  npy-file       Path to the Numpy array file to be created for the timeslice
+  field-name     Zero based index of the time slice to be extracted
+
+Optional arguments:
+  -h, --help     show this help message and exit
+  --dtype DTYPE  Numpy data type. If not provided a dtype compatible with the
+                 SEG Y data will be used.
+  --null NULL    Sample value to use for missing or short traces.
+
+Example:
+
+  extract_trace_header_field.py stack_final_int8.sgy slice_800.npy ensemble_num
+"""
+
+import argparse
+import os
+import sys
+import traceback
+
+import numpy as np
+
+from numpy import s_
+
+from segpy.reader import create_reader
+from segpy.trace_header import TraceHeaderRev1
+from segpy_numpy.dtypes import make_dtype
+from segpy_numpy.extract import extract_inline_3d, extract_trace_header_field_3d
+
+
+class DimensionalityError(Exception):
+    pass
+
+
+
+def extract_header_field(segy_filename, out_filename, field_name, null=None):
+    """Extract a timeslice from a 3D SEG Y file to a Numpy NPY file.
+
+    Args:
+        segy_filename: Filename of a SEG Y file.
+
+        out_filename: Filename of the NPY file.
+
+        inline_index: The zero-based index (increasing with depth) of the slice to be extracted.
+
+        null: Optional sample value to use for missing or short traces. Defaults to zero.
+    """
+    header_field = getattr(TraceHeaderRev1, field_name)
+
+    with open(segy_filename, 'rb') as segy_file:
+
+        segy_reader = create_reader(segy_file)
+        header_field_array = extract_trace_header_field_3d(segy_reader, header_field, null=null)
+        return header_field_array
+
+
+def nullable_int(s):
+    return None if not bool(s) else int(s)
+
+
+def main(argv=None):
+    parser = argparse.ArgumentParser()
+    parser.add_argument("segy_file", metavar="segy-file",
+                        help="Path to an existing SEG Y file of 3D seismic data")
+
+    parser.add_argument("npy_file", metavar="npy-file",
+                        help="Path to the Numpy array file to be created for the timeslice")
+
+    parser.add_argument("field_name", metavar="field-name", type=str,
+                        help="Name of the trace header field to be extracted", )
+
+    parser.add_argument("--null",  type=nullable_int, default="",
+                        help="Header value to use for missing or short traces.")
+
+    if argv is None:
+        argv = sys.argv[1:]
+
+    args = parser.parse_args(argv)
+
+    try:
+        extract_header_field(
+            segy_filename=args.segy_file,
+            out_filename=args.npy_file,
+            field_name=args.field_name,
+            null=args.null)
+    except (FileNotFoundError, IsADirectoryError) as e:
+        print(e, file=sys.stderr)
+        return os.EX_NOINPUT
+    except PermissionError as e:
+        print(e, file=sys.stderr)
+        return os.EX_NOPERM
+    except Exception as e:
+        traceback.print_exception(type(e), e, e.__traceback__, file=sys.stderr)
+        return os.EX_SOFTWARE
+    return os.EX_OK
+
+
+if __name__ == '__main__':
+    sys.exit(main())
+
+

examples/inline.py

@@ -22,7 +22,7 @@ Optional arguments:
 
 Example:
 
-  timeslice.py stack_final_int8.sgy slice_800.npy 800 --null=42.0 --dtype=f
+  inline.py stack_final_int8.sgy slice_800.npy 800 --null=42.0 --dtype=f
 """
 
 import argparse
@@ -58,8 +58,8 @@ def extract_inline(segy_filename, out_filename, inline_number, null=None):
     with open(segy_filename, 'rb') as segy_file:
 
         segy_reader = create_reader(segy_file)
-        inline_array = extract_inline_3d(segy_reader, inline_number, sample_numbers=s_[::2], null=null)
-        return inline_number
+        inline_array = extract_inline_3d(segy_reader, inline_number, null=null)
+        return inline_array
 
 
 def nullable_float(s):

segpy-ext/segpy-numpy/segpy_numpy/dtypes.py

@@ -9,7 +9,7 @@ NUMPY_DTYPES = {'ibm':     numpy.dtype('f4'),
                 'int8':    numpy.dtype('i1')}
 
 
-def make_dtype(data_sample_format):
+def make_dtype(data_sample_format): # TODO: What is the correct name for this arg?
     """Convert a SEG Y data sample format to a compatible numpy dtype.
 
     Note :

segpy-ext/segpy-numpy/segpy_numpy/extract.py

@@ -8,46 +8,127 @@ from segpy_numpy.dtypes import make_dtype
 class DimensionalityError:
     pass
 
+def extract_trace_header_field_3d(reader, field, null=None):
+    """Extract a single trace header field from all trace headers as an array.
 
-def extract_inline_3d(reader, inline_number, xline_numbers=None, sample_numbers=None, null=None):
+    Args:
+        reader: A SegYReader
+        field: A header field
+    """
+    shape = (reader.num_inlines(), reader.num_xlines())
+    dtype = make_dtype(field.value_type.SEG_Y_TYPE)
+    arr = _make_array(shape, dtype, null)
+    field_name = field.name
+    for inline_xline in reader.inline_xline_numbers():
+        inline_number, xline_number = inline_xline
+        trace_index = reader.trace_index(inline_xline)
+        trace_header = reader.trace_header(trace_index)
+        inline_index = reader.inline_numbers().index(inline_number)
+        xline_index = reader.xline_numbers().index(xline_number)
+        field_value = getattr(trace_header, field_name)
+        arr[inline_index, xline_index] = field_value
+    return arr
+
+
+def extract_trace(reader, trace_index, sample_numbers):
+    """Extract an single trace as a one-dimensional array.
+
+    Args:
+        reader: A SegYReader3D object.
+
+        trace_index: The index of the trace to be extracted.
+
+        sample_numbers: The sample numbers within each trace at which samples are to be extracted.
+            This argument can be specified in three ways:
+
+            None (the default) - All samples within the trace will be be extracted.
+
+            sequence - When a sequence, such as a range or a list is provided only those samples at
+                sample numbers corresponding to the items in the sequence will be extracted. The
+                samples will always be extracted in increasing numeric order and duplicate entries
+                will be ignored.  For example sample_numbers=range(100, 200, 2) will extract alternate
+                samples from sample number 100 to sample number 198 inclusive.
+
+            slice - When a slice object is provided the slice will be applied to the sequence of all
+                sample numbers. For example sample_numbers=slice(100, -100) will omit the first
+                one hundred and the last one hundred samples, irrespective of their numbers.
+
+    Returns:
+        A one-dimensional array.
+    """
+    if not reader.has_trace_index(trace_index):
+        raise ValueError("Inline number {} not present in {}".format(trace_index, reader))
+
+    sample_numbers = ensure_superset(range(0, reader.max_num_trace_samples()), sample_numbers)
+
+    trace_sample_start = sample_numbers[0]
+    trace_sample_stop = min(sample_numbers[-1] + 1, reader.num_trace_samples(trace_index))
+    trace_samples = reader.trace_samples(trace_index, trace_sample_start, trace_sample_stop)
+    arr = np.fromiter((trace_samples[sample_number - trace_sample_start] for sample_number in sample_numbers),
+                      make_dtype(reader.data_sample_format))
+    return arr
+
+
+def extract_inline_3d(reader_3d, inline_number, xline_numbers=None, sample_numbers=None, null=None):
     """Extract an inline as a two-dimensional array.
 
     Args:
-        reader:
+        reader: A SegYReader3D object.
 
         inline_number: The number of the inline to be extracted.
 
-        xline_numbers: Either a sequence of crossline numbers from which traces are
-            to be extracted or a slice object indicating that some slice of all
-            crossline numbers is to be used.  If None, traces will be extracted at
-            all crosslines.
+        xline_numbers: The crossline numbers within the inline at which traces are to be extracted.
+            This argument can be specified in three ways:
 
-        sample_numbers: An optional sorted sequence of samples numbers at which to
-            extract samples.  If not provided, samples will be extracted at all
-            depths.
+            None (the default) - All traces within the inline will be be extracted.
 
-        null: A null value
+            sequence - When a sequence, such as a range or a list is provided only those traces at
+                crossline numbers corresponding to the items in the sequence will be extracted. The
+                traces will always be extracted in increasing numeric order and duplicate entries
+                will be ignored.  For example xline_numbers=range(100, 200, 2) will extract alternate
+                traces from crossline number 100 to crossline number 198 inclusive.
+
+            slice - When a slice object is provided the slice will be applied to the sequence of all
+                crossline numbers. For example xline_numbers=slice(100, -100) will omit the first
+                one hundred and the last one hundred traces, irrespective of their numbers.
+
+        sample_numbers: The sample numbers within each trace at which samples are to be extracted.
+            This argument can be specified in three ways:
+
+            None (the default) - All samples within the trace will be be extracted.
+
+            sequence - When a sequence, such as a range or a list is provided only those samples at
+                sample numbers corresponding to the items in the sequence will be extracted. The
+                samples will always be extracted in increasing numeric order and duplicate entries
+                will be ignored.  For example sample_numbers=range(100, 200, 2) will extract alternate
+                samples from sample number 100 to sample number 198 inclusive.
+
+            slice - When a slice object is provided the slice will be applied to the sequence of all
+                sample numbers. For example sample_numbers=slice(100, -100) will omit the first
+                one hundred and the last one hundred samples, irrespective of their numbers.
+
+        null: A null value. When None is specified as the null value a masked array will be returned.
 
     Returns:
         A two-dimensional array. If null is None a masked array will be returned, otherwise
-        a simple array will be returned.. The first (slowest changing) index will correspond
+        a regular array will be returned. The first (slowest changing) index will correspond
         to the traces (index zero will correspond to the first crossline number). The
         second (fastest changing) index will correspond to the samples (index zero will
         correspond to the first sample number.
     """
-    if inline_number not in reader.inline_numbers():
-        raise ValueError("Inline number {} not present in {}".format(inline_number, reader))
+    if inline_number not in reader_3d.inline_numbers():
+        raise ValueError("Inline number {} not present in {}".format(inline_number, reader_3d))
 
-    xline_numbers = ensure_superset(reader.xline_numbers(), xline_numbers)
-    sample_numbers = ensure_superset(range(0, reader.max_num_trace_samples()), sample_numbers)
+    xline_numbers = ensure_superset(reader_3d.xline_numbers(), xline_numbers)
+    sample_numbers = ensure_superset(range(0, reader_3d.max_num_trace_samples()), sample_numbers)
     shape = (len(xline_numbers), len(sample_numbers))
-    dtype = make_dtype(reader.data_sample_format)
-    array = make_array(shape, dtype, null)
+    dtype = make_dtype(reader_3d.data_sample_format)
+    array = _make_array(shape, dtype, null)
 
     if isinstance(sample_numbers, range):
-        _populate_inline_array_over_sample_range(reader, inline_number, xline_numbers, sample_numbers, array)
+        _populate_inline_array_over_sample_range(reader_3d, inline_number, xline_numbers, sample_numbers, array)
     else:
-        _populate_inline_array_numbered_samples(reader, inline_number, xline_numbers, sample_numbers, array)
+        _populate_inline_array_numbered_samples(reader_3d, inline_number, xline_numbers, sample_numbers, array)
 
     return array
 
@@ -77,21 +158,10 @@ def _populate_inline_array_over_sample_range(reader, inline_number, xline_number
             array[xline_index, :] = trace_samples[source_slice]
 
 
-def make_array(shape, dtype, null=None):
+def _make_array(shape, dtype, null=None):
     """Make an array"""
     if null is None:
         return np.ma.masked_all(shape, dtype)
-
     array = np.empty(shape, dtype)
     array.fill(null)
     return array
-
-
-def start_and_stop(sequence):
-    """Obtain start and stop values from a sequence."""
-    sample_start = sequence[0]
-    try:
-        sample_stop = sequence.stop
-    except AttributeError:
-        sample_stop = sequence[-1] + 1
-    return sample_start, sample_stop

segpy-ext/segpy-numpy/test/__init__.py

@@ -0,0 +1 @@
+

segpy-ext/segpy-numpy/test/test_extract.py

@@ -0,0 +1,10 @@
+import unittest
+
+
+class MyTestCase(unittest.TestCase):
+    def test_something(self):
+        self.assertEqual(True, False)
+
+
+if __name__ == '__main__':
+    unittest.main()

segpy/util.py

@@ -3,7 +3,9 @@ import time
 import os
 import sys
 
+from collections.abc import Set
 from itertools import (islice, cycle, tee, chain, repeat)
+
 from segpy.sorted_set import SortedFrozenSet
 
 UNKNOWN_FILENAME = '<unknown>'
@@ -407,17 +409,21 @@ def is_range_superset_of_range(superset_range, subset_range):
     return True
 
 
-def is_subset(superset, subset):
-    """A more general version of set.issubset that is smart enough to work with ranges."""
+def is_superset(superset, subset):
+    """A more general version of set.issuperset that is smart enough to work with ranges."""
     if isinstance(subset, range) and isinstance(superset, range):
         return is_range_superset_of_range(superset, subset)
     if isinstance(superset, range):
         return all(item in superset for item in subset)
+    if isinstance(superset, set):
+        return superset.issuperset(subset)
+    if isinstance(subset, set):
+        return subset.issubset(superset)
     return set(superset).issuperset(subset)
 
 
 def ensure_superset(superset, subset):
-    """Obtains a subset of all items.
+    """Ensure that one collection is a subset of another.
 
     Args:
         all_items: A sequence containing all items.
@@ -437,7 +443,7 @@ def ensure_superset(superset, subset):
         return superset[subset]
     else:
         subset = make_sorted_distinct_sequence(subset)
-        if not is_subset(superset, subset):
+        if not is_superset(superset, subset):
             raise ValueError("subset_or_slice {!r} is not a subset of all_items {!r}"
                              .format(subset, superset))
         return subset