Merging Seismic Data Cubes¶

Often we receive seismic data cubes which we wish to merge that have different geometries. This workflow will guide you through the process of merging two such cubes.

In [2]:

from segysak import create3d_dataset
from scipy.optimize import curve_fit
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
from matplotlib.transforms import Affine2D
from segysak import open_seisnc, create_seismic_dataset

from dask.distributed import Client

# start a dask client for processing
client = Client()

from segysak import create3d_dataset from scipy.optimize import curve_fit import numpy as np import xarray as xr import matplotlib.pyplot as plt from matplotlib.transforms import Affine2D from segysak import open_seisnc, create_seismic_dataset from dask.distributed import Client # start a dask client for processing client = Client()

Creating some synthetic seismic surveys for us to play with. The online documentation uses a smaller example, but you can stress dask by commenting the dimension lines for smaller data and uncommenting the large desktop size volumes.

In [3]:

# large desktop size merge
# i1, i2, iN = 100, 750, 651
# x1, x2, xN = 300, 1250, 951
# t1, t2, tN = 0, 1848, 463
# trans_x, trans_y = 10000, 2000

# online example size merge
i1, i2, iN = 650, 750, 101
x1, x2, xN = 1150, 1250, 101
t1, t2, tN = 0, 200, 201
trans_x, trans_y = 9350, 1680

iline = np.linspace(i1, i2, iN, dtype=int)
xline = np.linspace(x1, x2, xN, dtype=int)
twt = np.linspace(t1, t2, tN, dtype=int)

affine_survey1 = Affine2D()
affine_survey1.rotate_deg(20).translate(trans_x, trans_y)

# create an empty dataset and empty coordinates
survey_1 = create_seismic_dataset(twt=twt, iline=iline, xline=xline)
survey_1["cdp_x"] = (("iline", "xline"), np.empty((iN, xN)))
survey_1["cdp_y"] = (("iline", "xline"), np.empty((iN, xN)))

stacked_iline = survey_1.iline.broadcast_like(survey_1.cdp_x).stack(
    {"ravel": (..., "xline")}
)
stacked_xline = survey_1.xline.broadcast_like(survey_1.cdp_x).stack(
    {"ravel": (..., "xline")}
)

points = np.dstack([stacked_iline, stacked_xline])
points_xy = affine_survey1.transform(points[0])
# put xy back in stacked and unstack into survey_1
survey_1["cdp_x"] = stacked_iline.copy(data=points_xy[:, 0]).unstack()
survey_1["cdp_y"] = stacked_iline.copy(data=points_xy[:, 1]).unstack()

# and the same for survey_2 but with a different affine to survey 1
affine_survey2 = Affine2D()
affine_survey2.rotate_deg(120).translate(11000, 3000)

survey_2 = create_seismic_dataset(twt=twt, iline=iline, xline=xline)
survey_2["cdp_x"] = (("iline", "xline"), np.empty((iN, xN)))
survey_2["cdp_y"] = (("iline", "xline"), np.empty((iN, xN)))

stacked_iline = survey_1.iline.broadcast_like(survey_1.cdp_x).stack(
    {"ravel": (..., "xline")}
)
stacked_xline = survey_1.xline.broadcast_like(survey_1.cdp_x).stack(
    {"ravel": (..., "xline")}
)

points = np.dstack([stacked_iline, stacked_xline])
points_xy = affine_survey2.transform(points[0])
# put xy back in stacked and unstack into survey_1
survey_2["cdp_x"] = stacked_iline.copy(data=points_xy[:, 0]).unstack()
survey_2["cdp_y"] = stacked_iline.copy(data=points_xy[:, 1]).unstack()

# large desktop size merge # i1, i2, iN = 100, 750, 651 # x1, x2, xN = 300, 1250, 951 # t1, t2, tN = 0, 1848, 463 # trans_x, trans_y = 10000, 2000 # online example size merge i1, i2, iN = 650, 750, 101 x1, x2, xN = 1150, 1250, 101 t1, t2, tN = 0, 200, 201 trans_x, trans_y = 9350, 1680 iline = np.linspace(i1, i2, iN, dtype=int) xline = np.linspace(x1, x2, xN, dtype=int) twt = np.linspace(t1, t2, tN, dtype=int) affine_survey1 = Affine2D() affine_survey1.rotate_deg(20).translate(trans_x, trans_y) # create an empty dataset and empty coordinates survey_1 = create_seismic_dataset(twt=twt, iline=iline, xline=xline) survey_1["cdp_x"] = (("iline", "xline"), np.empty((iN, xN))) survey_1["cdp_y"] = (("iline", "xline"), np.empty((iN, xN))) stacked_iline = survey_1.iline.broadcast_like(survey_1.cdp_x).stack( {"ravel": (..., "xline")} ) stacked_xline = survey_1.xline.broadcast_like(survey_1.cdp_x).stack( {"ravel": (..., "xline")} ) points = np.dstack([stacked_iline, stacked_xline]) points_xy = affine_survey1.transform(points[0]) # put xy back in stacked and unstack into survey_1 survey_1["cdp_x"] = stacked_iline.copy(data=points_xy[:, 0]).unstack() survey_1["cdp_y"] = stacked_iline.copy(data=points_xy[:, 1]).unstack() # and the same for survey_2 but with a different affine to survey 1 affine_survey2 = Affine2D() affine_survey2.rotate_deg(120).translate(11000, 3000) survey_2 = create_seismic_dataset(twt=twt, iline=iline, xline=xline) survey_2["cdp_x"] = (("iline", "xline"), np.empty((iN, xN))) survey_2["cdp_y"] = (("iline", "xline"), np.empty((iN, xN))) stacked_iline = survey_1.iline.broadcast_like(survey_1.cdp_x).stack( {"ravel": (..., "xline")} ) stacked_xline = survey_1.xline.broadcast_like(survey_1.cdp_x).stack( {"ravel": (..., "xline")} ) points = np.dstack([stacked_iline, stacked_xline]) points_xy = affine_survey2.transform(points[0]) # put xy back in stacked and unstack into survey_1 survey_2["cdp_x"] = stacked_iline.copy(data=points_xy[:, 0]).unstack() survey_2["cdp_y"] = stacked_iline.copy(data=points_xy[:, 1]).unstack()

Let's check that there is a bit of overlap between the two surveys.

In [4]:

# plotting every 10th line

plt.figure(figsize=(10, 10))
survey_1_plot = plt.plot(
    survey_1.cdp_x.values[::10, ::10], survey_1.cdp_y.values[::10, ::10], color="grey"
)
survey_2_plot = plt.plot(
    survey_2.cdp_x.values[::10, ::10], survey_2.cdp_y.values[::10, ::10], color="blue"
)
# plt.aspect("equal")
plt.legend([survey_1_plot[0], survey_2_plot[0]], ["survey_1", "survey_2"])

# plotting every 10th line plt.figure(figsize=(10, 10)) survey_1_plot = plt.plot( survey_1.cdp_x.values[::10, ::10], survey_1.cdp_y.values[::10, ::10], color="grey" ) survey_2_plot = plt.plot( survey_2.cdp_x.values[::10, ::10], survey_2.cdp_y.values[::10, ::10], color="blue" ) # plt.aspect("equal") plt.legend([survey_1_plot[0], survey_2_plot[0]], ["survey_1", "survey_2"])

Out[4]:

<matplotlib.legend.Legend at 0x7f350ca7a980>

Let's output these two datasets to disk so we can use lazy loading from seisnc. If the files are large, it's good practice to do this, because you will probably run out of memory.

We are going to fill survey one with a values of 1 and survey 2 with a value of 2, just for simplicity in this example. We also tell the dtype to be np.float32 to reduce memory usage.

In [5]:

# save out with new geometry
survey_1["data"] = (
    ("iline", "xline", "twt"),
    np.full((iN, xN, tN), 1, dtype=np.float32),
)
survey_1.seisio.to_netcdf("data/survey_1.seisnc")
del survey_1
survey_2["data"] = (
    ("iline", "xline", "twt"),
    np.full((iN, xN, tN), 2, dtype=np.float32),
)
survey_2.seisio.to_netcdf("data/survey_2.seisnc")
del survey_2

# save out with new geometry survey_1["data"] = ( ("iline", "xline", "twt"), np.full((iN, xN, tN), 1, dtype=np.float32), ) survey_1.seisio.to_netcdf("data/survey_1.seisnc") del survey_1 survey_2["data"] = ( ("iline", "xline", "twt"), np.full((iN, xN, tN), 2, dtype=np.float32), ) survey_2.seisio.to_netcdf("data/survey_2.seisnc") del survey_2

How to merge two different geometries¶

Let us reimport the surveys we created but in a chunked (lazy) way.

In [6]:

survey_1 = open_seisnc("data/survey_1.seisnc", chunks=dict(iline=10, xline=10, twt=100))
survey_2 = open_seisnc("data/survey_2.seisnc", chunks=dict(iline=10, xline=10, twt=100))

survey_1 = open_seisnc("data/survey_1.seisnc", chunks=dict(iline=10, xline=10, twt=100)) survey_2 = open_seisnc("data/survey_2.seisnc", chunks=dict(iline=10, xline=10, twt=100))

/tmp/ipykernel_2077/814101917.py:1: DeprecationWarning: open_seisnc will be removed in v0.6, please use the Xarray engine ds = xr.open_dataset(netcdf_file) method instead.
  survey_1 = open_seisnc("data/survey_1.seisnc", chunks=dict(iline=10, xline=10, twt=100))
/tmp/ipykernel_2077/814101917.py:2: DeprecationWarning: open_seisnc will be removed in v0.6, please use the Xarray engine ds = xr.open_dataset(netcdf_file) method instead.
  survey_2 = open_seisnc("data/survey_2.seisnc", chunks=dict(iline=10, xline=10, twt=100))

Check that the survey is chunked by looking at the printout for our datasets. Lazy and chunked data will have a dask.array<chunksize= where the values are usually displayed.

In [7]:

survey_1

survey_1

Out[7]:

<xarray.Dataset> Size: 8MB
Dimensions:  (iline: 101, xline: 101, twt: 201)
Coordinates:
  * iline    (iline) int64 808B 650 651 652 653 654 655 ... 746 747 748 749 750
  * xline    (xline) int64 808B 1150 1151 1152 1153 1154 ... 1247 1248 1249 1250
  * twt      (twt) int64 2kB 0 1 2 3 4 5 6 7 ... 193 194 195 196 197 198 199 200
    cdp_x    (iline, xline) float64 82kB dask.array<chunksize=(10, 10), meta=np.ndarray>
    cdp_y    (iline, xline) float64 82kB dask.array<chunksize=(10, 10), meta=np.ndarray>
Data variables:
    data     (iline, xline, twt) float32 8MB dask.array<chunksize=(10, 10, 100), meta=np.ndarray>
Attributes: (12/17)
    text:                
    ns:                  None
    sample_rate:         None
    measurement_system:  None
    d3_domain:           None
    epsg:                None
    ...                  ...
    percentiles:         None
    coord_scalar:        None
    coord_scaled:        None
    dimensions:          None
    vert_dimension:      None
    vert_domain:         None

xarray.Dataset

• Dimensions:
  • iline: 101
  • xline: 101
  • twt: 201
• Coordinates: (5)
  • iline
    (iline)
    int64
    650 651 652 653 ... 747 748 749 750

    array([650, 651, 652, 653, 654, 655, 656, 657, 658, 659, 660, 661, 662, 663,
           664, 665, 666, 667, 668, 669, 670, 671, 672, 673, 674, 675, 676, 677,
           678, 679, 680, 681, 682, 683, 684, 685, 686, 687, 688, 689, 690, 691,
           692, 693, 694, 695, 696, 697, 698, 699, 700, 701, 702, 703, 704, 705,
           706, 707, 708, 709, 710, 711, 712, 713, 714, 715, 716, 717, 718, 719,
           720, 721, 722, 723, 724, 725, 726, 727, 728, 729, 730, 731, 732, 733,
           734, 735, 736, 737, 738, 739, 740, 741, 742, 743, 744, 745, 746, 747,
           748, 749, 750])

  • xline
    (xline)
    int64
    1150 1151 1152 ... 1248 1249 1250

    array([1150, 1151, 1152, 1153, 1154, 1155, 1156, 1157, 1158, 1159, 1160, 1161,
           1162, 1163, 1164, 1165, 1166, 1167, 1168, 1169, 1170, 1171, 1172, 1173,
           1174, 1175, 1176, 1177, 1178, 1179, 1180, 1181, 1182, 1183, 1184, 1185,
           1186, 1187, 1188, 1189, 1190, 1191, 1192, 1193, 1194, 1195, 1196, 1197,
           1198, 1199, 1200, 1201, 1202, 1203, 1204, 1205, 1206, 1207, 1208, 1209,
           1210, 1211, 1212, 1213, 1214, 1215, 1216, 1217, 1218, 1219, 1220, 1221,
           1222, 1223, 1224, 1225, 1226, 1227, 1228, 1229, 1230, 1231, 1232, 1233,
           1234, 1235, 1236, 1237, 1238, 1239, 1240, 1241, 1242, 1243, 1244, 1245,
           1246, 1247, 1248, 1249, 1250])

  • twt
    (twt)
    int64
    0 1 2 3 4 5 ... 196 197 198 199 200

    array([  0,   1,   2, ..., 198, 199, 200])

  • cdp_x
    (iline, xline)
    float64
    dask.array<chunksize=(10, 10), meta=np.ndarray>

    Array Chunk Bytes 79.70 kiB 800 B Shape (101, 101) (10, 10) Dask graph 121 chunks in 2 graph layers Data type float64 numpy.ndarray

  • cdp_y
    (iline, xline)
    float64
    dask.array<chunksize=(10, 10), meta=np.ndarray>

    Array Chunk Bytes 79.70 kiB 800 B Shape (101, 101) (10, 10) Dask graph 121 chunks in 2 graph layers Data type float64 numpy.ndarray

• Data variables: (1)
  • data
    (iline, xline, twt)
    float32
    dask.array<chunksize=(10, 10, 100), meta=np.ndarray>

    Array Chunk Bytes 7.82 MiB 39.06 kiB Shape (101, 101, 201) (10, 10, 100) Dask graph 363 chunks in 2 graph layers Data type float32 numpy.ndarray

• Indexes: (3)
  • iline
    PandasIndex

    PandasIndex(Index([650, 651, 652, 653, 654, 655, 656, 657, 658, 659,
           ...
           741, 742, 743, 744, 745, 746, 747, 748, 749, 750],
          dtype='int64', name='iline', length=101))

  • xline
    PandasIndex

    PandasIndex(Index([1150, 1151, 1152, 1153, 1154, 1155, 1156, 1157, 1158, 1159,
           ...
           1241, 1242, 1243, 1244, 1245, 1246, 1247, 1248, 1249, 1250],
          dtype='int64', name='xline', length=101))

  • twt
    PandasIndex

    PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
           ...
           191, 192, 193, 194, 195, 196, 197, 198, 199, 200],
          dtype='int64', name='twt', length=201))

• Attributes: (17)

  text :

  ns :

  None

  sample_rate :

  None

  measurement_system :

  None

  d3_domain :

  None

  epsg :

  None

  corner_points :

  None

  corner_points_xy :

  None

  source_file :

  None

  srd :

  None

  datatype :

  None

  percentiles :

  None

  coord_scalar :

  None

  coord_scaled :

  None

  dimensions :

  None

  vert_dimension :

  None

  vert_domain :

  None

It will help with our other functions if we load the cdp_x and cdp_y locations into memory. This can be done using the load() method.

In [8]:

survey_1["cdp_x"] = survey_1["cdp_x"].load()
survey_1["cdp_y"] = survey_1["cdp_y"].load()
survey_2["cdp_x"] = survey_2["cdp_x"].load()
survey_2["cdp_y"] = survey_2["cdp_y"].load()

survey_1["cdp_x"] = survey_1["cdp_x"].load() survey_1["cdp_y"] = survey_1["cdp_y"].load() survey_2["cdp_x"] = survey_2["cdp_x"].load() survey_2["cdp_y"] = survey_2["cdp_y"].load()

In [9]:

survey_1

survey_1

Out[9]:

<xarray.Dataset> Size: 8MB
Dimensions:  (iline: 101, xline: 101, twt: 201)
Coordinates:
  * iline    (iline) int64 808B 650 651 652 653 654 655 ... 746 747 748 749 750
  * xline    (xline) int64 808B 1150 1151 1152 1153 1154 ... 1247 1248 1249 1250
  * twt      (twt) int64 2kB 0 1 2 3 4 5 6 7 ... 193 194 195 196 197 198 199 200
    cdp_x    (iline, xline) float64 82kB 9.567e+03 9.567e+03 ... 9.627e+03
    cdp_y    (iline, xline) float64 82kB 2.983e+03 2.984e+03 ... 3.111e+03
Data variables:
    data     (iline, xline, twt) float32 8MB dask.array<chunksize=(10, 10, 100), meta=np.ndarray>
Attributes: (12/17)
    text:                
    ns:                  None
    sample_rate:         None
    measurement_system:  None
    d3_domain:           None
    epsg:                None
    ...                  ...
    percentiles:         None
    coord_scalar:        None
    coord_scaled:        None
    dimensions:          None
    vert_dimension:      None
    vert_domain:         None

xarray.Dataset

• Dimensions:
  • iline: 101
  • xline: 101
  • twt: 201
• Coordinates: (5)
  • iline
    (iline)
    int64
    650 651 652 653 ... 747 748 749 750

    array([650, 651, 652, 653, 654, 655, 656, 657, 658, 659, 660, 661, 662, 663,
           664, 665, 666, 667, 668, 669, 670, 671, 672, 673, 674, 675, 676, 677,
           678, 679, 680, 681, 682, 683, 684, 685, 686, 687, 688, 689, 690, 691,
           692, 693, 694, 695, 696, 697, 698, 699, 700, 701, 702, 703, 704, 705,
           706, 707, 708, 709, 710, 711, 712, 713, 714, 715, 716, 717, 718, 719,
           720, 721, 722, 723, 724, 725, 726, 727, 728, 729, 730, 731, 732, 733,
           734, 735, 736, 737, 738, 739, 740, 741, 742, 743, 744, 745, 746, 747,
           748, 749, 750])

  • xline
    (xline)
    int64
    1150 1151 1152 ... 1248 1249 1250

    array([1150, 1151, 1152, 1153, 1154, 1155, 1156, 1157, 1158, 1159, 1160, 1161,
           1162, 1163, 1164, 1165, 1166, 1167, 1168, 1169, 1170, 1171, 1172, 1173,
           1174, 1175, 1176, 1177, 1178, 1179, 1180, 1181, 1182, 1183, 1184, 1185,
           1186, 1187, 1188, 1189, 1190, 1191, 1192, 1193, 1194, 1195, 1196, 1197,
           1198, 1199, 1200, 1201, 1202, 1203, 1204, 1205, 1206, 1207, 1208, 1209,
           1210, 1211, 1212, 1213, 1214, 1215, 1216, 1217, 1218, 1219, 1220, 1221,
           1222, 1223, 1224, 1225, 1226, 1227, 1228, 1229, 1230, 1231, 1232, 1233,
           1234, 1235, 1236, 1237, 1238, 1239, 1240, 1241, 1242, 1243, 1244, 1245,
           1246, 1247, 1248, 1249, 1250])

  • twt
    (twt)
    int64
    0 1 2 3 4 5 ... 196 197 198 199 200

    array([  0,   1,   2, ..., 198, 199, 200])

  • cdp_x
    (iline, xline)
    float64
    9.567e+03 9.567e+03 ... 9.627e+03

    array([[9567.47703869, 9567.13501854, 9566.7929984 , ..., 9533.95906464,
            9533.6170445 , 9533.27502435],
           [9568.41673131, 9568.07471116, 9567.73269102, ..., 9534.89875726,
            9534.55673712, 9534.21471697],
           [9569.35642393, 9569.01440378, 9568.67238364, ..., 9535.83844988,
            9535.49642974, 9535.1544096 ],
           ...,
           [9659.56691552, 9659.22489538, 9658.88287524, ..., 9626.04894148,
            9625.70692133, 9625.36490119],
           [9660.50660814, 9660.164588  , 9659.82256786, ..., 9626.9886341 ,
            9626.64661395, 9626.30459381],
           [9661.44630076, 9661.10428062, 9660.76226048, ..., 9627.92832672,
            9627.58630658, 9627.24428643]])

  • cdp_y
    (iline, xline)
    float64
    2.983e+03 2.984e+03 ... 3.111e+03

    array([[2982.95960707, 2983.89929969, 2984.83899231, ..., 3075.0494839 ,
            3075.98917652, 3076.92886914],
           [2983.30162721, 2984.24131983, 2985.18101245, ..., 3075.39150405,
            3076.33119667, 3077.27088929],
           [2983.64364735, 2984.58333997, 2985.52303259, ..., 3075.73352419,
            3076.67321681, 3077.61290943],
           ...,
           [3016.47758111, 3017.41727373, 3018.35696635, ..., 3108.56745795,
            3109.50715057, 3110.44684319],
           [3016.81960125, 3017.75929388, 3018.6989865 , ..., 3108.90947809,
            3109.84917071, 3110.78886333],
           [3017.1616214 , 3018.10131402, 3019.04100664, ..., 3109.25149824,
            3110.19119086, 3111.13088348]])

• Data variables: (1)
  • data
    (iline, xline, twt)
    float32
    dask.array<chunksize=(10, 10, 100), meta=np.ndarray>

    Array Chunk Bytes 7.82 MiB 39.06 kiB Shape (101, 101, 201) (10, 10, 100) Dask graph 363 chunks in 2 graph layers Data type float32 numpy.ndarray

• Indexes: (3)
  • iline
    PandasIndex

    PandasIndex(Index([650, 651, 652, 653, 654, 655, 656, 657, 658, 659,
           ...
           741, 742, 743, 744, 745, 746, 747, 748, 749, 750],
          dtype='int64', name='iline', length=101))

  • xline
    PandasIndex

    PandasIndex(Index([1150, 1151, 1152, 1153, 1154, 1155, 1156, 1157, 1158, 1159,
           ...
           1241, 1242, 1243, 1244, 1245, 1246, 1247, 1248, 1249, 1250],
          dtype='int64', name='xline', length=101))

  • twt
    PandasIndex

    PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
           ...
           191, 192, 193, 194, 195, 196, 197, 198, 199, 200],
          dtype='int64', name='twt', length=201))

• Attributes: (17)

  text :

  ns :

  None

  sample_rate :

  None

  measurement_system :

  None

  d3_domain :

  None

  epsg :

  None

  corner_points :

  None

  corner_points_xy :

  None

  source_file :

  None

  srd :

  None

  datatype :

  None

  percentiles :

  None

  coord_scalar :

  None

  coord_scaled :

  None

  dimensions :

  None

  vert_dimension :

  None

  vert_domain :

  None

We will also need an affine transform which converts from il/xl on survey 2 to il/xl on survey 1. This will relate the two grids to each other. Fortunately affine transforms can be added together which makes this pretty straight forward.

In [10]:

ilxl2_to_ilxl1 = (
    survey_2.seis.get_affine_transform()
    + survey_1.seis.get_affine_transform().inverted()
)

ilxl2_to_ilxl1 = ( survey_2.seis.get_affine_transform() + survey_1.seis.get_affine_transform().inverted() )

Let's check the transform to see if it makes sense. The colouring of our values by inline should be the same after the transform is applied. Cool.

In [11]:

# plotting every 50th line

n = 5

fig, axs = plt.subplots(ncols=2, figsize=(20, 10))
survey_1_plot = axs[0].scatter(
    survey_1.cdp_x.values[::n, ::n],
    survey_1.cdp_y.values[::n, ::n],
    c=survey_1.iline.broadcast_like(survey_1.cdp_x).values[::n, ::n],
    s=20,
)
survey_2_plot = axs[0].scatter(
    survey_2.cdp_x.values[::n, ::n],
    survey_2.cdp_y.values[::n, ::n],
    c=survey_2.iline.broadcast_like(survey_2.cdp_x).values[::n, ::n],
    s=20,
)
plt.colorbar(survey_1_plot, ax=axs[0])
axs[0].set_aspect("equal")

axs[0].set_title("inline numbering before transform")

new_ilxl = ilxl2_to_ilxl1.transform(
    np.dstack(
        [
            survey_2.iline.broadcast_like(survey_2.cdp_x).values.ravel(),
            survey_2.xline.broadcast_like(survey_2.cdp_x).values.ravel(),
        ]
    )[0]
).reshape((survey_2.iline.size, survey_2.xline.size, 2))

survey_2.seis.calc_corner_points()
s2_corner_points_in_s1 = ilxl2_to_ilxl1.transform(survey_2.attrs["corner_points"])

print("Min Iline:", s2_corner_points_in_s1[:, 0].min(), survey_1.iline.min().values)
min_iline_combined = min(
    s2_corner_points_in_s1[:, 0].min(), survey_1.iline.min().values
)
print("Max Iline:", s2_corner_points_in_s1[:, 0].max(), survey_1.iline.max().values)
max_iline_combined = max(
    s2_corner_points_in_s1[:, 0].max(), survey_1.iline.max().values
)

print("Min Xline:", s2_corner_points_in_s1[:, 1].min(), survey_1.xline.min().values)
min_xline_combined = min(
    s2_corner_points_in_s1[:, 1].min(), survey_1.xline.min().values
)
print("Max Xline:", s2_corner_points_in_s1[:, 1].max(), survey_1.xline.max().values)
max_xline_combined = max(
    s2_corner_points_in_s1[:, 1].max(), survey_1.xline.max().values
)
print(min_iline_combined, max_iline_combined, min_xline_combined, max_xline_combined)

survey_1_plot = axs[1].scatter(
    survey_1.cdp_x.values[::n, ::n],
    survey_1.cdp_y.values[::n, ::n],
    c=survey_1.iline.broadcast_like(survey_1.cdp_x).values[::n, ::n],
    vmin=min_iline_combined,
    vmax=max_iline_combined,
    s=20,
    cmap="hsv",
)


survey_2_plot = axs[1].scatter(
    survey_2.cdp_x.values[::n, ::n],
    survey_2.cdp_y.values[::n, ::n],
    c=new_ilxl[::n, ::n, 0],
    vmin=min_iline_combined,
    vmax=max_iline_combined,
    s=20,
    cmap="hsv",
)
plt.colorbar(survey_1_plot, ax=axs[1])

axs[1].set_title("inline numbering after transform")
axs[1].set_aspect("equal")

# plotting every 50th line n = 5 fig, axs = plt.subplots(ncols=2, figsize=(20, 10)) survey_1_plot = axs[0].scatter( survey_1.cdp_x.values[::n, ::n], survey_1.cdp_y.values[::n, ::n], c=survey_1.iline.broadcast_like(survey_1.cdp_x).values[::n, ::n], s=20, ) survey_2_plot = axs[0].scatter( survey_2.cdp_x.values[::n, ::n], survey_2.cdp_y.values[::n, ::n], c=survey_2.iline.broadcast_like(survey_2.cdp_x).values[::n, ::n], s=20, ) plt.colorbar(survey_1_plot, ax=axs[0]) axs[0].set_aspect("equal") axs[0].set_title("inline numbering before transform") new_ilxl = ilxl2_to_ilxl1.transform( np.dstack( [ survey_2.iline.broadcast_like(survey_2.cdp_x).values.ravel(), survey_2.xline.broadcast_like(survey_2.cdp_x).values.ravel(), ] )[0] ).reshape((survey_2.iline.size, survey_2.xline.size, 2)) survey_2.seis.calc_corner_points() s2_corner_points_in_s1 = ilxl2_to_ilxl1.transform(survey_2.attrs["corner_points"]) print("Min Iline:", s2_corner_points_in_s1[:, 0].min(), survey_1.iline.min().values) min_iline_combined = min( s2_corner_points_in_s1[:, 0].min(), survey_1.iline.min().values ) print("Max Iline:", s2_corner_points_in_s1[:, 0].max(), survey_1.iline.max().values) max_iline_combined = max( s2_corner_points_in_s1[:, 0].max(), survey_1.iline.max().values ) print("Min Xline:", s2_corner_points_in_s1[:, 1].min(), survey_1.xline.min().values) min_xline_combined = min( s2_corner_points_in_s1[:, 1].min(), survey_1.xline.min().values ) print("Max Xline:", s2_corner_points_in_s1[:, 1].max(), survey_1.xline.max().values) max_xline_combined = max( s2_corner_points_in_s1[:, 1].max(), survey_1.xline.max().values ) print(min_iline_combined, max_iline_combined, min_xline_combined, max_xline_combined) survey_1_plot = axs[1].scatter( survey_1.cdp_x.values[::n, ::n], survey_1.cdp_y.values[::n, ::n], c=survey_1.iline.broadcast_like(survey_1.cdp_x).values[::n, ::n], vmin=min_iline_combined, vmax=max_iline_combined, s=20, cmap="hsv", ) survey_2_plot = axs[1].scatter( survey_2.cdp_x.values[::n, ::n], survey_2.cdp_y.values[::n, ::n], c=new_ilxl[::n, ::n, 0], vmin=min_iline_combined, vmax=max_iline_combined, s=20, cmap="hsv", ) plt.colorbar(survey_1_plot, ax=axs[1]) axs[1].set_title("inline numbering after transform") axs[1].set_aspect("equal")

Min Iline: 640.7135889711899 650
Max Iline: 756.5591820391251 750
Min Xline: 1099.1258403243169 1150
Max Xline: 1214.9714333922225 1250
640.7135889711899 756.5591820391251 1099.1258403243169 1250

So we now have a transform that can convert ilxl from survey_2 to ilxl from survey_1 and we need to create a combined geometry for the merge. The dims need to cover the maximum range (inline, crossline) of two surveys. We are also going to use survey 1 as the base as we don't want to have to resample two cubes, although if you have a prefered geometry you could do that (custom affine transforms per survey required).

In [12]:

# create new dims
iline_step = survey_1.iline.diff("iline").mean().values
xline_step = survey_1.xline.diff("xline").mean().values

# create the new dimensions
new_iline_dim = np.arange(
    int(min_iline_combined // iline_step) * iline_step,
    int(max_iline_combined) + iline_step * 2,
    iline_step,
    dtype=np.int32,
)
new_xline_dim = np.arange(
    int(min_xline_combined // xline_step) * xline_step,
    int(max_xline_combined) + xline_step * 2,
    xline_step,
    dtype=np.int32,
)

# create a new empty dataset and a blank cdp_x for dims broadcasting
survey_comb = create_seismic_dataset(
    iline=new_iline_dim, xline=new_xline_dim, twt=survey_1.twt
)
survey_comb["cdp_x"] = (
    ("iline", "xline"),
    np.empty((new_iline_dim.size, new_xline_dim.size)),
)

# calculate the x and y using the survey 1 affine which is our base grid and reshape to the dataset grid
new_cdp_xy = affine_survey1.transform(
    np.dstack(
        [
            survey_comb.iline.broadcast_like(survey_comb.cdp_x).values.ravel(),
            survey_comb.xline.broadcast_like(survey_comb.cdp_x).values.ravel(),
        ]
    )[0]
)
new_cdp_xy_grid = new_cdp_xy.reshape((new_iline_dim.size, new_xline_dim.size, 2))

# plot to check
plt.figure(figsize=(10, 10))

survey_comb_plot = plt.plot(
    new_cdp_xy_grid[:, :, 0], new_cdp_xy_grid[:, :, 1], color="red", alpha=0.5
)
survey_1_plot = plt.plot(
    survey_1.cdp_x.values[::10, ::10], survey_1.cdp_y.values[::10, ::10], color="grey"
)
survey_2_plot = plt.plot(
    survey_2.cdp_x.values[::10, ::10], survey_2.cdp_y.values[::10, ::10], color="blue"
)

plt.legend(
    [survey_1_plot[0], survey_2_plot[0], survey_comb_plot[0]],
    ["survey_1", "survey_2", "survey_merged"],
)

# put the new x and y in the empty dataset
survey_comb["cdp_x"] = (("iline", "xline"), new_cdp_xy_grid[..., 0])
survey_comb["cdp_y"] = (("iline", "xline"), new_cdp_xy_grid[..., 1])
survey_comb = survey_comb.set_coords(("cdp_x", "cdp_y"))

# create new dims iline_step = survey_1.iline.diff("iline").mean().values xline_step = survey_1.xline.diff("xline").mean().values # create the new dimensions new_iline_dim = np.arange( int(min_iline_combined // iline_step) * iline_step, int(max_iline_combined) + iline_step * 2, iline_step, dtype=np.int32, ) new_xline_dim = np.arange( int(min_xline_combined // xline_step) * xline_step, int(max_xline_combined) + xline_step * 2, xline_step, dtype=np.int32, ) # create a new empty dataset and a blank cdp_x for dims broadcasting survey_comb = create_seismic_dataset( iline=new_iline_dim, xline=new_xline_dim, twt=survey_1.twt ) survey_comb["cdp_x"] = ( ("iline", "xline"), np.empty((new_iline_dim.size, new_xline_dim.size)), ) # calculate the x and y using the survey 1 affine which is our base grid and reshape to the dataset grid new_cdp_xy = affine_survey1.transform( np.dstack( [ survey_comb.iline.broadcast_like(survey_comb.cdp_x).values.ravel(), survey_comb.xline.broadcast_like(survey_comb.cdp_x).values.ravel(), ] )[0] ) new_cdp_xy_grid = new_cdp_xy.reshape((new_iline_dim.size, new_xline_dim.size, 2)) # plot to check plt.figure(figsize=(10, 10)) survey_comb_plot = plt.plot( new_cdp_xy_grid[:, :, 0], new_cdp_xy_grid[:, :, 1], color="red", alpha=0.5 ) survey_1_plot = plt.plot( survey_1.cdp_x.values[::10, ::10], survey_1.cdp_y.values[::10, ::10], color="grey" ) survey_2_plot = plt.plot( survey_2.cdp_x.values[::10, ::10], survey_2.cdp_y.values[::10, ::10], color="blue" ) plt.legend( [survey_1_plot[0], survey_2_plot[0], survey_comb_plot[0]], ["survey_1", "survey_2", "survey_merged"], ) # put the new x and y in the empty dataset survey_comb["cdp_x"] = (("iline", "xline"), new_cdp_xy_grid[..., 0]) survey_comb["cdp_y"] = (("iline", "xline"), new_cdp_xy_grid[..., 1]) survey_comb = survey_comb.set_coords(("cdp_x", "cdp_y"))

Resampling survey_2 to survey_1 geometry¶

Resampling one dataset to another requires us to tell Xarray/dask where we want the new traces to be. First we can convert the x and y coordinates of the combined survey to the il xl of survey 2 using the affine transform of survey 2. This transform works il/xl to x and y and therefore we need it inverted.

In [13]:

survey_2_new_ilxl_loc = affine_survey2.inverted().transform(new_cdp_xy)

survey_2_new_ilxl_loc = affine_survey2.inverted().transform(new_cdp_xy)

We then need to create a sampling DataArray. If we passed the iline and cross line locations from the previous cell unfortunately Xarray would broadcast the inline and xline values against each other. The simplest way is to create a new stacked flat dimension which we can unstack into a cube later.

We also need to give it unique names so xarray doesn't confuse the new dimensions with the old dimensions.

In [14]:

flat = (
    survey_comb.rename(  # use the combined survey
        {"iline": "new_iline", "xline": "new_xline"}
    )  # renaming inline and xline so they don't confuse xarray
    .stack(
        {"flat": ("new_iline", "new_xline")}
    )  # and flatten the iline and xline axes using stack
    .flat  # return just the flat coord object which is multi-index (so we can unstack later)
)
flat

flat = ( survey_comb.rename( # use the combined survey {"iline": "new_iline", "xline": "new_xline"} ) # renaming inline and xline so they don't confuse xarray .stack( {"flat": ("new_iline", "new_xline")} ) # and flatten the iline and xline axes using stack .flat # return just the flat coord object which is multi-index (so we can unstack later) ) flat

Out[14]:

<xarray.DataArray 'flat' (flat: 18054)> Size: 144kB
array([(640, 1099), (640, 1100), (640, 1101), ..., (757, 1249), (757, 1250),
       (757, 1251)], dtype=object)
Coordinates:
    cdp_x      (flat) float64 144kB 9.576e+03 9.575e+03 ... 9.634e+03 9.633e+03
    cdp_y      (flat) float64 144kB 2.932e+03 2.933e+03 ... 3.114e+03 3.114e+03
  * flat       (flat) object 144kB MultiIndex
  * new_iline  (flat) int32 72kB 640 640 640 640 640 640 ... 757 757 757 757 757
  * new_xline  (flat) int32 72kB 1099 1100 1101 1102 ... 1248 1249 1250 1251

xarray.DataArray

'flat'

• flat: 18054

• MultiIndex

  array([(640, 1099), (640, 1100), (640, 1101), ..., (757, 1249), (757, 1250),
         (757, 1251)], dtype=object)

• Coordinates: (5)
  • cdp_x
    (flat)
    float64
    9.576e+03 9.575e+03 ... 9.633e+03

    array([9575.52313979, 9575.18111964, 9574.8390995 , ..., 9634.16415492,
           9633.82213478, 9633.48011463])

  • cdp_y
    (flat)
    float64
    2.932e+03 2.933e+03 ... 3.114e+03

    array([2931.61508197, 2932.55477459, 2933.49446721, ..., 3112.58533186,
           3113.52502448, 3114.4647171 ])

  • flat
    (flat)
    object
    MultiIndex

    array([(640, 1099), (640, 1100), (640, 1101), ..., (757, 1249), (757, 1250),
           (757, 1251)], dtype=object)

  • new_iline
    (flat)
    int32
    640 640 640 640 ... 757 757 757 757

    array([640, 640, 640, ..., 757, 757, 757], dtype=int32)

  • new_xline
    (flat)
    int32
    1099 1100 1101 ... 1249 1250 1251

    array([1099, 1100, 1101, ..., 1249, 1250, 1251], dtype=int32)

• Indexes: (1)
  • flat
    new_iline
    new_xline
    PandasMultiIndex

    PandasIndex(MultiIndex([(640, 1099),
                (640, 1100),
                (640, 1101),
                (640, 1102),
                (640, 1103),
                (640, 1104),
                (640, 1105),
                (640, 1106),
                (640, 1107),
                (640, 1108),
                ...
                (757, 1242),
                (757, 1243),
                (757, 1244),
                (757, 1245),
                (757, 1246),
                (757, 1247),
                (757, 1248),
                (757, 1249),
                (757, 1250),
                (757, 1251)],
               name='flat', length=18054))

• Attributes: (0)

In [15]:

# create a dictionary of new coordinates to sample to. We use flat here because it will easily allow us to
# unstack the data after interpolation. It's necessary to use this flat dataset, otherwise xarray will try to
# interpolate on the cartesian product of iline and xline, which isn't really what we want.

resampling_data_arrays = dict(
    iline=xr.DataArray(survey_2_new_ilxl_loc[:, 0], dims="flat", coords={"flat": flat}),
    xline=xr.DataArray(survey_2_new_ilxl_loc[:, 1], dims="flat", coords={"flat": flat}),
)

resampling_data_arrays

# create a dictionary of new coordinates to sample to. We use flat here because it will easily allow us to # unstack the data after interpolation. It's necessary to use this flat dataset, otherwise xarray will try to # interpolate on the cartesian product of iline and xline, which isn't really what we want. resampling_data_arrays = dict( iline=xr.DataArray(survey_2_new_ilxl_loc[:, 0], dims="flat", coords={"flat": flat}), xline=xr.DataArray(survey_2_new_ilxl_loc[:, 1], dims="flat", coords={"flat": flat}), ) resampling_data_arrays

Out[15]:

{'iline': <xarray.DataArray (flat: 18054)> Size: 144kB
 array([653.01535386, 654.00016161, 654.98496936, ..., 780.41968002,
        781.40448778, 782.38929553])
 Coordinates:
   * flat       (flat) object 144kB MultiIndex
   * new_iline  (flat) int32 72kB 640 640 640 640 640 640 ... 757 757 757 757 757
   * new_xline  (flat) int32 72kB 1099 1100 1101 1102 ... 1248 1249 1250 1251,
 'xline': <xarray.DataArray (flat: 18054)> Size: 144kB
 array([1267.82560706, 1267.65195888, 1267.47831071, ..., 1126.55587331,
        1126.38222513, 1126.20857695])
 Coordinates:
   * flat       (flat) object 144kB MultiIndex
   * new_iline  (flat) int32 72kB 640 640 640 640 640 640 ... 757 757 757 757 757
   * new_xline  (flat) int32 72kB 1099 1100 1101 1102 ... 1248 1249 1250 1251}

The resampling can take a while with larger cubes, so it is good to use dask and to output the cube to disk at this step. When using dask for processing it is better to output to disk regularly, as this will improve how your code runs and the overall memory usage.

If you are executing this exmample locally, the task progress can be viewed by opening the dask client.

In [16]:

survey_2_resamp = survey_2.interp(**resampling_data_arrays)
survey_2_resamp_newgeom = (
    survey_2_resamp.drop_vars(("iline", "xline"))
    .unstack("flat")
    .rename({"new_iline": "iline", "new_xline": "xline"})
    .data
)
survey_2_resamp_newgeom.to_netcdf("data/survey_2_1.nc", compute=True, engine="h5netcdf")

survey_2_resamp = survey_2.interp(**resampling_data_arrays) survey_2_resamp_newgeom = ( survey_2_resamp.drop_vars(("iline", "xline")) .unstack("flat") .rename({"new_iline": "iline", "new_xline": "xline"}) .data ) survey_2_resamp_newgeom.to_netcdf("data/survey_2_1.nc", compute=True, engine="h5netcdf")

Combining Cubes¶

The last step is combining the cube is to load the two datasets to be combined and concatenate them together along a new axis. This will simplify reduction processes later.

In [17]:

survey_2_resamp_newgeom = xr.open_dataarray(
    "data/survey_2_1.nc", chunks=dict(iline=10, xline=10, twt=100), engine="h5netcdf"
)

survey_2_resamp_newgeom = xr.open_dataarray( "data/survey_2_1.nc", chunks=dict(iline=10, xline=10, twt=100), engine="h5netcdf" )

In [18]:

survey_2_resamp_newgeom.expand_dims({"survey": [2]})

survey_2_resamp_newgeom.expand_dims({"survey": [2]})

Out[18]:

<xarray.DataArray 'data' (survey: 1, twt: 201, iline: 118, xline: 153)> Size: 15MB
dask.array<broadcast_to, shape=(1, 201, 118, 153), dtype=float32, chunksize=(1, 100, 10, 10), chunktype=numpy.ndarray>
Coordinates:
  * survey   (survey) int64 8B 2
  * iline    (iline) int32 472B 640 641 642 643 644 645 ... 753 754 755 756 757
  * xline    (xline) int32 612B 1099 1100 1101 1102 1103 ... 1248 1249 1250 1251
  * twt      (twt) int64 2kB 0 1 2 3 4 5 6 7 ... 193 194 195 196 197 198 199 200
    cdp_x    (iline, xline) float64 144kB dask.array<chunksize=(10, 10), meta=np.ndarray>
    cdp_y    (iline, xline) float64 144kB dask.array<chunksize=(10, 10), meta=np.ndarray>

xarray.DataArray

'data'

• survey: 1
• twt: 201
• iline: 118
• xline: 153

• dask.array<chunksize=(1, 100, 10, 10), meta=np.ndarray>

  Array Chunk Bytes 13.84 MiB 39.06 kiB Shape (1, 201, 118, 153) (1, 100, 10, 10) Dask graph 576 chunks in 3 graph layers Data type float32 numpy.ndarray

• Coordinates: (6)
  • survey
    (survey)
    int64
    2

    array([2])

  • iline
    (iline)
    int32
    640 641 642 643 ... 754 755 756 757

    array([640, 641, 642, 643, 644, 645, 646, 647, 648, 649, 650, 651, 652, 653,
           654, 655, 656, 657, 658, 659, 660, 661, 662, 663, 664, 665, 666, 667,
           668, 669, 670, 671, 672, 673, 674, 675, 676, 677, 678, 679, 680, 681,
           682, 683, 684, 685, 686, 687, 688, 689, 690, 691, 692, 693, 694, 695,
           696, 697, 698, 699, 700, 701, 702, 703, 704, 705, 706, 707, 708, 709,
           710, 711, 712, 713, 714, 715, 716, 717, 718, 719, 720, 721, 722, 723,
           724, 725, 726, 727, 728, 729, 730, 731, 732, 733, 734, 735, 736, 737,
           738, 739, 740, 741, 742, 743, 744, 745, 746, 747, 748, 749, 750, 751,
           752, 753, 754, 755, 756, 757], dtype=int32)

  • xline
    (xline)
    int32
    1099 1100 1101 ... 1249 1250 1251

    array([1099, 1100, 1101, 1102, 1103, 1104, 1105, 1106, 1107, 1108, 1109, 1110,
           1111, 1112, 1113, 1114, 1115, 1116, 1117, 1118, 1119, 1120, 1121, 1122,
           1123, 1124, 1125, 1126, 1127, 1128, 1129, 1130, 1131, 1132, 1133, 1134,
           1135, 1136, 1137, 1138, 1139, 1140, 1141, 1142, 1143, 1144, 1145, 1146,
           1147, 1148, 1149, 1150, 1151, 1152, 1153, 1154, 1155, 1156, 1157, 1158,
           1159, 1160, 1161, 1162, 1163, 1164, 1165, 1166, 1167, 1168, 1169, 1170,
           1171, 1172, 1173, 1174, 1175, 1176, 1177, 1178, 1179, 1180, 1181, 1182,
           1183, 1184, 1185, 1186, 1187, 1188, 1189, 1190, 1191, 1192, 1193, 1194,
           1195, 1196, 1197, 1198, 1199, 1200, 1201, 1202, 1203, 1204, 1205, 1206,
           1207, 1208, 1209, 1210, 1211, 1212, 1213, 1214, 1215, 1216, 1217, 1218,
           1219, 1220, 1221, 1222, 1223, 1224, 1225, 1226, 1227, 1228, 1229, 1230,
           1231, 1232, 1233, 1234, 1235, 1236, 1237, 1238, 1239, 1240, 1241, 1242,
           1243, 1244, 1245, 1246, 1247, 1248, 1249, 1250, 1251], dtype=int32)

  • twt
    (twt)
    int64
    0 1 2 3 4 5 ... 196 197 198 199 200

    array([  0,   1,   2, ..., 198, 199, 200])

  • cdp_x
    (iline, xline)
    float64
    dask.array<chunksize=(10, 10), meta=np.ndarray>

    Array Chunk Bytes 141.05 kiB 800 B Shape (118, 153) (10, 10) Dask graph 192 chunks in 2 graph layers Data type float64 numpy.ndarray

  • cdp_y
    (iline, xline)
    float64
    dask.array<chunksize=(10, 10), meta=np.ndarray>

    Array Chunk Bytes 141.05 kiB 800 B Shape (118, 153) (10, 10) Dask graph 192 chunks in 2 graph layers Data type float64 numpy.ndarray

• Indexes: (4)
  • iline
    PandasIndex

    PandasIndex(Index([640, 641, 642, 643, 644, 645, 646, 647, 648, 649,
           ...
           748, 749, 750, 751, 752, 753, 754, 755, 756, 757],
          dtype='int32', name='iline', length=118))

  • xline
    PandasIndex

    PandasIndex(Index([1099, 1100, 1101, 1102, 1103, 1104, 1105, 1106, 1107, 1108,
           ...
           1242, 1243, 1244, 1245, 1246, 1247, 1248, 1249, 1250, 1251],
          dtype='int32', name='xline', length=153))

  • twt
    PandasIndex

    PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
           ...
           191, 192, 193, 194, 195, 196, 197, 198, 199, 200],
          dtype='int64', name='twt', length=201))

  • survey
    PandasIndex

    PandasIndex(Index([2], dtype='int64', name='survey'))

• Attributes: (0)

In [19]:

# concatenate survey with new dimension "survey".
survey_comb["data"] = xr.concat([survey_1.data, survey_2_resamp_newgeom], "survey")

# concatenate survey with new dimension "survey". survey_comb["data"] = xr.concat([survey_1.data, survey_2_resamp_newgeom], "survey")

/home/runner/.local/share/hatch/env/virtual/segysak/7aGXAoce/docs/lib/python3.10/site-packages/dask/array/core.py:4888: PerformanceWarning: Increasing number of chunks by factor of 14
  result = blockwise(

In [20]:

survey_comb

survey_comb

Out[20]:

<xarray.Dataset> Size: 29MB
Dimensions:  (iline: 118, xline: 153, twt: 201, survey: 2)
Coordinates:
  * iline    (iline) int64 944B 640 641 642 643 644 645 ... 753 754 755 756 757
  * xline    (xline) int64 1kB 1099 1100 1101 1102 1103 ... 1248 1249 1250 1251
  * twt      (twt) int64 2kB 0 1 2 3 4 5 6 7 ... 193 194 195 196 197 198 199 200
    cdp_x    (iline, xline) float64 144kB 9.576e+03 9.575e+03 ... 9.633e+03
    cdp_y    (iline, xline) float64 144kB 2.932e+03 2.933e+03 ... 3.114e+03
Dimensions without coordinates: survey
Data variables:
    data     (survey, iline, xline, twt) float32 29MB dask.array<chunksize=(1, 9, 9, 100), meta=np.ndarray>
Attributes: (12/17)
    ns:                  None
    sample_rate:         None
    text:                None
    measurement_system:  None
    d3_domain:           None
    epsg:                None
    ...                  ...
    percentiles:         None
    coord_scalar:        None
    coord_scaled:        None
    dimensions:          None
    vert_dimension:      None
    vert_domain:         None

xarray.Dataset

• Dimensions:
  • iline: 118
  • xline: 153
  • twt: 201
  • survey: 2
• Coordinates: (5)
  • iline
    (iline)
    int64
    640 641 642 643 ... 754 755 756 757

    array([640, 641, 642, 643, 644, 645, 646, 647, 648, 649, 650, 651, 652, 653,
           654, 655, 656, 657, 658, 659, 660, 661, 662, 663, 664, 665, 666, 667,
           668, 669, 670, 671, 672, 673, 674, 675, 676, 677, 678, 679, 680, 681,
           682, 683, 684, 685, 686, 687, 688, 689, 690, 691, 692, 693, 694, 695,
           696, 697, 698, 699, 700, 701, 702, 703, 704, 705, 706, 707, 708, 709,
           710, 711, 712, 713, 714, 715, 716, 717, 718, 719, 720, 721, 722, 723,
           724, 725, 726, 727, 728, 729, 730, 731, 732, 733, 734, 735, 736, 737,
           738, 739, 740, 741, 742, 743, 744, 745, 746, 747, 748, 749, 750, 751,
           752, 753, 754, 755, 756, 757])

  • xline
    (xline)
    int64
    1099 1100 1101 ... 1249 1250 1251

    array([1099, 1100, 1101, 1102, 1103, 1104, 1105, 1106, 1107, 1108, 1109, 1110,
           1111, 1112, 1113, 1114, 1115, 1116, 1117, 1118, 1119, 1120, 1121, 1122,
           1123, 1124, 1125, 1126, 1127, 1128, 1129, 1130, 1131, 1132, 1133, 1134,
           1135, 1136, 1137, 1138, 1139, 1140, 1141, 1142, 1143, 1144, 1145, 1146,
           1147, 1148, 1149, 1150, 1151, 1152, 1153, 1154, 1155, 1156, 1157, 1158,
           1159, 1160, 1161, 1162, 1163, 1164, 1165, 1166, 1167, 1168, 1169, 1170,
           1171, 1172, 1173, 1174, 1175, 1176, 1177, 1178, 1179, 1180, 1181, 1182,
           1183, 1184, 1185, 1186, 1187, 1188, 1189, 1190, 1191, 1192, 1193, 1194,
           1195, 1196, 1197, 1198, 1199, 1200, 1201, 1202, 1203, 1204, 1205, 1206,
           1207, 1208, 1209, 1210, 1211, 1212, 1213, 1214, 1215, 1216, 1217, 1218,
           1219, 1220, 1221, 1222, 1223, 1224, 1225, 1226, 1227, 1228, 1229, 1230,
           1231, 1232, 1233, 1234, 1235, 1236, 1237, 1238, 1239, 1240, 1241, 1242,
           1243, 1244, 1245, 1246, 1247, 1248, 1249, 1250, 1251])

  • twt
    (twt)
    int64
    0 1 2 3 4 5 ... 196 197 198 199 200

    array([  0,   1,   2, ..., 198, 199, 200])

  • cdp_x
    (iline, xline)
    float64
    9.576e+03 9.575e+03 ... 9.633e+03

    array([[9575.52313979, 9575.18111964, 9574.8390995 , ..., 9524.22011829,
            9523.87809815, 9523.536078  ],
           [9576.46283241, 9576.12081227, 9575.77879212, ..., 9525.15981091,
            9524.81779077, 9524.47577062],
           [9577.40252503, 9577.06050489, 9576.71848474, ..., 9526.09950353,
            9525.75748339, 9525.41546324],
           ...,
           [9683.58779118, 9683.24577104, 9682.90375089, ..., 9632.28476968,
            9631.94274954, 9631.60072939],
           [9684.5274838 , 9684.18546366, 9683.84344351, ..., 9633.2244623 ,
            9632.88244216, 9632.54042201],
           [9685.46717642, 9685.12515628, 9684.78313613, ..., 9634.16415492,
            9633.82213478, 9633.48011463]])

  • cdp_y
    (iline, xline)
    float64
    2.932e+03 2.933e+03 ... 3.114e+03

    array([[2931.61508197, 2932.55477459, 2933.49446721, ..., 3072.56897509,
            3073.50866771, 3074.44836033],
           [2931.95710212, 2932.89679474, 2933.83648736, ..., 3072.91099523,
            3073.85068785, 3074.79038047],
           [2932.29912226, 2933.23881488, 2934.1785075 , ..., 3073.25301538,
            3074.192708  , 3075.13240062],
           ...,
           [2970.94739845, 2971.88709108, 2972.8267837 , ..., 3111.90129157,
            3112.84098419, 3113.78067681],
           [2971.2894186 , 2972.22911122, 2973.16880384, ..., 3112.24331172,
            3113.18300434, 3114.12269696],
           [2971.63143874, 2972.57113136, 2973.51082398, ..., 3112.58533186,
            3113.52502448, 3114.4647171 ]])

• Data variables: (1)
  • data
    (survey, iline, xline, twt)
    float32
    dask.array<chunksize=(1, 9, 9, 100), meta=np.ndarray>

    Array Chunk Bytes 27.69 MiB 31.64 kiB Shape (2, 118, 153, 201) (1, 9, 9, 100) Dask graph 4464 chunks in 22 graph layers Data type float32 numpy.ndarray

• Indexes: (3)
  • iline
    PandasIndex

    PandasIndex(Index([640, 641, 642, 643, 644, 645, 646, 647, 648, 649,
           ...
           748, 749, 750, 751, 752, 753, 754, 755, 756, 757],
          dtype='int64', name='iline', length=118))

  • xline
    PandasIndex

    PandasIndex(Index([1099, 1100, 1101, 1102, 1103, 1104, 1105, 1106, 1107, 1108,
           ...
           1242, 1243, 1244, 1245, 1246, 1247, 1248, 1249, 1250, 1251],
          dtype='int64', name='xline', length=153))

  • twt
    PandasIndex

    PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
           ...
           191, 192, 193, 194, 195, 196, 197, 198, 199, 200],
          dtype='int64', name='twt', length=201))

• Attributes: (17)

  ns :

  None

  sample_rate :

  None

  text :

  None

  measurement_system :

  None

  d3_domain :

  None

  epsg :

  None

  corner_points :

  None

  corner_points_xy :

  None

  source_file :

  None

  srd :

  None

  datatype :

  None

  percentiles :

  None

  coord_scalar :

  None

  coord_scaled :

  None

  dimensions :

  None

  vert_dimension :

  None

  vert_domain :

  None

We can check the merged surveys by looking at some plots. If we select just the first survey the values will be 1. If we select just the second survey the values will be 2. And if we take the mean along the survey dimension, then where the surveys overlap the values will be 1.5.

For seismic data, a better form of conditioning and reduction might be required for merging traces together to ensure a smoother seam.

In [21]:

sel = survey_comb.sel(xline=1190)

fig, axs = plt.subplots(ncols=3, figsize=(30, 10))

sel.isel(survey=0).data.T.plot(ax=axs[0], yincrease=False, vmin=1, vmax=2)
sel.isel(survey=1).data.T.plot(ax=axs[1], yincrease=False, vmin=1, vmax=2)
sel.data.mean("survey").T.plot(ax=axs[2], yincrease=False, vmin=1, vmax=2)

sel = survey_comb.sel(xline=1190) fig, axs = plt.subplots(ncols=3, figsize=(30, 10)) sel.isel(survey=0).data.T.plot(ax=axs[0], yincrease=False, vmin=1, vmax=2) sel.isel(survey=1).data.T.plot(ax=axs[1], yincrease=False, vmin=1, vmax=2) sel.data.mean("survey").T.plot(ax=axs[2], yincrease=False, vmin=1, vmax=2)

Out[21]:

<matplotlib.collections.QuadMesh at 0x7f34f437f520>

In [22]:

survey_comb.twt

survey_comb.twt

Out[22]:

<xarray.DataArray 'twt' (twt: 201)> Size: 2kB
array([  0,   1,   2, ..., 198, 199, 200])
Coordinates:
  * twt      (twt) int64 2kB 0 1 2 3 4 5 6 7 ... 193 194 195 196 197 198 199 200

xarray.DataArray

'twt'

• twt: 201

• 0 1 2 3 4 5 6 7 8 9 10 ... 190 191 192 193 194 195 196 197 198 199 200

  array([  0,   1,   2, ..., 198, 199, 200])

• Coordinates: (1)
  • twt
    (twt)
    int64
    0 1 2 3 4 5 ... 196 197 198 199 200

    array([  0,   1,   2, ..., 198, 199, 200])

• Indexes: (1)
  • twt
    PandasIndex

    PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
           ...
           191, 192, 193, 194, 195, 196, 197, 198, 199, 200],
          dtype='int64', name='twt', length=201))

• Attributes: (0)

And if we inspect a time slice.

In [23]:

sel = survey_comb.sel(twt=100)

fig, axs = plt.subplots(ncols=3, figsize=(30, 10))

sel.isel(survey=0).data.T.plot(ax=axs[0], yincrease=False, vmin=1, vmax=2)
sel.isel(survey=1).data.T.plot(ax=axs[1], yincrease=False, vmin=1, vmax=2)
sel.data.mean("survey").T.plot(ax=axs[2], yincrease=False, vmin=1, vmax=2)

sel = survey_comb.sel(twt=100) fig, axs = plt.subplots(ncols=3, figsize=(30, 10)) sel.isel(survey=0).data.T.plot(ax=axs[0], yincrease=False, vmin=1, vmax=2) sel.isel(survey=1).data.T.plot(ax=axs[1], yincrease=False, vmin=1, vmax=2) sel.data.mean("survey").T.plot(ax=axs[2], yincrease=False, vmin=1, vmax=2)

Out[23]:

<matplotlib.collections.QuadMesh at 0x7f34dc1cf9a0>

In [ ]: