title: plotting seismic wiggle traces using matplotlib
body: <p><a href="https://i.stack.imgur.com/XMhoQ.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/XMhoQ.png" alt="enter image description here"></a></p>
  
  <p>I'm trying to recreate the above style of plotting using matplotlib.</p>
  
  <p>The raw data is stored in a 2D numpy array, where the fast axis is time.</p>
  
  <p>Plotting the lines is the easy bit. I'm trying to get the shaded-in areas efficiently.</p>
  
  <p>My current attempt looks something like:</p>
  
  <pre><code>import numpy as np
  from matplotlib import collections
  import matplotlib.pyplot as pylab
  
  #make some oscillating data
  panel = np.meshgrid(np.arange(1501), np.arange(284))[0]
  panel = np.sin(panel)
  
  #generate coordinate vectors.
  panel[:,-1] = np.nan #lazy prevents polygon wrapping 
  x = panel.ravel()
  y = np.meshgrid(np.arange(1501), np.arange(284))[0].ravel() 
  
  #find indexes of each zero crossing
  zero_crossings = np.where(np.diff(np.signbit(x)))[0]+1 
  
  #calculate scalars used to shift "traces" to plotting corrdinates
  trace_centers = np.linspace(1,284, panel.shape[-2]).reshape(-1,1) 
  gain = 0.5 #scale traces
  
  #shift traces to plotting coordinates
  x = ((panel*gain)+trace_centers).ravel()
  
  #split coordinate vectors at each zero crossing
  xpoly = np.split(x, zero_crossings)
  ypoly = np.split(y, zero_crossings)
  
  #we only want the polygons which outline positive values
  if x[0] &gt; 0:
      steps = range(0, len(xpoly),2)
  else:
      steps = range(1, len(xpoly),2)
  
  #turn vectors of polygon coordinates into lists of coordinate pairs
  polygons = [zip(xpoly[i], ypoly[i]) for i in steps if len(xpoly[i]) &gt; 2]
  
  #this is so we can plot the lines as well
  xlines = np.split(x, 284)
  ylines = np.split(y, 284)
  lines = [zip(xlines[a],ylines[a]) for a in range(len(xlines))]  
  
  #and plot
  fig = pylab.figure()
  ax = fig.add_subplot(111)
  col = collections.PolyCollection(polygons)
  col.set_color('k')
  ax.add_collection(col, autolim=True)
  col1 = collections.LineCollection(lines)
  col1.set_color('k')
  ax.add_collection(col1, autolim=True)
  ax.autoscale_view()
  pylab.xlim([0,284])
  pylab.ylim([0,1500])
  ax.set_ylim(ax.get_ylim()[::-1])
  pylab.tight_layout()
  pylab.show()
  </code></pre>
  
  <p>and the result is <a href="https://i.stack.imgur.com/HOLMG.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/HOLMG.png" alt="enter image description here"></a></p>
  
  <p>There are two issues:</p>
  
  <ol>
  <li><p>It does not fill perfectly because I am splitting on the array indexes closest to the zero crossings, not the exact zero crossings. I'm assuming that calculating each zero crossing will be a big computational hit.</p></li>
  <li><p>Performance. It's not that bad, given the size of the problem - around a second to render on my laptop, but i'd like to get it down to 100ms - 200ms.</p></li>
  </ol>
  
  <p>Because of the usage case I am limited to python with numpy/scipy/matplotlib. Any suggestions?</p>
  
  <p><strong>Followup:</strong></p>
  
  <p>Turns out linearly interpolating the zero crossings can be done with very little computational load.  By inserting the interpolated values into the data, setting negative values to nans, and using a single call to pyplot.fill, 500,000 odd samples can be plotted in around 300ms.</p>
  
  <p>For reference, Tom's method below on the same data took around 8 seconds.</p>
  
  <p>The following code assumes an input of a numpy recarray with a dtype that mimics a seismic unix header/trace definition.  </p>
  
  <pre><code>def wiggle(frame, scale=1.0):
          fig = pylab.figure()
          ax = fig.add_subplot(111)        
          ns = frame['ns'][0]
          nt = frame.size
          scalar = scale*frame.size/(frame.size*0.2) #scales the trace amplitudes relative to the number of traces
          frame['trace'][:,-1] = np.nan #set the very last value to nan. this is a lazy way to prevent wrapping
          vals = frame['trace'].ravel() #flat view of the 2d array.
          vect = np.arange(vals.size).astype(np.float) #flat index array, for correctly locating zero crossings in the flat view
          crossing = np.where(np.diff(np.signbit(vals)))[0] #index before zero crossing
          #use linear interpolation to find the zero crossing, i.e. y = mx + c. 
          x1=  vals[crossing]
          x2 =  vals[crossing+1]
          y1 = vect[crossing]
          y2 = vect[crossing+1]
          m = (y2 - y1)/(x2-x1)
          c = y1 - m*x1       
          #tack these values onto the end of the existing data
          x = np.hstack([vals, np.zeros_like(c)])
          y = np.hstack([vect, c])
          #resort the data
          order = np.argsort(y) 
          #shift from amplitudes to plotting coordinates
          x_shift, y = y[order].__divmod__(ns)
          ax.plot(x[order] *scalar + x_shift + 1, y, 'k')
          x[x&lt;0] = np.nan
          x = x[order] *scalar + x_shift + 1
          ax.fill(x,y, 'k', aa=True) 
          ax.set_xlim([0,nt])
          ax.set_ylim([ns,0])
          pylab.tight_layout()
          pylab.show()
  </code></pre>
  
  <p><a href="https://i.stack.imgur.com/vkyFU.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/vkyFU.png" alt="enter image description here"></a></p>
  
  <p>The full code is published at <a href="https://github.com/stuliveshere/PySeis" rel="nofollow noreferrer">https://github.com/stuliveshere/PySeis</a></p>