<jupyter_start><jupyter_text>Session 1: Seismic VolumesIn this session, you will learn how to **import, pre-process and view seismic volumes in python**. It also introduces you to the basics of using iPython notebooks. This sessions is suitable both for those **new to python programming** and those who are **new to the field of seismic interpretation**. However, even if you are already an expert with both, you might still learn something new or interesting from this session. Basics: the segyio library and  the .segy file formatThis iPython notebook will guide you through the session and the code examples. An iPython notebook allows you to interactlively write and execute python code within your browser.To begin, simply **run the code cell below**, by left-clicking on it and then clicking on the "Arrow" run button. Then wait for the spinning circle around the button to disappear, which means the cell has finished executing.<jupyter_code>import os
import sys
import time
import argparse
import numpy as np
import matplotlib.pyplot as plt<jupyter_output><empty_output><jupyter_text>In the cell above, we are importing some standard python libraries (os, sys, time and argparse), as well as the commonly used third party libraries **Numpy** and **Matplotlib**. Numpy is a powerful library that helps with all sort of numerical computations on matrices. Matplotlib, more specifically it's pyplot package, includes functionality for visualizing data. All these packages come pre-installed in google collab.The most commonly used data format for handling seismic volume data is **.segy**, though sometimes .sgy is used interchangeably. In order to open this binary file format in python, we are going to use a library called **segyio**. Since segyio is not included as a standard library in google collab, we have to download it using the python package installer (pip).Install segyio by **executing the cell below** and wait for it to finish, before continuing.<jupyter_code>!{sys.executable} -m pip install segyio
import segyio<jupyter_output><empty_output><jupyter_text>Using the code in the cell below, we will download a **sample segy file from github**. It is provided by **Equinor**, the authors of the segyio library. We are going to use this sample file for the rest of this notebook, so please **run the following cell**.<jupyter_code>!wget https://github.com/equinor/segyio-notebooks/blob/master/data/basic/F3_Dip_steered_median_subvolume_IL230-430_XL475-675_T1600-1800.sgy?raw=true -nv -O sample.sgy
segyfile = "sample.sgy"<jupyter_output><empty_output><jupyter_text>Opening and Accessing Seismic DataNow that we have a segy file, we can open and access it using segyio. Let's start by **defining a python function** that prints out information about the seismic volume included in our segyfile. **Run the following cell now**, in order to load the function into memory. After that, we'll inspect the code in detail.<jupyter_code>def print_info(segyfile):
    with segyio.open(segyfile, mode = 'r', iline=189, xline=193, strict=True, ignore_geometry=False) as segy:
        width = len(segy.ilines)
        length = len(segy.xlines)
        depth = len(segy.samples)
        dtype = segy.format
        print(f'Width: {width}\nLength: {length}\nDepth: {depth}\nNumber format: {dtype}')<jupyter_output><empty_output><jupyter_text>First we open our sample segy file using the **segyio.open()** function. There are a number of parameters to consider:* The first parameter is the **path to the segy file** in our filesystem.* We set mode to 'r' to open the file in **read-only mode**, since we do not want to change the file's content.* Inside segy files, seismic data is not organized as a 3D volume, but in a list of so called **traces**. Each of these traces is a single line of voxels reaching from the top the bottom of the survey. Therefore, the depth of our volume is equal to the length of a single trace. The width and length of the volume, howerver, is not directly inferrable from the data. Fortunately, most segy files include this information in their binary file headers. With the parameters **inline and xline**, we tell segyio where it has to look. The standard location is **byte 189** for width (iline) and **byte 193** for length (xline)* We set **strict to True** and **ignore_geometry to False**. In combination, this tells segyio that we expect our data to be a strict rectangular cuboid.Having opened the file, we use segyio to give us some information about it. With **segy.ilines**, **segy.xlines** and **segy.samples** we query the widths, lengths and depth of our 3D volume. With **segy.format**, we query the data type in which the numbers are saved.Run the cell below now, in order to pass our segyfile to the function and execute it.<jupyter_code>print_info(segyfile)<jupyter_output><empty_output><jupyter_text>We now have a better idea how our 3D volume looks.Apparently, it has a width of 201, a length of 201 and a depth of 51. The data consists of IEEE floating point numbers.In order to learn a bit more about the data, we will **define a function** that returns tha **mimimum and maxumim value** within our volume. We shall also use this opportunity to learn about accessing the data with segyio.Load the function into memory, by **running the cell below**.<jupyter_code>def find_value_range(segyfile):
    with segyio.open(segyfile, mode = 'r', iline=189, xline=193, strict=True, ignore_geometry=False) as segy:
        gmin, gmax = float('inf'), -1*float('inf')
        
        for f in segy.fast:
            min = np.min(f)
            if min < gmin:
                gmin = min
            max = np.max(f)
            if max > gmax:
                gmax = max
        return gmin, gmax<jupyter_output><empty_output><jupyter_text>Again, we use **segyio.open** to access our segyfile. We use a **standard algorithm** for finding the value range of an array,  which is to update global minimum and maxumim values while traversing the data in a loop.What is not standard is the way we access the data. Here, we used the  **"fast" attribute** of our opened segy file. Depending on the memory layout of the data, **one dimension is faster to traverse**. If we look at our 3D volume as a 2D array of 1D traces, then either all traces belonging to the same row are saved adjacent in memory, or all traces belonging to the same column. By calling the "fast" attribute on our opened segy file, we can directly access the data in the way that is fastest to traverse. Now **run the cell below** in order to execute the function on our segyfile.<jupyter_code>start=time.time()
min_val, max_val = find_value_range(segyfile)
finish=time.time()
print("Found value range: \nMin: {} \nMax: {} \nTook {:.2f} seconds".format(min_val, max_val, (finish-start)))<jupyter_output><empty_output><jupyter_text>Interactive ExerciseTry out the different modes of data access in segyio. * **Replace segy.fast with segy.slow** in the *find_value_range* function, then execute the last two cells again. Do you see a difference in computation time? * Now **try both segy.iline and segy.xline** instead. Which one is the fast and which one is the slow direction?* Now try accessing the file in trace mode using **segy.trace**. Is it faster or slower that the other modes?  Histograms and NormalizationIn this section, we are going to take a closer look at the data in our segy file. So far, we know about the value range, but we have no idea about the distribution of values.We start by defining a **function that visualizes the histogram** of a given array in a bar plot. Making use of the **numpy.unique()** function, we can access the unique values and their corresponfing count in an array. The **pyplot.bar()** function visualizes the histogram on bar graph. Run the cell below to load the method into memory.<jupyter_code>def plot_histogram(array):
    values, counts = np.unique(array, return_counts=True)
    plt.figure(figsize=(10,5))
    plt.bar(values, counts)
    plt.show()<jupyter_output><empty_output><jupyter_text>**Execute the cell below** in order to visualize the value histogram of the data in our sample segy file.<jupyter_code>with segyio.open(segyfile, mode = 'r', iline=189, xline=193, strict=True, ignore_geometry=False) as segy:
    array = np.array(segy.fast)
    plot_histogram(array)<jupyter_output><empty_output><jupyter_text>As you can see, the data roughly resembles a Gaussian Distribution, centered around zero. However, due to the large value range and the large number of unique values, the pyplot bar function has trouble properly visualizing the histogram.For the purpose of visualizing and processing our data, it might be desirable to **normalize our data to a smaller value range** with a smaller number of unique values. In the cell below, we are defining a method that normalizes the data in a segy file into a set number of bin values. First, we find the value range using our own **find_value_range()** function. Then we use the function **segyio.tools.cube()**. This method will convert the data from a list of traces into a 3D array, so we'll no longer have to worry about fast and slow access modes. We initialize a new empty output array of the same shape as the original volume using the **numpy.zeros()** function. We then traverse the segy data in a nested for loop, **normalize each value** to the new value range and **save the result to the output array**.**Run the cell below** in order to load the function into memory.<jupyter_code>def normalized_array(segyfile, bins):
    min_val, max_val = find_value_range(segyfile)
    segy = segyio.tools.cube(segyfile)
    width = segy.shape[0]
    length = segy.shape[1]
    depth = segy.shape[2]
    output_array = np.zeros(segy.shape,dtype=np.dtype(int))   
    bin_width = (max_val - min_val) / (bins-1)
    for x in range(width):
        for y in range(length):
            for z in range(depth):
                v = segy[x,y,z]
                v_norm = int((v - min_val) / bin_width)
                output_array[x,y,z]=v_norm
    return output_array<jupyter_output><empty_output><jupyter_text>In image processing, the most commonly used data are **grayscale images** with a **value range of 0 to 255**. This range will allow most image processing libraries to open, display and save seismic data as images.**Run the cell below** in oder to normalize our data to range of 0 to 255 and plot the resulting histogram.<jupyter_code>volume = normalized_array(segyfile, 256)
plot_histogram(volume)<jupyter_output><empty_output><jupyter_text>Note how the resulting histogram looks much smoother and more closely resembles a Gaussian Distribution. For some use cases it might be desirable to **center the data** and perform operations such as **histogram linearization** in order to improve the value disribution which in turn **enhances contrast** in the resulting seismic images. For the rest of this session, however,  we are going to leave it here and continue with the normalized array as-is.<jupyter_code>def get_oriented_slice(volume, orientation, index):
    try:
        if orientation == 'inline':
            return volume[index,:,:].T
        elif orientation == 'crossline':
            return volume[:,index,:].T
        elif orientation == 'timeslice':
            return volume[:,:,index]
        else:
            print("Error: Orientation must be one of 'inline', 'crossline, 'timeslice'")
            return np.array([])
    except IndexError:
        print(f"Index {index} is out of bounds for dimension {orientation}")
        return np.array([])

def plot_slice(slice, colormap):
    try:
        plt.figure(figsize=(8,8))
        plt.imshow(slice, cmap=colormap)
    except ValueError:
        print(f"{colormap} is not a valid colormap. Try 'gray' or 'seismic' for example.")<jupyter_output><empty_output><jupyter_text>Interactive Excercise**Execute the cell below** in order to display a single inline slice from our normalized seismic volume as a gray image, then try to play with different values for direction, index and colormap.* By **changing the direction** allows you to access the three dimensions of the volume.* By **changing the index**, you can access different slices in each dimenson. Make sure not o leave the bounds of each dimension!* The colormap changes in which colors the slice is displayed. **Try the 'seismic' colormap** for example. More valid entries for colormap with reference can be found found at https://matplotlib.org/stable/gallery/color/colormap_reference.html.<jupyter_code>#----Try Changing these values--------
direction = 'timeslice'
index = 25
colormap = 'gray'
#-------------------------------------

seismic_slice = get_oriented_slice(volume, direction,index)
if seismic_slice.any():
    plot_slice(seismic_slice, colormap)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08187271169276311, 'fraction_numerical': 0.006621484770585028, 'mean_word_length': 4.276613234689684, 'pattern_counts': {'":': 0, '<': 54, '<?xml version=': 0, '>': 54, 'https://': 2, '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': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '309097', 'license': 'MIT', 'n_tokens_mistral': 3627, 'n_tokens_neox': 3483, 'n_words': 1700, 'path': 'session_01.ipynb', 'repo_name': 'vrgeo/ml-tutorials', 'size': '12837'}