basically there for v 0.3

agilescientific · Dec 9, 2016 · 98f5cb9 · 98f5cb9
1 parent c0371b6
commit 98f5cb9
Show file tree

Hide file tree

Showing 8 changed files with 4,126 additions and 31 deletions.
diff --git a/.gitignore b/.gitignore
@@ -2,7 +2,7 @@ notebooks
 notebooks/*
 data/*.pdf
 data/*.png
-Untitled.ipynb
+.ipynb
 .ipynb_checkpoints
 
 # Byte-compiled / optimized / DLL files

diff --git a/Hacking.ipynb b/Hacking.ipynb
diff --git a/_version.py b/_version.py
@@ -0,0 +1,7 @@
+# -*- coding: utf-8 -*-
+"""
+Version.
+
+Doing it this way provides for access in setup.py and via __version__
+"""
+__version__ = "0.3"
diff --git a/config.yaml b/config.yaml
@@ -38,16 +38,16 @@ plot_height: 11
 # ...OR Inches per second (more -> taller plot)
 ips: 2
 
-# Min and max time in seconds. Set tmax: 0 or omit it for max time in file.
+# Min and max time in seconds.
 trange:
-  - 1
-  - 2
+  - 0  # Min 0 for top of data.
+  - 0  # Max 0 for bottom of data.
 
 # Wiggle skips (1 -> no skips, 2 -> every 2nd trace, etc.)
 skip: 1
 
 # Seismic data display parameters (wiggle, varden, both):
-display: wiggle
+display: vd
 
 # Multiplication factor for amplitudes 
 gain: 2

diff --git a/data/31_81_PR.png b/data/31_81_PR.png
diff --git a/plotter.py b/plotter.py
@@ -67,7 +67,7 @@ def watermark_seismic(ax, cfg):
     return ax
 
 
-def plot_header(head_ax, s, fs, cfg):
+def plot_header(head_ax, s, fs, cfg, version=''):
     """
     Plot EBCDIC or ASCII header.
     """
@@ -83,7 +83,7 @@ def plot_header(head_ax, s, fs, cfg):
     head_ax.set_yticks([])
     if cfg['credit']:
         head_ax.text(40, 42,
-                     'plot by github.com/agile-geoscience/seisplot',
+                     'plot by github.com/agile-geoscience/seisplot v{}'.format(version),
                      size=fs, color='lightgray',
                      ha='right', va='top'
                      )

diff --git a/seismic.py b/seismic.py
@@ -64,11 +64,18 @@ def __init__(self, data, dtype=float, params=None):
                 print(s.format(self.nxlines))
             self.data = self.data.reshape((self.ninlines, self.nxlines, self.data.shape[-1]))
 
+        # Make sure there are no singleton dimensions.
+        self.data = np.squeeze(self.data)
+
         if self.inlines is None:
             self.inlines = np.linspace(1, self.ninlines, self.ninlines)
         if self.xlines is None:
             self.xlines = np.linspace(1, self.nxlines, self.nxlines)
 
+        # Guarantee we're getting a unique list.
+        self.inlines = np.unique(self.inlines)
+        self.xlines = np.unique(self.xlines)
+
         self.tbasis = np.arange(0, self.nsamples * self.dt, self.dt)
         return
 
@@ -161,45 +168,62 @@ def from_seismic(cls, seismic, n=0, direction='inline'):
         if direction.lower()[0] == 'i':
             if n < 1:
                 n *= seismic.nxlines
-                n = int(n)
+                n = int(np.floor(n))
             data = seismic.data.copy()[n, ...]
             params['dimensions'] = ['i', 't']
-            nx = int(n * seismic.nxlines)
-            params['inlines'] = [seismic.inlines[nx + 1]] * seismic.nxlines
-            params['xlines'] = seismic.xlines[nx:int((n+1)*seismic.nxlines)]
+            params['inlines'] = seismic.inlines[n]
+            params['ninlines'] = 1
+            params['xlines'] = seismic.xlines
+            params['nxlines'] = len(seismic.xlines)
         elif direction.lower()[0] == 'x':
             if n < 1:
                 n *= seismic.ninlines
-                n = int(n)
+                n = int(np.floor(n))
             data = seismic.data.copy()[:, n, :]
             params['dimensions'] = ['x', 't']
-            nx = int(n*seismic.ninlines)
-            params['xlines'] = [seismic.xlines[nx + 1]] * seismic.ninlines
-            params['inlines'] = [seismic.inlines[(r*seismic.nxlines)+n]
-                                 for r in range(seismic.ninlines)]
+            params['inlines'] = seismic.inlines
+            params['ninlines'] = len(seismic.inlines)
+            params['xlines'] = seismic.xlines[n]
+            params['nxlines'] = 1
         elif direction.lower()[0] == 't':
             if n < 1:
                 n *= seismic.nsamples
-                n = int(n)
+                n = int(np.floor(n))
             data = seismic.data.copy()[..., n]
             params['dimensions'] = ['i', 'x']
         else:
             raise SeismicError("No corresponding data.")
         return cls(data, params=params)
 
     @property
-    def olines(self):
+    def inlineidx(self):
+        """
+        The inline names of every trace.
+        """
+        m = np.meshgrid(np.unique(self.inlines), np.unique(self.xlines))
+        return m[0].T.flatten()
+
+    @property
+    def xlineidx(self):
+        """
+        The inline names of every trace.
+        """
+        m = np.meshgrid(np.unique(self.inlines), np.unique(self.xlines))
+        return m[1].T.flatten()
+
+    @property
+    def olineidx(self):
         """
         The other-line numbers.
         """
-        return self.xlines if self.dimensions[0] == 'i' else self.inlines
+        return self.xlineidx if self.dimensions[0] == 'i' else self.inlineidx
 
     @property
-    def slines(self):
+    def slineidx(self):
         """
         The self-line numbers.
         """
-        return self.inlines if self.dimensions[0] == 'i' else self.xlines
+        return self.inlineidx if self.dimensions[0] == 'i' else self.xlineidx
 
     @property
     def slabel(self):
@@ -286,9 +310,11 @@ def plot_spectrum(self,
         f_peak = np.mean(peaks)
         f_min = np.amin(mis)
         f_max = np.amax(mas)
+        dt = 1000 // fs
+        f_nyquist = fs // 2
 
-        statstring = "\nMin: {:.1f} Hz\nPeak: {:.1f} Hz\nMax: {:.1f} Hz"
-        stats = statstring.format(f_min, f_peak, f_max)
+        statstring = "\nMin: {:.1f} Hz\nPeak: {:.1f} Hz\nMax: {:.1f} Hz\nNyquist ({} ms): {} Hz"
+        stats = statstring.format(f_min, f_peak, f_max, dt, f_nyquist)
 
         ax.plot(f, db, lw=0)  # Plot invisible line to get the min
         y_min = ax.get_yticks()[0]

diff --git a/seisplot.py b/seisplot.py
@@ -23,6 +23,7 @@
 from notice import Notice
 import utils
 import plotter
+from _version import __version__
 
 
 def main(target, cfg):
@@ -57,7 +58,13 @@ def main(target, cfg):
         direction = ['xline', 'inline']
 
     # Get the data.
-    ss = [Seismic.from_seismic(s, n=n, direction=d) for n, d in zip((n, xl), direction)]
+    try:
+        ss = [Seismic.from_seismic(s, n=n, direction=d) for n, d in zip((n, xl), direction)]
+    except IndexError as e:
+        # Perhaps misinterpreted 2D as 3D
+        s = Seismic.from_segy(target, params={'ndim': 2})
+        direction = ['inline']
+        ss = [Seismic.from_seismic(s, n=n, direction=d) for n, d in zip((n, xl), direction)]
 
     clip_val = np.percentile(s.data, cfg['percentile'])
 
@@ -178,7 +185,11 @@ def main(target, cfg):
     # Plot text header.
     start = (h - 1.5*mt - fhh) / h
     head_ax = fig.add_axes([ssl, start, wsl/w, fhh/h])
-    head_ax = plotter.plot_header(head_ax, s.header, fs=9, cfg=cfg)
+    head_ax = plotter.plot_header(head_ax,
+                                  s.header,
+                                  fs=9,
+                                  cfg=cfg,
+                                  version=__version__)
 
     # Plot histogram.
     # Params for histogram plot.
@@ -227,8 +238,8 @@ def main(target, cfg):
             im = ax.imshow(line.data.T,
                            cmap=cfg['cmap'],
                            clim=[-clip_val, clip_val],
-                           extent=[line.olines[0],
-                                   line.olines[-1],
+                           extent=[line.olineidx[0],
+                                   line.olineidx[-1],
                                    1000*line.tbasis[-1],
                                    line.tbasis[0]],
                            aspect='auto',
@@ -269,7 +280,7 @@ def main(target, cfg):
                         1000*cfg['trange'][0])
 
         # Crossing point. Will only work for non-arb lines.
-        ax.axvline(ss[i-1].slines[0],
+        ax.axvline(ss[i-1].slineidx[0],
                    c=utils.rgb_to_hex(cfg['highlight_colour']),
                    alpha=0.5
                    )
@@ -306,13 +317,13 @@ def main(target, cfg):
             ylim = ax.get_ylim()
             par1 = ax.twiny()
             par1.spines["top"].set_position(("axes", 1.0))
-            par1.plot(line.slines, np.zeros_like(line.slines), alpha=0)
+            par1.plot(line.slineidx, np.zeros_like(line.slineidx), alpha=0)
             par1.set_xlabel(utils.LABELS[line.slabel], fontsize=cfg['fontsize'])
             par1.set_ylim(ylim)
 
             # Adjust ticks
             tx = par1.get_xticks()
-            newtx = [line.slines[len(line.slines)*(i//len(tx))] for i, _ in enumerate(tx)]
+            newtx = [line.slineidx[len(line.slineidx)*(i//len(tx))] for i, _ in enumerate(tx)]
             par1.set_xticklabels(newtx, fontsize=cfg['fontsize']-2)
 
     t2 = time.time()
@@ -379,7 +390,7 @@ def main(target, cfg):
                         type=str,
                         nargs='?',
                         default='./*.[s,S]*[g,G][y,Y]',
-                        help='The path to one or more SEGY files. Uses Unix-style pathname expansion.')
+                        help='The path to one or more SEGY files. Uses Unix-style pathname expansion. Omit to find all SEGY files in current directory.')
     parser.add_argument('-o', '--out',
                         metavar='output file',
                         type=str,