As of yesterday, I no longer have a full-time day job.
I am looking for opportunities.
I’d love to hear about projects in geophysics, computational geoscience, data science, machine learning. Feel free to get in touch with me at matteo@mycarta.ca.
Thanks,
Matteo
There, I said it, and aloud: ‘EXCEL’!
I got this idea of making a modern (no, I am NOT kidding) educational tool to interactively construct and study Ricker wavelets after reading William Ashcroft’s A Petroleum Geologist’s Guide to Seismic Reflection. The book is a gem, in my opinion (I discovered it thanks to a blog post by Matt Hall a couple of years ago) and comes with some interesting tools on a DVD, which are unfortunately a bit outdated. So I went back to my old Geophysics Bible (Yilmaz’ book Seismic Data Analysis), and mashed a few ideas up; you can Download the Excel file here. Please feel free to use it, peruse it, and abuse it, in your classes, labs, nightmares, etcetera…
The tool is fairly simple. In Sheet 1 the user enters the dominant frequency of the desired Ricker wavelet, as shown in the middle of Figure 1. From that informatin the wavelet is constructed using the equation A = g^2 * 1/exp g^2 where g is the ration between frequency f (in increments of 5 Hz up to an arbitrary 125 Hz – but this could be easily changed!) and the dominant frequency f1 just entered. frequencies. The frequency spectrum of the wavelet is shown as a graph.
Figure 1
In Figure 2 I show the same Sheet, but for a wavelet of dominant frequency equal to 50 Hz.
Figure 2
A number of ancillary graphs, shown in Figure 3, display individual building blocks of the formula A = g^2 * 1/exp g^2.
Figure 3
In Sheet 2 the user can view a plot of the wavelet in the time domain, add Constant phase shifts and Linear phase shifts, and experiment with different combinations of them, as shown in Figures 4-7, below.
Figure 4
Figure 5
Figure 6
Figure 7
Finally, Sheet 3 displays a plot of the wavelet, and of the individual frequency components that make it up. Have fun!!!
Figure 8
A Petroleum Geologist’s Guide to Seismic Reflection, William Ashcroft, 2011, Wiley.
Seismic Data Analysis, Öz Yilmaz, 2001, Society of Exploration Geophysicists.
To plot a wavelet [in Python]. blog post by Evan Bianco.
In Visualization tips for geoscientists: MATLAB, Part III I showed there’s a qualitative correlation between occurrences of antimony mineralization in the southern Tuscany mining district and the distance from lineaments derived from the total horizontal derivative (also called maximum horizontal gradient).
Let’s take a look at the it below (distance from lineaments increases as the color goes from blue to green, yellow, and then red).
However, in a different map in the same post I showed that lineaments derived using the maxima of the hyperbolic tilt angle (Cooper and Cowan, 2006, Enhancing potential field data using filters based on the local phase) are offset systematically from those derived using the total horizontal derivative.
Let’s take a look at the it below: in this case Bouguer gravity values increase as the color goes from blue to green, yellow, and then red; white polygons are basement outcrops.
The lineaments from the total horizontal derivative are in black, those from the maxima of hyperbolic tilt angle are in gray. Which lineaments should be used?
The ideal way to map the location of density contrast edges (as a proxy for geological contacts) would be to gravity forward models, or even 3D gravity inversion, ideally constrained by all available independent data sources (magnetic or induced-polarization profiles, exploratory drilling data, reflection seismic interpretations, and so on).
The next best approach is to map edges using a number of independent gravity data enhancements, and then only use those that collocate.
Cooper and Cowan (same 2006 paper) demonstrate that no single-edge detector method is a perfect geologic-contact mapper. Citing Pilkington and Keating (2004, Contact mapping from gridded magnetic data – A comparison of techniques) they conclude that the best approach is to use “collocated solutions from different methods providing increased confidence in the reliability of a given contact location”.
I show an example of such a workflow in the image below. In the first column from the left is a map of the residual Bouguer gravity from a smaller area of interest in the southern Tuscany mining district (where measurements were made on a denser station grid). In the second column from the left are the lineaments extracted using three different (and independent) derivative-based data enhancements followed by skeletonization. The same lineaments are superimposed on the original data in the third column from the left. Finally, in the last column, the lineaments are combined into a single collocation map to increase confidence in the edge locations (I applied a mask so as to display edges only where at least two methods collocate).
If you want to learn more about this method, please read my note in the Geophysical tutorial column of The Leading Edge, which is available with open access here.
To run the open source Python code, download the iPython/Jupyter Notebook from GitHub.
With this notebook you will be able to:
1) create a full suite of derivative-based enhanced gravity maps;
2) extract and refine lineaments to map edges;
3) create a collocation map.
These technique can be easily adapted to collocate lineaments derived from seismic data, with which the same derivative-based enhancements are showing promising results (Russell and Ribordy, 2014, New edge detection methods for seismic interpretation.)
These below are two copies of a seismic horizon from the open source Penobscot 3D seismic survey coloured using two different colormaps (data from Hall, 2014).
Figure 1
Do you think either of them is ‘better’? If yes, can you explain why? If you are unsure and you want to learn how to answer such questions using perceptual principles and open source Python code, you can read my tutorial Evaluate and compare colormaps (Niccoli, 2014), one of the awesome Geophysical Tutorials from The Leading Edge. In the process you will learn many things, including how to calculate an RGB colormap’s intensity using a simple formula:
import numpy as np ntnst = 0.2989 * rgb[:,0] + 0.5870 * rgb[:,1] + 0.1140 * rgb[:,2] # get the intensity intensity = np.rint(ntnst) # rounds up to nearest integer
…and how to display the colormap as a colorbar with an overlay plot of the intensity as in Figure 2.
Figure 2
Hall, M. (2014) Smoothing surfaces and attributes. The Leading Edge 33, no. 2, 128–129. Open access at: https://github.com/seg/tutorials#february-2014
Niccoli, M. (2014) Evaluate and compare colormaps. The Leading Edge 33, no. 8., 910–912. Open access at: https://github.com/seg/tutorials#august-2014
I recently added to my Matlab File Exchange function, Perceptually improved colormaps, a colormap for periodic data like azimuth or phase. I am going to briefly showcase it using data from my degree thesis in geology, which I used before, for example in Visualization tips for geoscientists – Matlab. Figure 1, from that post, shows residual gravity anomalies in milligals.
Figure 1
Often we’re interested in characterizing these anomalies by calculating the direction of maximum dip at each point on the surface, and for that direction display the azimuth, or dip azimuth. I’ve done this for the surface of residual anomalies from Figure 1 and displayed the azimuth in Figure 2. Azimuth from 0 to 360 degrees are color-coded using Jet, Matlab’s standard colormap (until recently). Typically I do not trust azimuth values when the dip is close to zero because it is often contaminated by noise so I would use shading to de-saturate the colors where dip has the lowest values, but for ease of discussion I haven’t done so in this case.
Figure 2. Azimuth values color-coded with Jet.
There are two problems with Figure 2. First, the well-known problems with the jet colormap. For example, blue is too dark and blue areas appear as bands of constant colour. Yellow is much lighter than any other colour so we see artificial yellow edges that are not really present in the data. But there is an additional issue in Figure 2 because azimuths close in value to 0 and 360 degrees are colored with blue and red, respectively, instead of a single color as they should, causing an additional artificial edge.
In Figure 3 I recolored the map using a colormap that replicates those used in many geophysical software tools to display azimuth or phase data. This is better because it wraps around at 360 degrees but the perceptual issues are unresolved: in this case red, yellow and blue all appear as sharp perceptual edges.
Figure 3. Azimuth values color-coded with generic azimuth colormap.
Figure 4. Azimuth values color-coded with isoluminant azimuth colormap.
In Figure 4 I used my new colormap, called isoAZ (for isoluminant azimuth). This colormap is much better because not only does it wraps around at 360 degrees, but also lightness is held constant for all colors, which eliminates the perceptual anomalies. All the artificial yellow, red, and blue edges are gone, only real edges are left. This can be more easily appreciated in the figure below: if you hover with your mouse over it you are able to switch back and forth between Figure 3 and Figure 4.
From an interpretation point of view, azimuths 180 degrees apart are of opposing colours, which is ideal for dip azimuth data because it allows us to easily recognize folds where dips of opposite direction are juxtaposed at an edge. One example is the sharp edge in the northwest quadrant of Figure 4, where magenta is juxtaposed to green. If you look at Figure 1 you see that there’s a relative high in this area (the edge in Figure 4) with dips of opposite direction on either side (East and West, or 0 and 360 degrees).
The colormap was created in the Lightness-Chroma-Hue color space, a polar transform of the Lab color space, where lightness is the vertical axis and at each value of lightness, chroma is the radial coordinate and hue the polar angle. One limitation of this approach is that due to theirregular shape of the color gamut section at each lightness value, we can never exceed chroma values of about 38-40 (at lightness = 65 in Matlab; in Python, with extensive trial and error, I have not been able to go past 36 using the Scikit-image Color module), which make the resulting colors pale, pastely.
it creates For those that want to experiment with it further, I used just a few lines of code similar to the ones below:
radius = 38; % chroma
theta = linspace(0, 2*pi, 256)'; % hue
a = radius * cos(theta);
b = radius * sin(theta);
L = (ones(1, 256)*65)'; % lightness
Lab = [L, a, b];
RGB=colorspace('RGB<-Lab',Lab(end:-1:1,:));
This code is a modification from an example by Steve Eddins on a post on his Matlab Central blog. In Steve’s example the colormap cycles through the hues as lightness increases monotonically (which by the way is an excellent way to generate a perceptual rainbow). In this case lightness is kept constant and hue cycles through the entire 360 degrees and wraps around. Also, instead of using the Image Processing Toolbox, I used Colorspace, a free function from Matlab File Exchange, for the color space transformations.
For data like fracture orientation where azimuths 180 degrees apart are equivalent it is better to stack two of these isoluminant colormaps in a row. In this way we place opposing colors 90 degrees apart, whereas color 180 degrees apart are the same. You can do it using Matlab commands repmat or vertcat, as below:
radius = 38; % chroma
theta = linspace(0, 2*pi, 128)'; % hue
a = radius * cos(theta);
b = radius * sin(theta);
L = (ones(1, 128)*65)'; % lightness
Lab = [L, a, b];
rgb=colorspace('RGB<-Lab',Lab(end:-1:1,:));
RGB=vertcat(rgb,rgb);
Over the last couple of years I have looked at a number of apps of all sorts. Some were seismometers. Out of those, two had the desired (to me) 3-component recording and export capabilities: Seismometer by Yellowagents, which is for for iPhone, and Seismograph alpha by Calvico, which is for Android and it’s the one I am talking about in this post.
These are several screen captures from Seismograph’s download page.
seismograph alpha
This looked good, so I downloaded it, and played with it a bit.
You can change the color of the three components in the settings. I moved away from the default green, red, and blue (below, left) since the green and red can be confusing for color blind viewers, and used magenta, green, and black, which are less confusing (below, right). I also switched from black to white background and increased the saturation of the new colors.
You can also change the recording sensitivity on the fly with the pinch, see below.
Really cool. But after the initial excitement, I moved on. Because really, how many times can you tap the phone from every possible direction, export the data, load it in Excel, stare at it for a while, show it to your colleagues at work, and friends and family at home?
I recently had an opportunity to go back to this app and record some real data from a dynamite blast in a construction site!
Next to my office in Stavanger they have been building a new condo complex (personal note: as I publish this I am no longer in Stavanger, I moved back to Canada). Of course this being Norway, they soon had to deal with granite, and a whole pile of it.
OK, perhaps not so much rock and not as much as 25 tons of explosive at once, but you can see in the picture below (from Google Street view) that there’s a good amount of rock outcropping behind the parking lot and to the right of the car (and underneath it).
before
All that rock had to go, and so it is that for more than 3 weeks, there were a couple of charges shot every day. This is very good as once I set my mind on recording one shot, it took me a few experiments to get the setup right. But before I get to that, here is a photo I took a few weeks after the works had started, right after I recorded the shot I am using for this post.
after
The next two photos below show, respectively, one of the several excavation fronts, and the big pile of granite removed.
Looking at the excavation front
….
The big pile
As I said it was good there were so many shots as it took me some time to get the right setting and a good recording. The two biggest issues were:
– the positioning of the phone: I started in my office, first placing the phone on the window sill, but I soon realize the ‘coupling’ was not so good as the sill was aluminum and a bit shaky. The floor was a lot better, but I finally tried on cement in the parking lot outside of the building, which proved even better.
– the phone lock: I figured out that the recording stops when the phone locks and it is resumed when you unlock it. The first time this happened right before the shot, which I lost. Also, having to touch your phone to unlock it while recording will result in the recording of a very noisy section.
I finally got to get a good recording. Here’s the raw data for those readers that may want to play with it.
This is the part that took me a bit of work. The first thing to do is to generate a time vector from the recorded date column. This is a fairly common task when working with field data. In the code box below I added the first 50 rows of time data. The date, which I called calendar time, in the second column. The first column, with sample number, was not in the recorded data, I added it here for reference. The third column is the output time vector, which is cumulative time. To get to this result I copied the milliseconds only from the calendar time column. I used UltraEdit, my favourite text editor, which allows to work with columns and column portions, not just rows.
Then I used a MS Excel formula to generate the cumulative time vector. Once I got to this I recognized that the sample rate is variable. I added a column called dt, which is the difference between pairs of consecutive samples, to make it more clear. I am not too sure if THIS is a common occurrence or not. In my MSc thesis worked with ORION land seismometers equipped with three-component geophones and had to correct for drift and with dead time (a short interruption in recording every minute due to GPS clock updates) but not with uneven sample rate. Has anyone seen this before?
sample, calendar time, t(ms), dt(ms)
1, 013 06 27 12:19:36.972, 0, 11 2, 013 06 27 12:19:36.983, 11, 5 3, 013 06 27 12:19:36.988, 16, 2 4, 013 06 27 12:19:36.990, 18, 2 5, 013 06 27 12:19:36.992, 20, 3 6, 013 06 27 12:19:36.995, 23, 3 7, 013 06 27 12:19:36.998, 26, 1 8, 013 06 27 12:19:36.999, 27, 2 9, 013 06 27 12:19:37.001, 29, 2 10, 013 06 27 12:19:37.003, 31, 2 11, 013 06 27 12:19:37.005, 33, 2 12, 013 06 27 12:19:37.007, 35, 2 13, 013 06 27 12:19:37.009, 37, 9 14, 013 06 27 12:19:37.018, 46, 2 15, 013 06 27 12:19:37.020, 48, 17 16, 013 06 27 12:19:37.037, 65, 1 17, 013 06 27 12:19:37.038, 66, 2 18, 013 06 27 12:19:37.040, 68, 34 19, 013 06 27 12:19:37.074, 102, 2 20, 013 06 27 12:19:37.076, 104, 2 21, 013 06 27 12:19:37.078, 106, 2 22, 013 06 27 12:19:37.080, 108, 2 23, 013 06 27 12:19:37.082, 110, 2 24, 013 06 27 12:19:37.084, 112, 2 25, 013 06 27 12:19:37.086, 114, 5 26, 013 06 27 12:19:37.091, 119, 2 27, 013 06 27 12:19:37.093, 121, 4 28, 013 06 27 12:19:37.097, 125, 34 29, 013 06 27 12:19:37.131, 159, 5 30, 013 06 27 12:19:37.136, 164, 2 31, 013 06 27 12:19:37.138, 166, 2 32, 013 06 27 12:19:37.140, 168, 4 33, 013 06 27 12:19:37.144, 172, 2 34, 013 06 27 12:19:37.146, 174, 14 35, 013 06 27 12:19:37.160, 188, 5 36, 013 06 27 12:19:37.165, 193, 14 37, 013 06 27 12:19:37.179, 207, 7 38, 013 06 27 12:19:37.186, 214, 11 39, 013 06 27 12:19:37.197, 225, 12 40, 013 06 27 12:19:37.209, 237, 17 41, 013 06 27 12:19:37.226, 254, 2 42, 013 06 27 12:19:37.228, 256, 14 43, 013 06 27 12:19:37.242, 270, 7 44, 013 06 27 12:19:37.249, 277, 11 45, 013 06 27 12:19:37.260, 288, 16 46, 013 06 27 12:19:37.276, 304, 3 47, 013 06 27 12:19:37.279, 307, 14 48, 013 06 27 12:19:37.293, 321, 8 49, 013 06 27 12:19:37.301, 329, 11 50, 013 06 27 12:19:37.312, 340, 7
To get around this I decided to resample the data, which I did in Matlab.
I used a new sample rate of 10 ms after taking the average of all values in the dt column, which turned out to be almost exactly 10 ms. I suspect this may be the ‘nominal’ sample rate for the app, although I could not find confirmation of this anywhere.
Here’s the interpolated time vector (although it is not really interpolated), which I used to resample (this time there was actual interpolation) the x,y, and z components. For those interested, here is the final ASCII file with interpolated x,y,z,time.
Below is a plot of the data, which I generated in Matlab. There are two sections of high signal in all three recorded components. The first one is at the beginning of the recording, and it is due to the phone still moving for a couple of seconds after I touched it to star the recording. The second section, right after 2 x 10^4 ms, is the recorded shot.
That’s it for this post. In my next post or two I will take a closer look at the data: I will use a tap test to assess polarity, will run a spectral analysis, and try hodogram plots, and hopefully much more.
In a future posts I will take a look at some of the color palettes used for seismic amplitude display, and discuss ways we can design more perceptual and more efficient ones.
For now, I would like to ask readers to look at two sets of seismic images and answer the survey questions in each section. Far from being exhaustive sets, these are meant as a teaser to get a conversation started and exchange opinions and preferences.
The seismic line below is inline 424 from the F3 dataset, offshore Netherlands from the Open Seismic Repository (licensed CC-BY-SA).
I generated an animation, played at 0.5 frames/second, where 8 different color palette are alternated in sequence. Please click on the image to see a full resolution animation. I also generated a 0.25 frame/second version and a 1 frame/second version.
The images used to create the panel below are portions of seismic displays kindly provided by Steve Lynch of 3rd Science Solutions, generated using data released by PeruPetro. I am grateful to both.
Thanks to Matt Hall and Evan Bianco of Agile Geoscience for their suggestions.