{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Blending\nThis example shows how to use the :py:class:`pylops.waveeqprocessing.blending.Blending`\noperator to blend seismic data to mimic state-of-the-art simultaneous shooting\nacquisition systems.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\nimport numpy as np\nimport scipy as sp\n\nimport pylops\n\nplt.close(\"all\")\nnp.random.seed(0)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's start by considering a streamer seismic dataset and apply blending in\nso-called continuous blending mode\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "inputdata = np.load(\"../testdata/marchenko/input.npz\")\ndata = inputdata[\"R\"]\ndata = np.pad(data, ((0, 0), (0, 0), (0, 50)))\nwav = inputdata[\"wav\"]\nwav_c = np.argmax(wav)\nns, nr, nt = data.shape\n\n# time axis\ndt = 0.004\nt = np.arange(nt) * dt\n\n# convolve with wavelet\ndata = np.apply_along_axis(sp.signal.convolve, -1, data, wav, mode=\"full\")\ndata = data[..., wav_c:][..., :nt]\n\n# obc data\ndata_obc = data[:-1, :-1]\nns_obc, nr_obc, _ = data_obc.shape\n\n# streamer data\nnr_streamer = 21\nns_streamer = ns - nr_streamer\n\ndata_streamer = np.zeros((ns_streamer, nr_streamer, nt))\nfor isrc in range(ns_streamer):\n    data_streamer[isrc] = data[isrc, isrc : isrc + nr_streamer]\n\n# visualize\nisrcplot = [0, ns_obc // 2, ns_obc - 1]\nfig, axs = plt.subplots(1, 3, sharey=True, figsize=(12, 8))\nfig.suptitle(\"OBC data\")\nfor i, ax in enumerate(axs):\n    ax.imshow(\n        data_obc[isrcplot[i]].T,\n        cmap=\"gray\",\n        vmin=-0.1,\n        vmax=0.1,\n        extent=(0, nr, t[-1], 0),\n        interpolation=\"none\",\n    )\n    ax.set_title(f\"CSG {isrcplot[i]}\")\n    ax.set_xlabel(\"#Rec\")\n    ax.axis(\"tight\")\naxs[0].set_ylabel(\"t [s]\")\nplt.tight_layout()\n\nisrcplot = [0, ns_streamer // 2, ns_streamer - 1]\nfig, axs = plt.subplots(1, 3, sharey=True, figsize=(12, 8))\nfig.suptitle(\"Streamer data\")\nfor i, ax in enumerate(axs):\n    ax.imshow(\n        data_streamer[isrcplot[i]].T,\n        cmap=\"gray\",\n        vmin=-0.1,\n        vmax=0.1,\n        extent=(0, nr_streamer, t[-1], 0),\n        interpolation=\"none\",\n    )\n    ax.set_title(f\"CSG {isrcplot[i]}\")\n    ax.set_xlabel(\"#Rec\")\n    ax.axis(\"tight\")\naxs[0].set_ylabel(\"t [s]\")\nplt.tight_layout()\n\nirecplot = [0, nr_streamer // 2, nr_streamer - 1]\nfig, axs = plt.subplots(1, 3, sharey=True, figsize=(12, 8))\nfig.suptitle(\"Streamer data\")\nfor i, ax in enumerate(axs):\n    ax.imshow(\n        data_streamer[:, irecplot[i]].T,\n        cmap=\"gray\",\n        vmin=-0.1,\n        vmax=0.1,\n        extent=(0, ns_streamer, t[-1], 0),\n        interpolation=\"none\",\n    )\n    ax.set_title(f\"CRG {irecplot[i]}\")\n    ax.set_xlabel(\"#Src\")\n    ax.axis(\"tight\")\naxs[0].set_ylabel(\"t [s]\")\nplt.tight_layout()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can now consider the streamer seismic dataset and apply blending in\nso-called continuous blending mode\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "overlap = 0.5\nignition_times = np.random.normal(0, 0.6, ns_streamer)\nignition_times += (1 - overlap) * nt * dt\nignition_times[0] = 0.0\nignition_times = np.cumsum(ignition_times)\n\nplt.figure(figsize=(12, 4))\nplt.plot(ignition_times, \"k\")\nplt.title(\"Continuous blending times\")\n\nBop = pylops.waveeqprocessing.BlendingContinuous(\n    nt,\n    nr_streamer,\n    ns_streamer,\n    dt,\n    ignition_times,\n    dtype=\"complex128\",\n)\ndata_blended = Bop * data_streamer\ndata_pseudo = Bop.H * data_blended\n\nfig, ax = plt.subplots(1, 1, figsize=(4, 19))\nax.imshow(\n    data_blended.real.T,\n    cmap=\"gray\",\n    vmin=-0.1,\n    vmax=0.1,\n    extent=(0, ns_streamer, Bop.nttot * dt, 0),\n    interpolation=\"none\",\n)\nax.set_title(\"Blended CSG\")\nax.set_xlabel(\"#Rec\")\nax.set_ylabel(\"t [s]\")\nax.axis(\"tight\")\nax.set_ylim(10, 0)\nplt.tight_layout()\n\nfig, axs = plt.subplots(1, 2, sharey=True, figsize=(12, 8))\naxs[0].imshow(\n    data_streamer[:, 0].real.T,\n    cmap=\"gray\",\n    vmin=-0.01,\n    vmax=0.01,\n    extent=(0, ns_streamer, t[-1], 0),\n    interpolation=\"none\",\n)\naxs[0].set_title(\"Unblended CRG\")\naxs[0].set_xlabel(\"#Src\")\naxs[0].set_ylabel(\"t [s]\")\naxs[0].axis(\"tight\")\naxs[1].imshow(\n    data_pseudo[:, 0].real.T,\n    cmap=\"gray\",\n    vmin=-0.01,\n    vmax=0.01,\n    extent=(0, ns_streamer, t[-1], 0),\n    interpolation=\"none\",\n)\naxs[1].set_title(\"Pseudo-deblended CRG\")\naxs[1].set_xlabel(\"#Src\")\naxs[1].axis(\"tight\")\nplt.tight_layout()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Similarly we can consider the OBC data and apply both group and half blending\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Group\ngroup_size = 2\nn_groups = ns_obc // 2\nignition_times = np.abs(np.random.normal(0.2, 0.5, ns_obc))  # only positive shifts\nignition_times[0] = 0.0\n\nplt.figure(figsize=(12, 4))\nplt.plot(ignition_times.reshape(group_size, n_groups).T, \"k\")\nplt.title(\"Group blending times\")\n\nBop = pylops.waveeqprocessing.BlendingGroup(\n    nt,\n    nr_obc,\n    ns_obc,\n    dt,\n    ignition_times.reshape(group_size, n_groups),\n    group_size=group_size,\n    n_groups=n_groups,\n    dtype=\"complex128\",\n)\ndata_blended = Bop * data_obc\ndata_pseudo = Bop.H * data_blended\n\nfig, ax = plt.subplots(1, 1, figsize=(12, 8))\nax.imshow(\n    data_blended[n_groups // 2].real.T,\n    cmap=\"gray\",\n    vmin=-0.1,\n    vmax=0.1,\n    extent=(0, ns_streamer, t[-1], 0),\n    interpolation=\"none\",\n)\nax.set_title(\"Blended CSG\")\nax.set_xlabel(\"#Rec\")\nax.set_ylabel(\"t [s]\")\nax.axis(\"tight\")\nplt.tight_layout()\n\nfig, axs = plt.subplots(1, 2, sharey=True, figsize=(12, 8))\naxs[0].imshow(\n    data_obc[:, 10].real.T,\n    cmap=\"gray\",\n    vmin=-0.01,\n    vmax=0.01,\n    extent=(0, ns_streamer, t[-1], 0),\n    interpolation=\"none\",\n)\naxs[0].set_title(\"Unblended CRG\")\naxs[0].set_xlabel(\"#Src\")\naxs[0].set_ylabel(\"t [s]\")\naxs[0].axis(\"tight\")\naxs[1].imshow(\n    data_pseudo[:, 10].real.T,\n    cmap=\"gray\",\n    vmin=-0.01,\n    vmax=0.01,\n    extent=(0, ns_streamer, t[-1], 0),\n    interpolation=\"none\",\n)\naxs[1].set_title(\"Pseudo-deblended CRG\")\naxs[1].set_xlabel(\"#Src\")\naxs[1].axis(\"tight\")\nplt.tight_layout()\n\n# Half\ngroup_size = 2\nn_groups = ns_obc // 2\nignition_times = np.abs(np.random.normal(0.1, 0.5, ns_obc))  # only positive shifts\nignition_times[0] = 0.0\n\nplt.figure(figsize=(12, 4))\nplt.plot(ignition_times.reshape(group_size, n_groups).T, \"k\")\nplt.title(\"Half blending times\")\n\nBop = pylops.waveeqprocessing.BlendingHalf(\n    nt,\n    nr_obc,\n    ns_obc,\n    dt,\n    ignition_times.reshape(group_size, n_groups),\n    group_size=group_size,\n    n_groups=n_groups,\n    dtype=\"complex128\",\n    name=None,\n)\ndata_blended = Bop * data_obc\ndata_pseudo = Bop.H * data_blended\n\nfig, ax = plt.subplots(1, 1, figsize=(12, 8))\nax.imshow(\n    data_blended[n_groups // 2].real.T,\n    cmap=\"gray\",\n    vmin=-0.1,\n    vmax=0.1,\n    extent=(0, ns_streamer, t[-1], 0),\n    interpolation=\"none\",\n)\nax.set_title(\"Blended CSG\")\nax.set_xlabel(\"#Rec\")\nax.set_ylabel(\"t [s]\")\nax.axis(\"tight\")\nplt.tight_layout()\n\nfig, axs = plt.subplots(1, 2, sharey=True, figsize=(12, 8))\naxs[0].imshow(\n    data_obc[:, 10].real.T,\n    cmap=\"gray\",\n    vmin=-0.01,\n    vmax=0.01,\n    extent=(0, ns_streamer, t[-1], 0),\n    interpolation=\"none\",\n)\naxs[0].set_title(\"Unblended CRG\")\naxs[0].set_xlabel(\"#Src\")\naxs[0].set_ylabel(\"t [s]\")\naxs[0].axis(\"tight\")\naxs[1].imshow(\n    data_pseudo[:, 10].real.T,\n    cmap=\"gray\",\n    vmin=-0.01,\n    vmax=0.01,\n    extent=(0, ns_streamer, t[-1], 0),\n    interpolation=\"none\",\n)\naxs[1].set_title(\"Pseudo-deblended CRG\")\naxs[1].set_xlabel(\"#Src\")\naxs[1].axis(\"tight\")\nplt.tight_layout()"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.9.15"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}