Plotting - Githubissues

Another item I thought about in the 4000 mile drive is graphics. We need a way to provide some functionality without getting into the weeds of graphical libraries. Said another way we need a way to make it easy for someone to plot any mspass data object.

I suggest we use matplotlib as the initial engine as that is what obspy uses too. A first design choice is if the plotter should be all in one master plot engine or a set of functions? I think the answer is we group types of plots and have the plot engine accept any mspass data objects as inputs. e.g. I would suggest the matplotlib plots could be wrapped like this:

def wiggleplot(d, --set of parameters---)

where the function should handle d defined as a mspass TimeSeries, Seismogram, TimeSeriesEnsemble, or ThreeComponentEnsemble. The plot created for TimeSeries and Seismograms are simple matplotlib frame plots with 1 or 3 frames respectively. Ensembles need a different plot method and more options to work effectively. e.g. you will get in real problems if you naively plotted 1000 seismograms in 1000 subplots.

We might also consider a separate rasterplot/imageplot (name to be determined) for ensembles that will create matlab imagesc style plots from ensemble data. i.e. that plot function would only accept ensemble data. TimeSeriesEnsembles would be displayed as a single image while ThreeComponentEnsembles would be displayed with 3 panels of images.

I don't think any of these are hard to implement. Might be a good job for our student to get his feet wet.

Poking around I discovered this interesting package called Fatiando a terra. Seems to be an Italian open source project. Bad thing is it stated as working only with python 2.7. Came up searching for plotting with a wiggle trace variable error - standard black and white shaded seismic wiggle trace.

On the Fatiando github pages I think we could just steal (with appropriate citation) this code - search for seismic_wiggle. Immediately after is seismic_image which could be used as an option. I think we could use these with wrappers and minor mods to create a simple graphical plotter for ensembles.

Well that proved amazingly easy. Fatiando has a seismic_wiggles and seismic_image function that did exactly what we needed and had no dependencies on fatiando itself. Further, it appears they have depricated this code so we are pretty free to use it with no restrictions. We must, of course, acknowledge that the code we distribute was derived from fatiando.

My testing at this point was the unaltered seismic_wiggles and seismic_image functions. I wrote a small test python program to make sure they worked without any changes in python3 as I thought they should looking at the (fairly small) code. All these functions do really is exploit some features of pyplot I wasn't aware of. That was why it was easy to pull from the larger package and have it just work. An interesting feature I discovered in playing with these was that this code:

seismic_wiggle(d,0.005)
seismic_image(d,0.005,cmap='jet')

produced a filled wiggle trace with an image background using the jet colormap defined in pyplot. Produces a useful plot style not possible, for example, in seismic unix, but which I learned to really like when I had promax available.

Before I go very far on this, however, we should discuss the design of graphics further. Some key design questions are:

Do we do this in an object-oriented way and have a "plotter" (or some similar name) as a base class? The potential advantage of that is one could build a BasicSeismicPlotter as a base class that produces a simple plot with little to no decorations. It might be possible to create things like a picker or editor as children of BasicSeismicPlotter. I did something like that in C++ for seispp, but that is adventure land for me in the python world. Reading suggests it is routinely done though and might be good way to go.
Alternatively we could do this in a plain C procedural style with different functions created for different data types and types of plots. The pure procedural way is different names for each type object to be plotted and plot type. i.e. we'd need a set of names like wiggle,wiggle3c,image, and image3c to plot a scalar and 3c ensembles as wiggle trace and image forms respectively.
A better procedural model is generic names for plot types and have the function determine the data type and handle it appropriately. i.e. we would only have only two or three basic plotters called something like wiggleplot, wiggle_va_plot (for wiggle trace with filling and optional image overlay), and imageplot. (don't pick on the names - those just describe the three basic types I think we need for all data including singles. Well, image plots would be pretty useless for a single scalar or 3c seismogram. ) The point is each function creates a specific type of plot but figures out dynamically what the data type is. python does this better than C++ where RTTI requires the obnoxiously complicated typeid or a clean class heirarchy. Too bad, however, that python doesn't have a switch/case construct - that is a horrible deficiency in the language in my opinion but that is off topic.

Note also if we go the object oriented design we have a similar issue about how we handle different data types. i.e. do we use different names for different data types or have the code auto select and have names only for distinct plot types.

I recommend we produce version one using the procedural model based on item 3. We could easily adapt function produced that way to an oop form for a BasicSeismicPlot object.

Cool discovery! That is an interesting project that I don't really heard of, but I guess it make sense as it is deprecated. Anyway, I checked their license, which is BSD-3, and it is compatible with our license, so I think we should be able to use the code with proper citation and license included. I agree that an object-oriented interface is more clean here.

Made a fair amount of progress on this. I now have a single function I call wtvaplot that accepts any of our 4 ccore data object: TimeSeries, Seismogram, TimeSeriesEnsemble, and SeismogramEnsemble (python names not equal to C++ names for the ensembles). Seems to work fine, although this is a good example of how slow python can be if you have a loop running through lots of data. Just translating a 3C ensemble into 3 matrices for 20 seismograms 1000 points long creates a clear delay, although it is no doubt worse because I am running this under spyder. Irrelevant for now as I'll continue the 'make it work before you make it fast' mantra.

The point I want to bring up for discussion here is the design of the api we want to present to the user. I can either objectize the function I have and/or just use it as the implementation for the following proposed graphical engine. Functional features are:

Name: ReflectionStylePlotter. Kind of verbose but it describes what I have now produced. It can plot wiggle traces only, wiggle trace variable area with black shading or some other color choice, wiggle trace with image display underneath, or wiggle trace variable area with image display underneath. The display is properly called reflection style because time is a reversed y axis as used for seismic reflection sections. Some seismologists do no like that type display, but some do. I think we will want a different plotter for horizontal time displays. The way this stuff works it would be easier to start over than make that an option for this plotter engine.
I think the plotting object should be constructed to define a fixed style of plot (list above) and do nothing until a plot method is called. i.e. a bare bones api would start with something like this:
```
class ReflectionStylePlotter
def __init__(self,style='wtva',normalize=False,color_map=None):
  # actual constructor code 
def plot(d):
 # creates a plot for any valid input data d
```
Rest would be getters and putters for plot behavior. The purpose of them is to allow changing the plot style through putter methods rather than sort out a long list of arguments and side effects.
I think there are lots of ways to add functionality through the matplotlib.pyplot handle that is returned from pyplot with the gcf method. I think it is relatively simple to add functionality like a trace editor that could be implemented through inheritance of the simple plot object we can start with. The key, I think, is to make sure the object keeps a copy of the file handle returned by gcf(). matplotlib.pyplot has a very matlab approach to doing graphics that is a bit different from event driven full function gui libraries. Much simpler to implement from my experience, but likely with some fundamental limitations. I'm at the fringe of my knowledge here, but the approach seems feasible.

Anyway, I could use some feedback before I jump into the next stage here.

A different point from the previous one. Where should I put this module and what should the generic module be called? I think we discussed mspasspy.graphics so one could use this kind of construct:

from mspasspy.graphics import ReflectionStylePlotter

We could also call it something like "plotting", but that seems too restrictive. The module is likely to grow organically more than others as some people like nothing better than writing graphical displays.

I now have a bare bones object oriented implementation of the graphics engine discussed here. I ended up calling it SectionPlotter as a reasonable length name that describes the function well. Here is the API definition (cut and paste from the current implementation showing only the defs and intro comments):

class SectionPlotter:
    """
    This class was designed to plot data in a seismic reflection style 
    display.  That means time is the y axis but runs in a reverse direction
    and traces are equally spaced along the x axis. The data can be plotted
    as simple wiggle traces, wiggle trace variable area, an image format,
    and a wiggletrace variable area plot overlaying an image plot.  

    The object itself only defines the style of plot to be produced along
    with other assorted plot parameters.   The plot method can be called on 
    any mspass data object to produce a graphic display.   TimeSeries 
    data produce a frame with one seismogram plotted.  Seismogram data 
    will produce a one frame display with the 3 components arranged in 
    component order from left to right.  TimeSeriesEnsembles produce 
    conventional reflection style sections with the data displayed from 
    left to right in whatever order the ensemble is sorted to.  
    SeismogramEnsembles are the only type that create multiple windows.
    That is, each component is displayed in a different window.   The 
    display can be understood as a conversion of a SeismogramEnsemble to 
    a set of 3 TimeSeriesEnsembles that are displayed in 3 windows.   
    The component can be deduced only from the title string with 
    has :0, :1, and :2 appended to it to denote components 0, 1, and 2
    respectively.  

    The plots produced by this object are simple ones for display only. 
    The concept is that fancier tools for gui interaction like a trace 
    editor be created as children of this base class.  
    """
    def __init__(self,scale=1.0,normalize=False,title=None):
        """
        The constructor initializes all the internal variables that define 
        the possible behaviors of the plotting engine.  The constructor 
        currently does no sanity check on the compatibility of the parameters.
        The intent is overall behavior should be altered from the default 
        only through the change_style method.   The public attributes 
        (scale,normalize, and title) can be changed without hard.  Those 
        defined as private should be altered only indirectly through the 
        change_style method.  

        :param scale:  Set the amount by which the data should be scaled
          before plotting.   Note the assumption is the data amplitude is 
          of order 1.  Use normalize if that assumption is bad.
        :param normalize:   boolean used to determine if data need to be 
          autoscaled. If true a section scaling is produced from the minimum
          and maximum of the input data.  Relative amplitudes between 
          each scalar signal will be preserved.  To change scaling use 
          a mspass amplitude scaling function before calling the plotter.
        :param title: is a string written verbatim at the top of each 
          plot.  For 3C ensembles the string :0, :1, and :2 is added to 
          title for each of the 3 windows drawn for each of the 3 components
          of the parent ensemble.  
        """
...
def change_style(self,newstyle,fill_color='k',color_map='seismic'):
        """
        Use this method to change the plot style.   Options are described 
        below.   Some parameter combinations are illegal and will result in 
        a RuntimeError if used.  The illegal combos are listed below for 
        each style.

        :param newstyle:  This is the primary switch for this method and 
          is used to define the style of plot that will be produced. The 
          allowed values for this string are listed below along with 
          restrictions (when they exist) on related parameters.  

          wtva - wiggle trace variable area.  That means a wiggle traced with 
            the positive quadrants colored by the color defined by the fill_color
            parameter.   If fill_color is passed as None a RuntimeError will 
            be thrown.
          image - data are displayed as an image with the color_map parameter 
            defining the matplotlib color map used to map amplitudes into 
            a given color.   Note as with wtva the data are presumed to be 
            scaled to be of order 1.  color_map must not be None, which it 
            can be if the plot was changed from a previous style like wtva,
            or a RuntimeError exception will be thrown.
          colored_wtva is like wtva but the wtva plot is drawn on top of an 
            image plot.  Since this is the union of wtva and image a RuntimeError
            exception will be throw if either color_map or fill_color are null.
          wiggletrace will produce a standard line graph of each input as a 
            black line (that feature is currently frozen).   color_map and 
            fill_color are ignored for this style so no exceptions should 
            occur when the method is called with this value of newstyle. 
        """
...
    def plot(self,d):
        """
        Call this method to plot any data using the current style setup and any 
        details defined by public attributes.  

        :param d:  is the data to be plotted.   It can be any of the following:
            TimeSeries, Seismogram, TimeSeriesEnsemble, or SeismogramEnsemble.  
            If d is any other type the method will throw a RuntimeError exception.

        :Returns: a matplotlib.pylot plot handle.
        :rtype: The plot handle is what matplotlib.pyplot calls a gcf.  It 
          is the return of pyplot.gcf() and can be used to alter some properties
          of the figure.  See matplotlib documentation.
        """

I feel good about this API as it is simple to use. If one has a data object, d, of one of the supported types they want to plot the following construct is the minimal amount needed to create a useful graphics:

SectionPlotter plt()
plt.plot(d)

If the user found, for example, that the plot needed normalization (nearly guaranteed if they didn't do something beforehand) they could type

plt.plot(d,normalize=True)

Or if they wanted a simple wiggle trace plot

plt.change_style('wiggletrace')
plt.plot(d)

I think this is a good base class that fancier plotting could build on. A key feature of matplotlib is that following matlab the figure can be manipulated after creation with a file handle. For that reason the plot method returns that handle. I think, for example, that handle could be used to add functionality in children of SectionPlotter. That will take some research, however, and recommend we push that down the road for now. I think the basic api is sound and the implementation could change if we find that model doesn't work as I hope it might to produce something like a trace editor, which we will need.

Please comment on this design.

Next step is I think i will build a different base class for the plot most seismologists probably prefer with time running from left to right. (Note did you know the Soviet analog data had seismograms with time running from right to left? Their paper drums created records that were sort of inside-out compared to our paper records and analysis had to use time running right to left - weird and a pure digression by an old guy). It matters because so much of seismic graphics is visual pattern recognition so we should have support for that style plot. I considered having an all in one solution with the SectionPlotter becoming just a SeismicPlotter with a rotate option to produce the style of plot the current implementation produces. I thing it better to split them for two reasons. First, section plots will likely need to evolve to children with a different evolutionary path from what I'll call normal seismogram plots. It would be cleaner to keep them distinct. Second, I think i can improve the performance of SectionPlotter by building this one from scratch now that I fully understand how it works. i.e. after I get the "normal seismogram" plotter to work i may go back and rewrite the current SectionPlotter to make it perform better. I can already see it could get ponderously slow plotting a large ensemble.

I think that plotting API is pretty clear. I am also fine with decision on separating the two plotting style assuming SectionPlotter will be the only one to have more complex features such as picker and editor.

See you created a new repository for the tutorial material. That seems like a good plan. May be worth preserving the following monstrosity that is the test code I used to test the main functionality of the two plotting objects. Crude work so don't expect it to be clear, but you should be able to paste this file into a test.py and run the code to show you what the different plots styles look like. Not every piece is tested because what I'm called a "colored_wtva" plot also tests image plotting implicitly.

#import sys
#sys.path.append('/home/pavlis/src/mspass/python')
from mspasspy.graphics import wtva_raw
from mspasspy.graphics import image_raw
from mspasspy.graphics import wtvaplot
from mspasspy.graphics import SectionPlotter
from mspasspy.graphics import SeismicPlotter
import numpy as np

import matplotlib.pyplot as plt
import mspasspy.ccore as mspass

def rickerwave(f, dt):
    r"""
    Given a frequency and time sampling rate, outputs ricker function. The
    length of the function varies according to f and dt, in order for the
    ricker function starts and ends as zero. It is also considered that the
    functions is causal, what means it starts at time zero. To satisfy sampling
    and stability:

    .. math::

        f << \frac{1}{2 dt}.

    Here, we consider this as:

    .. math::

        f < 0.2 \frac{1}{2 dt}.

    Parameters:

    * f : float
        dominant frequency value in Hz
    * dt : float
        time sampling rate in seconds (usually it is in the order of ms)

    Returns:

    * ricker : float
        Ricker function for the given parameters.

    """
    assert f < 0.2*(1./(2.*dt)), "Frequency too high for the dt chosen."
    nw = 2.2/f/dt
    nw = 2*int(np.floor(nw/2)) + 1
    nc = int(np.floor(nw/2))
    ricker = np.zeros(nw)
    k = np.arange(1, nw+1)
    alpha = (nc-k+1)*f*dt*np.pi
    beta = alpha**2
    ricker = (1.-beta*2)*np.exp(-beta)
    return ricker

y=rickerwave(2.0,0.005)
ny=len(y)
#plt.plot(y)
#plt.show()
# create am empty 1000 x 20 matrix and then add ricker with 10 point 
# moveout per trace
d=np.zeros(shape=[1000,20])
nrows=1000
ncol=20
j0=50
dj=10
for i in range(ncol):
    for j in range(ny):
        jj=j0+j
        d[jj,i]=y[j]
    j0+=dj

wtva_raw(d,0.0,0.005)
image_raw(d,0.0,0.005,cmap='seismic')
plt.title('Test showing combined wtva and image')
plt.show()
# Verify this makes a new plot
wtva_raw(d,0.0,0.001)
plt.title('Test of wtva alone changing time scale')
plt.show()
image_raw(d,0.0,0.05,cmap='seismic')
plt.title('Test of image_raw alone')
plt.show()
wtva_raw(d,0.0,0.005,color=None)
image_raw(d,0.0,0.005,cmap='seismic')
plt.title('Test of wtva with only wiggles - no shading')
plt.show()
# Convert to a SeismogramEnsemble and test the ensemble plotters
# Each column of the matrix d will become one new Seismogram.  
tmpcore=mspass.CoreSeismogram(nrows)
tmp=mspass.Seismogram(tmpcore)
tmp.set_t0(0.0)
tmp.set_npts(nrows)
tmp.set_dt(0.005)
tmp.tref=mspass.TimeReferenceType.Relative
tmp.set_live()
# change only amplitudes for each component 
for k in range(3):
    scale=float(k+1)
    for i in range(nrows):
        tmp.u[k,i]=scale*d[i,0]

# Now push 20 copies to build an ensemble of 20 members
ens=mspass.SeismogramEnsemble()
for j in range(20):
    ens.member.append(tmp)

#print('trying to plot a Seismogram ensemble with 20 members')
#handle=wtvaplot(ens,cmap='seismic',title='test with shade and color')    
#print('repeat with only wiggles')
#handle2=wtvaplot(ens,wiggle_color=None,title='Test wiggles only',normalize=True)
x=ens.member[0]
print('Try plotting a single 3C seismogram')
handle3=wtvaplot(x,cmap='seismic',title='test plot of one 3C seismogram')
y=mspass.ExtractComponent(x,0)
print('Try plotting one TimeSeries object')
handle4=wtvaplot(y,cmap='seismic',title='test plot of one time series fancy')
print('Trying class wrapper for wtvaplot - calling constructor')
plotter=SectionPlotter(title='test SectionPlotter with ensemble')
plothandle=plotter.plot(ens)
print('trying to plot 3c seismogram')
plotter.title='test class wrapper with Seismogram object'
plothandle=plotter.plot(x)
print('trying same with wtva')
plotter.change_style('wtva')
plotter.title='test class wrapper with wtva'
plotter.plot(x)
print('trying with wiggle only')
plotter.change_style('wiggletrace')
plotter.title='test class wrapper as wiggletrace'
plotter.plot(x)
print('testing ensemble plotter with a linear moveout applied through t0 changes')
for j in range(20):
    ens.member[j].t0 += 0.25*j
plotter.change_style('colored_wtva',color_map='jet')
plotter.plot(ens)
print('testing image plot of same')
plotter.change_style('image',color_map='seismic')
plotter.plot(ens)
print('Trying to create SeismicPlotter object')
plt2=SeismicPlotter()
plt2.change_style('wtva')
plt2.title='test wtva with TimeSeries object'
plt2.plot(y)
plt.show()
plt2.title='test wtva with Seismogram object'
plt2.plot(x)
plt.show()
plt2.title='test wtva with SeismogramEnsemble object'
plt2.plot(ens)
plt.show()
tsens=mspass.EnsembleComponent(ens,0)
plt2.change_style('image')
plt2.title='test imageplot with TimeSeriesEnsemble object'
plt2.plot(tsens)
plt.show()
plt2.change_style('colored_wtva')
plt2.title='test colored_wtva with TimeSeriesEnsemble object'
plt2.plot(tsens)
plt.show()
plt2.title='test colored_wtva with SeismogramEnsemble object'
plt2.plot(ens)
plt.show()
plt2.title='test colored_wtva with TimeSeries object'
plt2.plot(y)
plt.show()
plt2.title='test colored_wtva with Seismogram object'
plt2.plot(x)
plt.show()

There are a few residual items I want to do on this plotting class before I move on to the tutorial.

Plotting requires a mix of scaling approaches to accomplish different tasks. First, because of the way the algorithm I borrowed here works the data must be scaled to order 1 or you get either straight lines (small amplitudes) a bar-code likeplot (amplitudes large). Second, sometimes you want true relative amplitude and sometimes you want to balance the amplitudes between signals to maximize dynamic range of the visual display. That requires a rich mix of amplitude scaling methods. I do not thing those belong in the plot api for use in combination with the normalize boolean. I suggest normalize does only the first - scale the data given uniformly to be of order 1. That is totally rational because with this type of display there is no scale for the individual signal amplitudes anyway. I do need to implement the normalize function with that model in the graphic code.
We need a generic scaling module. Interestingly enough I don't see anything in obspy that really addresses this issue. I think the reason is they think of data in a very seismic source centric way where there is a strong emphasis on preserving absolute amplitude in physical units of nm/s. We need to support that too, but we need more general methods to scale data. They key to doing so and not losing an absolute standard is to always update the "calib" attribute in the metadata when the data are scaled by a constant scale factor. (Note I'm talking here about a single scaling factor applied to all samples, not relative gain functions like the agc function produces or one could do with something like a exponential gain function - used frequently in GPR processing). I have started producing this module for ccore. Untested but so far I've pulled from seispp and implemented the following functions: peak amplitude scaling (Linfinity really), rms scaling (L2 scaling really), percentage clip level as used in seismic unix, and MAD (a variant of perf really with percentage 50% and mathematically L1). The scaling module will be a mix of python and C++ with this signature:
```
def scale(d,method='rms',starttime=None,endtime=None,perflevel=None,preserve_history=False):
```
Rather than write comment in that code, I'll describe the module here. The idea is d is any of 4 data types in ccore: TimeSeries, Seismogram, TimeSeriesEnsemble, and SeismogramEnsemble. method is rms, peak, perf, or mad. When perf is specified perflevel will need to also be define. startime and endtime are optional time window parameters defining the section of the data where the amplitude should be computed - default would be the entire waveform segment. An implementation detail is how to implement the history preservation option - i.e. in C++, python, or a mix (why is getting us in the weeks).

If you get a chance run the little test script above and give some thought to anything else you think should be added as a base level for making a simple plot with minimal decorations. i.e. are there other required decorations other than the title. Note - color scale is not the answer for the reason discussed above about amplitude scaling.

Finally I am able to test out the graphics code on my WSL setup - there was a nasty problem of X11 Forwarding that it took me a while to figure out. Anyway, I have run the above test code, and can see that the plot is clean and the API also makes sense.

One thing confused me when first using is that when wtvaplot plots 3C ensemble, it pops up one window at a time, and I need to close the previous window to see the next component. Is that done on purpose? Because I found later that they all show up at once with SeismicPlotter. I think the latter is more intuitive.

Another cosmetic issue is that the naming of different styles in plotter.change_style. It has a mix usage of abbreviation and long names e.g. wtva and wiggletrace. Since we are already naming things with wtvaplot, maybe we could always use abbreviations such as wtva, img_wtva, wt, and img. It shouldn't be more confusing than wtva alone.

Another different between SectionPlotter and SeismicPlotter is that the former pads zero on traces with no data and the latter leaves them empty. I don't have a preference between the two, but just wonder why they behave differently.

About the scaling function, I think you are correct that the generic scaling for data is different from scaling in the plot, and we should treat that as one of the algorithms. Also, this overlaps with the * and / operators you brought up in https://github.com/wangyinz/mspass/issues/60#issuecomment-684767298. It seems to me that the scaling algorithm is better done in C++ rather than Python. The scaling in the plot will be a pure Python thing and it has nothing to do with the scaling algorithm because there is no need to preserve the history for a plot.

Several things to respond to here. All helpful comments:

The 3 window thing and the zero padding differences was created by the change in implementation for the two types of plots. You appear to agree that the change was the right one, so I will take care of that. The first is trivial. The later might not happen for a while as it is more an esthetic issue. You may not have noted that the image plot has no difference if you use a common color map. Images always require data be inserted into a matrix, and if the start and end times are irregular require zero padding. A more subtle issue not apparent from that test is I also changed the sanity check to avoid creating a gigantic matrix. Version 1 uses a fixed size wall while the second uses a ratio test - largest single trace time/total time span < limit. The limit is something sensible for current generation bit maps of around 10000.
Good point about the style names. I'll aim to make those consistent and require a bit less typing.
I have spent the past couple days working on the scaling issue. I have implemented this with a mix of C++ and python and think what I've done will be a good model of how to adapt a C/C++ algorithm to be a full fledged mspass module. The entire set has a fairly lengthy set of complexities that cover a fair range of the issues one might face in any adaptation.

I fixed the two items listed in the previous comment as suggested. Item 3 is finished and was checked into master this morning. I am going to request the add_graphics branch be merged to master.

mspass-team / mspass

Plotting #65