Adds code for extraction of an inline.

This commit is contained in:
Robert Smallshire
2015-06-06 20:08:39 +02:00
parent 23ee609c50
commit 6d986a21a9
5 changed files with 216 additions and 2 deletions
+108
View File
@@ -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())
@@ -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
+2 -1
View File
@@ -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,
+4 -1
View File
@@ -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."""
+5
View File
@@ -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])