xarray
¶First, the usual imports.
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
The seismic data is too big for GitHub. The easiest thing to do is to download the HDF5 file:
and put it into the data
folder. Then you can just do this:
import h5py
h5f = h5py.File('data/Penobscot.h5','r')
seismic_data = h5f['amplitude'][:]
h5f.close()
seismic_data.shape
(251, 601, 481)
From the Terranubis site:
Corners, from Moreno and Galeano 2019 (I have not double-checked this, eg the CRS).
Corners | North coordinate | East coordinate |
---|---|---|
P1 | 4890130.78 | 731998.45 |
P2 | 4896722.05 | 728527.16 |
P3 | 4895689.47 | 742633.34 |
P4 | 4802330.74 | 739162.05 |
Can't figure out how to do the skewed grid, so I'm going to completely invent the coordinates:
import xarray as xr
t, x, i = seismic_data.shape
xs = 25.0 * np.arange(i) + 100000
ys = 12.5 * np.arange(x) + 200000
tslices = np.arange(t) * 0.004
seismic = xr.DataArray(seismic_data,
name='Amplitude',
coords=[tslices, ys, xs],
dims=['twt', 'y', 'x'],
)
seismic.loc[0.500].plot()
<matplotlib.collections.QuadMesh at 0x7fe6ed05b748>
del seismic_data
from shapely.geometry import LineString
line = LineString([(101000, 207200), (106000, 200200), (111000, 207200)])
line.length
17204.650534085253
Let's see how the line looks on the seismic data:
seismic.loc[0.500].plot()
plt.plot(*np.asarray(line.coords.xy))
[<matplotlib.lines.Line2D at 0x7fe7192eaa58>]
Define the distance between traces.
trace_distance = 100
n_traces = np.ceil(line.length / trace_distance)
trace_locs = trace_distance * np.arange(n_traces + 1)
print(f"Interpolating {n_traces} traces.")
Interpolating 173.0 traces.
from tqdm import tqdm
points = [line.interpolate(loc) for loc in trace_locs]
traces = np.stack([seismic.interp(x=p.x, y=p.y) for p in tqdm(points)])
100%|██████████| 174/174 [02:23<00:00, 1.25it/s]
Make a new xarray
object:
arbline = xr.DataArray(traces.T,
name='Amplitude',
coords=[tslices, trace_locs],
dims=['twt', 'offset'],
)
arbline.plot(figsize=(12, 5))
plt.gca().invert_yaxis()
Code and text © 2020 Agile Geoscience — CC-BY — Have fun!
Data: Penobscot, licensed CC-BY-SA, Nova Scotia Department of Energy https://opendtect.org/osr/.