diff --git a/GUI/demogtk.py b/GUI/demogtk.py
deleted file mode 100755
index cc62fc8..0000000
--- a/GUI/demogtk.py
+++ /dev/null
@@ -1,105 +0,0 @@
-#!/usr/bin/env python
-#Licence: GPLv2.0
-#Copyright: Dave Aitel
-
-import sys
-
-try:
- import pygtk
- #tell pyGTK, if possible, that we want GTKv2
- pygtk.require("2.0")
-except:
- print "You need to install pyGTK or GTKv2 or set your PYTHONPATH correctly"
- print "try: export PYTHONPATH=/usr/local/lib/python2.2/site-packages/"
- sys.exit(1)
-
-import gtk
-import gtk.glade
-import segypy
-
-#now we have both gtk and gtk.glade imported
-#Also, we know we are running GTKv2
-
-def insert_row(model,parent,firstcolumn,secondcolumn):
- myiter=model.insert_after(parent,None)
- model.set_value(myiter,0,firstcolumn)
- model.set_value(myiter,1,secondcolumn)
- return myiter
-
-class appgui:
- def __init__(self):
- """
- In this init we are going to display the main serverinfo window
- """
- gladefile="project1.glade"
- windowname="serverinfo"
- self.wTree=gtk.glade.XML (gladefile,windowname)
- #we only have one callback to register, but you could register
- #any number, or use a special class that
- #automatically registers all callbacks
- dic = { "on_button1_clicked" : self.button1_clicked,
- "on_serverinfo_destroy" : (gtk.mainquit)}
- self.wTree.signal_autoconnect (dic)
-
- self.logwindowview=self.wTree.get_widget("textview1")
- self.logwindow=gtk.TextBuffer(None)
- self.logwindowview.set_buffer(self.logwindow)
-
- import gobject
- self.treeview=self.wTree.get_widget("treeview1")
- self.treemodel=gtk.TreeStore(gobject.TYPE_STRING,
- gobject.TYPE_INT,
- gobject.TYPE_STRING)
- self.treeview.set_model(self.treemodel)
-
- self.treeview.set_headers_visible(True)
-
- renderer=gtk.CellRendererText()
- column=gtk.TreeViewColumn("Name",renderer, text=0)
- column.set_resizable(True)
- self.treeview.append_column(column)
-
- renderer=gtk.CellRendererText()
- column=gtk.TreeViewColumn("Description",renderer, text=1)
- column.set_resizable(True)
- self.treeview.append_column(column)
-
- self.treeview.show()
-
- model=self.treemodel
- treeSTH=insert_row(model,None,'Segy Trace Header', '')
- for key in segypy.TRACE_HEADER_DEF.keys():
- insert_row(model,treeSTH,key, 'An Empty Row')
- treeSH=insert_row(model,None,'Segy Header', 'SH')
- for key in segypy.HEADER_DEF.keys():
- insert_row(model,treeSH,key, 'An Empty Row')
-
-
-
- return
-
-
- #####CALLLBACKS
- def button1_clicked(self,widget):
- print "button clicked"
- host=self.wTree.get_widget("entry1").get_text()
- port=int(self.wTree.get_widget("spinbutton1").get_value())
- if host=="":
- return
- import urllib
- page=urllib.urlopen("http://"+host+":"+str(port)+"/")
- data=page.read()
- try:
- #this changed in a revision of pyGTK2
- self.logwindow.insert_at_cursor(data,len(data))
- except:
- self.logwindow.insert_at_cursor(data)
- #print data
- return
-
-
-app=appgui()
-gtk.main()
-
-
-
diff --git a/GUI/gui.py b/GUI/gui.py
deleted file mode 100644
index f2b6267..0000000
--- a/GUI/gui.py
+++ /dev/null
@@ -1,240 +0,0 @@
-#!/usr/bin/python2.4
-# This program is free software; you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# the Free Software Foundation; either version 2 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with this program; if not, write to the Free Software
-# Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
-#
-# Small test to demonstrate glade.XML.signal_autoconnect on an instance
-#
-
-import pygtk
-pygtk.require('2.0')
-
-
-import matplotlib
-matplotlib.use('GTK')
-from matplotlib.figure import Figure
-from matplotlib.axes import Subplot
-from matplotlib.backends.backend_gtk import FigureCanvasGTK, NavigationToolbar
-from matplotlib.numerix import arange, sin, pi
-
-import Numeric as numpy
-
-import gtk, gtk.glade
-
-import segypy
-
-
-def insert_row(model,parent,firstcolumn,secondcolumn,thirdcolumn):
- myiter=model.insert_after(parent,None)
- model.set_value(myiter,0,firstcolumn)
- model.set_value(myiter,1,secondcolumn)
- model.set_value(myiter,2,thirdcolumn)
-# print secondcolumn
-# return myiter
-
-def update_segyplot(self):
- self.figure = Figure(figsize=(6,4), dpi=72)
- self.axis = self.figure.add_subplot(111)
- self.axis.set_xlabel('Trace Header')
- self.axis.set_ylabel('Time [s]')
- self.axis.set_title(self.segy[1]['filename'])
- self.axis.grid(True)
- self.axis.imshow(self.segy[0])
- self.canvas = FigureCanvasGTK(self.figure) # a gtk.DrawingArea
- self.canvas.show()
-
- self.graphview = self.xml.get_widget("vbox4")
- self.graphview.pack_start(self.canvas, True, True)
-
-def init_treeviews(self):
- import gobject
-
- ###############################################################
- # SEGYHEADER VIEW
- self.treeview1 = self.xml.get_widget("treeview1")
- self.treemodel1 = gtk.TreeStore(gobject.TYPE_STRING,
- gobject.TYPE_STRING,
- gobject.TYPE_STRING)
- self.treeview1.set_model(self.treemodel1)
- self.treeview1.set_headers_visible(True)
-
- # COLUMN 1
- renderer= gtk.CellRendererText()
- self.tvcolumn1 = gtk.TreeViewColumn('Name', renderer, text=0)
- self.tvcolumn1.set_resizable(True)
- self.treeview1.append_column(self.tvcolumn1)
-
- # COLUMN 2
- renderer= gtk.CellRendererText()
- self.tvcolumn2 = gtk.TreeViewColumn('Meaning', renderer, text=1)
- self.tvcolumn2.set_resizable(True)
- self.treeview1.append_column(self.tvcolumn2)
-
- # COLUMN 3
- renderer= gtk.CellRendererText()
- self.tvcolumn3 = gtk.TreeViewColumn('Value', renderer, text=2)
- self.tvcolumn3.set_resizable(True)
- self.treeview1.append_column(self.tvcolumn3)
-
- ###############################################################
- # SEGYHEADER VIEW
-
- self.treeview2 = self.xml.get_widget("treeview2")
- self.treemodel2 = gtk.TreeStore(gobject.TYPE_STRING,
- gobject.TYPE_STRING,
- gobject.TYPE_STRING)
- self.treeview2.set_model(self.treemodel2)
- self.treeview2.set_headers_visible(True)
-
- # COLUMN 1
- renderer= gtk.CellRendererText()
- self.tvcolumn1 = gtk.TreeViewColumn('Name', renderer, text=0)
- self.tvcolumn1.set_resizable(True)
- self.treeview2.append_column(self.tvcolumn1)
-
- # COLUMN 2
- renderer= gtk.CellRendererText()
- self.tvcolumn2 = gtk.TreeViewColumn('Meaning', renderer, text=1)
- self.tvcolumn2.set_resizable(True)
- self.treeview2.append_column(self.tvcolumn2)
-
- # COLUMN 3
- renderer= gtk.CellRendererText()
- self.tvcolumn3 = gtk.TreeViewColumn('Value', renderer, text=2)
- self.tvcolumn3.set_resizable(True)
- self.treeview2.append_column(self.tvcolumn3)
-
-
-def update_segyheader(self):
- import gobject
-
-
- # INSERT INTO TREEMODEL
- self.SHtree= {"init": 1}
- for key in segypy.HEADER_DEF.keys():
- SHkey=segypy.HEADER_DEF[key]
- if (SHkey.has_key('descr')):
- descr = segypy.HEADER_DEF[key]['descr'][0][self.segy[1][key]]
- else:
- descr = ''
- insert_row(self.treemodel1,None,key,descr,self.segy[1][key])
-
- self.treeview1.show()
-
-def update_segytraceheader(self,itrace=1):
- import gobject
-
- # INSERT INTO TREEMODEL
- self.STHtree= {"init": 1}
- for key in segypy.TRACE_HEADER_DEF.keys():
- STHkey=segypy.TRACE_HEADER_DEF[key]
- if (STHkey.has_key('descr')):
- try:
- descr = segypy.TRACE_HEADER_DEF[key]['descr'][0][self.segy[2][key][itrace-1]]
- except:
- descr='Not Defined'
- else:
- descr = ''
- insert_row(self.treemodel2,None,key,descr,self.segy[2][key][itrace-1])
-
- self.treeview2.show()
-
-
-class SimpleTest:
- def __init__(self):
- self.xml = gtk.glade.XML('segygui/segygui.glade')
- self.xml.signal_autoconnect(self)
- init_treeviews(self)
- self.window1 = self.xml.get_widget("window1")
- self.window1.resize(600,400)
-
- def on_cut1_activate(self, button):
- self.hpaned1= self.xml.get_widget("hpaned1")
- self.vbox3= self.xml.get_widget("vbox3")
- pos = self.hpaned1.get_position()
- if (pos==1):
- pos=100
- self.vbox3.visible = 1
- else:
- pos=1
- self.vbox3.visible = 0
- pos = self.hpaned1.set_position(pos)
- print self.vbox3.visible
-
- #width,height = self.vbox3.get_size()
- #print width,height
- # self.window1.resize(400,400)
-
-
- def on_about1_activate(self, button):
- dialog = gtk.AboutDialog()
- dialog.set_name("SegyPY GUI")
- dialog.set_version(segypy.version.__str__())
- dialog.set_copyright("\302\251 Copyright 200x Thomas Mejer Hansen")
- dialog.set_website("http://segymat.sourceforge.net/segypy/")
- authors=["Thomas Mejer Hansen"]
- dialog.set_authors(authors)
- try:
- f=open('LICENSE')
- license = f.read(100000)
- except IOError:
- license = 'Could not read license file....';
- dialog.set_license(license)
- ## Close dialog on user response
- dialog.connect ("response", lambda d, r: d.destroy())
- dialog.show()
-
- def on_open1_activate(self, button):
- dialog = gtk.FileChooserDialog("Open..",
- None,
- gtk.FILE_CHOOSER_ACTION_OPEN,
- (gtk.STOCK_CANCEL, gtk.RESPONSE_CANCEL,
- gtk.STOCK_OPEN, gtk.RESPONSE_OK))
- dialog.set_default_response(gtk.RESPONSE_OK)
-
- filter = gtk.FileFilter()
- filter.set_name("All files")
- filter.add_pattern("*")
-
- dialog.add_filter(filter)
- filter = gtk.FileFilter()
- filter.set_name("SEG-Y files")
- filter.add_pattern("*.segy")
- filter.add_pattern("*.sgy")
- filter.add_pattern("*.seg-y")
- dialog.add_filter(filter)
- filter = gtk.FileFilter()
- filter.set_name("SU files")
- filter.add_pattern("*.su")
- dialog.add_filter(filter)
- response = dialog.run()
- if response == gtk.RESPONSE_OK:
- print dialog.get_filename(), 'selected'
- filename = dialog.get_filename()
- dialog.destroy()
- self.segy = segypy.read_segy(filename)
-
- update_segyheader(self)
- update_segytraceheader(self)
- update_segyplot(self)
-
- def on_new1_activate(self, button):
- self.segy = segypy.read_segy('../data_4byteINT.segy')
- update_segyheader(self)
- update_segytraceheader(self,1)
- update_segyplot(self)
-
-
-test = SimpleTest()
-
-gtk.main()
diff --git a/GUI/project1.glade b/GUI/project1.glade
deleted file mode 100644
index 97f8fcf..0000000
--- a/GUI/project1.glade
+++ /dev/null
@@ -1,192 +0,0 @@
-
-
-
-
-
-
- True
- Server Info
- GTK_WINDOW_TOPLEVEL
- GTK_WIN_POS_NONE
- False
- True
- False
- True
- False
- False
- GDK_WINDOW_TYPE_HINT_NORMAL
- GDK_GRAVITY_NORTH_WEST
- True
-
-
-
-
- True
- False
- 0
-
-
-
- True
- False
- 8
-
-
-
- True
- Host:
- False
- False
- GTK_JUSTIFY_LEFT
- False
- False
- 0.5
- 0.5
- 0
- 0
- PANGO_ELLIPSIZE_NONE
- -1
- False
- 0
-
-
- 0
- False
- False
-
-
-
-
-
- True
- True
- True
- True
- 0
-
- True
- *
- False
-
-
- 0
- True
- True
-
-
-
-
- 0
- True
- True
-
-
-
-
-
- True
- True
- 1
- 0
- False
- GTK_UPDATE_ALWAYS
- False
- False
- 80 1 65535 1 10 10
-
-
- 0
- False
- False
-
-
-
-
-
- True
- True
- GO!
- True
- GTK_RELIEF_NORMAL
- True
-
-
-
- 0
- False
- False
-
-
-
-
-
- True
- True
- GTK_POLICY_ALWAYS
- GTK_POLICY_ALWAYS
- GTK_SHADOW_NONE
- GTK_CORNER_TOP_LEFT
-
-
-
- True
- True
- True
- False
- True
- GTK_JUSTIFY_LEFT
- GTK_WRAP_NONE
- True
- 0
- 0
- 0
- 0
- 0
- 0
-
-
-
-
-
- 0
- True
- False
-
-
-
-
-
- True
- True
- GTK_POLICY_ALWAYS
- GTK_POLICY_ALWAYS
- GTK_SHADOW_NONE
- GTK_CORNER_TOP_LEFT
-
-
-
- 130
- 226
- True
- True
- True
- False
- False
- True
- False
- False
- False
-
-
-
-
- 0
- True
- True
-
-
-
-
-
-
-
diff --git a/GUI/segygui/segygui.glade b/GUI/segygui/segygui.glade
deleted file mode 100644
index 0a7d679..0000000
--- a/GUI/segygui/segygui.glade
+++ /dev/null
@@ -1,497 +0,0 @@
-
-
-
-
-
-
- True
- Segy Viewer
- GTK_WINDOW_TOPLEVEL
- GTK_WIN_POS_CENTER
- False
- True
- False
- stock_chart-reorganize
- True
- False
- False
- GDK_WINDOW_TYPE_HINT_NORMAL
- GDK_GRAVITY_NORTH_WEST
- True
- False
-
-
-
- True
- False
- 0
-
-
-
-
- 0
- False
- False
-
-
-
-
-
- 385
- True
- True
-
-
-
- True
- False
- 0
-
-
-
- True
- True
- True
- True
- GTK_POS_TOP
- False
- False
-
-
-
-
-
-
-
- True
- Plot
- False
- False
- GTK_JUSTIFY_LEFT
- False
- False
- 0.5
- 0.5
- 0
- 0
- PANGO_ELLIPSIZE_NONE
- -1
- False
- 0
-
-
- tab
-
-
-
-
-
- True
- True
- GTK_POLICY_ALWAYS
- GTK_POLICY_ALWAYS
- GTK_SHADOW_IN
- GTK_CORNER_TOP_LEFT
-
-
-
- True
- True
- True
- False
- False
- True
- False
- False
- False
-
-
-
-
- False
- True
-
-
-
-
-
- True
- Header
- False
- False
- GTK_JUSTIFY_LEFT
- False
- False
- 0.5
- 0.5
- 0
- 0
- PANGO_ELLIPSIZE_NONE
- -1
- False
- 0
-
-
- tab
-
-
-
-
-
- True
- False
- 0
-
-
-
- True
- False
- 0
-
-
-
- 15
- True
- GTK_ARROW_LEFT
- GTK_SHADOW_OUT
- 0.5
- 0.5
- 0
- 0
-
-
-
- 0
- True
- True
-
-
-
-
-
- True
- True
- False
- GTK_POS_TOP
- 0
- GTK_UPDATE_DISCONTINUOUS
- False
- 0 0 100 0 0 0
-
-
- 0
- True
- True
-
-
-
-
-
- 15
- True
- GTK_ARROW_RIGHT
- GTK_SHADOW_OUT
- 0.5
- 0.5
- 0
- 0
-
-
-
- 0
- True
- True
-
-
-
-
- 0
- False
- False
-
-
-
-
-
- True
- True
- GTK_POLICY_ALWAYS
- GTK_POLICY_ALWAYS
- GTK_SHADOW_IN
- GTK_CORNER_TOP_LEFT
-
-
-
- True
- True
- True
- False
- False
- True
- False
- False
- False
-
-
-
-
- 0
- True
- True
-
-
-
-
- False
- True
-
-
-
-
-
- True
- Trace
- False
- False
- GTK_JUSTIFY_LEFT
- False
- False
- 0.5
- 0.5
- 0
- 0
- PANGO_ELLIPSIZE_NONE
- -1
- False
- 0
-
-
- tab
-
-
-
-
- 0
- True
- True
-
-
-
-
- True
- False
-
-
-
-
-
- True
- False
- 0
-
-
-
-
-
-
-
-
-
-
-
-
-
-
- True
- True
-
-
-
-
- 0
- True
- True
-
-
-
-
-
- True
- True
-
-
- 0
- False
- False
-
-
-
-
-
-
-
diff --git a/GUI/test/test.glade b/GUI/test/test.glade
deleted file mode 100644
index 71c67da..0000000
--- a/GUI/test/test.glade
+++ /dev/null
@@ -1,251 +0,0 @@
-
-
-
-
-
-
- 800
- 600
- True
- SegyPY GUI
- GTK_WINDOW_TOPLEVEL
- GTK_WIN_POS_NONE
- False
- True
- False
- True
- False
- False
- GDK_WINDOW_TYPE_HINT_NORMAL
- GDK_GRAVITY_NORTH_WEST
- True
- False
-
-
-
- 307
- True
- True
- True
- 3
- 1
- False
- 0
- 0
-
-
-
-
- 0
- 1
- 0
- 1
- fill
-
-
-
-
-
-
- True
- GTK_ORIENTATION_HORIZONTAL
- GTK_TOOLBAR_BOTH
- True
- True
-
-
-
-
-
-
-
-
-
-
-
-
-
-
- 0
- 1
- 1
- 2
- fill
-
-
-
-
-
-
- True
- False
- 0
-
-
-
-
-
-
-
-
-
-
- 0
- 1
- 2
- 3
- fill
-
-
-
-
-
-
-
diff --git a/GUI/test/test.glade.bak b/GUI/test/test.glade.bak
deleted file mode 100644
index 7ba9dd7..0000000
--- a/GUI/test/test.glade.bak
+++ /dev/null
@@ -1,250 +0,0 @@
-
-
-
-
-
-
- 800
- 600
- True
- SegyPY GUI
- GTK_WINDOW_TOPLEVEL
- GTK_WIN_POS_NONE
- False
- True
- False
- True
- False
- False
- GDK_WINDOW_TYPE_HINT_NORMAL
- GDK_GRAVITY_NORTH_WEST
- True
- False
-
-
-
- True
- True
- True
- 3
- 1
- False
- 0
- 0
-
-
-
-
- 0
- 1
- 0
- 1
- fill
-
-
-
-
-
-
- True
- GTK_ORIENTATION_HORIZONTAL
- GTK_TOOLBAR_BOTH
- True
- True
-
-
-
-
-
-
-
-
-
-
-
-
-
-
- 0
- 1
- 1
- 2
- fill
-
-
-
-
-
-
- True
- False
- 0
-
-
-
-
-
-
-
-
-
-
- 0
- 1
- 2
- 3
- fill
-
-
-
-
-
-
-
diff --git a/GUI/test/test.gladep b/GUI/test/test.gladep
deleted file mode 100644
index 0d232b1..0000000
--- a/GUI/test/test.gladep
+++ /dev/null
@@ -1,8 +0,0 @@
-
-
-
-
- Test
- test
- FALSE
-
diff --git a/GUI/test/test.gladep.bak b/GUI/test/test.gladep.bak
deleted file mode 100644
index 0d232b1..0000000
--- a/GUI/test/test.gladep.bak
+++ /dev/null
@@ -1,8 +0,0 @@
-
-
-
-
- Test
- test
- FALSE
-
diff --git a/catalog.py b/catalog.py
new file mode 100644
index 0000000..16ec6f8
--- /dev/null
+++ b/catalog.py
@@ -0,0 +1,357 @@
+from collections import Mapping, Sequence, OrderedDict
+from fractions import Fraction
+import repr
+from util import contains_duplicates, measure_stride, minmax
+
+
+class CatalogBuilder:
+ """Use a catalog builder to construct optimised, immutable mappings.
+
+ A CatalogBuilder is useful when, depending on the particular keys and
+ values used, a more compact or efficient representation of the mapping
+ is possible than, say, a regular dictionary. The CatalogBuilder
+ accumulates values and then, once all values have been added, analyzes
+ the keys and values to produce an more optimized representation of the
+ mapping.
+ """
+
+ def __init__(self, mapping=None):
+ """Initialize a Catalog Builder.
+
+ Args:
+ mapping: An optional mapping (such as a dictionary) of items.
+ """
+ self._catalog = []
+ if mapping is not None:
+ for key, value in mapping.items():
+ self.add(key, value)
+
+ def add(self, index, value):
+ """Add an item.
+
+ Each index must be unique if create() is to be subsequently called successfully,
+ although duplicate index values will be accepted by this call without complaint.
+ """
+ self._catalog.append((index, value))
+
+ def create(self):
+ """Create a possibly more optimized representation of the mapping.
+
+ In this worst case, this method returns an object which is essentially an immutable dictionary. In the best
+ case, the space savings can be vast.
+
+ Returns:
+ A mapping, if a unique mapping from indexes to values is possible, otherwise None.
+ """
+
+ # This method examines the contents of the mapping using various heuristics to come up with
+ # a better representation.
+
+ if len(self._catalog) < 2:
+ return DictionaryCatalog(self._catalog)
+
+ # In-place sort by index
+ self._catalog.sort(key=lambda index_value: index_value[0])
+
+ if contains_duplicates(index for index, value in self._catalog):
+ return None
+
+ if all(isinstance(index, Sequence) and (len(index) == 2) for index, value in self._catalog):
+ return self._create_catalog_2()
+
+ return self._create_catalog_1()
+
+ def _create_catalog_1(self):
+ """Create a catalog for one-dimensional integer keys (i.e. scalars)
+ """
+ index_min = self._catalog[0][0]
+ index_max = self._catalog[-1][0]
+ index_stride = measure_stride(index for index, value in self._catalog)
+
+ if index_stride is None:
+ # Dictionary strategy - arbitrary keys and values
+ return DictionaryCatalog(self._catalog)
+
+ self._catalog.sort(key=lambda index_value: index_value[1])
+ value_min = self._catalog[0][1]
+ value_max = self._catalog[-1][1]
+ value_stride = measure_stride(value for index, value in self._catalog)
+
+ if index_stride is not None and value_stride is None:
+ # Regular index - regular keys and arbitrary values
+ return RegularCatalog(index_min, index_max, index_stride, (value for index, value in self._catalog))
+
+ assert (index_stride is not None) and (value_stride is not None)
+ catalog = LinearRegularCatalog(index_min, index_max, index_stride, value_min, value_max, value_stride)
+ return catalog
+
+ def _create_catalog_2(self):
+ """Create a catalog for two-dimensional integer keys.
+
+ Each key must be a two-element sequence.
+ """
+ i_min, i_max = minmax(i for (i, j), value in self._catalog)
+ j_min, j_max = minmax(j for (i, j), value in self._catalog)
+
+ is_rm, diff = self._is_row_major(i_min, j_min, j_max)
+ if is_rm:
+ return RowMajorCatalog(i_min, i_max, j_min, j_max, diff)
+
+ def _is_row_major(self, i_min, j_min, j_max):
+ """Does row major ordering predict values from keys?
+
+ Args:
+ i_min: The minimum i value.
+ j_min: The minimum j value.
+ j_max: The maximum j value.
+
+ Returns:
+ A 2-tuple containing, in the first element True if the values can
+ be predicted from the keys by assuming a row-major ordering,
+ otherwise False. If True, the second element will be a constant
+ offset, otherwise it can be ignored.
+ ordering
+ """
+ diff = None
+ for (i, j), actual_value in self._catalog:
+ proposed_value = (i - i_min) * j_max + (j - j_min)
+ current_diff = actual_value - proposed_value
+ if diff is None:
+ diff = current_diff
+ if current_diff != diff:
+ return False, None
+ return True, diff
+
+ return DictionaryCatalog(self._catalog)
+
+
+class RowMajorCatalog(Mapping):
+ """A mapping which assumes a row-major ordering of a two-dimensional matrix.
+
+ This is the ordering of items in a two-dimensional matrix where in the (i, j)
+ key tuple the j value changes fastest when iterating through the items in order.
+
+ A RowMajorCatalog predicts the value v from the key (i, j) according to the
+ following formula:
+
+ v = (i - i_min) * j_max + (j - j_min) + c
+
+ for
+ i_min <= i <= i_max
+ j_min <= j <= j_max
+
+ and where c is an integer constant to allow zero- or one-based indexing.
+ """
+
+ def __init__(self, i_min, i_max, j_min, j_max, c):
+ """Initialize a RowMajorCatalog.
+
+ Args:
+ i_min (int): The minimum i value.
+ i_max (int): The maximum i value.
+ j_min (int): The minimum j value.
+ j_max (int): The maximum j value.
+ c (int): The constant offset
+ """
+ self._i_min = i_min
+ self._i_max = i_max
+ self._j_min = j_min
+ self._j_max = j_max
+ self._c = c
+
+ @property
+ def i_min(self):
+ """Minimum i value"""
+ return self._i_min
+
+ @property
+ def i_max(self):
+ """Maximum i value"""
+ return self._i_max
+
+ @property
+ def j_min(self):
+ """Minimum j value"""
+ return self._j_min
+
+ @property
+ def j_max(self):
+ """Maximum j value"""
+ return self._j_max
+
+ def min(self):
+ """Minimum (i, j) value"""
+ return self._i_min, self._j_min
+
+ def max(self):
+ """Maximum (i, j) value"""
+ return self._j_min, self._j_max
+
+ def __getitem__(self, key):
+ i, j = key
+ if not (self._i_min <= i <= self._i_max) and (self._j_min <= j <= self._j_max):
+ raise KeyError("{!r} key {!r} out of range".format(self, key))
+ value = (i - self._i_min) * self._j_max + (j - self._j_min) + self._c
+ return value
+
+ def __contains__(self, key):
+ i, j = key
+ return (self._i_min <= i <= self._i_max) and (self._j_min <= j <= self._j_max)
+
+ def __len__(self):
+ return (self._i_max - self._i_min) * (self._j_max - self._j_min)
+
+ def __iter__(self):
+ for i in range(self._i_min, self._i_max + 1):
+ for j in range(self._j_min, self._j_max + 1):
+ yield (i, j)
+
+ def __repr__(self):
+ return '{}({}, {}, {}, {}, {})'.format(self.__class__.__name__,
+ self._i_min, self._i_max, self._j_min, self._j_max, self._c)
+
+
+class DictionaryCatalog(Mapping):
+ """An immutable, ordered, dictionary mapping.
+ """
+
+ def __init__(self, items):
+ self._items = OrderedDict(items)
+
+ def __getitem__(self, key):
+ return self._items[key]
+
+ def __iter__(self):
+ return iter(self._items)
+
+ def __len__(self):
+ return len(self._items)
+
+ def __contains__(self, item):
+ return item in self._items
+
+ def __repr__(self):
+ return '{}({})'.format(self.__class__.__name__, repr.repr(self._items.items()))
+
+
+class RegularCatalog(Mapping):
+ """Mapping with keys ordered with regular spacing along the number line.
+
+ The values associated with the keys are arbitrary.
+ """
+
+ def __init__(self, key_min, key_max, key_stride, values):
+ """Initialize a RegularCatalog.
+
+ The catalog is initialized by provided a description of how the keys
+ along the number line, and an iterable series of corresponding values.
+
+ Args:
+ key_min: The minimum key.
+ key_max: The maximum key.
+ key_stride: The difference between successive keys.
+ values: An iterable series of values corresponding to the keys.
+ """
+ key_range = key_max - key_min
+ if key_range % key_stride != 0:
+ raise ("RegularIndex key range {!r} is not a multiple of stride {!r}".format(key_stride, key_range))
+ self._key_min = key_min
+ self._key_max = key_max
+ self._key_stride = key_stride
+ self._values = list(values)
+ num_keys = key_range // key_stride
+ if num_keys != len(self._values):
+ raise ValueError("RegularIndex key range and values inconsistent")
+
+ def __getitem__(self, key):
+ if not (self._key_min <= key <= self._key_max):
+ raise KeyError("{!r} key {!r} out of range".format(self, key))
+ offset = key - self._key_min
+ if offset % self._key_stride != 0:
+ raise KeyError("{!r} does not contain key {!r}".format(self, key))
+ index = offset // self._key_stride
+ return self._values[index]
+
+ def __len__(self):
+ return len(self._values)
+
+ def __contains__(self, key):
+ return (self._key_min <= key <= self._key_max) and ((key - self._key_min) % self._key_stride == 0)
+
+ def __iter__(self):
+ return iter(range(self._key_min, self._key_max + 1, self._index_stride))
+
+ def __repr__(self):
+ return '{}({}, {}, {}, {})'.format(self.__class__.__name, self._key_min, self._key_max, self._key_stride,
+ repr.repr(self._values))
+
+
+class LinearRegularCatalog(Mapping):
+ """A mapping which assumes a linear relationship between keys and values.
+
+ This is the ordering of items in a two-dimensional matrix where in the (i, j)
+ key tuple the j value changes fastest when iterating through the items in order.
+
+ A LinearRegularCatalog predicts the value v from the key according to the
+ following formula:
+
+ v = (value_max - value_min) / (key_max - key_min) * (key - key_min) + value_min
+ """
+
+ def __init__(self, key_min, key_max, key_stride, value_min, value_max, value_stride):
+ """Initialize a LinearRegularCatalog.
+
+ Args:
+ key_min: The minimum key.
+ key_max: The maximum key.
+ key_stride: The difference between successive keys.
+ value_min: The minimum value.
+ value_max: The maximum value.
+ """
+ key_range = key_max - key_min
+ if key_range % key_stride != 0:
+ raise ("{} key range {!r} is not a multiple of key stride {!r}".format(
+ self.__class__.__name__, key_stride, key_range))
+ self._key_min = key_min
+ self._key_max = key_max
+ self._key_stride = key_stride
+
+ value_range = value_max - value_min
+ if value_range % value_stride != 0:
+ raise ("{} value range {!r} is not a multiple of value stride {!r}".format(
+ self.__class__.__name__, value_stride, value_range))
+ self._value_min = value_min
+ self._value_max = value_max
+ self._value_stride = value_stride
+
+ num_keys = (self._key_max - self._key_min) // self._key_stride
+ num_values = (self._value_max - self._value_min) // self._value_stride
+ if num_keys != num_values:
+ raise ("{} inconsistent number of keys {} and values {}".format(
+ self.__class__.__name__, num_keys, num_values))
+
+ self._m = Fraction(self._value_max - self._value_min, self._key_max - self._key_min)
+
+ def __getitem__(self, key):
+ if not (self._key_min <= key <= self._key_max):
+ raise KeyError("{!r} key {!r} out of range".format(self, key))
+ offset = key - self._key_min
+ if offset % self._key_stride != 0:
+ raise KeyError("{!r} does not contain key {!r}".format(self, key))
+
+ v = self._m * (key - self._key_min) + self._value_min
+ assert v.denominator == 1
+ return v.numerator
+
+ def __len__(self):
+ return 1 + (self._key_max - self._key_min) // self._key_stride
+
+ def __contains__(self, key):
+ return (self._key_min <= key <= self._key_max) and ((key - self._key_min) % self._key_stride == 0)
+
+ def __iter__(self):
+ return iter(range(self._key_min, self._key_max + 1, self._key_stride))
+
+ def __repr__(self):
+ return '{}({}, {}, {}, {}, {}, {})'.format(self.__class__.__name__,
+ self._key_min, self._key_max, self._key_stride, self._value_min, self._value_max, self._value_stride)
\ No newline at end of file
diff --git a/datatypes.py b/datatypes.py
new file mode 100644
index 0000000..ad94441
--- /dev/null
+++ b/datatypes.py
@@ -0,0 +1,59 @@
+DATA_SAMPLE_FORMAT = {1: 'ibm',
+ 2: 'l',
+ 3: 'h',
+ 5: 'f',
+ 8: 'b'}
+
+# A mapping from SEG Y data types to format characters used by the struct module,
+# known a 'ctypes'
+CTYPES = {'l': 'l',
+ 'long': 'l',
+ 'int32': 'l',
+
+ 'L': 'L',
+ 'ulong': 'L',
+ 'uint32': 'L',
+
+ 'h': 'h',
+ 'short': 'h',
+ 'int16': 'h',
+
+ 'H': 'H',
+ 'ushort': 'H',
+ 'uint16': 'H',
+
+ 'c': 'b',
+ 'char': 'b',
+ 'b': 'b',
+
+ 'B': 'B',
+ 'uchar': 'B',
+
+ 'f': 'f',
+ 'float': 'f',
+
+ 'ibm': 'ibm'}
+
+# TODO This is redundant with data in the SH_def below
+CTYPE_DESCRIPTION = {'ibm': 'IBM float',
+ 'l': '32 bit signed integer',
+ 'L': '32 bit unsigned integer',
+ 'h': '16 bit signed integer',
+ 'H': '16 bit unsigned integer',
+ 'f': 'IEEE float32',
+ 'b': '8 bit signed char',
+ 'B': '8 bit unsigned char'}
+
+
+SIZES = dict(l=4,
+ L=4,
+ h=2,
+ H=2,
+ b=1,
+ B=1,
+ f=4,
+ ibm=4)
+
+
+def size_in_bytes(ctype):
+ return SIZES[ctype]
diff --git a/ibm_float.py b/ibm_float.py
index 3fa7d7b..828c1cd 100644
--- a/ibm_float.py
+++ b/ibm_float.py
@@ -1,22 +1,49 @@
-import struct
+import sys
+_P24 = float(pow(2, 24))
-def ibm2ieee2(ibm_float):
- """
- ibm2ieee2(ibm_float)
- Used by permission
- (C) Secchi Angelo
- with thanks to Howard Lightstone and Anton Vredegoor.
- """
- dividend = float(16 ** 6)
+if sys.version_info >= (3, 0):
- if ibm_float == 0:
- return 0.0
- istic, a, b, c = struct.unpack('>BBBB', ibm_float)
- if istic >= 128:
- sign = -1.0
- istic -= 128
- else:
- sign = 1.0
- mant = float(a << 16) + float(b << 8) + float(c)
- return sign * 16 ** (istic - 64) * (mant / dividend)
+ def ibm2ieee(big_endian_bytes):
+ """Interpret a bytes object as a big-endian IBM float.
+
+ Args:
+ big_endian_bytes (bytes): A string containing at least four bytes.
+
+ Returns:
+ The floating point value.
+ """
+ a, b, c, d = big_endian_bytes
+
+ if a == b == c == c == 0:
+ return 0.0
+
+ sign = 1 if (a & 0x80) else -1
+ exponent = a & 0x7f
+ mantissa = ((b << 16) | (c << 8) | d) / _P24
+ value = sign * mantissa * pow(16, exponent - 64)
+ return value
+else:
+
+ def ibm2ieee(big_endian_bytes):
+ """Interpret a byte string as a big-endian IBM float.
+
+ Args:
+ big_endian_bytes (str): A string containing at least four bytes.
+
+ Returns:
+ The floating point value.
+ """
+ a = ord(big_endian_bytes[0])
+ b = ord(big_endian_bytes[1])
+ c = ord(big_endian_bytes[2])
+ d = ord(big_endian_bytes[3])
+
+ if a == b == c == c == 0:
+ return 0.0
+
+ sign = 1 if (a & 0x80) else -1
+ exponent = a & 0x7f
+ mantissa = ((b << 16) | (c << 8) | d) / _P24
+ value = sign * mantissa * pow(16, exponent - 64)
+ return value
\ No newline at end of file
diff --git a/makerelease.sh b/makerelease.sh
deleted file mode 100755
index dcf37d9..0000000
--- a/makerelease.sh
+++ /dev/null
@@ -1,15 +0,0 @@
-#!/bin/sh
-VERSION=0.3
-
-DIR=segypy-$VERSION
-
-mkdir $DIR
-cp segypy.py $DIR/.
-cp testsegy.py $DIR/.
-cp LICENSE $DIR/.
-cp README $DIR/.
-
-zip -r $DIR.zip $DIR
-tar cvfz $DIR.tgz $DIR
-
-rm -fr $DIR
\ No newline at end of file
diff --git a/plotting.py b/plotting.py
deleted file mode 100644
index ba5b2bf..0000000
--- a/plotting.py
+++ /dev/null
@@ -1,34 +0,0 @@
-def image_segy(Data):
- """
- imageSegy(Data)
- Image segy Data
- """
- import pylab
- pylab.imshow(Data)
- pylab.title('pymat test')
- pylab.grid(True)
- pylab.show()
-
-
-def wiggle(Data, SH, skipt=1, maxval=8, lwidth=0.1):
- """
- wiggle(Data, SH)
- """
- import pylab
-
- t = range(SH['ns'])
-
- for i in range(0, SH['ntraces'], skipt):
-
- trace = Data[:, i]
- trace[0] = 0
- trace[SH['ns'] - 1] = 0
- pylab.plot(i + trace / maxval, t, color='black', linewidth=lwidth)
- for a in range(len(trace)):
- if trace[a] < 0:
- trace[a] = 0
-
- pylab.fill(i + Data[:, i] / maxval, t, 'k', linewidth=0)
- pylab.title(SH['filename'])
- pylab.grid(True)
- pylab.show()
diff --git a/portability.py b/portability.py
new file mode 100644
index 0000000..83ca77f
--- /dev/null
+++ b/portability.py
@@ -0,0 +1,25 @@
+import os
+from util import file_length
+
+
+def seekable(fh):
+ """Determine whether a file-like object supports seeking.
+
+ Args:
+ fh: The file-like-object to be tested.
+
+ Returns:
+ True if the file supports seeking, otherwise False.
+ """
+ try:
+ return fh.seekable()
+ except AttributeError:
+ try:
+ pos = fh.tell()
+ try:
+ fh.seek(0, os.SEEK_END)
+ finally:
+ fh.seek(pos)
+ except AttributeError:
+ return False
+ return True
\ No newline at end of file
diff --git a/reader.py b/reader.py
new file mode 100644
index 0000000..88ab259
--- /dev/null
+++ b/reader.py
@@ -0,0 +1,383 @@
+from __future__ import print_function
+
+from portability import seekable
+from util import file_length, filename_from_handle
+from datatypes import DATA_SAMPLE_FORMAT, CTYPE_DESCRIPTION, CTYPES
+from toolkit import (extract_revision, bytes_per_sample, read_reel_header, catalog_traces, read_binary_values,
+ compile_trace_header_format, TraceHeader, REEL_HEADER_NUM_BYTES, TRACE_HEADER_NUM_BYTES)
+
+
+def create_reader(fh, endian='>'):
+ """Create a SegYReader (or one of its subclasses) based on performing a scan of SEG Y data.
+
+ This function is the preferred method for creating SegYReader objects. It reads basic header
+ information and attempts to build indexes for traces, CDP numbers (for 2D surveys), and inline
+ and cross line co-ordinates (for 3D surveys) to facilitate subsequent random-access to traces.
+
+ Args:
+ fh: A file-like-object open in binary mode positioned such that the
+ beginning of the reel header will be the next byte to be read. For disk-based
+ SEG Y files, this is the beginning of the file.
+
+ endian: '>' for big-endian data (the standard and default), '<' for
+ little-endian (non-standard)
+
+ Returns:
+ A SegYReader object. Depending on the exact type of the SegYReader returned different capabilities may be
+ available. Inspect the returned object to determine these capabilities, or be prepared for capabilities not
+ defined in the SegYReader base class to be unavailable. The underlying file-like object must remain open
+ for the duration of use of the returned reader object. It is the caller's responsibility to close the
+ underlying file.
+
+ Example:
+
+ with open('my_seismic_data.sgy', 'rb') as fh:
+ reader = create_reader(fh)
+ print(reader.num_traces())
+
+ """
+ if fh.encoding is not None:
+ raise TypeError("SegYReader must be provided with a binary mode file object")
+
+ if not seekable(fh):
+ raise TypeError("SegYReader must be provided with a seekable file object")
+
+ if fh.closed:
+ raise ValueError("SegYReader must be provided with an open file object")
+
+ num_file_bytes = file_length(fh)
+ if num_file_bytes < REEL_HEADER_NUM_BYTES:
+ raise ValueError("SEG Y file {!r} of {} bytes is too short".format(filename_from_handle(fh),
+ num_file_bytes))
+ if endian not in ('<', '>'):
+ raise ValueError("Unrecognised endian value {!r}".format(endian))
+
+ reel_header = read_reel_header(fh, endian)
+ revision = extract_revision(reel_header)
+ bps = bytes_per_sample(reel_header, revision)
+
+ trace_catalog, cdp_catalog, line_catalog = catalog_traces(fh, bps, endian)
+
+ if cdp_catalog is not None and line_catalog is None:
+ return SegYReader2D(fh, reel_header, trace_catalog, cdp_catalog, endian)
+
+ if cdp_catalog is None and line_catalog is not None:
+ return SegYReader3D(fh, reel_header, trace_catalog, line_catalog, endian)
+
+ return SegYReader(fh, reel_header, trace_catalog, endian)
+
+
+class SegYReader(object):
+ """A basic SEG Y reader.
+
+ Use to obtain read the reel header, the trace headers or trace values. Traces can be accessed
+ only by trace index.
+ """
+
+ def __init__(self, fh, reel_header, trace_catalog, endian='>'):
+ """Initialize a SegYReader around a file-like-object.
+
+ Note:
+ Usually a SegYReader is most easily constructed using the create_reader() function.
+
+ Args:
+ fh: A file-like object, which must support seeking and support binary reading.
+
+ reel_header: A dictionary containing reel header data.
+
+ trace_catalog: A mapping from zero-based trace index to the byte-offset to
+ individual traces within the file.
+
+ endian: '>' for big-endian data (the standard and default), '<' for
+ little-endian (non-standard)
+ """
+ self._fh = fh
+ self._endian = endian
+ self._trace_header_format = compile_trace_header_format(self._endian)
+
+ self._reel_header = reel_header
+ self._trace_catalog = trace_catalog
+
+ self._revision = extract_revision(self._reel_header)
+ self._bytes_per_sample = bytes_per_sample(self._reel_header, self.revision)
+
+ def trace_indexes(self):
+ """An iterator over zero-based trace indexes.
+
+ Returns:
+ An iterator which yields integers in the range zero to num_traces() - 1
+ """
+ return iter(self._trace_catalog)
+
+ def num_traces(self):
+ """The number of traces"""
+ return len(self._trace_catalog)
+
+ def read_trace(self, trace_index):
+ """Read a specific trace.
+
+ Args:
+ trace_index: An integer in the range zero to num_traces() - 1
+
+ Returns:
+ A 2-tuple containing a TraceHeader as the first item and a sequence of numeric trace samples
+ as the second item.
+
+ Example:
+
+ first_trace_header, first_trace_samples = segy_reader.read_trace(0)
+ """
+ if not (0 <= trace_index < self.num_traces()):
+ raise ValueError("Trace index out of range.")
+ trace_header = self.read_trace_header(trace_index)
+ num_samples = trace_header.ns
+ dsf = self._reel_header['DataSampleFormat']
+ ctype = DATA_SAMPLE_FORMAT[dsf]
+ pos = self._trace_catalog[trace_index] + TRACE_HEADER_NUM_BYTES
+ trace_values = read_binary_values(self._fh, pos, ctype, num_samples, self._endian)
+ return trace_header, trace_values
+
+ def read_trace_header(self, trace_index):
+ """Read a specific trace.
+
+ Args:
+ trace_index: An integer in the range zero to num_traces() - 1
+
+ Returns:
+ A TraceHeader corresponding to the requested trace.
+
+ Example:
+
+ first_trace_header, first_trace_samples = segy_reader.read_trace(0)
+ """
+ if not (0 <= trace_index < self.num_traces()):
+ raise ValueError("Trace index out of range.")
+ pos = self._trace_catalog[trace_index]
+ self._fh.seek(pos)
+ data = self._fh.read(TRACE_HEADER_NUM_BYTES)
+ trace_header = TraceHeader._make(self._trace_header_format.unpack(data))
+ return trace_header
+
+ @property
+ def dimensionality(self):
+ """The spatial dimensionality of the data.
+
+ Returns:
+ 3 for 3D seismic volumes, 2 for 2D seismic lines, 1 for a single trace, otherwise 0.
+ """
+ return self._dimensionality()
+
+ def _dimensionality(self):
+ return 1 if self.num_traces() == 1 else 0
+
+ @property
+ def reel_header(self):
+ """The reel header, sometimes known as the binary header.
+
+ Returns:
+ A dictionary containing data from the reel header.
+ """
+ return self._reel_header
+
+ @property
+ def filename(self):
+ """The filename.
+
+ Returns:
+ The filename if it could be determined, otherwise ''
+ """
+ return filename_from_handle(self._fh)
+
+ @property
+ def revision(self):
+ """The SEG Y revision.
+
+ Returns:
+ Either datatypes.SEGY_REVISION_0 or datatypes.SEGY_REVISION_1
+ """
+ return self._revision
+
+ @property
+ def bytes_per_sample(self):
+ """The number of bytes per trace sample.
+ """
+ return self._bytes_per_sample
+
+ @property
+ def data_sample_format(self):
+ """The data type of the samples in machine-readable form.
+
+ Returns:
+ One of the values from datatypes.DATA_SAMPLE_FORMAT
+ """
+ return DATA_SAMPLE_FORMAT[self._reel_header['DataSampleFormat']]
+
+ @property
+ def data_sample_format_description(self):
+ """A descriptive human-readable description of the data sample format
+ """
+ return CTYPE_DESCRIPTION[CTYPES[self.data_sample_format]]
+
+
+class SegYReader3D(SegYReader):
+ """A reader for 3D seismic data.
+
+ In addition to the capabilities provided by the SegYReader base class, this reader
+ provides an index to facilitate random access to individual traces via crossline and
+ inline co-ordinates.
+ """
+
+ def __init__(self, fh, reel_header, trace_catalog, line_catalog, endian='>'):
+ """Initialize a SegYReader3D around a file-like-object.
+
+ Note:
+ Usually a SegYReader is most easily constructed using the create_reader() function.
+
+ Args:
+ fh: A file-like object, which must support seeking and support binary reading.
+
+ reel_header: A dictionary containing reel header data.
+
+ trace_catalog: A mapping from zero-based trace indexes to the byte-offset to
+ individual traces within the file.
+
+ line_catalog: A mapping from (xline, inline) tuples to trace_indexes.
+
+ endian: '>' for big-endian data (the standard and default), '<' for
+ little-endian (non-standard)
+ """
+ super(SegYReader3D, self).__init__(fh, reel_header, trace_catalog, endian)
+ self._line_catalog = line_catalog
+
+ def _dimensionality(self):
+ return 3
+
+ def num_inlines(self):
+ """The number of distinct inlines in the survey
+ """
+ try:
+ return self._line_catalog.j_max - self._line_catalog.j_min
+ except AttributeError:
+ # TODO: Memoize
+ return len(set(j for i, j in self._line_catalog))
+
+ def num_xlines(self):
+ """The number of distinct crosslines in the survey
+ """
+ try:
+ return self._line_catalog.i_max - self._line_catalog.i_min
+ except AttributeError:
+ # TODO: Memoize
+ return len(set(i for i, j in self._line_catalog))
+
+ def inline_xlines(self):
+ """An iterator over all (xline_number, inline_number) tuples corresponding to traces.
+ """
+ return iter(self._line_catalog)
+
+ def trace_index(self, xline, inline):
+ """Obtain the trace index given an xline and a inline.
+
+ Note:
+ Do not assume that all combinations of crossline and inline co-ordinates are valid.
+ The volume may not be rectangular. Valid values can be obtained from the
+ inline_xlines() iterator.
+
+ Furthermore, inline and crossline numbers should not be relied upon to be zero- or one-based
+ indexes (although they may be).
+
+ Args:
+ xline: A crossline number.
+ inline: An inline number.
+
+ Returns:
+ A trace index which can be used with read_trace().
+ """
+ return self._line_catalog[(xline, inline)]
+
+
+class SegYReader2D(SegYReader):
+
+ def __init__(self, fh, reel_header, trace_catalog, cdp_catalog, endian='>'):
+ """Initialize a SegYReader2D around a file-like-object.
+
+ Note:
+ Usually a SegYReader is most easily constructed using the create_reader() function.
+
+ Args:
+ fh: A file-like object, which must support seeking and support binary reading.
+
+ reel_header: A dictionary containing reel header data.
+
+ trace_catalog: A mapping from zero-based trace index to the byte-offset to
+ individual traces within the file.
+
+ cdp_catalog: A mapping from CDP numbers to trace_indexes.
+
+ endian: '>' for big-endian data (the standard and default), '<' for
+ little-endian (non-standard)
+ """
+ super(SegYReader2D, self).__init__(fh, reel_header, trace_catalog, endian)
+ self._cdp_catalog = cdp_catalog
+
+ def _dimensionality(self):
+ return 2
+
+ def cdp_numbers(self):
+ """An iterator over all cdp numbers corresponding to traces.
+ """
+ return iter(self._cdp_catalog)
+
+ def num_cdps(self):
+ return len(self._cdp_catalog)
+
+ def trace_index(self, cdp_number):
+ """Obtain the trace index given an xline and a inline.
+
+ Note:
+ Do not assume that all combinations of crossline and inline co-ordinates are valid.
+ The volume may not be rectangular. Valid values can be obtained from the
+ inline_xlines() iterator.
+
+ Furthermore, inline and crossline numbers should not be relied upon to be zero- or one-based
+ indexes (although they may be).
+
+ Args:
+ xline: A crossline number.
+ inline: An inline number.
+
+ Returns:
+ A trace index which can be used with read_trace().
+ """
+ return self._cdp_catalog[cdp_number]
+
+
+def main(argv=None):
+ import sys
+ if argv is None:
+ argv = sys.argv[1:]
+
+ filename = argv[0]
+
+ with open(filename, 'rb') as segy_file:
+ segy_reader = create_reader(segy_file)
+ print("Filename: ", segy_reader.filename)
+ print("SEG Y revision: ", segy_reader.revision)
+ print("Number of traces: ", segy_reader.num_traces())
+ print("Data format: ", segy_reader.data_sample_format_description)
+ print("Dimensionality: ", segy_reader.dimensionality)
+
+ try:
+ print("Number of CDPs: ", segy_reader.num_cdps())
+ except AttributeError:
+ pass
+
+ try:
+ print("Number of inlines: ", segy_reader.num_inlines())
+ print("Number of crosslines: ", segy_reader.num_xlines())
+ except AttributeError:
+ pass
+
+if __name__ == '__main__':
+ main()
+
diff --git a/header_definition.py b/reel_header_definition.py
similarity index 100%
rename from header_definition.py
rename to reel_header_definition.py
diff --git a/revisions.py b/revisions.py
index bf19be9..e91fb2b 100644
--- a/revisions.py
+++ b/revisions.py
@@ -1,15 +1,19 @@
-SEGY_REVISION_0 = 0x0000
-SEGY_REVISION_1 = 0x0100
+SEGY_REVISION_0 = 0
+SEGY_REVISION_1 = 1
VARIANTS = {SEGY_REVISION_0: SEGY_REVISION_0,
SEGY_REVISION_1: SEGY_REVISION_1,
+ 0: SEGY_REVISION_0,
+ 1: SEGY_REVISION_1,
100: SEGY_REVISION_1,
256: SEGY_REVISION_1}
+
class SegYRevisionError(Exception):
pass
+
def canonicalize_revision(revision):
"""Canonicalize a SEG Y revision.
diff --git a/segypy.py b/segypy.py
deleted file mode 100644
index ba23890..0000000
--- a/segypy.py
+++ /dev/null
@@ -1,565 +0,0 @@
-"""
-A Python module for reading/writing/manipulating
-SEG-Y formatted files
-
-segy.readSegy : Read SEGY file
-segy.read_reel_header : Get SEGY header
-segy.read_trace_header : Get SEGY Trace header
-segy.read_all_trace_headers : Get all SEGY Trace headers
-segy.getSegyTrace : Get SEGY Trace header and trace data for one trace
-
-segy.writeSegy : Write a data to a SEGY file
-segy.writeSegyStructure : Writes a segpy data structure to a SEGY file
-
-segy.read_binary_value : Get a value from a binary string
-segy.ibm2ieee : Convert IBM floats to IEEE
-
-segy.version : The version of SegyPY
-segy.verbose : Amount of verbose information to the screen
-"""
-#
-# segpy : A Python module for reading and writing SEG-Y formatted data
-#
-# Forked by Robert Smallshire from the original segypy by
-#
-# (C) Thomas Mejer Hansen, 2005-2006
-#
-# with contributions from Pete Forman and Andrew Squelch 2007
-from collections import namedtuple
-import os
-
-import sys
-import struct
-import logging
-
-from numpy import (transpose, reshape, zeros, arange)
-
-from revisions import canonicalize_revision
-from header_definition import HEADER_DEF
-from trace_header_definition import TRACE_HEADER_DEF
-from ibm_float import ibm2ieee2
-
-FORMAT = '%(asctime)-15s %(levelname)s %(message)s'
-logging.basicConfig(level=logging.INFO, format=FORMAT)
-logger = logging.getLogger('segpy.segypy')
-
-version = '0.3.1'
-
-REEL_HEADER_NUM_BYTES = 3600
-TRACE_HEADER_NUM_BYTES = 240
-
-l_char = struct.calcsize('c')
-l_uchar = struct.calcsize('B')
-l_float = struct.calcsize('f')
-
-
-DATA_SAMPLE_FORMAT = {1: 'ibm',
- 2: 'l',
- 3: 'h',
- 5: 'f',
- 8: 'B'}
-
-CTYPES = {'l': 'l', 'long': 'l', 'int32': 'l',
- 'L': 'L', 'ulong': 'L', 'uint32': 'L',
- 'h': 'h', 'short': 'h', 'int16': 'h',
- 'H': 'H', 'ushort': 'H', 'uint16': 'H',
- 'c': 'c', 'char': 'c',
- 'B': 'B', 'uchar': 'B',
- 'f': 'f', 'float': 'f',
- 'ibm': 'ibm'}
-
-# TODO This is redundant with data in the SH_def below
-CTYPE_DESCRIPTION = {'ibm': 'IBM float',
- 'l': '32 bit integer',
- 'h': '16 bit integer',
- 'f': 'IEEE float',
- 'B': '8 bit char'}
-
-
-def size_in_bytes(ctype):
- if ctype == 'l' and struct.calcsize(ctype) == 8:
- return 4 # 64-bit issue?
- return struct.calcsize(ctype) if ctype != 'ibm' else struct.calcsize('f')
-
-
-def get_default_segy_header(ntraces=100, ns=100):
- """
- header = getDefaultSegyHeader()
- """
- # TraceSequenceLine
- header = {'Job': {'pos': 3200, 'type': 'int32', 'def': 0}}
-
- for key in HEADER_DEF:
- header[key] = HEADER_DEF[key].get('def', 0)
-
- header['ntraces'] = ntraces
- header['ns'] = ns
-
- return header
-
-
-def get_default_segy_trace_headers(ntraces=100, ns=100, dt=1000):
- """
- SH = getDefaultSegyTraceHeader()
- """
- # INITIALIZE DICTIONARY
- trace_header = {'TraceSequenceLine': {'pos': 0, 'type': 'int32'}}
-
- for key in TRACE_HEADER_DEF:
- trace_header[key] = zeros(ntraces)
-
- for a in range(ntraces):
- trace_header['TraceSequenceLine'][a] = a + 1
- trace_header['TraceSequenceFile'][a] = a + 1
- trace_header['FieldRecord'][a] = 1000
- trace_header['TraceNumber'][a] = a + 1
- trace_header['ns'][a] = ns
- trace_header['dt'][a] = dt
- return trace_header
-
-
-def read_trace_header(f, reel_header, trace_header_name='cdp', endian='>'):
- """Read a particular trace header entry for all traces.
-
- Args:
- f: A seekable file-linke-object containing the binary SEG Y data.
-
- reel_header: A dictionary obtained by parsing the reel-header using the read_reel_header() function. The
- dictionary must contain integer values for both the 'ntraces' (number of traces) and 'ns' (number of
- samples per trace) keys.
-
- trace_header_name: The name of the entry to read from the trace header.
-
- endian: '>' or '<' for big-endian and little-endian respectively.
-
- Returns:
- A numpy.ndarray containing the requested header entry values for each trace, in SEG Y file order.
- """
- bps = get_byte_per_sample(reel_header)
-
- # MAKE SOME LOOKUP TABLE THAT HOLDS THE LOCATION OF HEADERS
- trace_header_pos = TRACE_HEADER_DEF[trace_header_name]['pos']
-
- # TODO: Be consistent between 'type' and 'format' here.
- trace_header_format = TRACE_HEADER_DEF[trace_header_name]['type']
- ntraces = reel_header['ntraces']
- trace_header_values = zeros(ntraces)
- binary_reader = create_binary_reader(f, trace_header_format, endian)
- start_pos = trace_header_pos + REEL_HEADER_NUM_BYTES
- stride = reel_header['ns'] * bps + TRACE_HEADER_NUM_BYTES
- end_pos = start_pos + (ntraces - 1) * stride + 1
- for i, pos in enumerate(xrange(start_pos, end_pos, stride)):
- trace_header_values[i] = binary_reader(pos)
- return trace_header_values
-
-
-def read_all_trace_headers(f, reel_header, endian='>'):
- """Read all trace headers.
-
- Args:
- f: A seekable file-linke-object containing the binary SEG Y data.
-
- reel_header: A dictionary obtained by parsing the reel-header using the read_reel_header() function. The
- dictionary must contain integer values for both the 'ntraces' (number of traces) and 'ns' (number of
- samples per trace) keys.
-
- endian: '>' or '<' for big-endian and little-endian respectively.
-
- Returns:
- A dictionary mapping string trace header names to numpy.ndarray trace header values, in SEG Y file order.
- """
- trace_headers = {'filename': reel_header['filename']}
-
- logger.debug('read_all_trace_headers : '
- 'trying to get all segy trace headers')
-
- for i, key in enumerate(TRACE_HEADER_DEF):
- trace_header = read_trace_header(f, reel_header, key, endian)
- trace_headers[key] = trace_header
- logger.info("read_all_trace_headers : {!r} {} of {} ".format(key, i, len(TRACE_HEADER_DEF)))
-
- return trace_headers
-
-
-TraceAttributeSpec = namedtuple('Record', ['name', 'pos', 'type'])
-
-
-def compile_trace_header_format(endian='>'):
-
- record_specs = sorted([TraceAttributeSpec(name, TRACE_HEADER_DEF[name]['pos'], TRACE_HEADER_DEF[name]['type']) for name in TRACE_HEADER_DEF],
- key=lambda r : r.pos)
-
- fmt = [endian]
- length = 0
- for record_spec in record_specs:
-
- shortfall = length - record_spec.pos
- if shortfall:
- fmt.append(str(shortfall) + 'x') # Ignore bytes
- length += shortfall
-
- ctype = CTYPES[record_spec.type]
- fmt.append(ctype)
- length += size_in_bytes(ctype)
-
- assert length == TRACE_HEADER_NUM_BYTES
-
- return struct.Struct(''.join(fmt))
-
-
-def file_length(f):
- pos = f.tell()
- f.seek(0, os.SEEK_END)
- file_size = f.tell()
- f.seek(pos, os.SEEK_SET)
- return file_size
-
-
-def _filename(f):
- return f.name if hasattr(f, 'name') else ''
-
-
-def read_segy(f, endian='>'):
- """
- data, header, trace_headers = read_reel_header(f)
- """
-
- file_size = file_length(f)
-
- logger.debug("read_segy : Length of data : {0}".format(file_size))
-
- reel_header = read_reel_header(f, endian) # modified by A Squelch
-
- # GET TRACE
- index = REEL_HEADER_NUM_BYTES
- bytes_per_sample = get_byte_per_sample(reel_header)
- num_data = (file_size - REEL_HEADER_NUM_BYTES) / bytes_per_sample
-
- data, reel_header, trace_headers = read_traces(f,
- reel_header,
- num_data,
- bytes_per_sample,
- index,
- endian)
-
- logger.debug("read_segy : Read segy data") # modified by A Squelch
-
- return data, reel_header, trace_headers
-
-
-def read_traces(f,
- reel_header,
- num_data,
- bytes_per_sample,
- index,
- endian='>'): # added by A Squelch
- """Read the trace data.
-
- values, SegyHeader, SegyTraceHeaders = read_traces(data,
- reel_header,
- num_data,
- bytes_per_sample,
- index)
- """
-
- # Calculate number of dummy samples needed to account for Trace Headers
- num_dummy_samples = TRACE_HEADER_NUM_BYTES / bytes_per_sample
- logger.debug("read_traces : num_dummy_samples = " + str(num_dummy_samples))
-
- # READ ALL SEGY TRACE HEADERS
- trace_headers = read_all_trace_headers(f, reel_header, endian)
-
- logger.info("read_traces : Reading segy data")
-
- dsf = reel_header['DataSampleFormat']
- ctype = DATA_SAMPLE_FORMAT[dsf]
- description = CTYPE_DESCRIPTION[ctype]
- logger.debug("read_traces : Assuming DSF = {0}, {1}".format(
- dsf, description))
- values, _ = read_binary_value(f, index, ctype, endian, num_data)
-
- logger.debug("read_traces : - reshaping")
- values = reshape(values,
- (reel_header['ntraces'],
- reel_header['ns'] + num_dummy_samples))
- logger.debug("read_traces : - stripping header dummy data")
- values = values[:, num_dummy_samples:
- (reel_header['ns'] + num_dummy_samples)]
- logger.debug("read_traces : - transposing")
- values = transpose(values)
-
- # SOMEONE NEEDS TO IMPLEMENT A NICER WAY DO DEAL WITH DSF = 8
- if reel_header['DataSampleFormat'] == 8:
- for i in arange(reel_header['ntraces']):
- for j in arange(reel_header['ns']):
- if values[i][j] > 128:
- values[i][j] = values[i][j] - 256
-
- logger.debug("read_traces : Finished reading segy data")
-
- return values, reel_header, trace_headers
-
-
-def read_reel_header(f, endian='>'):
- """
- reel_header = read_reel_header(file_handle)
- """
- filename = _filename(f)
- reel_header = {'filename': filename}
- for key in HEADER_DEF:
- pos = HEADER_DEF[key]['pos']
- format = HEADER_DEF[key]['type']
-
- reel_header[key], index = read_binary_value(f, pos, format, endian)
-
- logger.debug(str(pos) + " " +
- str(format) +
- " Reading " + key +
- "=" + str(reel_header[key]))
-
- # SET NUMBER OF BYTES PER DATA SAMPLE
- bps = get_byte_per_sample(reel_header)
-
- file_size = file_length(f)
- ntraces = (file_size - REEL_HEADER_NUM_BYTES) / \
- (reel_header['ns'] * bps + TRACE_HEADER_NUM_BYTES)
- reel_header['ntraces'] = ntraces
-
- logger.debug('read_reel_header : successfully read ' + filename)
-
- return reel_header
-
-
-def write_segy(filename, data, dt=1000, trace_header_in=None, header_in=None):
- """
- write_segy(filename, data, dt)
-
- Write SEGY
-
- See also read_segy
-
- (c) 2005, Thomas Mejer Hansen
-
- MAKE OPTIONAL INPUT FOR ALL SEGYHTRACEHEADER VALUES
-
- """
- if header_in is None:
- header_in = {}
- if trace_header_in is None:
- trace_header_in = {}
-
- logger.debug("write_segy : Trying to write " + filename)
-
- shape = data.shape
- ns = shape[0]
- ntraces = shape[1]
- print ntraces, ns
-
- header = get_default_segy_header(ntraces, ns)
- trace_header = get_default_segy_trace_headers(ntraces, ns, dt)
-
- # Add trace_header_in, if exists...
- for key in trace_header_in.keys():
- print key
- for a in range(ntraces):
- trace_header[key] = trace_header_in[key][a]
-
- # Add header_in, if exists...
- for key in header_in.keys():
- print key
- header[key] = header_in[key]
-
- write_segy_structure(filename, data, header, trace_header)
-
-
-def write_segy_structure(filename,
- data,
- header,
- trace_header,
- endian='>'): # modified by A Squelch
- """
- writeSegyStructure(filename, data, header, trace_header)
-
- Write SEGY file using SegyPy data structures
-
- See also readSegy
-
- (c) 2005, Thomas Mejer Hansen
-
- """
-
- logger.debug("writeSegyStructure : Trying to write " + filename)
-
- f = open(filename, 'wb')
-
- # VERBOSE INF
- revision = canonicalize_revision(header['SegyFormatRevisionNumber'])
- dsf = header['DataSampleFormat']
-
- try: # block added by A Squelch
- data_descriptor = HEADER_DEF['DataSampleFormat']['descr'][revision][dsf]
- except KeyError:
- logging.critical(" An error has occurred interpreting a SEGY binary"
- "header key")
- logging.critical(" Please check the Endian setting for this "
- "file: {0}".format(header["filename"]))
- sys.exit()
-
- logger.debug("writeSegyStructure : SEG-Y revision = " + str(revision))
- logger.debug("writeSegyStructure : DataSampleFormat = " +
- str(dsf) +
- "(" + data_descriptor + ")")
-
- # WRITE SEGY HEADER
-
- for key in HEADER_DEF.keys():
- pos = HEADER_DEF[key]["pos"]
- format = HEADER_DEF[key]["type"]
- value = header[key]
- put_value(value, f, pos, format, endian)
-
- # SEGY TRACES
- ctype = HEADER_DEF['DataSampleFormat']['datatype'][revision][dsf]
- bps = HEADER_DEF['DataSampleFormat']['bps'][revision][dsf]
-
- sizeT = TRACE_HEADER_NUM_BYTES + header['ns'] * bps
-
- for itrace in range(header['ntraces']):
- index = REEL_HEADER_NUM_BYTES + itrace * sizeT
- logger.debug('Writing Trace #' +
- str(itrace + 1) +
- '/' + str(header['ntraces']))
- # WRITE SEGY TRACE HEADER
- for key in TRACE_HEADER_DEF.keys():
- pos = index + TRACE_HEADER_DEF[key]['pos']
- format = TRACE_HEADER_DEF[key]['type']
- value = trace_header[key][itrace]
- logger.debug(str(pos) + " " +
- str(format) +
- " Writing " + key +
- "=" + str(value))
- put_value(value, f, pos, format, endian)
-
- # Write Data
- cformat = endian + ctype
- for s in range(header['ns']):
- strVal = struct.pack(cformat, data[s, itrace])
- f.seek(index +
- TRACE_HEADER_NUM_BYTES +
- s * struct.calcsize(cformat))
- f.write(strVal)
-
- f.close()
-
-
-def put_value(value, fileid, index, ctype='l', endian='>', number=1):
- """
- putValue(data, index, ctype, endian, number)
- """
- ctype = CTYPES[ctype]
-
- cformat = endian + ctype*number
-
- logger.debug('putValue : cformat : ' + cformat + ' ctype = ' + ctype)
-
- str_val = struct.pack(cformat, value)
- fileid.seek(index)
- fileid.write(str_val)
-
- return 1
-
-
-def create_binary_reader(f, ctype='l', endian='>'):
- """Create a unary callable which reads a given binary data type from a file.
- """
- ctype = CTYPES[ctype]
- size = size_in_bytes(ctype)
-
- cformat = endian + ctype
-
- def reader(index):
- f.seek(index, os.SEEK_SET)
- data = f.read(size)
- # TODO: Check the content of data before proceeding
- value = struct.unpack(cformat, data)
- return value[0]
-
- return reader
-
-
-def read_binary_value(f, index, ctype='l', endian='>', number=1):
- """
- read_binary_value(data, index, ctype, endian, number)
- """
-
- ctype = CTYPES[ctype]
-
- size = size_in_bytes(ctype)
-
- cformat = endian + ctype * number
-
- logger.debug('read_binary_value : cformat : ' + cformat)
-
- index_end = index + size * number
-
- f.seek(index, os.SEEK_SET)
- data = f.read(size * number)
- if ctype == 'ibm':
- # ASSUME IBM FLOAT DATA
- value = range(number)
- for i in arange(number):
- index_ibm = i * 4
- value[i] = ibm2ieee2(data[index_ibm: index_ibm + 4])
- # this returns an array as opposed to a tuple
- else:
- # TODO: Check the content of data before proceeding
- value = struct.unpack(cformat, data)
-
- if ctype == 'B':
- logger.warning('read_binary_value : '
- 'Inefficient use of 1 byte Integer...', 1)
-
- logger.debug('read_binary_value : ' +
- 'start = ' + str(index) +
- ' size = ' + str(size) +
- ' number = ' + str(number) +
- ' value = ' + str(value) +
- ' cformat = ' + str(cformat))
-
- if number == 1:
- return value[0], index_end
- else:
- return value, index_end
-
-
-def get_byte_per_sample(header):
- revision = canonicalize_revision(header["SegyFormatRevisionNumber"])
- dsf = header["DataSampleFormat"]
-
- try: # block added by A Squelch
- bps = HEADER_DEF["DataSampleFormat"]["bps"][revision][dsf]
- except KeyError:
- # TODO: This should not be a critical failure - should just convert
- # exception
- logging.critical(" An error has occurred interpreting a SEGY "
- "binary header key")
- logging.critical("Please check the Endian setting for "
- "this file: {0}".format(header["filename"]))
- sys.exit()
-
- logger.debug("getBytePerSample : bps = " + str(bps))
-
- return bps
-
-
-def main():
- s = compile_trace_header_format()
- pass
-
- #filename = 'data/stack_final_scaled50_int8.sgy'
- #with open(filename, 'rb') as segy:
- # data, header, trace_headers = read_segy(segy)
-
-
-if __name__ == '__main__':
- main()
diff --git a/testsegy.py b/testsegy.py
deleted file mode 100644
index 936cbfc..0000000
--- a/testsegy.py
+++ /dev/null
@@ -1,80 +0,0 @@
-#!/usr/bin/env python
-#
-# testsegy.py
-#
-
-import segypy
-
-
-def test_read(filename):
- # Read Segy File
- with open(filename, 'rb') as f:
- return segypy.read_segy(f, filename)
-
-
-def test_write(filename_out, data, sh, sth):
- sh['DataSampleFormat'] = 5
- sh['SegyFormatRevisionNumber'] = 100
- segypy.write_segy_structure(filename_out, data, sh, sth)
-
- segypy.wiggle(data, sh, 2, .1, .1)
-
- f_ieee = 'data_IEEE.segy'
- f_ibm = 'data_IBM_REV1.segy'
- d_ieee, sh, sth = segypy.read_segy(f_ieee)
- d_ibm, sh, sth = segypy.read_segy(f_ibm)
-
- return d_ieee, d_ibm, data
-
-
-def test_render(d_ieee, d_ibm, data):
- import pylab
-
- # imshow(Data)
- pylab.figure(1)
- pylab.imshow(d_ieee)
- pylab.title('ieee')
- pylab.show()
-
- pylab.figure(2)
- pylab.imshow(d_ibm)
- pylab.title('IBM')
- pylab.show()
-
- pylab.figure(3)
- pylab.imshow(data)
- pylab.title('TEST')
- pylab.show()
-
-
-def parse_args():
- import argparse
-
- parser = argparse.ArgumentParser()
-
- parser.add_argument('filename',
- help='The SEGY input file')
-
- parser.add_argument('--output', '-o',
- dest='outfile',
- default='',
- help='The output file (optional).')
-
- parser.add_argument('--render_test', '-r',
- dest='render_test',
- action='store_true',
- help='Whether to test rendering (only applies if '
- 'output is produced.)')
-
- return parser.parse_args()
-
-if __name__ == '__main__':
- args = parse_args()
-
- read_results = test_read(args.filename)
-
- if args.outfile:
- write_results = test_write(args.outfile, *read_results)
-
- if args.render_test:
- test_render(*write_results)
diff --git a/toolkit.py b/toolkit.py
new file mode 100644
index 0000000..e53cea8
--- /dev/null
+++ b/toolkit.py
@@ -0,0 +1,261 @@
+from array import array
+from collections import namedtuple
+import itertools
+import os
+import struct
+
+from catalog import CatalogBuilder
+from datatypes import CTYPES, size_in_bytes
+from reel_header_definition import HEADER_DEF
+from ibm_float import ibm2ieee
+from revisions import canonicalize_revision
+from trace_header_definition import TRACE_HEADER_DEF
+
+
+REEL_HEADER_NUM_BYTES = 3600
+TRACE_HEADER_NUM_BYTES = 240
+
+
+def extract_revision(reel_header):
+ """Obtain the SEG Y revision from the reel header.
+
+ Args:
+ reel_header: A dictionary containing a reel header, such as obtained
+ from read_reel_header()
+
+ Returns:
+ One of the constants revisions.SEGY_REVISION_0 or
+ revisions.SEGY_REVISION_1
+ """
+ raw_revision = reel_header['SegyFormatRevisionNumber']
+ return canonicalize_revision(raw_revision)
+
+
+def bytes_per_sample(reel_header, revision):
+ """Determine the number of bytes per sample from the reel header.
+
+ Args:
+ reel_header: A dictionary containing a reel header, such as obtained
+ from read_reel_header()
+
+ revision: One of the constants revisions.SEGY_REVISION_0 or
+ revisions.SEGY_REVISION_1
+
+ Returns:
+ An integer number of bytes per sample.
+ """
+ dsf = reel_header['DataSampleFormat']
+ bps = HEADER_DEF["DataSampleFormat"]["bps"][revision][dsf]
+ return bps
+
+
+def samples_per_trace(reel_header):
+ """Determine the number of samples per trace from the reel header.
+
+ Note: There is no requirement for all traces to be of the same length,
+ so this value should be considered indicative only, and as such is
+ mostly useful in the absence of other information. The actual number
+ of samples for a specific trace should be retrieved from individual
+ trace headers.
+
+ Args:
+ reel_header: A dictionary containing a reel header, such as obtained
+ from read_reel_header()
+
+ Returns:
+ An integer number of samples per trace
+ """
+ return reel_header['ns']
+
+
+def trace_length_bytes(reel_header, bps):
+ """Determine the trace length in bytes from the reel header.
+
+ Note: There is no requirement for all traces to be of the same length,
+ so this value should be considered indicative only, and as such is
+ mostly useful in the absence of other information. The actual number
+ of samples for a specific trace should be retrieved from individual
+ trace headers.
+
+ Args:
+ reel_header: A dictionary containing a reel header, such as obtained
+ from read_reel_header()
+
+ bps: The number of bytes per sample, such as obtained from a call to
+ bytes_per_sample()
+
+ """
+ return samples_per_trace(reel_header) * bps + TRACE_HEADER_NUM_BYTES
+
+
+def read_reel_header(fh, endian='>'):
+ """Read the SEG Y reel header, also known as the binary header.
+
+ Args:
+ fh: A file-like-object open in binary mode positioned such that the
+ beginning of the reel header will be the next byte to be read.
+
+ endian: '>' for big-endian data (the standard and default), '<' for
+ little-endian (non-standard)
+ """
+ reel_header = {}
+ for key in HEADER_DEF:
+ pos = HEADER_DEF[key]['pos']
+ ctype = HEADER_DEF[key]['type']
+ values = tuple(read_binary_values(fh, pos, ctype, 1, endian))
+ reel_header[key] = values[0]
+ return reel_header
+
+
+def catalog_traces(fh, bps, endian='>'):
+ """Determine the file offsets of each trace in the SEG Y file.
+
+ Args:
+ fh: A file-like-object open in binary mode.
+
+ bps: The number of bytes per sample, such as obtained by a call to
+ bytes_per_sample()
+
+ endian: '>' for big-endian data (the standard and default), '<' for
+ little-endian (non-standard)
+
+ Returns:
+ An immutable sequence containing byte offsets to the beginning of each trace.
+ """
+ trace_header_format = compile_trace_header_format(endian)
+
+ pos_begin = REEL_HEADER_NUM_BYTES
+
+ trace_catalog_builder = CatalogBuilder()
+ line_catalog = CatalogBuilder()
+ cdp_catalog = CatalogBuilder()
+
+ for trace_number in itertools.count():
+ fh.seek(pos_begin)
+ data = fh.read(TRACE_HEADER_NUM_BYTES)
+ if len(data) < TRACE_HEADER_NUM_BYTES:
+ break
+ trace_header = TraceHeader._make(trace_header_format.unpack(data))
+ num_samples = trace_header.ns
+ samples_bytes = num_samples * bps
+ trace_catalog_builder.add(trace_number, pos_begin)
+ # Should we check the data actually exists?
+ line_catalog.add((trace_header.Inline3D, trace_header.Crossline3D), trace_number)
+ cdp_catalog.add(trace_header.cdp, trace_number)
+ pos_end = pos_begin + TRACE_HEADER_NUM_BYTES + samples_bytes
+ pos_begin = pos_end
+
+ return (trace_catalog_builder.create(),
+ cdp_catalog.create(),
+ line_catalog.create())
+
+
+def read_binary_values(fh, pos, ctype='l', count=1, endian='>'):
+ """Read a series of values from a binary file.
+
+ Args:
+ fh: A file-like-object open in binary mode.
+
+ pos: The file offset in bytes from the beginning from which the data
+ is to be read.
+
+ ctype: The SEG Y data type.
+
+ number: The number of items to be read.
+ Returns:
+ A sequence containing count items.
+ """
+ fmt = CTYPES[ctype]
+ item_size = size_in_bytes(fmt)
+ block_size = item_size * count
+
+ fh.seek(pos, os.SEEK_SET)
+ buf = fh.read(block_size)
+
+ if len(buf) < block_size:
+ raise EOFError("{} bytes requested but only {} available".format(block_size, len(buf)))
+
+ values = unpack_ibm_floats(buf, count) if fmt == 'ibm' else unpack_values(buf, count, item_size, fmt)
+ assert len(values) == count
+ return values
+
+
+def unpack_ibm_floats(data, count):
+ """Unpack a series of binary-encoded big-endian single-precision IBM floats.
+
+ Args:
+ data: A sequence of bytes. (Python 2 - a str object, Python 3 - a bytes object)
+
+ count: The number of floats to be read.
+
+ Returns:
+ A sequence of floats.
+ """
+ return array('f', (ibm2ieee(data[i: i+4]) for i in range(0, count * 4, 4)))
+
+
+def unpack_values(buf, count, item_size, fmt, endian='>'):
+ """Unpack a series items from a byte string.
+
+ Args:
+ data: A sequence of bytes. (Python 2 - a str object, Python 3 - a bytes object)
+
+ count: The number of floats to be read.
+
+ fmt: A format code (one of the values in the datatype.CTYPES dictionary)
+
+ Returns:
+ A sequence of objects with type corresponding to the format code.
+ """
+ c_format = '{}{}{}'.format(endian, count, fmt)
+ return struct.unpack(c_format, buf)
+ # We could use array.fromfile() here. On the one hand it's likely to be faster and more compact,
+ # On the other, it only works on "real" files, not arbitrary file-like-objects and it would require us
+ # to handle endian byte swapping ourselves.
+
+
+_TraceAttributeSpec = namedtuple('Record', ['name', 'pos', 'type'])
+
+
+def compile_trace_header_format(endian='>'):
+ """Compile a format string for use with the struct module from the trace header definition.
+
+ Args:
+ endian: '>' for big-endian data (the standard and default), '<' for
+ little-endian (non-standard)
+
+ Returns:
+ A string which can be used with the struct module for parsing trace headers.
+ """
+
+ record_specs = sorted([_TraceAttributeSpec(name, TRACE_HEADER_DEF[name]['pos'], TRACE_HEADER_DEF[name]['type']) for name in TRACE_HEADER_DEF],
+ key=lambda r : r.pos)
+
+ fmt = [endian]
+ length = 0
+ for record_spec in record_specs:
+
+ shortfall = length - record_spec.pos
+ if shortfall:
+ fmt.append(str(shortfall) + 'x') # Ignore bytes
+ length += shortfall
+
+ ctype = CTYPES[record_spec.type]
+ fmt.append(ctype)
+ length += size_in_bytes(ctype)
+
+ assert length == TRACE_HEADER_NUM_BYTES
+
+ return struct.Struct(''.join(fmt))
+
+
+def _compile_trace_header_record():
+ """Build a TraceHeader namedtuple from the trace header definition"""
+ record_specs = sorted([_TraceAttributeSpec(name, TRACE_HEADER_DEF[name]['pos'], TRACE_HEADER_DEF[name]['type']) for name in TRACE_HEADER_DEF],
+ key=lambda r : r.pos)
+ return namedtuple('TraceHeader', (record_spec.name for record_spec in record_specs))
+
+
+TraceHeader = _compile_trace_header_record()
+
+
diff --git a/util.py b/util.py
new file mode 100644
index 0000000..bae04c8
--- /dev/null
+++ b/util.py
@@ -0,0 +1,110 @@
+import itertools
+import os
+
+
+def pairwise(iterable):
+ """Pairwise iteration.
+
+ Args:
+ iterable: An iterable series.
+
+ Returns:
+ An iterator over 2-tuples.
+ """
+ a, b = itertools.tee(iterable)
+ next(b, None)
+ return itertools.izip(a, b)
+
+
+def contains_duplicates(sorted_iterable):
+ """Determine in an iterable series contains duplicates.
+
+ Args:
+ sorted_iterable: Any iterable series which must be sorted in either
+ ascending or descending order.
+
+ Returns:
+ True if sorted_iterable contains duplicates, otherwise False.
+ """
+ for a, b in pairwise(sorted_iterable):
+ if a == b:
+ return True
+ return False
+
+
+def measure_stride(iterable):
+ """Determine whether successive numeric items differ by a contant amount.
+
+ Args:
+ iterable: An iterable series of numeric values.
+
+ Returns:
+ The difference between successive values (e.g. item[1] - item[0]) if
+ that difference is the same between all successive pairs, otherwise
+ None.
+ """
+ stride = None
+ for a, b in pairwise(iterable):
+ new_stride = b - a
+ if stride is None:
+ stride = new_stride
+ elif stride != new_stride:
+ return None
+ return stride
+
+def minmax(iterable):
+ """Return the minimum and maximum of an iterable series.
+
+ This function requires only a single pass over the data.
+
+ Args:
+ iterable: An iterable series for which to determine the minimum and
+ maximum values.
+
+ Returns:
+ A 2-tuple containing the minimum and maximum values.
+ """
+ iterator = iter(iterable)
+ try:
+ first = next(iterator)
+ except StopIteration:
+ raise ValueError("minmax() arg is an empty iterable series")
+ minimum = first
+ maximum = first
+ for item in iterator:
+ minimum = min(minimum, item)
+ maximum = max(maximum, item)
+ return minimum, maximum
+
+
+def file_length(fh):
+ """Determine the length of a file-like object in bytes.
+
+ Args:
+ fh: A seekable file-like-object.
+
+ Returns:
+ An integer length in bytes.
+ """
+ pos = fh.tell()
+ try:
+ fh.seek(0, os.SEEK_END)
+ return fh.tell()
+ finally:
+ fh.seek(pos, os.SEEK_SET)
+
+
+def filename_from_handle(fh):
+ """Determine the name of the file underlying a file-like object.
+
+ Args:
+ fh: A file-like object.
+
+ Returns:
+ A string containing the file name, or '' if it could not
+ be determined.
+ """
+ try:
+ return fh.name
+ except AttributeError:
+ return ''