Import and Process FMI DLIS Using dlisio

# Install a pip package in the current Jupyter kernel
import sys
!{sys.executable} -m pip install dlisio

# Import libraries (Jupyter knows about these libraries so a special import is not required)
from dlisio import dlis
from tqdm import tqdm
import math
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import ipywidgets as widgets
from IPython.display import display
from scipy import interpolate
from scipy import ndimage

# Ensure that we can see pandas tables of up to 1000 rows as DLIS parameter lists can be very long
pd.set_option('display.max_rows', 1000)

filenames = []
filenames.append("Pasca A4(AD-1) Schlumberger FMI (Main Pass 8,131'-8,577') DLIS.dlis")
filenames.append("Pasca_A4(AD-1)_SLB_S2R8_FMI_MainPass_Static_Processed_Image_data.dlis")
filenames.append("Pasca_A4(AD-1)_SLB_S2R8_FMI_MainPass_Dynamic_Processed_Image_data.dlis")
import os.path
for fn in filenames:
    print(fn)
    print(os.path.exists(fn))

for fn in filenames:
    with dlis.load(fn) as files:
        print('Data contained in file: "', fn, '"\n')
        for f in files:
            print(f.describe())

with dlis.load(filenames[0]) as files:
    for f in files:
        for o in f.origins:
            print(o.describe())

# Load our two processed FMI DLIS logs
file = []
f, *tail = dlis.load(filenames[1])
lf, *tail = dlis.load(filenames[1])
file.append(lf)
lf, *tail = dlis.load(filenames[2])
file.append(lf)

desc_parameters = ['CN', 'WN', 'FN',                                     # Well ID
                   'LOND', 'LATD',                                       # Well location
                   'DLAB', 'TLAB']                                       # Time and date of well logging

# Create a table out of the parameters in desc_parameters
parameter_table = summarize(file[0].parameters, name='Name', long_name='Long name', values='Value(s)')
desc_table = parameter_table.loc[parameter_table['Name'].isin(desc_parameters)].copy()

# Sort the table in the same order as desc_parameters
categorical_sorter = pd.Categorical(desc_parameters, desc_parameters, ordered=True)
desc_table['Name'] = desc_table['Name'].astype(categorical_sorter.dtype)
desc_table.sort_values(by='Name', inplace=True)

# Display the table
display(desc_table)

for f in file:
    for frame in f.frames:

        # Search through the channels for the index and obtain the units
        for channel in frame.channels:
            if channel.name == frame.index:
                depth_units = channel.units

        print(f'Frame Name: \t\t {frame.name}')
        print(f'Index Type: \t\t {frame.index_type}')
        print(f'Depth Interval: \t {frame.index_min} - {frame.index_max} {depth_units}')
        print(f'Depth Spacing: \t\t {frame.spacing} {depth_units}')
        print(f'Direction: \t\t {frame.direction}')
        print(f'Num of Channels: \t {len(frame.channels)}')
        print(f'Channel Names: \t\t {str(frame.channels)}')
        print('\n')

# Static FMI channels
channel_table = summarize(file[0].channels, name='Channel Name', units='Units', dimension='Dimension', frame='Frame')
channel_table.sort_values('Channel Name')

def get_frame(logical_file, name):
    """
    Return a frame with a given name within a logical file
    """
    [frame] = [x for x in logical_file.frames if x.name == name]
    return frame

def index_of(frame):
    """
    Return the index channel of the frame
    """
    return next(ch for ch in frame.channels if ch.name == frame.index)

def get_channel(frame, name):
    """
    Get a channel with a given name from a given frame; fail if the frame does not have exactly one such channel
    """
    [channel] = [x for x in frame.channels if x.name == name]
    return channel

# Specify frames and channels to lookup
lookup_frame = ['B34276', 'B34276']
lookup_channel = ['FMI_STA', 'FMI_DYN']

# Set arrays to hold output
fmi = []
depth = []
fmi_image = []

# Iterate over the static and dynamic files
for i in range(2):

    # Get frame and its index channel
    fmi_frame = get_frame(file[i], lookup_frame[i])
    fmi_index = index_of(fmi_frame)

    # Get the metadata
    fmi.append(get_channel(fmi_frame, lookup_channel[i]))
    fmi_curves = fmi_frame.curves()
    depth.append(fmi_curves[fmi_frame.index])
    fmi_image.append(fmi_curves[lookup_channel[i]])

    # Print out the properties of the channel/curve
    print('Shape of depth index curve array: \t', depth[i].shape)
    print('Shape of FMI curve array:\t', fmi_image[i].shape)
    print(f'Name: \t\t{fmi[i].name}')
    print(f'Long Name: \t{fmi[i].long_name}')
    print(f'Units: \t\t{fmi[i].units}')
    print(f'Properties: \t{fmi[i].properties}')
    print(f'Dimension: \t{fmi[i].dimension}') #if >1, then data is an array
    print('\n')

def paint_channel(ax, curve, y_axis, x_axis, **kwargs):
    """
    Plot an image channel into a matplotlib axes using an index channel for the y-axis
    
    Parameters
    ----------
    
        ax : matplotlib.axes
        
        curve : numpy array
            The curve to be plotted
        
        y_axis : numpy array 
            An array of values that defines the y-axis
                
        x_axis : numpy array 
            An array of values that defines the x-axis
            
        **kwargs : dict 
            Keyword arguments to be passed on to ax.imshow()
    """
    # Determine the extent of the image so that the pixel centres correspond with the correct axis values
    dx = np.mean(x_axis[1:] - x_axis[:-1])
    dy = np.mean(y_axis[1:] - y_axis[:-1])
    extent = (x_axis[0] - dx/2, x_axis[-1] + dx/2, y_axis[0] - dy/2, y_axis[-1] + dy/2)
    
    # Determine the correct orientation of the image
    if y_axis[1] < y_axis[0]:   # Frame recorded from the bottom to the top of the well
        origin = 'lower'
    else:                       # Frame recorded from the top to the bottom of the well
        origin = 'upper'
    
    return ax.imshow(curve, aspect='auto', origin=origin, extent=extent, **kwargs)

def fmi_plot(MinDepth, MaxDepth, StaticFMI, DynamicFMI):   

    # Determine the maximum absolute value so that we can balance the colormap maximum value
    fmi_static_max = np.percentile(np.abs(fmi_image[0]), StaticFMI)
    fmi_dynamic_max = np.percentile(np.abs(fmi_image[1]), DynamicFMI)
    
    # Parameters for plotting the waveforms
    fmi_static_pltargs = {
        'cmap': 'YlOrBr',
        'vmin': 0,
        'vmax': fmi_static_max,
    }
    fmi_dynamic_pltargs = {
        'cmap': 'YlOrBr',
        'vmin': 0,
        'vmax': fmi_dynamic_max,
    }    
    
    fmi_samples = np.arange(fmi[0].dimension[0]) # x values to use in plotting

    # Create figure and axes
    fig, axes = plt.subplots(ncols=2, figsize=(10, 20), constrained_layout=True)

    # Plot static data as a colourmapped image
    ax = axes[0]
    im = paint_channel(ax, fmi_image[0], depth[0], fmi_samples, **fmi_static_pltargs)
    cbar = fig.colorbar(im, ax=ax, location='top')
    cbar.set_label(f'{fmi[0].name}: {fmi[0].long_name}')
    ax.set(ylim=[MaxDepth, MinDepth], ylabel='Depth $z$ [ft]')

    # Plot dynamic data as a colourmapped image
    ax = axes[1]
    im = paint_channel(ax, fmi_image[1], depth[1], fmi_samples, **fmi_dynamic_pltargs)
    cbar = fig.colorbar(im, ax=ax, location='top')
    cbar.set_label(f'{fmi[1].name}: {fmi[1].long_name}')
    ax.set(ylim=[MaxDepth, MinDepth])
    
    for ax in axes:
        ax.set(xlim=[0,360], xticks=np.arange(0, 360, 45), xlabel='Image Orientation [deg]')
        ax.grid(True)

fmi_plot(depth[0].min(), depth[0].max(), 50, 50)

# Set max depth for testing cropped image of delta_d == 0.4 ft
MaxDepth = 8570.5
MinDepth = MaxDepth - 0.4

# Get our cropped section for testing and show the range of index values
dint = 0.008333333333333226
d1 = int((8578.533333333327 - MaxDepth) / dint)
d2 = d1 + 49
d3 = d1 + 24
print('d1: ', d1, 'd2: ', d2, 'd3 (target row): ', d3)
fmi_orig_crop = np.copy(fmi_image[0][d1:d2, :])

# Display
nx = fmi_orig_crop[0].size
ny = len(fmi_orig_crop)
fig, axes = plt.subplots(figsize=(7.2, 3), constrained_layout=True)
im = paint_channel(axes, fmi_orig_crop, np.arange(ny)[::-1], np.arange(nx), **fmi_static_fixpltargs)
cbar = fig.colorbar(im, ax=axes, location='top')
cbar.set_label('fmi_orig_crop: original image')

# Interpolate using linear method after replacing negative values with NaN (this is to get low frequencies)
fmi_linear = np.copy(fmi_orig_crop)
for j in range(ny):
    r = 0
    line = np.copy(fmi_orig_crop[j])
    while line[0] < 0 or line[nx-1] < 0:
        r -= 1
        line = np.roll(fmi_orig_crop[j], r) # roll the values so the start and end of each row are not blank
    count = 0
    for i in range(nx):
        if line[i] < 0:
            count += 1
            line[i] = np.nan
    missing_data_df = pd.DataFrame(line, columns=['value'])
    missing_data_df.interpolate(method='linear', fill_value='extrapolate', limit_direction='both', inplace=True)
    fmi_linear[j] = missing_data_df.to_numpy().flatten()
    fmi_linear[j] = np.roll(fmi_linear[j], -r) # roll back if necessary

# Display code ommitted from snippet

# Interpolate using Akima method after replacing negative values with NaN (this is to get low frequencies)
fmi_akima = np.copy(fmi_orig_crop)
for j in range(ny):
    r = 0
    line = np.copy(fmi_orig_crop[j])
    while line[0] < 0 or line[nx-1] < 0:
        r -= 1
        line = np.roll(fmi_orig_crop[j], r) # roll the values so the start and end of each row are not blank
    count = 0
    for i in range(nx):
        if line[i] < 0:
            count += 1
            line[i] = np.nan
    missing_data_df = pd.DataFrame(line, columns=['value'])
    missing_data_df.interpolate(method='akima', fill_value='extrapolate', limit_direction='both', inplace=True)
    fmi_akima[j] = missing_data_df.to_numpy().flatten()
    fmi_akima[j] = np.roll(fmi_akima[j], -r) # roll back if necessary

# Display code ommitted from snippet

# Compress the image by omitting negative values
fmi_compress = []
max_len = 0
for j in range(ny):
    fmi_compress_line = []
    for val in fmi_orig_crop[j]:
        if val >= 0:
            fmi_compress_line.append(val)
    fmi_compress.append(fmi_compress_line)
    line_len = len(fmi_compress_line)
    if line_len > max_len:
        max_len = line_len
        
# Now resample every line to make it the same length
for j in range (ny):
    resamp_f = interpolate.interp1d(np.arange(len(fmi_compress[j])), fmi_compress[j], fill_value='extrapolate')
    fmi_compress[j] = resamp_f(np.arange(max_len))
    
# Display code ommitted from snippet

# Fill in gaps in the original cropped image using a mirroring approach
fmi_fill = []
for j in range(ny):
    r = 0
    fmi_fill_line = np.copy(fmi_orig_crop[j])
    while fmi_fill_line[0] < 0 or fmi_fill_line[nx-1] < 0:
        fmi_fill_line = np.roll(fmi_orig_crop[j], r)
        r -= 1
    i = 0
    ii = 0
    for val in fmi_fill_line:
        if val >= 0:
            i += 1
            ii = i
        else:
            idx = ii - 1
            if (idx < 0):
                idx += 360
            val = fmi_fill_line[idx]
            fmi_fill_line[i] = val
            i += 1
            ii -= 1
    fmi_fill_line = np.roll(fmi_fill_line, -r)
    fmi_fill.append(fmi_fill_line)
    
# Display code ommitted from snippet

def heal_fmi(fmi_orig):
    nx = fmi_orig.size
    
    # Unadjusted FFT of chosen row
    t = np.arange(nx)
    sp = np.fft.fft(fmi_orig)
    freq = np.fft.fftfreq(t.shape[-1])

    # Interpolate using linear method after replacing negative values with NaN (this is to get low frequencies)
    r = 0
    line = np.copy(fmi_orig)
    while line[0] < 0 or line[nx-1] < 0:
        r -= 1
        if (r < -nx):
            return fmi_orig
        line = np.roll(fmi_orig, r) # roll the values so the start and end of each row are not blank
    count = 0
    for i in range(nx):
        if line[i] < 0:
            count += 1
            line[i] = np.nan
    if count == nx:
        return fmi_orig
    missing_data_df = pd.DataFrame(line, columns=['value'])
    missing_data_df.interpolate(method='linear', fill_value='extrapolate', limit_direction='both', inplace=True)
    fmi_interp = missing_data_df.to_numpy().flatten()
    fmi_interp = np.roll(fmi_interp, -r) # roll back if necessary    
    
    # FFT to get low frequencies
    sp_interp = np.fft.fft(fmi_interp)

    # Remove negative values which indicate there are no data (this is to get high frequencies)
    fmi_compress = []
    for val in fmi_orig:
        if val >= 0:
            fmi_compress.append(val)

    # FFT to get high frequencies and then resample
    t_compress = np.arange(len(fmi_compress))
    freq_compress = np.fft.fftfreq(t_compress.shape[-1])
    sp_compress = np.fft.fft(fmi_compress)
    resamp_f = interpolate.interp1d(freq_compress, sp_compress, fill_value='extrapolate')
    sp_compress_resamp = resamp_f(freq)

    # Combine low and high frequencies
    sp_healed = np.copy(sp_interp)
    for i in range(sp_interp.size):
        if abs(freq[i]) > 0.05 and abs(freq[i] < 0.5):
            sp_healed[i] = sp_compress_resamp[i]
    fmi_healed = np.fft.ifft(sp_healed)
    
    # Splice original and healed data (we don't want to overwrite the original values)
    fmi_splice = np.copy(fmi_orig)
    for i in range(nx):
        if fmi_orig[i] < 0:
            fmi_splice[i] = max(fmi_healed[i].real, 0)
            
    # Return our healed FMI
    return fmi_splice

test_row = 15812
print('d: ', test_row)

# Original row without negative values
fmi_nan = np.copy(fmi_image[0][test_row])
nx = fmi_nan.size
t = np.arange(nx)
for i in range(fmi_image[0][test_row].size):
    if fmi_image[0][test_row][i] < 0:
        fmi_nan[i] = np.nan

# Heal using 1D FFT
fmi_heal = heal_fmi(fmi_image[0][test_row])
plt.plot(t, fmi_heal, t, fmi_nan)
plt.show()

test_row = 15812
d1 = test_row - 24
d2 = test_row + 25
print('d1: ', d1, 'd2: ', d2)
fmi_orig_crop = np.copy(fmi_image[0][d1:d2, :])

# Original row without negative values
fmi_nan = np.copy(fmi_image[0][test_row])
nx = fmi_nan.size
t = np.arange(nx)
for i in range(fmi_image[0][test_row].size):
    if fmi_image[0][test_row][i] < 0:
        fmi_nan[i] = np.nan

# Heal using 2D FFT
fmi_heal2 = heal_fmi2(fmi_orig_crop)
plt.plot(t, fmi_heal2[24], t, fmi_nan)
plt.show()

fmi_static = np.copy(fmi_image[0])
ny = len(fmi_image[0])
for row in tqdm(range(ny), 'Healing FMI Image (1D FFT)'):
    fmi_static[row] = heal_fmi(fmi_image[0][row])

fmi_static2 = np.copy(fmi_image[0])
ny = len(fmi_image[0])
for row in tqdm(range(ny), 'Healing FMI Image (2D FFT)'):
    if (row < 24):
        d1 = 0
        d2 = 49
        d3 = row
    elif (row > (ny - 24)):
        d1 = ny - 49
        d2 = ny
        d3 = 49 - (ny - row)
    else:
        d1 = row - 24
        d2 = d1 + 49
        d3 = 24
    fmi_orig_crop = np.copy(fmi_image[0][d1:d2, :])
    fmi_heal_crop = heal_fmi2(fmi_orig_crop)
    fmi_static2[row] = fmi_heal_crop[d3]

def heal_fmi_final(fmi_orig):
    nx = fmi_orig.size
    
    # Unadjusted FFT of chosen row
    t = np.arange(nx)
    sp = np.fft.fft(fmi_orig)
    freq = np.fft.fftfreq(t.shape[-1])

    # Interpolate using linear method after replacing negative values with NaN (this is to get low frequencies)
    r = 0
    line = np.copy(fmi_orig)
    while line[0] < 0 or line[nx-1] < 0:
        r -= 1
        if (r < -nx):
            return fmi_orig
        line = np.roll(fmi_orig, r) # roll the values so the start and end of each row are not blank
    count = 0
    for i in range(nx):
        if line[i] < 0:
            count += 1
            line[i] = np.nan
    if count == nx:
        return fmi_orig
    missing_data_df = pd.DataFrame(line, columns=['value'])
    missing_data_df.interpolate(method='linear', fill_value='extrapolate', limit_direction='both', inplace=True)
    fmi_linear = missing_data_df.to_numpy().flatten()
    fmi_linear = np.roll(fmi_linear, -r) # roll back if necessary    
    
    # Interpolate using Akima method after replacing negative values with NaN (this is to get low frequencies)
    r = 0
    line = np.copy(fmi_orig)
    while line[0] < 0 or line[nx-1] < 0:
        r -= 1
        if (r < -nx):
            return fmi_orig
        line = np.roll(fmi_orig, r) # roll the values so the start and end of each row are not blank
    count = 0
    for i in range(nx):
        if line[i] < 0:
            count += 1
            line[i] = np.nan
    if count == nx:
        return fmi_orig
    missing_data_df = pd.DataFrame(line, columns=['value'])
    missing_data_df.interpolate(method='akima', fill_value='extrapolate', limit_direction='both', inplace=True)
    fmi_akima = missing_data_df.to_numpy().flatten()
    fmi_akima = np.roll(fmi_akima, -r) # roll back if necessary    

    # Combine the linear and Akima interpolations
    fmi_interp = np.copy(fmi_linear)
    for i in range(fmi_linear.size):
        fmi_interp[i] *= 0.8
        fmi_interp[i] += 0.2 * fmi_akima[i]
    
    # FFT to get low frequencies
    sp_interp = np.fft.fft(fmi_interp)

    # Remove negative values which indicate there are no data (this is to get high frequencies)
    fmi_compress = []
    for val in fmi_orig:
        if val >= 0:
            fmi_compress.append(val)

    # FFT to get high frequencies and then resample
    t_compress = np.arange(len(fmi_compress))
    freq_compress = np.fft.fftfreq(t_compress.shape[-1])
    sp_compress = np.fft.fft(fmi_compress)
    resamp_f = interpolate.interp1d(freq_compress, sp_compress, fill_value='extrapolate')
    sp_compress_resamp = resamp_f(freq)
    
    # Fill in gaps in the original cropped image using a mirroring approach
    r = 0
    fmi_fill_line = np.copy(fmi_orig)
    while fmi_fill_line[0] < 0 or fmi_fill_line[nx-1] < 0:
        fmi_fill_line = np.roll(fmi_orig, r)
        r -= 1
    i = 0
    ii = 0
    for val in fmi_fill_line:
        if val >= 0:
            i += 1
            ii = i
        else:
            idx = ii - 1
            if (idx < 0):
                idx += 360
            val = fmi_fill_line[idx]
            if (val < 0):
                val = fmi_fill_line[i - 1]
            fmi_fill_line[i] = val
            i += 1
            ii -= 1
    fmi_fill_line = np.roll(fmi_fill_line, -r)  
    
    # FFT to get high frequencies for mirror filled line and blend with compressed high frequencies
    sp_texture = np.fft.fft(fmi_fill_line)
    for i in range(sp_texture.size):
        sp_texture[i] *= 0.8
        sp_texture[i] += 0.2 * sp_compress_resamp[i]
        
    # Combine low and high frequencies
    sp_healed = np.copy(sp_interp)
    for i in range(sp_interp.size):
        if abs(freq[i]) > 0.05 and abs(freq[i] < 0.5):
            sp_healed[i] = sp_texture[i]
    fmi_healed = np.fft.ifft(sp_healed)
    
    # Splice original and healed data (we don't want to overwrite the original values)
    fmi_splice = np.copy(fmi_orig)
    for i in range(nx):
        if fmi_orig[i] < 0:
            fmi_splice[i] = max(fmi_healed[i].real, 0)
            
    # Return our healed FMI
    return fmi_splice

fmi_static = np.copy(fmi_image[0])
ny = len(fmi_image[0])
for row in tqdm(range(ny), 'Healing FMI Image (Final 1D FFT)'):
    fmi_static[row] = heal_fmi_final(fmi_image[0][row])