From dcac37022db7553cd7641187c2e99e067de9630b Mon Sep 17 00:00:00 2001 From: Robert Smallshire Date: Mon, 1 Sep 2014 13:58:34 +0200 Subject: [PATCH] General cleanup. --- GUI/demogtk.py | 105 ---- GUI/gui.py | 240 -------- GUI/project1.glade | 192 ------ GUI/segygui/segygui.glade | 497 --------------- GUI/test/test.glade | 251 -------- GUI/test/test.glade.bak | 250 -------- GUI/test/test.gladep | 8 - GUI/test/test.gladep.bak | 8 - catalog.py | 357 +++++++++++ datatypes.py | 59 ++ ibm_float.py | 65 +- makerelease.sh | 15 - plotting.py | 34 -- portability.py | 25 + reader.py | 383 ++++++++++++ ...definition.py => reel_header_definition.py | 0 revisions.py | 8 +- segypy.py | 565 ------------------ testsegy.py | 80 --- toolkit.py | 261 ++++++++ util.py | 110 ++++ 21 files changed, 1247 insertions(+), 2266 deletions(-) delete mode 100755 GUI/demogtk.py delete mode 100644 GUI/gui.py delete mode 100644 GUI/project1.glade delete mode 100644 GUI/segygui/segygui.glade delete mode 100644 GUI/test/test.glade delete mode 100644 GUI/test/test.glade.bak delete mode 100644 GUI/test/test.gladep delete mode 100644 GUI/test/test.gladep.bak create mode 100644 catalog.py create mode 100644 datatypes.py delete mode 100755 makerelease.sh delete mode 100644 plotting.py create mode 100644 portability.py create mode 100644 reader.py rename header_definition.py => reel_header_definition.py (100%) delete mode 100644 segypy.py delete mode 100644 testsegy.py create mode 100644 toolkit.py create mode 100644 util.py 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 - - - - True - GTK_PACK_DIRECTION_LTR - GTK_PACK_DIRECTION_LTR - - - - True - _File - True - - - - - - - True - gtk-new - True - - - - - - - True - gtk-open - True - - - - - - - True - gtk-save - True - - - - - - - True - gtk-save-as - True - - - - - - - True - - - - - - True - gtk-quit - True - - - - - - - - - - - True - _Edit - True - - - - - - - True - Segy Header - True - - - - - - - True - Trace Header - True - - - - - - - True - gtk-cut - True - - - - - - - True - gtk-copy - True - - - - - - - True - gtk-paste - True - - - - - - - True - gtk-delete - True - - - - - - - - - - - True - _View - True - - - - - - True - _Help - True - - - - - - - True - _About - True - - - - - - - - - - 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 - - - - True - GTK_PACK_DIRECTION_LTR - GTK_PACK_DIRECTION_LTR - - - - True - _File - True - - - - - - - True - gtk-new - True - - - - - - - True - gtk-open - True - - - - - - - True - gtk-save - True - - - - - - - True - gtk-save-as - True - - - - - - - True - - - - - - True - gtk-quit - True - - - - - - - - - - - True - _Edit - True - - - - - - - True - gtk-cut - True - - - - - - - True - gtk-copy - True - - - - - - - True - gtk-paste - True - - - - - - - True - gtk-delete - True - - - - - - - - - - - True - _View - True - - - - - - True - _Help - True - - - - - - - True - _About - True - - - - - - - - - - 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 - - - - True - GTK_PACK_DIRECTION_LTR - GTK_PACK_DIRECTION_LTR - - - - True - _File - True - - - - - - - True - gtk-new - True - - - - - - - True - gtk-open - True - - - - - - - True - gtk-save - True - - - - - - - True - gtk-save-as - True - - - - - - - True - - - - - - True - gtk-quit - True - - - - - - - - - - - True - _Edit - True - - - - - - - True - gtk-cut - True - - - - - - - True - gtk-copy - True - - - - - - - True - gtk-paste - True - - - - - - - True - gtk-delete - True - - - - - - - - - - - True - _View - True - - - - - - True - _Help - True - - - - - - - True - _About - True - - - - - - - - - - 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 ''