Category Archives: Uncategorized

Using Python to calculate northern hemisphere’s surface land coverage

Posted on by
2

Yesterday during my lunch break I was rather bored; it is unseasonably cold for the fall, even in Calgary, and a bit foggy too.
For something to do I browsed the Earth Science beta on Stack Exchange looking for interesting questions (as an aside, I encourage readers to look at the unanswered questions).
There was one that piqued my curiosity, “In the northern hemisphere only, what percentage of the surface is land?”.
It occurred to me that I could get together an answer using an equal area projection map and a few lines of Python code; and indeed in 15 minutes I whipped-up this workflow:

  • Invert and import this B/W image of equal area projection (Peters) for the Northern hemisphere (land = white pixels).
Peters_projection,_black_north

Source of original image (full globe): Wikimedia Commons

  • Store the image as a Numpy array.
  • Calculate the total number of pixels in the image array (black + white).
  • Calculate the total number of white pixels (1s) by summing the entire array. Black pixels (0s) will not contribute.
  • Calculate percentage of white pixels.

The result I got is 40.44%. Here’s the code:

# import libraries
import numpy as np
from skimage import io
from matplotlib import pyplot as plt

# import image
url = 'https://mycartablog.com/wp-content/uploads/2018/09/peters_projection_black_north.png'
north_equal_area = io.imread(url, as_grey=True)

# check the image
fig = plt.figure(figsize=(20, 10))
ax = fig.add_subplot(1, 1, 1)
ax.set_xticks([])
ax.set_yticks([])
plt.imshow(north_equal_area, cmap = 'gray');

# Do the calculations
r, c = np.shape(north_equal_area)
sz =  r*c
s = np.sum(north_equal_area)
print(np.round(s/sz*100, decimals=2))
>>> 40.44

As suggested in a comment to my initial answer, I run the same Python script for the entire globe and got the expected 30% land coverage:

# import image 
url = 'https://mycartablog.com/wp-content/uploads/2018/09/peters_projection_black_full.png'
equal_area = io.imread(url1, as_grey=True)

# Do the calculations 
r1, c1= np.shape(equal_area)
sz1 =  r1*c1
s1 = np.sum(equal_area)
print(np.round(s1/sz1*100, decimals=2))
>>> 30.08

 

What is acquisition footprint noise in seismic data?

Posted on by
7

Acquisition footprint is a noise field that appears on 3D seismic amplitude slices or horizons as an interwoven linear crosshatching parallel to the source line and receiver line directions. It is for the most part an expression of inadequate acquisition geometry, resulting in insufficient sampling of the seismic wave field (aliasing) and irregularities in the offset and azimuth distribution, particularly in the cross line direction.

Sometimes source-generated noise and incorrect processing (for example residual NMO due to erroneous velocity picks, incomplete migration, or other systematic errors) can accentuate the footprint.

This noise can interfere with the mapping of stratigraphic features and fault patterns, posing a challenge to seismic interpreters working in both exploration and development settings.

To demonstrate the relevance of the phenomenon I show below a gallery of examples from the literature of severe footprint in land data: an amplitude time slice (Figure 1a) and a vertical section (Figure 1b) from a Saudi Arabian case study, some seismic attributes (Figures 2, 3, 4, and 5), and also some modeled streamer data (Figure 6).

Bannagi combo

Figure 1. Amplitude time slice (top, time = 0.44 s) showing footprint in both inline and crossline direction, and amplitude section (bottom) highlighting the effect in the vertical direction. From Al-Bannagi et al. Copyrighted material.

Penobscop_sobel

Figure 2. Edge detection (Sobel filter) on the Penobscot 3D horizon (average time ~= 0.98 s) displaying N-S footprint. From Hall.

F3_shallow_sobel

Figure 3. Edge detection (Sobel filter) on a shallow horizon (average time ~= 0.44 s)  from the F3 Netherlands 3D survey displaying E-W footprint.

Davogustto and Marfurt

Figure 4. Similarity attribute (top , time = 0.6 s), and most positive curvature (bottom, time = 1.3 s), both showing footprint. From Davogustto and Marfurt. Copyrighted material.

Chopra-Larsen

Figure 5. Amplitude time slice (top, time = 1.32 s) the corresponding  coherence section  (bottom) both showing footprint. From Chopra and Larsen. Copyrighted material.

Long et al

Figure 6. Acquisition footprint in the form of low fold striation due to dip streamer acquisition. From Long et al. Copyrighted material.

In my next post I will review (with more examples form literature) some strategies available to either prevent or minimize the footprint with better acquisition parameters and modeling of the stack response; I will also discuss some ways the footprint can be attenuated after the acquisition of the data (with bin regularization/interpolation, dip-steered median filters, and kx ky filters, from simple low-pass to more sophisticated ones) when the above mentioned strategies are not available, due to time/cost constraint or because the interpreter is working with legacy data.

In subsequent posts I will illustrate a workflow to model synthetic acquisition footprint using Python, and how to automatically remove it in the Fourier domain with frequency filters, and then how to remove it from real data.

References

Al-Bannagi et al. 2005 – Acquisition footprint suppression via the truncated SVD technique: Case studies from Saudi Arabia: The Leading Edge, SEG, 24, 832– 834.

Chopra and Larsen,  2000 – Acquisition Footprint, Its Detection and Removal: CSEG Recorder, 25 (8).

Davogusto and Martfurt, 2011 – Footprint Suppression Applied to Legacy Seismic Data Volumes: 31st Annual GCSSEPM Foundation Bob F Perkins Research Conference 2011.

F3 Netherlands open access 3D:  info on SEG Wiki

Hall, 2014 –  Sobel filtering horizons (open source Jupyter Notebook on GitHub).

Long et al., 2004 – On the issue of strike or dip streamer shooting for 3D multi-streamer acquisition: Exploration Geophysics, 35(2), 105-110.

Penobscot open access 3D:  info on SEG Wiki

A study of Ricker wavelets in MS Excel

Posted on by
Reply

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.

sheet 1 input and amplitude spectrum 30 hz

Figure 1

In Figure 2 I show the same Sheet, but for a wavelet of dominant frequency equal to 50 Hz.

sheet 1 input and amplitude spectrum 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

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.

sheet 2 zero phase

Figure 4

sheet 2 linear phase shift 180

Figure 5

sheet 2 constant phase shift 90

Figure 6

sheet 2 constant 90 and linear 180 phase shifts

Figure 7

Finally, Sheet 3 displays a plot of the wavelet, and of the individual frequency components that make it up. Have fun!!!

sheet 3 component sinusoids

Figure 8

References and further playing

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.

Geology in pictures – meanders and oxbow lakes

Posted on by
2

I love meandering rivers. This is a short post with some (I think) great images  of meandering rivers and oxbow lakes.

In my previous post Google Earth and a 5 minute book review: Geology illustrated I reviewed one of my favorite coffee table books: Geology Illustrated by John S. Shelton.

The last view in that post was as close a perspective replica of Figure 137 in the book as I could get using Google Earth, showing the meander belt of the Animas River a few miles from Durango, Colorado.

What a great opportunity to create a short time lapse: a repeat snapshots of the same landscape nearly 50 years apart. This is priceless: 50 years are in many cases next to nothing in geological terms, and yet there are already some significant differences in the meanders in the two images, which I have included together below.

Shelton_137_copyright

Meander belt and floodplain of the Animas River a few miles from Durango, Colorado. From Geology Illustrated, John S. Shelton, 1966. Copyright of the University of Washington Libraries, Special Collections (John Shelton Collection, Shelton 2679).

Meander belt and floodplain of the Animas River as above, as viewed on Google Earth in 2013.

For example, it looks like the meander cutoff in the lower left portion of the image had likely ‘just’ happened in the 60s, whereas at the time the imagery used by Google Earth was acquired, the remnant oxbow lake seems more clearly defined. Another oxbow lake in the center has nearly disappeared, perhaps in part due to human activity.

Below are two pictures of the Fraser River in the Robson Valley that I took in the summer of 2013 during a hike to McBride Peak. This is an awesome example of oxbow lake. We can still see where the river used to run previous to the jump.

horseshoe_panorama

horseshoe_closeup

And this is a great illustration (from Infographic illustrations) showing how these oxbow lakes are created. I love the hands-on feel…ox-bow-lake-infographic

The last is an image of tight meander loop from a previous post: Some photos of Northern British Columbia wildlife and geology.

tight_loop

 

A good divergent color palette for Matlab

Posted on by
16

INTRODUCTION

Before starting my series on perceptual color palettes I thought it was worth mentioning an excellent function I found some time ago on the Matlab File Exchange. The function is called Light and Bartlein Color Maps. It was a Matlab Pick of the week, and it can be used to create four color palettes discussed in the EOS paper by Light and Bartlein. Each of these palettes is suited for a specific task, and the authors claim they are non confusing for viewers with color vision deficiencies.

In the remainder of this post I will showcase one of the palettes, called orange-white-purple, as it is good divergent scheme [1]. With the code below I am going to load the World Topography Matlab demo data, create  the palette and use it to display the data.

%% load World Topography Matlab demo
load topo;

%% create Light Bartlein orange-white-purple diverging scheme
LB=flipud(lbmap(256,'BrownBlue')); % flip it so blue is for negative(ocean)
                                   % and green for positive (land)

%% plot map
fig2 = figure;
imagesc(flipud(topo));
axis equal
axis tight
axis off
set(fig2,'Position',[720 400 980 580]);
title(' Non-symmetric divergent orange-white-purple palette','Color',...
    'k','FontSize',12,'FontWeight','demi');
colormap(LB);
colorbar;

And here is the result below. I like this color scheme better than many othera for divergent data. One only issue in the figure, although not inherently due to the palette itself [2], is that the centre of the palette is not at the zero. This is a problem since the zero is such an important element in ratio data, in this case representing sea level.

MAKING THE PALETTE SYMMETRIC AROUND THE ZERO

The problem fortunately can be easily fixed by clipping the data limit to a symmetric range. In Matlab this has to be done programmatically, and rather than going about it with trial and error I like to do it automatically with the code below:

Continue reading →