Why Whole-Well Completeness Numbers Can Hide Critical Data Gaps

Visualising Well Log Data Availability by Formation Across Multiple Wells

In my previous article, I walked through transforming a default matplotlib heatmap into something publication-ready. The result was a clean figure showing which curves were available across a set of wells.

But there is a problem with that approach.

A curve showing 75% completeness across an entire well sounds acceptable.

Until you realise the missing 25% sits entirely inside your reservoir.

Whole-well completeness numbers can hide the most important gaps in your dataset.

You can see that in my previous article here: *From Default to Publication-Ready: Transforming Matplotlib Heatmaps in 8 Steps.*

In this article, we build on that heatmap and take it one step further by filtering completeness by formation interval and displaying the results across multiple formations in a single figure.

By the end, you will have a visual tool that shows not just what data exists, but where in the well it actually exists.

Setting Up the Dataset

Before we start building, we need data that includes depth. Unlike the previous article where we only needed binary present or missing flag per well, here we need depth-indexed measurements so we can filter by formation interval.

For this example we will use a synthetic dataset. If you already have data you can skip ahead to the code to just before step 1. This will allow you to create the completeness matrix and carry on.

The synthetic data we are creating represents 8 wells with as realistic as possible depth coverage and curve availability patterns. Each curve has its own missingness behaviour. For example some curves tend to be less consistently acquired in shallower wells, PEF and RHOB often co-acquired, so their missingness can correlate, and SP is frequently absent in more modern wells. This is a result of different priorities when designing the logging programme.

import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
import numpy as np
import pandas as pd

np.random.seed(42)

wells = [
    'Ardmore-1',  'Ardmore-4',
    'Brindal-1',  'Brindal-3',  'Brindal-5',  'Brindal-8',
    'Caxton-2',   'Caxton-6',
]

curves = ['GR', 'RDEP', 'RSHA', 'RHOB', 'NPHI', 'DTC', 'DTS', 'PEF', 'CALI', 'SP']

# Field-average top and base depth for each formation (metres)
formation_avg = {
    'Thornfield':  (800,  1000),
    'Veldwick':    (1000, 1400),
    'Oskarby':     (1400, 1900),
    'Dravenport':  (1900, 2400),
}

formation_names = list(formation_avg.keys())

In a real project, each well has its own formation top and base picks rather than a single shared depth range. We model this by applying a small random depth shift to the field-average formation depths for each well, keeping formations internally consistent within a well while reflecting the variation you would see across a field.

rng = np.random.default_rng(42)

def make_well_formation_tops(wells, formation_avg, rng, jitter=50):
    well_tops = {}
    for well in wells:
        shift = rng.integers(-jitter, jitter)
        well_tops[well] = {
            fm: (top + shift, base + shift)
            for fm, (top, base) in formation_avg.items()
        }
    return well_tops

well_formation_tops = make_well_formation_tops(wells, formation_avg, rng)

We also need a function that generates depth-indexed availability data for each well. In a real workflow, this would come from your LAS files. You would load each well, check for non-null values at each depth sample, and calculate completeness within each formation interval.

def generate_well_data(well_idx):
    depth = np.arange(500, 3500, 0.5)
    n = len(depth)
    data = {}

    for curve in curves:
        presence = np.ones(n)
        if curve == 'DTS':
            presence[depth < 1600] = np.where(
                np.random.random(np.sum(depth < 1600)) > 0.4, 1, 0)
        elif curve == 'DTC':
            presence[depth < 1200] = np.where(
                np.random.random(np.sum(depth < 1200)) > 0.6, 1, 0)
        elif curve == 'PEF':
            presence[depth < 1500] = np.where(
                np.random.random(np.sum(depth < 1500)) > 0.5, 1, 0)
        elif curve == 'SP':
            presence = np.where(np.random.random(n) > 0.45, 1, 0)
        elif curve == 'RHOB':
            presence[depth < 1100] = np.where(
                np.random.random(np.sum(depth < 1100)) > 0.7, 1, 0)
        if well_idx % 3 == 0 and curve in ['NPHI', 'RHOB']:
            mid = len(depth) // 2
            presence[mid - 200:mid] = 0
        if well_idx % 4 == 0 and curve == 'RDEP':
            presence[:500] = 0
        data[curve] = presence
    return pd.DataFrame(data, index=depth)

Now we can calculate completeness per formation per well, using each well’s own formation tops rather than a single shared depth range.

results = {}
for formation in formation_names:
    comp_matrix = np.zeros((len(wells), len(curves)))
    for w_idx, well in enumerate(wells):
        top, base = well_formation_tops[well][formation]
        well_data = generate_well_data(w_idx)
        interval = well_data[
            (well_data.index >= top) & (well_data.index < base)
        ]
        if len(interval) > 0:
            comp_matrix[w_idx] = (interval.sum() / len(interval) * 100).values
    results[formation] = pd.DataFrame(comp_matrix, index=wells, columns=curves)

Step 1: Plotting the Raw Data with imshow()

Before adding any customisation, it helps to see what matplotlib gives you out of the box. We pass the completeness matrix (results) for a single formation directly to imshow() with no other arguments. As we have multiple formations we can select one if them and call upon .values.

fig, ax = plt.subplots(figsize=(10, 5))
ax.imshow(results['Veldwick'].values)
plt.show()

Default imshow output with no labels, wrong colourmap, unreadable

We get back a colour plot, but it has the same problems as before. No labels, a colourmap designed for continuous data rather than percentages, and no way to read individual cell values. This is our starting point.

Step 2: Applying a Colourmap That Communicates Quality

Before touching the colourmap, it is worth pausing on what the data actually represents. In the previous article, each cell held a simple 0 or 1: curve present or missing. Here, each cell holds a value between 0 and 100, representing the percentage of depth samples within that formation interval where the curve has a valid measurement. A cell showing 60% means you have data for just over half the interval, which may or may not be acceptable depending on where in the formation that gap sits.

In this context, completeness simply means the percentage of depth samples within the formation interval that contain a valid value for a given curve. In a real workflow you would normally replace LAS null values (for example -999.25) before calculating this percentage. It is important to note that completeness measures coverage, not data quality. A curve can be present at every depth sample and still contain poor measurements caused by washouts, tool sticking, or acquisition issues.

That shift in what the data represents also changes what the colourmap needs to do. The default viridis colourmap treats every value as equally interesting, which makes sense for continuous scientific data, but not for completeness percentages. Here, what matters is immediately separating the problem cells from the acceptable ones.

The RdYlGn diverging colourmap maps directly to how a petrophysicist reads a quality matrix: red flags a problem, yellow signals something worth checking, and green confirms the data is usable. We also set an explicit normalisation from 0 to 100 so the colour scale always represents the full percentage range, regardless of what values happen to appear in the data.

cmap = plt.cm.RdYlGn
norm = mcolors.Normalize(vmin=0, vmax=100)

fig, ax = plt.subplots(figsize=(10, 5))
ax.imshow(results['Veldwick'].values, cmap=cmap, norm=norm)
plt.show()

RdYlGn colourmap applied to show red/yellow/green tp communicate quality.

Step 3: Adding Percentage Labels to Each Cell

The colourmap gives you a quick visual triage, but colour alone is imprecise. Two cells that both look yellow could be 55% and 72%, which may mean very different things depending on your minimum data requirements for a given evaluation.

Adding text annotations to each cell makes the figure precise as well as readable. We also adjust the text colour dynamically: white text on dark red or dark green cells, black text on the lighter yellows and mid-greens in between. The threshold values of 35 and 80 were chosen to match the points where the background colour is dark enough to make black text hard to read.

fig, ax = plt.subplots(figsize=(10, 5))
im = ax.imshow(results['Veldwick'].values, cmap=cmap, norm=norm, aspect='auto')

for i in range(len(wells)):
    for j in range(len(curves)):
        val = results['Veldwick'].values[i, j]
        text_colour = 'white' if val < 35 or val > 80 else 'black'
        ax.text(j, i, f'{val:.0f}%',
                ha='center', va='center',
                fontsize=7, fontfamily='monospace',
                color=text_colour, fontweight='bold')

plt.show()

Cell annotations added to show exact percentages in each cell

Step 4: Adding Axis Labels for Wells and Curves

Without axis labels, the heatmap is uninterpretable outside the context of the code that generated it. We set tick positions to match the row and column indices of the matrix, then assign the well names and curve mnemonics as labels. The x-axis labels are rotated 45 degrees to prevent overlap, and both axes use a monospace font to keep the label alignment clean.

ax.set_xticks(range(len(curves)))
ax.set_xticklabels(curves, rotation=45, ha='right',
                   fontsize=9, fontfamily='monospace')
ax.set_yticks(range(len(wells)))
ax.set_yticklabels(wells, fontsize=8, fontfamily='monospace')

Axis labels added and both wells and curves can be easily identified

Step 5: Adding Cell Borders with Minor Grid Lines

At this point the cells run together visually, making it harder to isolate individual well-curve combinations. Adding white grid lines on the minor tick positions draws a clear border around each cell without introducing heavy chart furniture. We also remove the outer spines, which gives the figure a cleaner, more modern appearance and lets the grid lines define the structure instead.

ax.set_xticks(np.arange(-0.5, len(curves), 1), minor=True)
ax.set_yticks(np.arange(-0.5, len(wells), 1), minor=True)
ax.grid(which='minor', color='white', linewidth=1.5)
ax.tick_params(which='minor', bottom=False, left=False)

for spine in ax.spines.values():
    spine.set_visible(False)

Gridlines and clean spines help clear up the cells and separates them.

Step 6: Building the Multi-Formation Panel

This is the step that makes the figure genuinely useful for formation-level analysis. Instead of showing completeness for a single formation, we use plt.subplots() to create a stacked panel with one heatmap per formation. Each panel shares the same colour scale and layout, so comparisons across formations are direct.

The x-axis labels are only shown on the bottom panel to avoid repetition, and each panel gets a formation name label on the left using ax.text with transform=ax.transAxes, which positions it relative to the axes rather than the data. The hspace parameter in subplots_adjust is kept tight so the panels read as a single connected figure rather than separate charts, and left is set wide enough to give the formation labels room outside the plot area.

You can now scan down any column and see immediately whether a curve’s completeness holds across all formations or collapses in a specific zone.

n_formations = len(formation_names)
fig, axes = plt.subplots(n_formations, 1, figsize=(12, 14), facecolor='white')

for idx, (ax, formation) in enumerate(zip(axes, formation_names)):
    df = results[formation]
    ax.imshow(df.values, cmap=cmap, norm=norm, aspect='auto')
    ax.set_xticks(np.arange(-0.5, len(curves), 1), minor=True)
    ax.set_yticks(np.arange(-0.5, len(wells), 1), minor=True)
    ax.grid(which='minor', color='white', linewidth=1.5)
    ax.tick_params(which='minor', bottom=False, left=False)
    ax.set_yticks(range(len(wells)))
    ax.set_yticklabels(wells, fontsize=8, fontfamily='monospace')
    ax.set_xticks(range(len(curves)))
    if idx == n_formations - 1:
        ax.set_xticklabels(curves, rotation=45, ha='right',
                           fontsize=9, fontfamily='monospace')
    else:
        ax.set_xticklabels([])
    for i in range(len(wells)):
        for j in range(len(curves)):
            val = df.values[i, j]
            text_colour = 'white' if val < 35 or val > 80 else 'black'
            ax.text(j, i, f'{val:.0f}%', ha='center', va='center',
                    fontsize=7, fontfamily='monospace',
                    color=text_colour, fontweight='bold')
    ax.text(-0.15, 0.5, formation,
            transform=ax.transAxes,
            fontsize=10, fontweight='bold',
            rotation=0, va='center', ha='center')
    for spine in ax.spines.values():
        spine.set_visible(False)

plt.subplots_adjust(hspace=0.08, left=0.18, right=0.9, top=0.98)

Faceted multi-formation view with all four formations stacked vertically

Step 7: Adding a Colourbar and Figure Title

The final additions are a shared colourbar and a figure title. The colourbar is placed manually using fig.add_axes() so it sits outside the subplot grid without disturbing the panel layout. It maps directly to the same normalisation used across all panels, making it valid for the entire figure rather than any individual formation.

The suptitle sits above all panels and names the figure as a whole. Together, these elements turn the panel into something you could drop into a report or present to a team without additional explanation.

cbar_ax = fig.add_axes([0.92, 0.15, 0.015, 0.7])
cb = fig.colorbar(
    plt.cm.ScalarMappable(norm=norm, cmap=cmap),
    cax=cbar_ax
)
cb.set_label('Completeness (%)', fontsize=9)
cb.ax.tick_params(labelsize=8)

fig.suptitle('Well Log Data Completeness by Formation',
             fontsize=14, fontweight='bold', y=1.01)

Final polished figure — colourbar, title, and formation labels complete

Reading the Figure

Looking at the final output, a few patterns become immediately obvious that would have been hidden in a whole-well completeness view.

DTC and DTS completeness degrades significantly in the shallower Thornfield and Veldwick formations across most wells. If your workflow depends on acoustic data, you now know exactly which wells and which formations have the problem.

PEF shows a similar pattern, dropping off in shallower intervals across several wells. Since PEF and RHOB are often acquired on the same logging run, their missingness tends to track together.

GR, RDEP, and CALI are consistently green across all formations and wells, confirming these are reliably acquired regardless of depth.

This figure tells you which formation has a completeness problem. If you want to see exactly where within a single well those gaps sit, I have covered that in a previous article on Substack. Link at the end.

What This Figure Cannot Tell You

This figure shows completeness within predefined depth intervals. It does not tell you whether data quality within those intervals is good. A curve can be present at every depth sample but still be affected by washouts, tool stick, or other acquisition issues.

It also assumes consistent formation tops across all wells. In reality, formation depths vary between wells and the top and base values should come from your well-by-well picks rather than fixed values. The structure of the code makes that straightforward to implement. Instead of a single depth range per formation, you would pass a dictionary of per-well tops and bases.

Wrapping Up

Starting from a default imshow() call and ending with a faceted multi-formation completeness panel, each step in this article changed one thing and moved the figure closer to something genuinely useful.

The key shift from the previous article is moving from binary present/missing to percentage completeness within specific intervals. That change transforms the figure from a data inventory check into a tool that informs decisions.

Whilst the examples here use well log curves and geological formations, the same approach applies to any dataset where you need to assess completeness across two categorical dimensions within defined intervals. If you have depth-indexed measurements, time-series data split into phases, or any other grouped structure, the core pattern of filtering by interval and visualising completeness as a heatmap transfers well.

Knowing which formation has a data gap is the first step. The next question is where in the well that gap actually sits. A formation-level percentage tells you there is a problem, but it does not show you whether the missing data is clustered at the top of the interval, scattered throughout, or concentrated right at the zone of interest. Over on Substack I built a depth-based coverage heatmap that runs directly in the terminal using lasio and rich, so you can see the shape of your data before interpretation starts.

*→ Visualising Well Log Data Availability in the Terminal with Python*

Full Code

import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
import numpy as np
import pandas as pd

np.random.seed(42)

wells = [
    'Ardmore-1',  'Ardmore-4',
    'Brindal-1',  'Brindal-3',  'Brindal-5',  'Brindal-8',
    'Caxton-2',   'Caxton-6',
]
curves = ['GR', 'RDEP', 'RSHA', 'RHOB', 'NPHI', 'DTC', 'DTS', 'PEF', 'CALI', 'SP']

formation_avg = {
    'Thornfield':  (800,  1000),
    'Veldwick':    (1000, 1400),
    'Oskarby':     (1400, 1900),
    'Dravenport':  (1900, 2400),
}
formation_names = list(formation_avg.keys())

rng = np.random.default_rng(42)

def make_well_formation_tops(wells, formation_avg, rng, jitter=50):
    well_tops = {}
    for well in wells:
        shift = rng.integers(-jitter, jitter)
        well_tops[well] = {
            fm: (top + shift, base + shift)
            for fm, (top, base) in formation_avg.items()
        }
    return well_tops

well_formation_tops = make_well_formation_tops(wells, formation_avg, rng)

def generate_well_data(well_idx):
    depth = np.arange(500, 3500, 0.5)
    n = len(depth)
    data = {}
    for curve in curves:
        presence = np.ones(n)
        if curve == 'DTS':
            presence[depth < 1600] = np.where(
                np.random.random(np.sum(depth < 1600)) > 0.4, 1, 0)
        elif curve == 'DTC':
            presence[depth < 1200] = np.where(
                np.random.random(np.sum(depth < 1200)) > 0.6, 1, 0)
        elif curve == 'PEF':
            presence[depth < 1500] = np.where(
                np.random.random(np.sum(depth < 1500)) > 0.5, 1, 0)
        elif curve == 'SP':
            presence = np.where(np.random.random(n) > 0.45, 1, 0)
        elif curve == 'RHOB':
            presence[depth < 1100] = np.where(
                np.random.random(np.sum(depth < 1100)) > 0.7, 1, 0)
        if well_idx % 3 == 0 and curve in ['NPHI', 'RHOB']:
            mid = len(depth) // 2
            presence[mid - 200:mid] = 0
        if well_idx % 4 == 0 and curve == 'RDEP':
            presence[:500] = 0
        data[curve] = presence
    return pd.DataFrame(data, index=depth)

results = {}
for formation in formation_names:
    comp_matrix = np.zeros((len(wells), len(curves)))
    for w_idx, well in enumerate(wells):
        top, base = well_formation_tops[well][formation]
        well_data = generate_well_data(w_idx)
        interval = well_data[
            (well_data.index >= top) & (well_data.index < base)
        ]
        if len(interval) > 0:
            comp_matrix[w_idx] = (interval.sum() / len(interval) * 100).values
    results[formation] = pd.DataFrame(comp_matrix, index=wells, columns=curves)

cmap = plt.cm.RdYlGn
norm = mcolors.Normalize(vmin=0, vmax=100)

n_formations = len(formation_names)
fig, axes = plt.subplots(n_formations, 1, figsize=(12, 14), facecolor='white')

for idx, (ax, formation) in enumerate(zip(axes, formation_names)):
    df = results[formation]
    ax.imshow(df.values, cmap=cmap, norm=norm, aspect='auto')
    ax.set_xticks(np.arange(-0.5, len(curves), 1), minor=True)
    ax.set_yticks(np.arange(-0.5, len(wells), 1), minor=True)
    ax.grid(which='minor', color='white', linewidth=1.5)
    ax.tick_params(which='minor', bottom=False, left=False)
    ax.set_yticks(range(len(wells)))
    ax.set_yticklabels(wells, fontsize=8, fontfamily='monospace')
    ax.set_xticks(range(len(curves)))
    if idx == n_formations - 1:
        ax.set_xticklabels(curves, rotation=45, ha='right',
                           fontsize=9, fontfamily='monospace')
    else:
        ax.set_xticklabels([])
    for i in range(len(wells)):
        for j in range(len(curves)):
            val = df.values[i, j]
            text_colour = 'white' if val < 35 or val > 80 else 'black'
            ax.text(j, i, f'{val:.0f}%', ha='center', va='center',
                    fontsize=7, fontfamily='monospace',
                    color=text_colour, fontweight='bold')
    ax.text(-0.15, 0.5, formation,
            transform=ax.transAxes,
            fontsize=10, fontweight='bold',
            rotation=0, va='center', ha='center')
    for spine in ax.spines.values():
        spine.set_visible(False)

cbar_ax = fig.add_axes([0.92, 0.15, 0.015, 0.7])
cb = fig.colorbar(
    plt.cm.ScalarMappable(norm=norm, cmap=cmap),
    cax=cbar_ax
)
cb.set_label('Completeness (%)', fontsize=9)
cb.ax.tick_params(labelsize=8)
fig.suptitle('Well Log Data Completeness by Formation',
             fontsize=14, fontweight='bold')
plt.subplots_adjust(hspace=0.08, left=0.18, right=0.9, top=0.97)
plt.savefig('formation_completeness.png', dpi=150,
            bbox_inches='tight', facecolor='white')
plt.show()

The dataset used in this article is entirely synthetic. Well names, formation names, and depth values are fictional and do not correspond to any real field or geological sequence.