diff --git a/examples/inline.py b/examples/inline.py new file mode 100644 index 0000000..ac02931 --- /dev/null +++ b/examples/inline.py @@ -0,0 +1,108 @@ +#!/usr/bin/env python3 + +"""Extract an inline from a 3D seismic volume to a 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: inline.py [-h] [--dtype DTYPE] [--null NULL] + segy-file npy-file inline-number + +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 + slice-index 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: + + timeslice.py stack_final_int8.sgy slice_800.npy 800 --null=42.0 --dtype=f +""" + +import argparse +import os +import sys +import traceback + +import numpy as np + +from numpy import s_ + +from segpy.reader import create_reader +from segpy_numpy.dtypes import make_dtype +from segpy_numpy.extract import extract_inline_3d + + +class DimensionalityError(Exception): + pass + + +def extract_inline(segy_filename, out_filename, inline_number, 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. + """ + 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 + + +def nullable_float(s): + return None if not bool(s) else float(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("inline_number", metavar="inline-number", type=int, + help="Zero based index of the inline to be extracted", ) + + parser.add_argument("--null", type=nullable_float, default="", + help="Sample value to use for missing or short traces.") + + if argv is None: + argv = sys.argv[1:] + + args = parser.parse_args(argv) + + try: + extract_inline( + segy_filename=args.segy_file, + out_filename=args.npy_file, + inline_number=args.inline_number, + 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()) + diff --git a/segpy-ext/segpy-numpy/segpy_numpy/extract.py b/segpy-ext/segpy-numpy/segpy_numpy/extract.py new file mode 100644 index 0000000..052f6d8 --- /dev/null +++ b/segpy-ext/segpy-numpy/segpy_numpy/extract.py @@ -0,0 +1,97 @@ +"""Tools for interoperability between Segpy and Numpy arrays.""" +from functools import singledispatch +import numpy as np +from segpy.util import make_sorted_distinct_sequence + +from segpy_numpy.dtypes import make_dtype + + +class DimensionalityError: + pass + + +def extract_inline_3d(reader, inline_number, xline_numbers=None, sample_numbers=None, null=None, packed=False): + """Extract an inline as a two-dimensional array. + + Args: + reader: + + inline_number: The number of the inline to be extracted. + + xline_numbers: An optional sorted sequence of crossline numbers at which to + extract samples. If not provided, samples will be extracted at all + crosslines. + + sample_numbers: An optional sorted sequence of samples numbers at which to + extract samples. If not provided, samples will be extracted at all + depths. + + null: A null value + + 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 + 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 xline_numbers is None: + xline_numbers = reader.xline_numbers() + elif isinstance(xline_numbers, slice): + xline_numbers = reader.xline_numbers()[xline_numbers] + else: + xline_numbers = make_sorted_distinct_sequence(xline_numbers) + + if sample_numbers is None: + sample_numbers = range(0, reader.max_num_trace_samples()) + elif isinstance(sample_numbers, slice): + sample_numbers = range(0, reader.max_num_trace_samples())[sample_numbers] + else: + sample_numbers = make_sorted_distinct_sequence(sample_numbers) + + shape = (len(xline_numbers), len(sample_numbers)) + dtype = make_dtype(reader.data_sample_format) + + if null is None: + array = np.ma.masked_all(shape, dtype) + else: + array = np.empty(shape, dtype) + array.fill(null) + + for xline_index, xline_number in enumerate(xline_numbers): + inline_xline_number = (inline_number, xline_number) + if reader.has_trace_index(inline_xline_number): + # Read the trace + trace_index = reader.trace_index(inline_xline_number) + + sample_start = sample_numbers[0] + try: + sample_stop = sample_numbers.stop + except AttributeError: + sample_stop = sample_numbers[-1] + 1 + + num_trace_samples = reader.num_trace_samples(trace_index) + sample_stop = max(sample_stop, num_trace_samples) + + trace_samples = reader.trace_samples(trace_index, sample_start, sample_stop) + + try: + sample_step = sample_numbers.step + except AttributeError: + sample_step = None + + if sample_step is not None: + # Assign to a slice of the target array + source_slice = slice( + sample_numbers.start - sample_start, + sample_numbers.stop - sample_start, + sample_step) + array[xline_index, :] = trace_samples[source_slice] + else: + # Assign element by element + for sample_index, sample_number in enumerate(sample_numbers): + array[xline_number, sample_index] = trace_samples[sample_number - sample_start] + return array \ No newline at end of file diff --git a/segpy/catalog.py b/segpy/catalog.py index ab0eb08..56459d7 100644 --- a/segpy/catalog.py +++ b/segpy/catalog.py @@ -511,7 +511,8 @@ class LinearRegularCatalog(Mapping): key_max: The maximum key. key_stride: The difference between successive keys. value_start: The value corresponding to the minimum key. - value_max: The value corresponding to the maximum key. + value_stop: The value corresponding to the maximum key. + value_stride: Raises: ValueError: There is any inconsistency in the keys, strides, diff --git a/segpy/reader.py b/segpy/reader.py index 9eee0d8..ce154d9 100644 --- a/segpy/reader.py +++ b/segpy/reader.py @@ -292,6 +292,7 @@ class SegYReader(object): self._revision = extract_revision(self._binary_reel_header) self._bytes_per_sample = bytes_per_sample( self._binary_reel_header, self.revision) + self._max_num_trace_samples = None def __getstate__(self): """Copy the reader's state to a pickleable dictionary. @@ -363,7 +364,9 @@ class SegYReader(object): def max_num_trace_samples(self): """The number of samples in the trace_samples with the most samples.""" - return max(self._trace_length_catalog.values()) + if self._max_num_trace_samples is None: + self._max_num_trace_samples = max(self._trace_length_catalog.values()) + return self._max_num_trace_samples def num_trace_samples(self, trace_index): """The number of samples in the specified trace_samples.""" diff --git a/segpy/util.py b/segpy/util.py index 01dbac2..f30975a 100644 --- a/segpy/util.py +++ b/segpy/util.py @@ -346,6 +346,11 @@ def make_sorted_distinct_sequence(iterable): An immutable collection which supports the Sized, Iterable, Container and Sequence protocols. """ + if isinstance(iterable, range): + if iterable.step > 0: + return iterable + else: + return reversed(iterable) sorted_set = SortedFrozenSet(iterable) if len(sorted_set) == 1: return single_item_range(sorted_set[0])