.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "tutorials/nc_CCPStack_SingleStation.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_tutorials_nc_CCPStack_SingleStation.py>`
        to download the full example code

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_tutorials_nc_CCPStack_SingleStation.py:


Common Conversion Point Stack of a Single Station
=================================================

The following notebook carries you through how to get from 
a set of receiver functions of a single station to a 3D 
common conversion point (CCP) stack. 
As an example as in the previous notebooks,
we again use IU-HRV as example station. 

.. note::

    It is assumed here that you have successfully computed the
    receiver functions from the `DataCollection{SAC | HDF5}.py`.


Load the Receiver functions
+++++++++++++++++++++++++++

.. GENERATED FROM PYTHON SOURCE LINES 22-25

.. code-block:: default

    # sphinx_gallery_thumbnail_number = 1
    # sphinx_gallery_dummy_images = 1








.. GENERATED FROM PYTHON SOURCE LINES 26-27

So, first load the receiver functions into a `RFStream`.

.. GENERATED FROM PYTHON SOURCE LINES 27-50

.. code-block:: default

    import os
    import matplotlib.pyplot as plt
    from pyglimer.rf.create import read_rf
    from pyglimer.plot.plot_utils import set_mpl_params
    set_mpl_params()

    # Define the location of the database
    databaselocation = os.path.join('static_data', 'database_sac')  # or "database_hdf5"

    try: db_base_path = ipynb_path
    except NameError:
        try: db_base_path = os.path.dirname(os.path.realpath(__file__))
        except NameError: db_base_path = os.getcwd()
    databaselocation = os.path.join(db_base_path, databaselocation)

    rfst = read_rf(os.path.join(
        databaselocation, 'waveforms', 'RF', 'P', 'IU', 'HRV', '*.sac'))

    # Check traces
    print("Number of loaded RFs: ", len(rfst))
    rfst[0].plot()
    plt.show(block=False)




.. image-sg:: /tutorials/images/sphx_glr_nc_CCPStack_SingleStation_001.png
   :alt: nc CCPStack SingleStation
   :srcset: /tutorials/images/sphx_glr_nc_CCPStack_SingleStation_001.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Number of loaded RFs:  7




.. GENERATED FROM PYTHON SOURCE LINES 51-56

Compute Common Conversion Point Stack
+++++++++++++++++++++++++++++++++++++

This is similar to the single station stacks.


.. GENERATED FROM PYTHON SOURCE LINES 56-82

.. code-block:: default


    from pyglimer.ccp import init_ccp
    import os
    import numpy as np

    inter_bin_distance = 0.1
    velocity_model = 'iasp91.dat'

    ccp_init_dict = {
        'spacing': inter_bin_distance,
        'vel_model': velocity_model,
        "binrad": np.cos(np.radians(30)),
        "phase": 'P',
        "preproloc": os.path.join(databaselocation, "waveforms", "preprocessed"),
        "rfloc": os.path.join(databaselocation, "waveforms", "RF"),
        "network": "IU",
        "station": "HRV",
        "compute_stack": True,
        "save": 'ccp_IU_HRV.pkl',
        'format': 'sac',
        'mc_backend': 'joblib',
    }

    # Initialize bins
    ccpstack = init_ccp(**ccp_init_dict)





.. rst-class:: sphx-glr-script-out

 .. code-block:: none


      0%|          | 0/2 [00:00<?, ?it/s]
     50%|#####     | 1/2 [00:00<00:00,  6.41it/s]
    100%|##########| 2/2 [00:00<00:00, 12.69it/s]




.. GENERATED FROM PYTHON SOURCE LINES 83-84

Finalizing the CCP Stack checking whether locations are on-land or not.

.. GENERATED FROM PYTHON SOURCE LINES 84-87

.. code-block:: default


    ccpstack.conclude_ccp(keep_water=True)








.. GENERATED FROM PYTHON SOURCE LINES 88-90

Plot Bins
+++++++++

.. GENERATED FROM PYTHON SOURCE LINES 90-94

.. code-block:: default


    import matplotlib.pyplot as plt
    ccpstack.plot_bins()




.. image-sg:: /tutorials/images/sphx_glr_nc_CCPStack_SingleStation_002.png
   :alt: nc CCPStack SingleStation
   :srcset: /tutorials/images/sphx_glr_nc_CCPStack_SingleStation_002.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 95-108

Use the ``CCPStack`` object to image the subsurface
+++++++++++++++++++++++++++++++++++++++++++++++++++

Given a ``CCPStack`` object there are multiple ways of getting an image of 
the subsurface. You could either compute a three-dimensional volume and plot 
slices with respect to the volume or you can create cross sections slices.

Compute a volume and slice
--------------------------

To do that we first need to read the created ``CCPStack``. Second, we create
a discretization for the are of interest, and third, we compute the volume and
some slices.

.. GENERATED FROM PYTHON SOURCE LINES 108-124

.. code-block:: default


    import numpy as np
    from pyglimer.ccp.ccp import read_ccp

    # Reading the stack
    ccpstack = read_ccp(filename='ccp_IU_HRV.pkl', fmt=None)

    # Create discretization
    lats = np.arange(41, 43.5, 0.05)
    lons = np.arange(-72.7, -69.5, 0.05)
    z = np.linspace(-10, 200, 211)

    # Plotting the volume and slices
    vplot = ccpstack.plot_volume_sections(
        lons, lats, zmax=211, lonsl=-71.45, latsl=42.5, zsl=23)




.. image-sg:: /tutorials/images/sphx_glr_nc_CCPStack_SingleStation_003.png
   :alt: nc CCPStack SingleStation
   :srcset: /tutorials/images/sphx_glr_nc_CCPStack_SingleStation_003.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Creating Spherical Nearest Neighbour class ... done.
    Creating interpolators ... done.
    KDTree interpolation ... (intensive part)
    KDTree interpolation ... (intensive part)--->  001/311
    KDTree interpolation ... (intensive part)--->  032/311
    KDTree interpolation ... (intensive part)--->  063/311
    KDTree interpolation ... (intensive part)--->  094/311
    KDTree interpolation ... (intensive part)--->  125/311
    KDTree interpolation ... (intensive part)--->  156/311
    KDTree interpolation ... (intensive part)--->  187/311
    KDTree interpolation ... (intensive part)--->  218/311
    KDTree interpolation ... (intensive part)--->  249/311
    KDTree interpolation ... (intensive part)--->  280/311
    KDTree interpolation ... (intensive part)--->  311/311
    ... done.
    Regular Grid interpolation ... (less intensive part) done.




.. GENERATED FROM PYTHON SOURCE LINES 125-133

Slice ``CCPStack`` directly using
---------------------------------

Here we first define two slices with their waypoints, and then plot a cross 
section in the back end, both RF sections and Illumination section are 
computed using a KDTree algorithm taking as input the ``CCPStack`` locations.
For the illumination, this sometimes leads to weird results; see profile A


.. GENERATED FROM PYTHON SOURCE LINES 133-159

.. code-block:: default


    # Create points waypoints for the cross section
    lat0 = np.array([42.5, 42.5])
    lon0 = np.array([-72.7, -69.5])
    lat1 = np.array([41.5, 43.5])
    lon1 = np.array([-71.45, -71.45])

    # Set RF boundaries
    mapextent=[-73.5, -69, 41, 44]
    depthextent=[0, 200]
    vmin = -0.1
    vmax =  0.1
 
    # Plot cross sections
    ax1, geoax = ccpstack.plot_cross_section(
        lat0, lon0, ddeg=0.01, z0=23, vmin=vmin, vmax=vmax,
        mapplot=True, minillum=1, label="A",
        depthextent=depthextent,
        mapextent=mapextent)
    ax2, _ = ccpstack.plot_cross_section(
        lat1, lon1, ddeg=0.01, vmin=vmin, vmax=vmax,
        geoax=geoax, mapplot=True,  
        minillum=1, label="B",
        depthextent=depthextent)
    plt.show(block=False)




.. rst-class:: sphx-glr-horizontal


    *

      .. image-sg:: /tutorials/images/sphx_glr_nc_CCPStack_SingleStation_004.png
         :alt: nc CCPStack SingleStation
         :srcset: /tutorials/images/sphx_glr_nc_CCPStack_SingleStation_004.png
         :class: sphx-glr-multi-img

    *

      .. image-sg:: /tutorials/images/sphx_glr_nc_CCPStack_SingleStation_005.png
         :alt: nc CCPStack SingleStation
         :srcset: /tutorials/images/sphx_glr_nc_CCPStack_SingleStation_005.png
         :class: sphx-glr-multi-img

    *

      .. image-sg:: /tutorials/images/sphx_glr_nc_CCPStack_SingleStation_006.png
         :alt: nc CCPStack SingleStation
         :srcset: /tutorials/images/sphx_glr_nc_CCPStack_SingleStation_006.png
         :class: sphx-glr-multi-img





.. GENERATED FROM PYTHON SOURCE LINES 160-181

Note that here we are only running the example only with a very limited
dataset, for which CCP Stacking makes limited sense. It still gives you
an idea on how you could do this for a larger area.
``z0`` here is the depth of the illumination map and ``minillum`` the minimum 
number of hitpoints per bin to be counted.

.. note::

    The waypoints do not have to be just start and endpoint, they can be a 
    series of points.

Compute a "dirty" global stack
------------------------------

In short, we are assuming that latitudes and longitudes are cartesian
entities, with a small correction for the area of each bin that depends on the
change in metric width of a degree of longitude as a function of latitude.
This really is a dirty way of interpolation because the spherical nature of
the Earth is disregarded.

We first read the RFs needed for the computation

.. GENERATED FROM PYTHON SOURCE LINES 181-188

.. code-block:: default


    from pyglimer.rf.create import read_rf
    from pyglimer.plot.plot_utils import set_mpl_params
    set_mpl_params()

    rfst = read_rf(f"{databaselocation}/waveforms/RF/P/IU/HRV/*.sac")








.. GENERATED FROM PYTHON SOURCE LINES 189-191

Then, we migrate the receiver functions, i.e., move-out correct the RFs wrt.
to a given model

.. GENERATED FROM PYTHON SOURCE LINES 191-194

.. code-block:: default


    z, rfz = rfst.moveout(vmodel='iasp91.dat')








.. GENERATED FROM PYTHON SOURCE LINES 195-202

and, compute the CCP volume like a 3D histogram by providing the extent 
of the region of interest, and ``dlon``, which describes the longitudinal
spacing.

.. note::

    The latitudinal spacing is automatically computed.

.. GENERATED FROM PYTHON SOURCE LINES 202-209

.. code-block:: default


    # Set extent=[minlon, maxlon, minlat, maxlat]
    extent = [-72.7, -69.5, 41, 43.5]

    # Compute Dirty Stack
    latc, lonc, zc, ccp, illum = rfz.dirty_ccp_stack(dlon=0.075, extent=extent)








.. GENERATED FROM PYTHON SOURCE LINES 210-212

Then, we can use the 3D histogram to plot the stack using the same tool 
as the 'clean' stacking.

.. GENERATED FROM PYTHON SOURCE LINES 212-216

.. code-block:: default


    from pyglimer.plot.plot_volume import VolumePlot

    vp = VolumePlot(lonc, latc, zc[:211], ccp[:,:,:211], xl=-71.45, yl=42.5, zl=23)
    plt.show(block=False)


.. image-sg:: /tutorials/images/sphx_glr_nc_CCPStack_SingleStation_007.png
   :alt: nc CCPStack SingleStation
   :srcset: /tutorials/images/sphx_glr_nc_CCPStack_SingleStation_007.png
   :class: sphx-glr-single-img






.. rst-class:: sphx-glr-timing

   **Total running time of the script:** ( 0 minutes  35.851 seconds)


.. _sphx_glr_download_tutorials_nc_CCPStack_SingleStation.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example




    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: nc_CCPStack_SingleStation.py <nc_CCPStack_SingleStation.py>`

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: nc_CCPStack_SingleStation.ipynb <nc_CCPStack_SingleStation.ipynb>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_