How to store (serialize) shear-wave splitting results?

[adigitoleo]

Hi,

I'm very interested in doing shear-wave splitting in Julia, and I found this package some time ago but haven't had the time to play around much yet. Currently, I am processing some results from measurements taken in SplitLab. While SplitLab is very feature-rich, there are many reasons that I want to leave it behind (not least that I much prefer Julia to Matlab).

Over the last couple of years, I toyed with some implementations in Python but found the performance to be unreasonably slow (but I'm certainly no NumPy or SciPy guru). Then I discovered Julia, and now I'm just starting to try a few ideas.

I have a lot of notes floating around that I still need to organize, but my first gripe with the SplitLab style implementation is the results have too many "fields", in Julia-speak. In some cases it's obvious: clearly we don't need to store both an array of timestamp components AND a string representation of the date... but in other cases there are some tricky problems, and I'm still on the fence about many things. For example, I've just uploaded some snippets (no splitting yet) with my idea for a more minimalist Result type. So far, I've only used it to process my SplitLab results in Julia.

From your experience of using shear-wave splitting analysis, do we really need to store:

• the eigenvalues of the grid search matrices
• the cross-correlation matrix

in the final result?

By including references to the original traces in Result it also becomes harder to serialize, I think (although I do have potential references to SeisIO.SeisChannel in mine, I haven't used them so far). On the other hand, information about the time window, frequency band of the filter, resolution of the search (in both dt and phi) are all important and removing them could make the result harder to reproduce.

The question is: where should we draw the line? What information should persist in files/databases? How much could be factored out into a "configuration file" (e.g. resolution parameters?). I think it's worth thinking about, as more shear-wave splitting packages arise, e.g.:

• SplitPy: https://paudetseis.github.io/SplitPy

  Teleseismic SKS splitting in Python using routines based on SplitLab. Issue about persistent result storage.

• Pytheas: https://github.com/ispingos/pytheas-splitting

  Local SKS splitting in Python, with Cluster Analysis and a Qt GUI.

• SplitWavePy: http://splitwavepy.readthedocs.io

  CLI shear-wave splitting routines in Python, built on ObsPy.

Sorry that this was a bit long-winded, and I'm still getting around to the SeisModels PR (I'll look into what exactly I need for the ray tracing first).

[anowacki]

Thanks for your thoughts on this, Leo. I am open to discussion and thanks for the pointers to other projects I am not familiar with. Sorry for the length of my own reply!

Contents of Result type

To begin with, Result is just a way of returning several things from the splitting function cleanly, so exactly what is in it I think should just be a case of what is returned from the calculation to the user. This in my mind is separate from serialisation or later analysis of the 'raw' result.

Eigenvalue surface

The eigenvalue surface has to be constructed during the search for optimum splitting parameters and determination of uncertainties, and so regardless of whether we return it, that memory is allocated and the array filled. So returning a reference to that array is no issue in terms of performance.

Beyond that, looking at the eigenvalue (and RC if you like) surface alongside the traces is perhaps the only way to really diagnose the quality of the result, and so the user should always have the option of seeing it. Likewise, plotting a result requires this. So I don't think that removing it from the output is very useful to the user. Optionally not returning it would be okay if you plan to hold millions of Results in memory and want the GC to release those matrices and avoid running out of memory.

Your proposal

I understand why you have included things like 'quality' and 'phase' in your proposed Result type, but I think in many ways I consider these to be derived quantities that could go in a ProcessedResult or QCedResult type which contains a Result (which should probably be called a ShearWaveSplittingResult instead), or maybe just a DataFrame. Things like quality or whether something is a null are either subjective, or require further processing of the raw result with variable parameters, or with a choice or combination of measures. Conversely, there is no way from your Result to later construct a different measure of uncertainty or 'nullness' because there aren't details on the analysis window, grid search or eigenvalue/RC surfaces. I also couldn't plot your Result and look to see how good the splitting result was.

Something SeisSplit aims to support is all the geometries supported by Seis itself, so assuming we know longitude and latitude and aren't worried about elevation or depth is not ideal. Some users will care about eastings, northings and elevations; some longitude, latitude and not so much altitude.

Note that currently measures of confidence/uncertainty in the result are not returned from splitting nor stored in Result, but rather calculated from a Result by other functions (e.g., quality). I quite like this but would be open to persuasion that they are so useful that they should be automatically calculated by splitting.

Inclusion of references to Traces

The reason I added references to the input Traces into the structure is because that way, it the user can simply do plot(split) to get a diagnostic plot of the splitting, including the waveforms. It is also then easy to query a result for other things known about the data, such as station name, location, the event associated with it, and so on. If you have the traces in memory, returning another reference to them is no problem, though perhaps again this could be optional.

Serialisation

At the moment SeisSplit doesn't have any tools to store results to disk, and I don't think that any serialisation should aim to completely preserve the Result structure for the reasons above. I think a plain-text format for serialisation would be nice, and perhaps something compatible with say SHEBA. That said, SHEBA's output formats aren't designed for precise coordinates, Cartesian geometries or high-frequency (>100 Hz) data so it would have to be a variation on that theme.

To be reproducible, in theory the only things to keep when serialising are:

• Two file names (or two files themselves, perhaps cut as short as required)
• Start and end window
• Grid search parameters

But in practice it's also useful to store:

• Optimum parameters and uncertainties
• Recovered source polarisation and uncertainty
• Information about the station beyond the file name like name, coordinates, etc.
• Other QC information such as SNR and Wuestefeld et al.'s Q

On top of that, there are other pieces of metadata that you might like to store:

• Processing applied to the data (e.g., filters, instrument response, etc.)
• Data acquisition information (operator, project name, etc.)

I think the place for this last category is elsewhere. Whilst Seis.Traces and SeisIO.SeisChannels (or maybe SeisIO.SeisDatas) contain fields to put arbitrary stuff into them (.meta for Traces), eventually this stuff needs to be organised in a way appropriate to the user and I don't think the result of a shear-wave splitting calculation is the place for that. What if this result is sits alongside the picked arrival time for the P- and S-waves, the event location, moment tensor, and so on? This is more a tabular data structure at this point.

Serialisation I think should optionally include the eigenvalue surface and maybe cross-correlation surface. This leaves open the possibility of later interrogation of the result, and in fact these surfaces are used directly by some codes (e.g., James Wookey's shear wave splitting tomography).

I don't think serialisation of shear wave splitting results needs to include the traces themselves, though optionally leaving around the windowed data may be useful for the future if the original data are lost.

---

Apologies again for the long response, but I hope there are some useful points in there. Please do let me know what you still think should be done differently in SeisSplit!

[adigitoleo]

Thanks for the information! I'm quite inexperienced in this area, so it's very helpful to read through your reasoning. It looks like a well thought-out structure, I think it fits better with what I really wanted to do :) Taking the quality control out of the result type is a better approach than removing the eigenvalue/correlation matrix.

I think you're also right that I was conflating the serialization format with the Result type. I think serialization is what I was really getting at, but since I hadn't really developed splitting routines yet I didn't really understand that I was trying to serialize the wrong stuff! The main motivation for my Result was in fact a need for some "container" to store result data for easier analysis (e.g. plotting on maps). The SplitLab format proved to be a bit frustrating for progammatic data access (highly unpredictable Types, etc). I'll change the issue title, because it's probably best to narrow this down a bit.

I was also thinking plain text, but then (optionally) storing the eigenvalue/correlation matrices could be tricky, since it probably can't be just a flat format...

---

I'll definitely try to use SeisSplit more when I get a chance. I am pretty sure it will be able to replace my feeble attempts which were mostly written as an excuse to learn some Julia. One feature that I will need is flexibility in the choice of seismic phase (all of the projects linked in the OP are limited to SKS/SKKS) and ability to do source-side measurements (hence my "phase" field). Hower, some of that is probably best discussed by email. Would you mind if we discussed some more ideas via email?

[anowacki]

No problem—I'm happy to discuss things generally via email and we can raise issues again at the repo if there are concrete things to do.

I was also thinking plain text, but then (optionally) storing the eigenvalue/correlation matrices could be tricky, since it probably can't be just a flat format...

SHEBA solves this by writing plan-text 'matrices' (i.e., it puts the eigenvalues (λ₂ actually) in columns and rows); these work fine and are self-describing. You could imagine using NetCDF or HDF5 files but it seems overkill to me. Having the eigenvalues in a separate file to the other information about the result makes sense to me, as would only writing it optionally. Coordination with other splitting packages would be cool and allow for easy movement between packages, but designing a SeisSplit convention in the meantime is okay as well.

One feature that I will need is flexibility in the choice of seismic phase (all of the projects linked in the OP are limited to SKS/SKKS) and ability to do source-side measurements

This is definitely a goal of SeisSplit—it has already been used for splitting along downhole geophones and on seismic nodes. There isn't a simple receiver-side correction option yet (it can be done manually), but it is an easy thing to add, and while we're at it, we should add source-side corrections too (see #3).

[adigitoleo]

I somehow missed the write_fasst_* functions in src/io.jl last time. I'm not familiar with the FASST program, do you have a link to that code? It could be easier to use writedlm from DelimitedFiles for this instead of @sprintf.

[anowacki]

FASST is yet another splitting code, which I wrote in modern Fortran years ago. Since I don't really plan on maintaining it, I haven't made it public.

It could be easier to use writedlm from DelimitedFiles for this instead of @sprintf

True, but there is something appealing about having plain-text formats with columns that line up, which is why the IO for FASST files is as it is. However, it gets a bit hard to generalise the precision you need for IO, so I agree just printing the full precision is more reliable.

[adigitoleo]

After finally using CSV.jl today after stubbornly restricting myself to DelimitedFiles for recent weeks, I've found it to be a significant and justifiable improvement. There are several issues about readdlm that can be found on the julilalang tracker, with some ideas about removing these functions from Base altogether. I think CSV.jl is probably a safer bet for I/O of simple text (data) files, despite not being "core Julia". I'm sure there are some edge cases I haven't encountered yet, so I'll keep using it for a bit more before perhaps starting some work on this.

[anowacki]

Okay, thanks for letting me know of your progress in looking into this. What are the issues with DelimitedFiles that would cause problems? From my perspective, the advantages of using a stdlib are quite large—no extra dependencies and functions that are in the system image so fast to run the first time. We'd need to be using the additional functionality of CSV.jl to make it worth adding it as a dependency.

[adigitoleo]

It's mainly CSV.read vs readdlm: the former can parse mixed-type tables for free (without copying columns in most cases, although I've yet to thoroughly test that). I've also found the missingstring parameter to be quite useful, i.e. the ability to parse some string (or "") as missing data during file reading. The ability to just drop the table directly into a "sink", like a DataFrame or TypedTable is also really nice. We might want multi-line headers, which are also easier...

I'm not sure if there are any issues with readdlm that would outright prevent using it, so I can still explore that option as well. It's true that the latency of loading CSV.jl is a bit annoying, and minimal dependencies is a good goal.