Using Plotly 3D Surface Plots to Visualise Geological Surfaces | by Andy McDonald | Jun, 2023


Visualising the Subsurface using Python Data Visualisation Libraries

3D Surface Plot of the Hugin Formation. Image by the author.

Within geoscience, it is essential to have a full understanding of the geological surfaces present within the subsurface. By knowing the exact location of the formation and its geometry, it becomes easier to identify potential new targets for oil and gas exploration as well as potential carbon capture and storage locations. There are a variety of methods we can use to refine these surfaces ranging from seismic data to well log derived formation tops. Most often, these techniques are used in combination with each other to refine the final surface.

This article continues from my previous articles, which focused on the process of extrapolating well log measurements across a region to understand and visualise geospatial variations. In this article, we will look at how we can create 3D surfaces using interactive Plotly charts.

As modelling geological surfaces is a complex process and often involves multiple iterations and refinement, this article demonstrates a very simple example of how we can visualise this data with Python.

To see how we can use Plotly to visualise our geological formation top across an area we will be using two sets of data.

The first data set is derived from 28 intersections with the formation derived from wellbore picks, which are used as input for kriging to produce a low-resolution surface. In contrast, the second set of data is derived from interpreted seismic data, providing a much higher resolution surface.

Both sets of data are from the Equinor Volve dataset. Details of which are at the bottom of this article.

You can also check out the following articles within this mini-series at the links below:

Before attempting to do any work with the data, we first need to import the libraries we will need. These are:

  • pandas — to read our data, which is in csv format
  • matplotlib to create our visualisation
  • pykrige to carry out the kriging
  • numpy for some numerical calculations
  • plotly.graph_objects to visualise surface in 3D
import pandas as pd
import matplotlib.pyplot as plt
import plotly.graph_objects as go

from pykrige import OrdinaryKriging
import numpy as np

Next, we can load the data using pd.read_csv().

As this data contains information about geological surfaces from all of the wells within the Volve field, we can use query() to extract the data for the formation we need. In this case, we will be looking at the Hugin Formation.

df = pd.read_csv('Data/Volve/Well_picks_Volve_v1 copy.csv')
df_hugin = df.query('SURFACE == "Hugin Fm. VOLVE Top"')
df_hugin

When we run the above code, we get back the following table. You may notice that some wells have encountered the Hugin Formation multiple times. This is likely due to the wells penetrating the formation multiple times either due to wellbore or formation geometry.

Pandas dataframe containing information about formation picks for the Hugin Formation in the Volve field. Image by the author.

In my previous articles, I have focused on how we can use a process known as kriging to “fill in the gaps” between the measurement points. We won’t be covering the details of this process in this article; however, you can check this article for more information.

Once our data has been loaded, we can run the kriging process by calling upon pykrige’s OrdinaryKriging method.

Within this call, we pass in our x and y data which represents the Easting and Northing position of where the wellbore encountered the formation within the subsurface.

As we want to generate a surface of the Hugin Formation, we need to use the TVDSS — True Vertical Depth Subsea — measurement. This gives a true reflection of how deep down the surface is below the selected datum.

OK = OrdinaryKriging(x=df_hugin['Easting'], 
y=df_hugin['Northing'],
z=df_hugin['TVDSS'],
variogram_model='linear',
verbose=True, enable_plotting=True)
Ordinary Kriging results for the Hugin Formation. Image by the author.

Once the model has been generated, we can apply it to two arrays that cover the entire range of the well/penetration points.

This allows us to generate a grid of values which are then passed into the OrdinaryKriging object we generated above.

gridx = np.arange(433986, 438607, 50, dtype='float64')
gridy = np.arange(6477539, 6479393, 50,dtype='float64')

zstar, ss = OK.execute('grid', gridx, gridy)

To finish things off, we can generate a simple 2D map view of our surface using matplotlib’s imshow method.

fig, ax = plt.subplots(figsize=(15,5))

# Create a 2D image plot of the data in 'zstar'
# The 'extent' parameter sets the bounds of the image in data coordinates
# 'origin' parameter sets the part of the image that corresponds to the origin of the axes
image = ax.imshow(zstar, extent=(433986, 438607, 6477539, 6479393),
origin='lower')

# Set the labels for the x-axis and y-axis
ax.set_xlabel('X Location (m)', fontsize=14, fontweight='bold')
ax.set_ylabel('Y Location (m)', fontsize=14, fontweight='bold')

# Add contours
contours = ax.contour(gridx, gridy, zstar, colors='black')

colorbar = fig.colorbar(image)

colorbar.set_label('DTC (us/ft)', fontsize=14, fontweight='bold')

# Display the plot
plt.show()

2D surface and contour map of the Hugin Formation after applying kriging. Image by the author.

To convert our 2D surface to 3D, we need to use Plotly. We could use matplotlib to do this; however, from my experience, it is easier, smoother and more interactive to generate the 3D visualisations with Plotly.

In the code below, we first create our grid of coordinates. To do this, we use numpy’s linspace function. This function will create a set of evenly spaced numbers over a specified range. For our dataset and example, the range extends from the minimum to the maximum of xgrid_extent and ygrid_extent.

The total number of values used within this range will be equivalent to the number of x and y elements present in the zstar grid we created above.

Once our grid is formed, we then call upon the Plotly.

First, we create our figure object and then use fig.add_trace to add our 3D surface plot to the figure.

Once this has been added, we need to tweak the layout of our plot so that we have axis labels, a defined width and height, and some padding.

xgrid_extent = [433986, 438607]
ygrid_extent = [6477539, 6479393]

x = np.linspace(xgrid_extent[0], xgrid_extent[1], zstar.shape[1])
y = np.linspace(ygrid_extent[0], ygrid_extent[1], zstar.shape[0])

fig = go.Figure()

fig.add_trace(go.Surface(z=zstar, x=x, y=y))

fig.update_layout(scene = dict(
xaxis_title='X Location',
yaxis_title='Y Location',
zaxis_title='Depth'),
width=1000,
height=800,
margin=dict(r=20, l=10, b=10, t=10))

fig.show()

When we run the code above, we get back the following interactive plot showing the geological surface of the Hugin formation based on the multiple encounters from the drilled wellbores.

3D Surface plot of the Hugin Formation generated using Plotly. Image by the author.

The Volve dataset has a number of fully interpreted surfaces that have been generated from geological interpretations, including seismic data.

This data contains the x and y locations of data points across the field, as well as our TVDSS data (z).

The data supplied on the Volve data portal is in the form of a .dat file, however, with a bit of manipulation within a text editor, it can easily be converted to a CSV file and loaded using pandas.

hugin_formation_surface = pd.read_csv('Data/Volve/Hugin_Fm_Top+ST10010ZC11_Near_190314_adj2_2760_EasyDC+STAT+DEPTH.csv')
X, Y and Z locations of the Hugin Formation. Image by the author.

Once we have the data loaded, we can make things easier for ourselves by extracting the x, y and z data to variables.

x = hugin_formation_surface['x']
y = hugin_formation_surface['y']
z = hugin_formation_surface['z']

We then need to create a regularly spaced grid between our min and max positions within the x and y data locations. This can be done using numpy’s meshgrid.

xi = np.linspace(x.min(), x.max(), 100)
yi = np.linspace(y.min(), y.max(), 100)
xi, yi = np.meshgrid(xi, yi)

There are several ways to interpolate between a series of data points. The method chosen will depend on the form your data is in (regularly sampled data points vs irregularly sampled), data size and computing power.

If we have a large dataset like the one here, it will be much more computationally expensive with some of the methods such as the Radial Basis Function.

For this example, we will use the LinearNDInterpolator method within scipy to build our model, and then apply it to our z (TVDSS) variable.

In order for us to interpolate between points, we need to reshape xi, yi to 1D arrays for interpolation, as LinearNDInterpolator expects 1D array.

xir = xi.ravel()
yir = yi.ravel()

interp = LinearNDInterpolator((x, y), z)
zi = interp(xir, yir)

Once that has been computed, we can move on to creating our 3D Surface with Plotly Graph Objects.

fig = go.Figure()

fig.add_trace(go.Surface(z=zi, x=xi, y=yi, colorscale='Viridis'))

fig.update_layout(scene = dict(
xaxis_title='Easting (m)',
yaxis_title='Northing (m)',
zaxis_title='Depth',
zaxis=dict(autorange='reversed')),
width=1000,
height=800,
margin=dict(r=20, l=10, b=10, t=10))

fig.update_traces(contours_z=dict(show=True,
usecolormap=True,
project_z=True,
highlightcolor="white"))

fig.show()

When we run the above code, we get back the following 3D surface plot of the Hugin Formation.

Fully interpreted geological surface of the Hugin Formation. Image by the author.

When we compare this plot to the one generated from wellbore measurements, we can definitely see similarities in the overall shape with the valley in the middle. However, the seismic-derived surface provides much greater detail than the well-derived formation tops one.



Source link

Leave a Comment