Skip to content

Commit

Permalink
Merge pull request #4 from cgre-aachen/dev_pybore
Browse files Browse the repository at this point in the history 
Dev pybore
  • Loading branch information
@AlexanderJuestel
AlexanderJuestel authored Sep 27, 2023
2 parents 259135c + 60f8ae5 commit d61d4ca
Show file tree
Hide file tree
Showing 4 changed files with 804 additions and 391 deletions.
14 changes: 12 additions & 2 deletions README.md
View file
Original file line number Diff line number Diff line change
@@ -1,2 +1,12 @@
# pybore
Processing and visualizing borehole data
# PyBore - Processing and visualizing borehole data

`pybore` is a Python library for processing and visualizing borehole data. It is similar to the `welly`package building
upon several libraries such as `lasio` and `wellpathpy` to process borehole data.

## Functionality
The functionality of `pybore` includes:

* Reading and visualizing well logs using `lasio`
* Reading and visualizing well deviation paths and desurveying the well using `wellpathpy`
* Parsing XML files written in [BoreholeML](https://www.bgr.bund.de/Infogeo/DE/Home/BoreholeML/boreholeml_node.html)

950 changes: 581 additions & 369 deletions notebooks/borehole.ipynb
View file

Large diffs are not rendered by default.

191 changes: 181 additions & 10 deletions pybore/borehole.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@
from shapely.geometry import Point
from pyproj import CRS
import matplotlib.pyplot as plt
from typing import Union


class Borehole:
Expand Down Expand Up @@ -80,11 +81,6 @@ def add_well_logs(self,
# Creating well logs
self.logs = Logs(path=path)

#def get_mesh_along_borehole_path(self,
# log):



# def read_boreholeml(self):
# def write_boreholeml(self):

Expand Down Expand Up @@ -215,9 +211,9 @@ def lines_from_points(points):
poly.lines = cells
return poly

spline = lines_from_points(np.c_[self.easting_rel+x,
self.northing_rel+y,
self.tvd])
spline = lines_from_points(np.c_[self.easting_rel + x,
self.northing_rel + y,
-self.tvd])

tube = spline.tube(radius=radius)

Expand Down Expand Up @@ -285,5 +281,180 @@ def plot_well_logs(self,

return fig, ax

# def get_log_along_well_path(self):
# see GemGIS function
def plot_well_log_along_path(self,
log,
coordinates,
spacing=0.5,
radius_factor=75):

# Importing pyvista
try:
import pyvista as pv
except ModuleNotFoundError:
ModuleNotFoundError('PyVista package not installed')

if not {'Northing', 'Easting', 'True Vertical Depth Below Sea Level'}.issubset(coordinates.columns):
raise ValueError('The coordinates DataFrame must contain a northing, easting and true vertical depth '
'below sea level column')

coordinates = coordinates[['Easting', 'Northing', 'True Vertical Depth Below Sea Level']].to_numpy()

logs = self.df.reset_index()[['MD', log]]

points = resample_between_well_deviation_points(coordinates=coordinates,
spacing=spacing)

# polyline_well_path = polyline_from_points(points=coordinates)

polyline_well_path_resampled = pv.Spline(points)

points_along_spline = get_points_along_spline(spline=polyline_well_path_resampled,
dist=logs['MD'].values)

polyline_along_spline = polyline_from_points(points=points_along_spline)

polyline_along_spline['values'] = logs[log].values

tube_along_spline = polyline_along_spline.tube(scalars='values',
radius_factor=radius_factor)

return tube_along_spline


def resample_between_well_deviation_points(coordinates: np.ndarray,
spacing: Union[float, int]) -> np.ndarray:
"""Resampling between points that define the path of a well
Parameters
__________
coordinates: np.ndarray
Nx3 Numpy array containing the X, Y, and Z coordinates that define the path of a well
Returns
_______
points_resampled: np.ndarray
Resampled points along a well
.. versionadded:: 1.0.x
"""

# Checking that the coordinates are provided as np.ndarray
if not isinstance(coordinates, np.ndarray):
raise TypeError('Coordinates must be provided as NumPy Array')

# Checking that three coordinates are provided for each point
if coordinates.shape[1] != 3:
raise ValueError('Three coordinates X, Y, and Z must be provided for each point')

# Creating list for storing points
list_points = []

# Iterating over points and creating additional points between all other points
for i in range(len(coordinates) - 1):
dist = np.linalg.norm(coordinates[i] - coordinates[i + 1])
num_points = int(dist // spacing)
points = np.linspace(coordinates[i], coordinates[i + 1], num_points + 1)
list_points.append(points)

# Converting lists of points into np.ndarray
points_resampled = np.array([item for sublist in list_points for item in sublist])

return points_resampled


def polyline_from_points(points: np.ndarray):
"""Creating PyVista PolyLine from points
Parameters
__________
points: np.ndarray
Points defining the PolyLine
Return
______
poly: pv.core.pointset.PolyData
.. versionadded:: 1.0.x
"""

# Importing pyvista
try:
import pyvista as pv
except ModuleNotFoundError:
ModuleNotFoundError('PyVista package not installed')

# Checking that the points are of type PolyData Pointset
if not isinstance(points, np.ndarray):
raise TypeError('The points must be provided as NumPy Array')

# Creating PolyData Object
poly = pv.PolyData()

# Assigning points
poly.points = points

# Creating line values
the_cell = np.arange(0, len(points), dtype=np.int_)
the_cell = np.insert(the_cell, 0, len(points))

# Assigning values to PolyData
poly.lines = the_cell

return poly


def get_points_along_spline(spline,
dist: np.ndarray):
"""Returning the closest point on the spline a given a length along a spline.
Parameters
__________
spline: pv.core.pointset.PolyData
Spline with the resampled vertices
dist: np.ndarray
np.ndarray containing the measured depths (MD) of values along the well path
Return
______
spline.points[idx_list]: pv.core.pyvista_ndarray.pyvista_ndarray
PyVista Array containing the selected points
.. versionadded:: 1.0.x
"""

# Importing pyvista
try:
import pyvista as pv
except ModuleNotFoundError:
ModuleNotFoundError('PyVista package not installed')

# Checking that the spline is a PyVista PolyData Pointset
if not isinstance(spline, pv.core.pointset.PolyData):
raise TypeError('The well path/the spline must be provided as PyVista PolyData Pointset')

# Checking that the distances are provided as np.ndarray
if not isinstance(dist, np.ndarray):
raise TypeError('The distances must be provided as np.ndarray')

# Creating list for storing indices
idx_list = []

# Getting index of spline that match with a measured value and append index to list of indices
for distance in dist:
idx = np.argmin(np.abs(spline.point_data['arc_length'] - distance))
idx_list.append(idx)

points = spline.points[idx_list]

return points
40 changes: 30 additions & 10 deletions pybore/boreholeml.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -414,6 +414,15 @@ def __init__(self,
for layer in child1:
if layer.tag == '{http://www.infogeo.de/boreholeml/3.0}layer':
for interval in layer:
lithology_rockname = []
lithology_percentage = []
lithology_rockColor = []
sublayer_rockcode = []
sublayer_rocknametext = []
sublayer_rockname = []
sublayer_rockColor = []
sublayer_from = []
sublayer_to = []
if interval.tag == '{http://www.infogeo.de/boreholeml/3.0}Interval':
for item in interval:
if item.tag == '{http://www.infogeo.de/boreholeml/3.0}from':
Expand All @@ -432,11 +441,11 @@ def __init__(self,
if subitem.tag == '{http://www.infogeo.de/boreholeml/3.0}Lithology':
for subsubitem in subitem:
if subsubitem.tag == '{http://www.infogeo.de/boreholeml/3.0}rockName':
rockname.append(subsubitem.text)
lithology_rockname.append(subsubitem.text)
if subsubitem.tag == '{http://www.infogeo.de/boreholeml/3.0}percentage':
percentage.append(subsubitem.text)
lithology_percentage.append(subsubitem.text)
if subsubitem.tag == '{http://www.infogeo.de/boreholeml/3.0}rockColor':
rockColor.append(subsubitem.text)
lithology_rockColor.append(subsubitem.text)

if item.tag == '{http://www.infogeo.de/boreholeml/3.0}stratigraphy':
for subitem in item:
Expand All @@ -448,33 +457,44 @@ def __init__(self,
for child in subsubitem:
if child.tag == '{http://www.isotc211.org/2005/gmd}LocalisedCharacterString':
lithostratigraphy.append(child.text)
else:
lithostratigraphy.append('')

if item.tag == '{http://www.infogeo.de/boreholeml/3.0}sublayer':
for child1 in item:
if child1.tag == '{http://www.infogeo.de/boreholeml/3.0}Component':
for child2 in child1:
if child2.tag == '{http://www.infogeo.de/boreholeml/3.0}rockCode':
rockcode_sublayer.append(child2.text)
sublayer_rockcode.append(child2.text)
if child2.tag == '{http://www.infogeo.de/boreholeml/3.0}rockNameText':
for child3 in child2:
if child3.tag == '{http://www.isotc211.org/2005/gmd}LocalisedCharacterString':
rocknametext_sublayer.append(
sublayer_rocknametext.append(
child3.text)

if child2.tag == '{http://www.infogeo.de/boreholeml/3.0}lithology':
for child3 in child2:
if child3.tag == '{http://www.infogeo.de/boreholeml/3.0}Lithology':
for subitem in child3:
if subitem.tag == '{http://www.infogeo.de/boreholeml/3.0}rockName':
rockname_sublayer.append(
sublayer_rockname.append(
subitem.text)
if subitem.tag == '{http://www.infogeo.de/boreholeml/3.0}rockColor':
rockColor_sublayer.append(
sublayer_rockColor.append(
subitem.text)
if child2.tag == '{http://www.infogeo.de/boreholeml/3.0}from':
from_sublayer.append(child2.text)
sublayer_from.append(child2.text)
if child2.tag == '{http://www.infogeo.de/boreholeml/3.0}to':
to_sublayer.append(child2.text)

sublayer_to.append(child2.text)
rockColor_sublayer.append(sublayer_rockColor)
from_sublayer.append(sublayer_from)
to_sublayer.append(sublayer_to)
rockname_sublayer.append(sublayer_rockname)
rocknametext_sublayer.append(sublayer_rocknametext)
rockcode_sublayer.append(sublayer_rockcode)
rockname.append(lithology_rockname)
percentage.append(lithology_percentage)
rockColor.append(lithology_rockColor)
if level4.tag == '{http://www.infogeo.de/boreholeml/3.0}drillingProcess':
for child1 in level4:
if child1.tag == '{http://www.infogeo.de/boreholeml/3.0}DrillingProcess':
Expand Down

0 comments on commit d61d4ca

Please sign in to comment.