<jupyter_start><jupyter_text>Query Data and Plot Together***This is a duplicate of the `local_01_getAllData` notebook. Just added a cell at the beginning and ran on a SageMaker instance instead!******Things run a little faster here since we are closer to the S3 buckets***In this notebook we will explain how to query well, horizon, and seismic data.We will explain how to extract seismic data around the wellbore using `real-simple-seismic`. Same is doable with`OpenVDS` files as well using the `get_minicube` convenience function. Since `OpenVDS` has bricked indexing as wellit will bring chunks of data much faster than `real-simple-seismic`. If you are using a deep neural networkwe recommend you stick to `OpenVDS`, if not `real-simple-seismic` easier to use.Finally, we will show nice plots that show all data together!Let's start by importing some of the tools we need.<jupyter_code>%%capture
%pip install s3fs==0.5.0
%pip install real-simple-seismic
from math import floor, ceil
from os.path import join

from rss.client import rssFromS3
from pandas import concat, IndexSlice, read_csv, read_json<jupyter_output><empty_output><jupyter_text>Loading AWS Credentials**THIS ONLY APPLIES TO ON-PREM (LOCAL) ACCESS TO S3 DATA. IF YOU ARE ON SAGEMAKER DON'T DO THIS.**AWS credentials must be stored in environment variables.To access competition data on-premise, you will need these environment variables set. They expire occasionally, so please update your variables when you get permission errors. See instructions PDF on how to obtain them.Variables: `AWS_ACCESS_KEY_ID`, `AWS_SECRET_ACCESS_KEY`, and `AWS_SESSION_TOKEN`If these are not set, you will get a "Forbidden" error. Configuring this depends on your operating system. It needs to be set in your global environment or within your virtual environment **before** launching Jupyter server or scripts. It is usually set using terminal. AWS SSO logins provides a convenient way to copy these for your account into your terminal or command-prompt. It will look something like this:| **Linux / UNIX / MacOS**                                | **Windows**                                        ||:--------------------------------------------------------|---------------------------------------------------:|| `export AWS_ACCESS_KEY_ID="your_access_key"`            | `SET AWS_ACCESS_KEY_ID=your_access_key`            || `export AWS_SECRET_ACCESS_KEY="your_secret_access_key"` | `SET AWS_SECRET_ACCESS_KEY=your_secret_access_key` || `export AWS_SESSION_TOKEN="your_long_session_token"`    | `SET AWS_SESSION_TOKEN=your_long_session_token`    |After credentials are set, we can read data from the S3 buckets just like we are on SageMaker environment. Reading the wells from S3.<jupyter_code>%%time
well_bucket = 's3://sagemaker-gitc2021/poseidon/wells/'
well_file = 'poseidon_geoml_training_wells.json.gz'

well_df = read_json(
    path_or_buf=join(well_bucket, well_file),
    compression='gzip',
)

well_df.set_index(['well_id', 'twt'], inplace=True)

well_df<jupyter_output>CPU times: user 170 ms, sys: 18.1 ms, total: 188 ms
Wall time: 256 ms<jupyter_text>Get well names and separate `DataFrame` per well, and put in a dictionary.We also build inline / crossline ranges (min/max) for each well here. This will give us the range toquery from the seismic data.<jupyter_code>well_names = well_df.index.levels[0].to_list()
num_wells = len(well_names)

wells = {}
il_ranges = {}
xl_ranges = {}
for well_name in well_names:
    well = well_df.loc[well_name]
    well_ils_xls = well[['inline', 'xline']]
    il_xl_min = well_ils_xls.min()
    il_xl_max = well_ils_xls.max()

    wells[well_name] = well
    il_ranges[well_name] = floor(il_xl_min[0]), ceil(il_xl_max[0])
    xl_ranges[well_name] = floor(il_xl_min[1]), ceil(il_xl_max[1])<jupyter_output><empty_output><jupyter_text>Pretty print well names and il/xl ranges<jupyter_code>for well in wells.keys():
    print(f"{well}:\tIL Range: {il_ranges[well]}\t\tXL Range: {xl_ranges[well]}")<jupyter_output>well_01:	IL Range: (2992, 2994)		XL Range: (2202, 2210)
well_11:	IL Range: (3523, 3524)		XL Range: (2397, 2401)
well_21:	IL Range: (2765, 2770)		XL Range: (1490, 1493)
well_25:	IL Range: (1193, 1196)		XL Range: (2512, 2518)
well_27:	IL Range: (2856, 2859)		XL Range: (2285, 2291)
well_33:	IL Range: (3179, 3185)		XL Range: (2337, 2348)<jupyter_text>Now, we need to extract seismic values around the wellbore. Let's start withmounting `real-simple-seismic` volumes of interest so we can query later.**Poseidon data vertical sampling is between 0 and 6,000 milliseconds.****Seismic is sampled every 4 ms, and velocities are sampled every 10 ms.**<jupyter_code>seismic_bucket = r's3://sagemaker-gitc2021/poseidon/seismic/rss/'

seismic_files = [
    'psdn11_TbsdmF_Near_Nov_11_32bit',  # near angle stack (AVO)
    'psdn11_TbsdmF_Mid_Nov_11_32bit',  # mid angle stack (AVO)
    'psdn11_TbsdmF_Far_Nov_11_32bit',  # far angle stack (AVO)
    'Final_PSDM_intVel_gridded_D2T_10ms',  # PSDM interval velocities in time
]

# We loop through all files, and make connections to all four.
# This takes about a minute.
# The handles are later used for querying.
rss_handles = []
for seismic_file in seismic_files:
    current_rss = rssFromS3(filename=seismic_bucket + seismic_file, client_kwargs={})

    rss_handles.append(current_rss)<jupyter_output>Connection complete.<jupyter_text>Now we can query only the traces around the well (makes a miniature 3d cube, which is usually called a minicube).We do this for each well and for every cube.The dimensions of the minicubes are:1. inlines2. crosslines3. z_valuesSo if you want to slice the first crossline from the minicube you would do `cube_data[:, 0, :]` or if you want tolook at the 3rd inline, you would do `cube_data[2, :, :]`<jupyter_code>from itertools import product
import numpy as np

# let's have nicer names for seismic data
seismic_names = ('near', 'mid', 'far', 'velocity')
# seismic_names = ('near', 'mid', 'far')


well_seismic = {}
for well in wells.keys():
    il_min, il_max = il_ranges[well]
    xl_min, xl_max = xl_ranges[well]

    ils = range(il_min, il_max + 1, 1)  # inline increment is 1
    xls = range(xl_min, xl_max + 1, 1)  # xline increment is 1

    # Use product to have all il/xl pairs within our ranges
    ilxl_pairs = list(product(ils, xls))

    cube_data = {}
    # In rss, for minicube extraction we query trace by trace.
    for rss_handle, seismic_name in zip(rss_handles, seismic_names):
        tmp = []
        for trace_ilxl in ilxl_pairs:
            tmp.append(rss_handle.trace(*trace_ilxl)[0])  # rss returns values and a live mask, we don't need the mask now

        # Traces aren't in a "cube" shape, they come sequential. So we reshape!
        cube_data[seismic_name] = np.asarray(tmp).reshape(len(ils), len(xls), -1)

    well_seismic[well] = cube_data<jupyter_output><empty_output><jupyter_text>Now we have all the traces around the well. Next, lets load up the horizons so we can display all of them together.<jupyter_code>%%capture

horizon_bucket = 's3://sagemaker-gitc2021/poseidon/horizons/'
horizon_files = ['water_bottom', 'top_heywood', 'top_johnson', 'top_jameison', 'top_near_plover']

horizons = []
for horizon in horizon_files:
    horizons.append(
        read_csv(filepath_or_buffer=horizon_bucket + f'{horizon}.csv.gz',
                 index_col=[0, 1],
                 compression='gzip'))

horizons = concat(horizons, axis=1)
horizons<jupyter_output><empty_output><jupyter_text>Now that we loaded all the horizons, wells, and seismic minicubes, let's plot all of them together.*Note some of the missing logs. To augment or not to augment.* Some key comments1. We have different ranges for seismic vs. velocity (vmin/vmax).2. We filter the horizon `DataFrame` in the minicube range using a `pd.IndexSlice` instance,which is a convenient way to slice a `DataFrame`.3. We are showing the **center** inline and all crosslines for the wellbore.4. Deviations look exaggerated because of the small lateral scale.5. We can now extract seismic traces **along** the wellbores (linear interpolator.or nearest neighbor interpolation + averaging).6. Typical regression models will map from seismic -> well data using values extracted at seismic(s).7. A CNN (deep neural network) based approach would use larger cubes around well log points for convolutions.8. Sorry, for some repetition in the code for labeling etc!<jupyter_code>import matplotlib.pyplot as plt

colors = ['gray_r']*3 + ['plasma']
vmins = [-39000]*3 + [1500]
vmaxs = [39000]*3 + [5500]

for well in wells.keys():
    il_min, il_max = il_ranges[well]
    xl_min, xl_max = xl_ranges[well]
    ils = range(il_min, il_max + 1, 1)  # inline increment is 1
    xls = range(xl_min, xl_max + 1, 1)  # inline increment is 1

    # We just want to plot the "center" inline and all crosslines
    num_il = il_max - il_min + 1
    plot_il_idx = num_il // 2

    fig, axes = plt.subplots(1, 7, sharey='all', figsize=(7, 10))

    # Let's filter horizon to the il/xl range we need

    minicube_horizons = horizons.loc[IndexSlice[il_min:il_max, xl_min:xl_max], :]

    for idx, (axis, seismic_name) in enumerate(zip(axes[:6], seismic_names)):
        seismic = well_seismic[well][seismic_name]
        # mean, std = np.nanmean(seismic), np.nanstd(seismic)
        axis.imshow(seismic[plot_il_idx, :, :].T, aspect='auto',
                    extent=[xl_min, xl_max, 6000, 0], interpolation='bicubic',
                    vmin=vmins[idx], vmax=vmaxs[idx], cmap=colors[idx])

        axis.set_xlabel('Xline')
        axis.plot(wells[well]['xline'], wells[well].index.values, 'r')
        axis.set_title(seismic_name)
        axis.tick_params(axis='x', rotation=30)
        hrz_lines = axis.plot(minicube_horizons.loc[ils[plot_il_idx]], lw=3)

    axes[4].plot(wells[well]['dtc'], wells[well].index.values, 'k', lw=1)
    axes[4].set_title('DTC'), axes[4].set_xlabel('DTC')
    axes[5].plot(wells[well]['dts'], wells[well].index.values, 'k', lw=1)
    axes[5].set_title('DTS'), axes[5].set_xlabel('DTS')
    axes[6].plot(wells[well]['rhob'], wells[well].index.values, 'k', lw=1)
    axes[6].set_title('RHOB'), axes[6].set_xlabel('RHOB')

    hrz_colors = [hrz_line.get_color() for hrz_line in hrz_lines]
    hrz_labels = horizons.keys()
    for idx, value in enumerate(minicube_horizons.loc[ils[plot_il_idx], round(np.mean(xls))].values):
        axes[4].axhline(value, lw=3, color=hrz_colors[idx])
        axes[4].tick_params(axis='x', rotation=30)
        axes[5].axhline(value, lw=3, color=hrz_colors[idx])
        axes[5].tick_params(axis='x', rotation=30)
        axes[6].axhline(value, lw=3, color=hrz_colors[idx], label=hrz_labels[idx])
        axes[6].tick_params(axis='x', rotation=30)

    handles, labels = axes[-1].get_legend_handles_labels()
    fig.legend(handles, labels, loc='lower center', ncol=len(hrz_labels),
               fontsize='x-small', markerscale=0.1, frameon=False)

    axes[0].set_ylabel('Two-way-time (ms)')
    axes[0].set_ylim([3500, 1500])
    fig.suptitle(well + f': inline {ils[plot_il_idx]}')
    fig.tight_layout()

    for axis in axes:
        axis.grid('minor')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12017512508934954, 'fraction_numerical': 0.025911365260900644, 'mean_word_length': 4.1819444444444445, 'pattern_counts': {'":': 0, '<': 30, '<?xml version=': 0, '>': 31, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '203524', 'license': 'Apache-2.0', 'n_tokens_mistral': 4091, 'n_tokens_neox': 3740, 'n_words': 1130, 'path': 'seismic_inversion_2021/tutorials/cloud_01_getAllData.ipynb', 'repo_name': 'jayden-mao/geoml_competitions', 'size': '11192'}