Notebook

Wellbore Paths in PyVista¶

I just learned about PyVista a few weeks ago, and I finally had some time to sit down with it and see some of the things it can do. So for an example I decided to plot up a few well pads with horizontal wells from the Denver Basin. I spent some time cleaning and organizing the data so it would be easier to deal with and I could focus on learning PyVista and PVGeo. To install PyVista, PVGeo, and panel run the first cell

In [ ]:

#!pip install pyvista
#!pip install panel
#!pip install PVGeo

Now we are going to import pandas for reading in our well data, numpy for making stratigraphic surfaces, pyvista and PVGeo for visualization. Also we are going to set the panel.extension to 'vtk' and set the pyvista theme to `document

In [ ]:

import pandas as pd
import numpy as np
import pyvista as pv
import panel
import PVGeo

panel.extension("vtk")
pv.set_plot_theme("document")
import warnings

warnings.filterwarnings("ignore")

Let's read in the data using pandas and get a quick look at it

In [ ]:

data = pd.read_csv("wellbore_and_gamma.csv", index_col=[0])
data.head()

Let's get the names of the wells from the data and assign them to a variable

In [ ]:

wells = data.name.unique()

Now we want to go through the data and build a dict for each well

In [ ]:

# split out the wells into a dictionary
well_dict = {}
for well in wells:
    well_dict["{0}".format(well)] = data[data.name == well]

Next let's take each well in the well_dict and get the points referenced to one another so they all plot up in real space relative to one another. We are then creating a numpy array with the 'easting' (x-value), the 'northing' (y-value), and the depth which is our true vertical depth. We now have each entry in the points_dict with an x, y, and z value.

In [ ]:

# now turn the wells into points
points_dict = {}
for points in well_dict:
    points_dict["{0}".format(points)] = np.array(
        list(
            zip(
                well_dict[points].easting + well_dict[points].x,
                well_dict[points].northing + well_dict[points].y,
                -well_dict[points].depth,
            )
        )
    )

From the points we want to create PolyData to plot in PyVista. To create this PolyData we use PVGeo.points_to_poly_data() on each entry in the points_dict. Then we use PVGeo.filters.AddCellConnToPoints to create cells connecting the points. At this point we have everything we need to plot up wellbore trajectories. But let's have some fun!

In [ ]:

# now turn points into polydata
lines_dict = {}
for lines in points_dict:
    poly = PVGeo.points_to_poly_data(points_dict[lines])
    lines_dict[
        "{0}".format(lines)
    ] = PVGeo.filters.AddCellConnToPoints().apply(poly)

To color and size the PolyData by a value we use the GR column from the original dataset. We assigne the GR gamma-ray values to each PolyData object meaning we now have DataArrays for each well that we can use to color our wellbore. Finally, we use .tube() to create tubes at each point with radius 10 that scale with the GR gamma-ray values.

In [ ]:

for path in lines_dict:
    lines_dict[path]["GR"] = data[data.name == path].GR.values
    lines_dict[path].tube(radius=10, scalars="GR", inplace=True)

We could plot up the wellbores at this point and have a great plot, but let's make some stratigraphic surfaces. I grabbed some tops from the COGCC website for one of these wells. To make surfaces we need x and y values spanning the area we are interested in. So I use the minimum and maximum from the dataset and then add 1,000 feet in the x and y direction for the surfaces. Next we use meshgrid and create our grid. Then for each surface I create an empty array using np.empty() and .fill it with the depth to each top. For the last step for each surface we use StructuredGrid to create a structured PolyData grid. Here we have the Niobrara, Shannon, Sussex, Parkman, and Pierre formation tops.

In [ ]:

x = np.arange(
    data.x.min() - data.easting.min() - 1000,
    data.x.max() + data.easting.max() + 1000,
    200,
)
y = np.arange(
    data.y.min() - data.northing.min() - 1000,
    data.y.max() + data.northing.max() + 1000,
    200,
)
x, y = np.meshgrid(x, y)
r = np.sqrt(x ** 2 + y ** 2)

nio = np.empty(r.shape)
nio.fill(-5930)
niobraragrid = pv.StructuredGrid(x, y, nio)

shn = np.empty(r.shape)
shn.fill(-4454)
shannongrid = pv.StructuredGrid(x, y, shn)

ssx = np.empty(r.shape)
ssx.fill(-4009)
sussexgrid = pv.StructuredGrid(x, y, ssx)

prk = np.empty(r.shape)
prk.fill(-3332)
parkmangrid = pv.StructuredGrid(x, y, prk)

pie = np.empty(r.shape)
pie.fill(-2267)
pierregrid = pv.StructuredGrid(x, y, pie)

All that work and now we get so see how it turned out! We start by creating our Plotter and then add_mesh for each pice of PolyData that we have created so far. Lastly, we call .show() and we have an interactive set of wells!

In [ ]:

plotter = pv.Plotter()
for well in lines_dict:
    plotter.add_mesh(lines_dict[well])
plotter.add_mesh(niobraragrid, opacity=0.5, color="white")
plotter.add_mesh(shannongrid, opacity=0.5, color="black")
plotter.add_mesh(sussexgrid, opacity=0.5, color="purple")
plotter.add_mesh(parkmangrid, opacity=0.5, color="orange")
plotter.add_mesh(pierregrid, opacity=0.5, color="green")

plotter.show()

This notebook is licensed as CC-BY, use and share to your hearts content.

In [ ]: