From 077e11cfed9acb45757c4ed502f1966abd1afa38 Mon Sep 17 00:00:00 2001 From: jonnymaserati Date: Sat, 6 Jan 2024 15:25:17 +0100 Subject: [PATCH 1/4] Update to utilize to make callbacks easier. Bump version. Update unit tests to automate in vscode. Update volve example to use new . --- examples/volve_wells.py | 8 +- tests/test_clearance_iscwsa.py | 131 +++--- tests/test_connector.py | 266 ++++++----- ...ror_mwdrev5_iscwsa_validation_results.xlsx | Bin 21243 -> 21244 bytes tests/test_fluid.py | 46 +- tests/test_iscwsa_mwd_error.py | 125 +++--- tests/test_survey_interpolate.py | 89 ++-- tests/test_survey_parameters.py | 61 +-- tests/test_utils.py | 148 +++---- welleng/version.py | 2 +- welleng/visual.py | 419 ++++++++++-------- 11 files changed, 667 insertions(+), 628 deletions(-) diff --git a/examples/volve_wells.py b/examples/volve_wells.py index f2483ae..9ee844f 100644 --- a/examples/volve_wells.py +++ b/examples/volve_wells.py @@ -92,8 +92,12 @@ # create a trimesh scene and plot with welleng plotter print("Making a scene and plotting...") -scene = we.mesh.make_trimesh_scene(data) -we.visual.plot(scene) +plt = we.visual.Plotter() +for well in data.values(): + plt.add(well) +plt.show() +# scene = we.mesh.make_trimesh_scene(data) +# we.visual.plot(scene) ########################################################################## # if you wanted to export a transformed scene so that you can, for example diff --git a/tests/test_clearance_iscwsa.py b/tests/test_clearance_iscwsa.py index 069f3f9..94b7035 100644 --- a/tests/test_clearance_iscwsa.py +++ b/tests/test_clearance_iscwsa.py @@ -1,3 +1,4 @@ +import unittest from welleng.survey import Survey, make_survey_header from welleng.clearance import IscwsaClearance import numpy as np @@ -16,7 +17,7 @@ data = json.load(open(filename)) -def generate_surveys(data): +def generate_surveys(self, data=data): # Generate surveys for imported wells surveys = {} @@ -56,74 +57,74 @@ def generate_surveys(data): return surveys -def test_minimize_sf(data=data): - surveys = generate_surveys(data) - reference = surveys["Reference well"] - offset = surveys["09 - well"] - - result = IscwsaClearance(reference, offset, minimize_sf=False) - result_min = IscwsaClearance(reference, offset, minimize_sf=True) - - idx = np.where(result_min.ref.interpolated == False) - - # Check that interpolated survey is not corrupted - for attr in [ - 'azi_grid_rad', 'azi_mag_rad', 'azi_true_rad', 'cov_hla', 'cov_nev', - 'pos_nev', 'pos_xyz', 'md', 'radius' - ]: - assert np.allclose( - getattr(result.ref, attr), getattr(result_min.ref, attr)[idx] - ) - - pass - - for attr in [ - 'Rr', 'calc_hole', 'distance_cc', 'eou_boundary', - 'eou_separation', 'hoz_bearing', 'idx', 'masd', 'off_cov_hla', - 'off_cov_nev', 'off_delta_hlas', 'off_delta_nevs', 'off_pcr', - 'ref_cov_hla', 'ref_cov_nev', 'ref_delta_hlas', 'ref_delta_nevs', - 'ref_nevs', 'ref_pcr', 'sf', 'wellbore_separation' - ]: - # `toolface_bearing` and `trav_cyl_azi_deg` are a bit unstable when - # well paths are parallel. - - assert np.allclose( - getattr(result, attr), getattr(result_min, attr)[idx], - rtol=1e-01, atol=1e-02 - ) - - pass - - -def test_clearance_iscwsa(data=data, rtol=1e-02, atol=1e-03): - surveys = generate_surveys(data) - reference = surveys["Reference well"] - - # Perform clearance checks for each survey - for well in surveys: - if well != "09 - well": - continue - if well == "Reference well": - continue - else: - offset = surveys[well] - # skip well 10 - if well in ["10 - well"]: +class TestClearanceIscwsa(unittest.TestCase): + def test_minimize_sf(self, data=data): + surveys = generate_surveys(data) + reference = surveys["Reference well"] + offset = surveys["09 - well"] + + result = IscwsaClearance(reference, offset, minimize_sf=False) + result_min = IscwsaClearance(reference, offset, minimize_sf=True) + + idx = np.where(result_min.ref.interpolated == False) + + # Check that interpolated survey is not corrupted + for attr in [ + 'azi_grid_rad', 'azi_mag_rad', 'azi_true_rad', 'cov_hla', 'cov_nev', + 'pos_nev', 'pos_xyz', 'md', 'radius' + ]: + assert np.allclose( + getattr(result.ref, attr), getattr(result_min.ref, attr)[idx] + ) + + pass + + for attr in [ + 'Rr', 'calc_hole', 'distance_cc', 'eou_boundary', + 'eou_separation', 'hoz_bearing', 'idx', 'masd', 'off_cov_hla', + 'off_cov_nev', 'off_delta_hlas', 'off_delta_nevs', 'off_pcr', + 'ref_cov_hla', 'ref_cov_nev', 'ref_delta_hlas', 'ref_delta_nevs', + 'ref_nevs', 'ref_pcr', 'sf', 'wellbore_separation' + ]: + # `toolface_bearing` and `trav_cyl_azi_deg` are a bit unstable when + # well paths are parallel. + + assert np.allclose( + getattr(result, attr), getattr(result_min, attr)[idx], + rtol=1e-01, atol=1e-02 + ) + + pass + + def test_clearance_iscwsa(self, data=data, rtol=1e-02, atol=1e-03): + surveys = generate_surveys(data) + reference = surveys["Reference well"] + + # Perform clearance checks for each survey + for well in surveys: + if well != "09 - well": + continue + if well == "Reference well": continue else: - for b in [False, True]: - result = IscwsaClearance(reference, offset, minimize_sf=b) + offset = surveys[well] + # skip well 10 + if well in ["10 - well"]: + continue + else: + for b in [False, True]: + result = IscwsaClearance(reference, offset, minimize_sf=b) + assert np.allclose( + result.sf[np.where(result.ref.interpolated == False)], + np.array(data["wells"][well]["SF"]), + rtol=rtol, atol=atol + ) - assert np.allclose( - result.sf[np.where(result.ref.interpolated == False)], - np.array(data["wells"][well]["SF"]), - rtol=rtol, atol=atol - ) - - pass + pass # make above test runnanble separately if __name__ == '__main__': - test_minimize_sf(data=data) - test_clearance_iscwsa(data=data) + unittest.main() + # test_minimize_sf(data=data) + # test_clearance_iscwsa(data=data) diff --git a/tests/test_connector.py b/tests/test_connector.py index 2418249..b451689 100644 --- a/tests/test_connector.py +++ b/tests/test_connector.py @@ -1,145 +1,139 @@ import inspect import sys +import unittest + +import numpy as np + from welleng.connector import Connector from welleng.survey import Survey, from_connections -import numpy as np -def test_md_hold(): - # test hold with only md provided - c = Connector( - vec1=[0, 0, 1], - md2=500, - ) - assert ( - c.inc_target == c.inc1 - and c.azi_target == c.azi1 - and c.pos_target[2] == c.md_target - ), "Failed c1" - assert c.method == 'hold', "Unexpected method" - - assert isinstance(from_connections(c), Survey) - - -def test_md_and_vec(): - # test with md2 and vec2 provided (minimum curvature) - c = Connector( - vec1=[0, 0, 1], - md2=1000, - vec2=[0, 1, 0] - ) - assert c.method == 'min_curve' - - -def test_pos(): - # test with pos2 provided (minimum distance) - c = Connector( - vec1=[0, 0, 1], - pos2=[100, 100, 1000], - ) - assert c.md_target > c.pos1[2], "Failed c3" - - -def test_pos_and_dls(): - # test with pos2 needing more aggressive dls (minimum curvature) - c = Connector( - vec1=[0, 0, 1], - pos2=[200, 400, 200] - ) - assert c.method == 'min_curve_to_target' - - -def test_pos_and_vec(): - # test with pos2 and vec2 provided - vec1 = [-1, -1, 1] - vec2 = [1, -1, 0] - c = Connector( - pos1=[0., 0., 0], - vec1=vec1 / np.linalg.norm(vec1), - pos2=[0., 1000., 500.], - vec2=vec2 / np.linalg.norm(vec2), - ) - assert c.method == 'curve_hold_curve' - - # test if interpolator and survey functions are working - assert isinstance(from_connections(c, step=30), Survey) - - -def test_pos_inc_azi(): - # test with pos2, inc1 and azi1 provided - c = Connector( - pos1=[0., 0., 0], - inc1=0., - azi1=90, - pos2=[1000., 1000., 1000.], - vec2=[0., 0., 1.], - ) - assert c.method == 'curve_hold_curve' - - -def test_dls2(): - # test with different dls for second curve section - c = Connector( - pos1=[0., 0., 0], - vec1=[0., 0., 1.], - pos2=[0., 100., 1000.], - vec2=[0., 0., 1.], - dls_design2=5 - ) - assert c.radius_design2 < c.radius_design - - -def test_radius_critical(): - # test with dls_critical requirement (actual dls < dls_design) - c = Connector( - pos1=[0., 0., 0], - vec1=[0., 0., 1.], - pos2=[0., 100., 100.], - vec2=[0., 0., 1.], - ) - assert c.radius_critical < c.radius_design - - -def test_min_curve(): - # test min_curve (inc2 provided) - c = Connector( - pos1=[0., 0., 0], - vec1=[0., 0., 1.], - inc2=30, - ) - assert c.method == 'min_curve' - - -def test_radius_critical_with_min_curve(): - # test min_curve with md less than required radius - c = Connector( - pos1=[0., 0., 0], - inc1=0, - azi1=0, - md2=500, - inc2=90, - azi2=0, - ) - assert c.radius_critical < c.radius_design - - -def one_function_to_run_them_all(): - """ - Function to gather the test functions so that they can be tested by - running this module. - - https://stackoverflow.com/questions/18907712/python-get-list-of-all- - functions-in-current-module-inspecting-current-module - """ - test_functions = [ - obj for name, obj in inspect.getmembers(sys.modules[__name__]) - if (inspect.isfunction(obj) - and name.startswith('test') - and name != 'all') - ] - - [f() for f in test_functions] +class ConnectorTest(unittest.TestCase): + def test_md_hold(self): + # test hold with only md provided + c = Connector( + vec1=[0, 0, 1], + md2=500, + ) + assert ( + c.inc_target == c.inc1 + and c.azi_target == c.azi1 + and c.pos_target[2] == c.md_target + ), "Failed c1" + assert c.method == 'hold', "Unexpected method" + + assert isinstance(from_connections(c), Survey) + + def test_md_and_vec(self): + # test with md2 and vec2 provided (minimum curvature) + c = Connector( + vec1=[0, 0, 1], + md2=1000, + vec2=[0, 1, 0] + ) + assert c.method == 'min_curve' + + def test_pos(self): + # test with pos2 provided (minimum distance) + c = Connector( + vec1=[0, 0, 1], + pos2=[100, 100, 1000], + ) + assert c.md_target > c.pos1[2], "Failed c3" + + def test_pos_and_dls(self): + # test with pos2 needing more aggressive dls (minimum curvature) + c = Connector( + vec1=[0, 0, 1], + pos2=[200, 400, 200] + ) + assert c.method == 'min_curve_to_target' + + def test_pos_and_vec(self): + # test with pos2 and vec2 provided + vec1 = [-1, -1, 1] + vec2 = [1, -1, 0] + c = Connector( + pos1=[0., 0., 0], + vec1=vec1 / np.linalg.norm(vec1), + pos2=[0., 1000., 500.], + vec2=vec2 / np.linalg.norm(vec2), + ) + assert c.method == 'curve_hold_curve' + + # test if interpolator and survey functions are working + assert isinstance(from_connections(c, step=30), Survey) + + def test_pos_inc_azi(self): + # test with pos2, inc1 and azi1 provided + c = Connector( + pos1=[0., 0., 0], + inc1=0., + azi1=90, + pos2=[1000., 1000., 1000.], + vec2=[0., 0., 1.], + ) + assert c.method == 'curve_hold_curve' + + def test_dls2(self): + # test with different dls for second curve section + c = Connector( + pos1=[0., 0., 0], + vec1=[0., 0., 1.], + pos2=[0., 100., 1000.], + vec2=[0., 0., 1.], + dls_design2=5 + ) + assert c.radius_design2 < c.radius_design + + def test_radius_critical(self): + # test with dls_critical requirement (actual dls < dls_design) + c = Connector( + pos1=[0., 0., 0], + vec1=[0., 0., 1.], + pos2=[0., 100., 100.], + vec2=[0., 0., 1.], + ) + assert c.radius_critical < c.radius_design + + def test_min_curve(self): + # test min_curve (inc2 provided) + c = Connector( + pos1=[0., 0., 0], + vec1=[0., 0., 1.], + inc2=30, + ) + assert c.method == 'min_curve' + + def test_radius_critical_with_min_curve(self): + # test min_curve with md less than required radius + c = Connector( + pos1=[0., 0., 0], + inc1=0, + azi1=0, + md2=500, + inc2=90, + azi2=0, + ) + assert c.radius_critical < c.radius_design +# def one_function_to_run_them_all(): +# """ +# Function to gather the test functions so that they can be tested by +# running this module. + +# https://stackoverflow.com/questions/18907712/python-get-list-of-all- +# functions-in-current-module-inspecting-current-module +# """ +# test_functions = [ +# obj for name, obj in inspect.getmembers(sys.modules[__name__]) +# if (inspect.isfunction(obj) +# and name.startswith('test') +# and name != 'all') +# ] + +# [f() for f in test_functions] if __name__ == '__main__': - one_function_to_run_them_all() + unittest.main() + # one_function_to_run_them_all() diff --git a/tests/test_data/error_mwdrev5_iscwsa_validation_results.xlsx b/tests/test_data/error_mwdrev5_iscwsa_validation_results.xlsx index 846d9ce8659a869a10b5db0a48814006ea2b4263..4e134980b19c0ad526d4de01ba6445c9799023af 100644 GIT binary patch delta 607 zcmeypl=06}M!o=VW)=|!1_lm>Yn5sd6ZuvPf@r}wg`n3!(bJ3!41$vd0!8YVPCV#$ z*np=kKK0e!(#un27wL%$s85(MgXNK!BTLP1=ynE|;`dZ_Iwkh6!>i@mb;H~Pc zU9kK@HCK-2I=Aa(fs5{^xu@wfD#@LE8KJVs`{6!*zv5kO%oV-we>nt9eaa@Y(QKBf z&!Ic|$KI?l*`e{*JtMu*altg1QW-}Pruua{!37`Bw^>G>Fuc_&)wWRLn(=fqp%1@) ze6Z1)4#sG`_3*n zH1cYt=BJ{~!drco@2i|0WB+Tt({<_H+y5KHv8-y0|G&(R5j9L^#>`O`*nFGOf{g_f zG@Hda0v#c&9M5x1VAfB!vh)|2)Lx@nOzZX-uJai}j3NVug@xbP@ zF)++gXJ8NphB*v0FrJ+3AE*RY)7ZyoBJaq+Ah(>6K?0@#NH;WcO`Z^_4zY1}pp;F# zCj&!8j(%}TWln0bUPW$BfHxzP2s0vvx4XaYIsnuu#>BwDjcgL5;bg%eDTqn>K~juT rlf8oEAlw`vw|sJYkS5qIv)i0nq=D(;Iy1U^V(nVO@D3kzSQDdsKh0tFjKf{?kb-WOPenjvrnC`XjyP}x1CbKv1<<|Zn?T7 zQD|DCp_uqm({7`hBV}R{%r_Pt?-RNb%3E(f&EVIU^uu#>p3J!8#Cte^H|g0YH`6`u zww?d{Y3XaBm*)}|N@XwUT(Po0YnR5H`O`VB36`(VUG1 z6f~P9I079ZtX$7?Okmbzd!MHeAr0S1Muf~{|+#d2k{^l zFfhzfXJ8NphB*v0FrJtk5U2!J)7ZyoBJaq+Ah(>6K?0@#NH;WcPM#R34zY1hpj1?x zCj&!8j(%}TWln0bUPW$BfHxzP2s0vvx4XaYIsnuu#>BwDjcgL5fe!;id{JsnvA!ON zMAtakAV`W)a NDArray: ).reshape((-1, 2, 4)) -def test_annular_volume(): - av = annular_volume( - od=ureg('12.25 inch').to('meter'), - id=ureg(f'{9+5/8} inch').to('meter'), - length=ureg('1000 meter') - ) - - assert av.m == 3.491531223156194 - assert str(av.u) == 'meter ** 3' - - pass - - -def test_decimal2dms(): - degrees, minutes, seconds, direction = decimal2dms( - (LAT[0] + LAT[1] / 60 + LAT[2] / 3600, LAT[3]) - ) - assert (degrees, minutes, round(seconds, 4), direction) == LAT - - dms = decimal2dms(np.array([ - (LAT[0] + LAT[1] / 60 + LAT[2] / 3600, LAT[3]), - (LON[0] + LON[1] / 60 + LON[2] / 3600, LON[3]) - ]), ndigits=4) - assert np.all(np.equal( - dms, - np.array((np.array(LAT, dtype=object), np.array(LON, dtype=object))) - )) - - -def test_dms2decimal(): - decimal = dms2decimal((-LAT[0], LAT[1], LAT[2], LAT[3])) # check it handles westerly - assert np.all(np.equal( - decimal, - np.array([ - -(LAT[0] + LAT[1] / 60 + LAT[2] / 3600), - LAT[3] - ], dtype=object) - )) - - decimals = dms2decimal((LAT, LON)) - assert np.all(np.equal( - decimals, - np.array((dms2decimal(LAT), dms2decimal(LON))) - )) - - decimals = dms2decimal(((LAT, LON), (LON, LAT))) - assert np.all(np.equal( - decimals, - np.array(( - (dms2decimal(LAT), dms2decimal(LON)), - (dms2decimal(LON), dms2decimal(LAT)) +class UtilsTest(unittest.TestCase): + def test_annular_volume(self): + av = annular_volume( + od=ureg('12.25 inch').to('meter'), + id=ureg(f'{9+5/8} inch').to('meter'), + length=ureg('1000 meter') + ) + + assert av.m == 3.491531223156194 + assert str(av.u) == 'meter ** 3' + + pass + + def test_decimal2dms(self): + degrees, minutes, seconds, direction = decimal2dms( + (LAT[0] + LAT[1] / 60 + LAT[2] / 3600, LAT[3]) + ) + assert (degrees, minutes, round(seconds, 4), direction) == LAT + + dms = decimal2dms(np.array([ + (LAT[0] + LAT[1] / 60 + LAT[2] / 3600, LAT[3]), + (LON[0] + LON[1] / 60 + LON[2] / 3600, LON[3]) + ]), ndigits=4) + assert np.all(np.equal( + dms, + np.array((np.array(LAT, dtype=object), np.array(LON, dtype=object))) + )) + + def test_dms2decimal(self): + decimal = dms2decimal((-LAT[0], LAT[1], LAT[2], LAT[3])) # check it handles westerly + assert np.all(np.equal( + decimal, + np.array([ + -(LAT[0] + LAT[1] / 60 + LAT[2] / 3600), + LAT[3] + ], dtype=object) )) - )) + decimals = dms2decimal((LAT, LON)) + assert np.all(np.equal( + decimals, + np.array((dms2decimal(LAT), dms2decimal(LON))) + )) -def test_dms2decimal2dms(): - _dms = _generate_random_dms(int(1e3), 8) - decimal = dms2decimal(_dms) - dms = decimal2dms(decimal, 8) + decimals = dms2decimal(((LAT, LON), (LON, LAT))) + assert np.all(np.equal( + decimals, + np.array(( + (dms2decimal(LAT), dms2decimal(LON)), + (dms2decimal(LON), dms2decimal(LAT)) + )) + )) - assert np.all(np.equal(_dms, dms)) + def test_dms2decimal2dms(self): + _dms = _generate_random_dms(int(1e3), 8) + decimal = dms2decimal(_dms) + dms = decimal2dms(decimal, 8) + assert np.all(np.equal(_dms, dms)) -def one_function_to_run_them_all(): - """ - Function to gather the test functions so that they can be tested by - running this module. - https://stackoverflow.com/questions/18907712/python-get-list-of-all- - functions-in-current-module-inspecting-current-module - """ - test_functions = [ - obj for name, obj in inspect.getmembers(sys.modules[__name__]) - if (inspect.isfunction(obj) - and name.startswith('test') - and name != 'all') - ] +# def one_function_to_run_them_all(): +# """ +# Function to gather the test functions so that they can be tested by +# running this module. - for f in test_functions: - f() +# https://stackoverflow.com/questions/18907712/python-get-list-of-all- +# functions-in-current-module-inspecting-current-module +# """ +# test_functions = [ +# obj for name, obj in inspect.getmembers(sys.modules[__name__]) +# if (inspect.isfunction(obj) +# and name.startswith('test') +# and name != 'all') +# ] - pass +# for f in test_functions: +# f() + +# pass if __name__ == '__main__': - one_function_to_run_them_all() + unittest.main() + # one_function_to_run_them_all() diff --git a/welleng/version.py b/welleng/version.py index a62e53d..32a90a3 100644 --- a/welleng/version.py +++ b/welleng/version.py @@ -1 +1 @@ -__version__ = '0.7.8' +__version__ = '0.8.0' diff --git a/welleng/visual.py b/welleng/visual.py index 182bde4..8b2aa21 100644 --- a/welleng/visual.py +++ b/welleng/visual.py @@ -4,11 +4,17 @@ except ImportError: TRIMESH = False try: - from vedo import trimesh2vedo, Lines, Sphere + import vedo + from vedo import Lines, Mesh + from vedo import Plotter as VedoPlotter + from vedo import Point, Text2D, mag, plotter_instance + from vedo.addons import Icon, compute_visible_bounds + from vedo.utils import buildPolyData VEDO = True except ImportError: VEDO = False import numpy as np +from scipy.spatial import KDTree try: import plotly.graph_objects as go @@ -17,220 +23,226 @@ except ImportError: PLOTLY = False -from vtk import ( - vtkCubeAxesActor, vtkNamedColors, vtkInteractorStyleTerrain -) -from vtkmodules.vtkRenderingCore import ( - vtkRenderWindow, vtkRenderWindowInteractor, vtkRenderer -) -from vtkmodules.vtkInteractionWidgets import ( - vtkCameraOrientationWidget, - vtkOrientationMarkerWidget -) +from vtk import vtkAxesActor, vtkCubeAxesActor, vtkNamedColors +from . import mesh from .version import __version__ as VERSION # VEDO = False -class Plotter(vtkRenderer): - def __init__(self, data=None, **kwargs): - super().__init__() +class Plotter(VedoPlotter): + def __init__(self, *args, **kwargs): """ - A vtk wrapper for quickly visualizing well trajectories for QAQC - purposes. Initiates a vtkRenderer instance and populates it with - mesh data and a suitable axes for viewing well trajectory data. - - Parameters - ---------- - data: obj or list(obj) - A trimesh.Trimesh object or a list of trimesh.Trimesh - objects or a trimesh.scene object - names: list of strings (default: None) - A list of names, index aligned to the list of well meshes. - colors: list of strings (default: None) - A list of color or colors. If a single color is listed then this is - applied to all meshes in data, otherwise the list of colors is - indexed to the list of meshes. + Notes + ----- + On account of Z or TVD pointing down in the drilling world, the + coordinate system is right handed. In order to map coordinates in the + NEV (or North, East, Vertical) reference correctly, North coordinates + are plotted on the X axis and East on the Y axis. Be mindful of this + adding objects to a scene. """ - assert all((VEDO, TRIMESH)), \ - "ImportError: try pip install welleng[easy]" - - if data is not None: - names = kwargs.get('names') - - if isinstance(data, trimesh.scene.scene.Scene): - meshes = [v for k, v in data.geometry.items()] - if names is None: - names = list(data.geometry.keys()) - - # handle a single mesh being passed - elif isinstance(data, trimesh.Trimesh): - meshes = [data] - - else: - meshes = data - if names is not None: - assert len(names) == len(data), \ - "Names must be length of meshes list else None" - - colors = kwargs.get('colors') - if colors is not None: - if len(colors) == 1: - colors = colors * len(meshes) - else: - assert len(colors) == len(meshes), \ - "Colors must be length of meshes list, 1 else None" - - points = kwargs.get('points') - if points is not None: - points = [ - Sphere(p, r=30, c='grey') - for p in points - ] - - meshes_vedo = [] - for i, mesh in enumerate(meshes): - if i == 0: - vertices = np.array(mesh.vertices) - start_locations = np.array([mesh.vertices[0]]) - else: - vertices = np.concatenate( - (vertices, np.array(mesh.vertices)), - axis=0 - ) - start_locations = np.concatenate( - (start_locations, np.array([mesh.vertices[0]])), - axis=0 - ) - - # convert to vedo mesh - m_vedo = trimesh2vedo(mesh) - if colors is not None: - m_vedo.c(colors[i]) - if names is not None: - m_vedo.name = names[i] - # m_vedo.flag() - meshes_vedo.append(m_vedo) - - for mesh in meshes_vedo: - if mesh is None: - continue - self.AddActor(mesh) - - for obj in kwargs.values(): - if isinstance(obj, list): - for item in obj: - try: - self.AddActor(item) - except TypeError: - pass - else: - try: - self.AddActor(obj) - except TypeError: - pass - - self.namedColors = vtkNamedColors() - - self.colors = {} - self.colors['background'] = kwargs.get('background', 'LightGrey') - self.colors['background2'] = kwargs.get('background', 'Lavender') + super().__init__(*args, **kwargs) - def add_axes(self, **kwargs): - axes = CubeAxes(self, **kwargs) - self.AddActor(axes) + self.wells = [] - def get_center(self): - min, max = np.array(self.ComputeVisiblePropBounds()).reshape(-1, 3) - center = min + (max - min) / 2 - - return tuple(center) + pass - def show(self, add_axes=True, orientation_marker=False, **kwargs): + def add(self, obj, *args, **kwargs) -> None: + """Modified method to support direct plotting of + ``welleng.mesh.WellMesh`` instances and for processing the callback + to print well data when the pointer is hovered of a well trajectory. + + If the ``obj`` is a ``welleng.mesh.WellMesh`` instance, then the args + and kwargs will be passed to the `vedo.Mesh` instance to facilate e.g. + color options etc. + + Notes + ----- + ``welleng.mesh.WellMesh`` creates ``trimesh.Mesh`` instances, a legacy + of using the ``trimesh`` library for detecting mesh collisions when + developing automated well trajectory planning. Therefore, to visualize + the meshes with ``vedo`` and ``vtk``, the meshes need to be converted. + + Meshes in ``welleng`` typically reference an 'NEV' coordinate system, + which is [North, East, Vertical]. To map correctly to ``vtk``, North + needs to be mapped to X and East to Y on account of Z pointing down. """ - Convenient method for opening a window to view the rendered scene. - """ - if add_axes: - self.add_axes(**kwargs) - - self.GetActiveCamera().Azimuth(kwargs.get('azimuth', 30)) - self.GetActiveCamera().Elevation(kwargs.get('elevation', 30)) - self.GetActiveCamera().SetViewUp(0, 0, -1) - # self.GetActiveCamera().SetPosition(tuple(pos_new)) - self.GetActiveCamera().SetFocalPoint(self.get_center()) - - self.ResetCamera() - self.SetBackground( - self.namedColors.GetColor3d(self.colors['background']) - ) - self.SetBackground2( - self.namedColors.GetColor3d(self.colors['background2']) - ) - - setSize = kwargs.get('setSize', (1200, 900)) - - renderWindow = vtkRenderWindow() - - renderWindow.AddRenderer(self) - renderWindow.SetSize(*(setSize)) - renderWindow.SetWindowName(f'welleng {VERSION}') - - renderWindowInteractor = vtkRenderWindowInteractor() - renderWindowInteractor.SetRenderWindow(renderWindow) - renderWindowInteractor.SetInteractorStyle(vtkInteractorStyleTerrain()) - - renderWindow.Render() - self.GetActiveCamera().Zoom(0.8) + if isinstance(obj, mesh.WellMesh): + poly = buildPolyData(obj.mesh.vertices, obj.mesh.faces) + vedo_mesh = Mesh(poly, *args, *kwargs) + setattr(obj, 'vedo_mesh', vedo_mesh) + self.wells.append(obj) + super().add(obj.vedo_mesh) + else: + super().add(obj) - interactive = kwargs.get('interactive', True) + pass - if orientation_marker: - cow = vtkCameraOrientationWidget() - cow.SetParentRenderer(self) - # cow.SetDefaultRenderer(self) + def _initiate_axes_actor(self): + plt = vedo.plotter_instance + r = plt.renderers.index(plt.renderer) + + axact = vtkAxesActor() + axact.SetShaftTypeToCylinder() + axact.SetCylinderRadius(0.03) + axact.SetXAxisLabelText("N") + axact.SetYAxisLabelText("E") + axact.SetZAxisLabelText("V") + axact.GetXAxisShaftProperty().SetColor(1, 0, 0) + axact.GetYAxisShaftProperty().SetColor(0, 1, 0) + axact.GetZAxisShaftProperty().SetColor(0, 0, 1) + axact.GetXAxisTipProperty().SetColor(1, 0, 0) + axact.GetYAxisTipProperty().SetColor(0, 1, 0) + axact.GetZAxisTipProperty().SetColor(0, 0, 1) + bc = np.array(plt.renderer.GetBackground()) + if np.sum(bc) < 1.5: + lc = (1, 1, 1) + else: + lc = (0, 0, 0) + axact.GetXAxisCaptionActor2D().GetCaptionTextProperty().BoldOff() + axact.GetYAxisCaptionActor2D().GetCaptionTextProperty().BoldOff() + axact.GetZAxisCaptionActor2D().GetCaptionTextProperty().BoldOff() + axact.GetXAxisCaptionActor2D().GetCaptionTextProperty().ItalicOff() + axact.GetYAxisCaptionActor2D().GetCaptionTextProperty().ItalicOff() + axact.GetZAxisCaptionActor2D().GetCaptionTextProperty().ItalicOff() + axact.GetXAxisCaptionActor2D().GetCaptionTextProperty().ShadowOff() + axact.GetYAxisCaptionActor2D().GetCaptionTextProperty().ShadowOff() + axact.GetZAxisCaptionActor2D().GetCaptionTextProperty().ShadowOff() + axact.GetXAxisCaptionActor2D().GetCaptionTextProperty().SetColor(lc) + axact.GetYAxisCaptionActor2D().GetCaptionTextProperty().SetColor(lc) + axact.GetZAxisCaptionActor2D().GetCaptionTextProperty().SetColor(lc) + axact.PickableOff() + icn = Icon(axact, size=0.1) + plt.axes_instances[r] = icn + icn.SetInteractor(plt.interactor) + icn.EnabledOn() + icn.InteractiveOff() + plt.widgets.append(icn) + + def _pointer_callback(self, event): + i = event.at + pt2d = event.picked2d + objs = self.at(i).objects + pt3d = self.at(i).compute_world_coordinate(pt2d, objs=objs) + if mag(pt3d) < 0.01: + if self.pointer is None: + return + # if self.pointer in self.at(i).actors: + self.at(i).remove(self.pointer) + self.pointer = None + self.pointer_text.text('') + self.render() + return + if self.pointer is None: + self.pointer = Point().color('red').pos(pt3d) + else: + self.pointer.pos(pt3d) + self.at(i).add(self.pointer) - cow.Off() - cow.EnabledOn() + well_data = self._get_closest_well(pt3d, objs) + if well_data is None: + self.pointer_text.text(f'point coordinates: {np.round(pt3d, 3)}') + else: + survey = well_data.get('well').s + idx = well_data.get('idx_survey') + name = survey.header.name + md = survey.md[idx] + inc = survey.inc_deg[idx] + azi_grid = survey.azi_grid_deg[idx] + self.pointer_text.text(f''' + well name: {name}\n + md: {md:.2f}\t inc: {inc:.2f}\t azi: {azi_grid:.2f}\n + point coordinates: {np.round(pt3d, 3)} + ''') + self.render() + + def _well_vedo_meshes(self): + return [well.vedo_mesh for well in self.wells] + + def _get_closest_well(self, pos, objs) -> dict: + wells = [ + well for well in self.wells + if well.vedo_mesh in objs + ] + if not bool(wells): + return + + results = np.zeros((len(wells), 3)) + for i, well in enumerate(wells): + tree = KDTree(well.vertices.reshape(-1, 3)) + distance, idx_vertices = tree.query(pos) + results[i] = np.array([distance, idx_vertices, well.n_verts]) + + winner = np.argmin(results[:, 0]) + distance, idx_vertices, n_verts = results[winner] + + return { + 'well': wells[winner], + 'distance': distance, + 'idx_vertices': int(idx_vertices), + 'n_verts': int(n_verts), + 'idx_survey': int(idx_vertices // n_verts) + } + + def show(self, axes=None, *args, **kwargs): + # check if there's an axes and if so remove them + if self.axes is not None: + self.remove(self.axes) + + self._initiate_axes_actor() + + self.add_callback('mouse move', self._pointer_callback) + self.pointer_text = Text2D("", pos='bottom-right', s=0.5, c='black') + self.add(self.pointer_text) + self.pointer = None + + self = super().show( + viewup=[0, 0, -1], mode=8, + axes=CubeAxes() if axes is None else axes, + title=f'welleng {VERSION}', + *args, **kwargs + ) - if interactive: - renderWindowInteractor.Start() + return self class CubeAxes(vtkCubeAxesActor): - def __init__(self, renderer, **kwargs): + def __init__( + self, + **kwargs + ): super().__init__() - # Determine the bounds from the meshes/actors being plotted rounded - # up to the nearest 100 units. - bounds = np.array(renderer.ComputeVisiblePropBounds()) + # # Determine the bounds from the meshes/actors being plotted rounded + # # up to the nearest 100 units. + # bounds = np.array(renderer.ComputeVisiblePropBounds()) + + # with np.errstate(divide='ignore', invalid='ignore'): + # self.bounds = tuple(np.nan_to_num( + # ( + # np.ceil(np.abs(bounds) / 100) * 100 + # ) * (bounds / np.abs(bounds)) + # )) - with np.errstate(divide='ignore', invalid='ignore'): - self.bounds = tuple(np.nan_to_num( - ( - np.ceil(np.abs(bounds) / 100) * 100 - ) * (bounds / np.abs(bounds)) - )) + plt = vedo.plotter_instance + r = plt.renderers.index(plt.renderer) + vbb = compute_visible_bounds()[0] + self.SetBounds(vbb) + self.SetCamera(plt.renderer.GetActiveCamera()) namedColors = vtkNamedColors() self.colors = {} for n in range(1, 4): self.colors[f'axis{n}Color'] = namedColors.GetColor3d( - kwargs.get('axis1Color', 'DarkGrey') + kwargs.get(f'axis{n}Color', 'Black') ) - self.SetLabelScaling(0, 0, 0, 0) - - # Try and prevent scientific numbering on axes - self.SetXLabelFormat("%.0f") - self.SetYLabelFormat("%.0f") - self.SetZLabelFormat("%.0f") - - self.SetUseTextActor3D(1) + # self.SetUseTextActor3D(1) - self.SetBounds(self.bounds) - self.SetCamera(renderer.GetActiveCamera()) + # self.SetBounds(self.bounds) + # self.SetCamera(renderer.GetActiveCamera()) self.GetTitleTextProperty(0).SetColor(self.colors['axis1Color']) self.GetLabelTextProperty(0).SetColor(self.colors['axis1Color']) @@ -244,26 +256,45 @@ def __init__(self, renderer, **kwargs): self.GetLabelTextProperty(2).SetColor(self.colors['axis3Color']) self.GetLabelTextProperty(2).SetOrientation(45.0) - self.DrawXGridlinesOn() - self.DrawYGridlinesOn() - self.DrawZGridlinesOn() self.SetGridLineLocation(self.VTK_GRID_LINES_FURTHEST) + for a in ('X', 'Y', 'Z'): + getattr(self, f'Get{a}AxesLinesProperty')().SetColor( + namedColors.GetColor3d('Black') + ) + getattr(self, f'SetDraw{a}Gridlines')(1) + getattr(self, f'Get{a}AxesGridlinesProperty')().SetColor( + namedColors.GetColor3d('Grey') + ) + getattr(self, f'{a}AxisMinorTickVisibilityOff')() - self.XAxisMinorTickVisibilityOff() - self.YAxisMinorTickVisibilityOff() - self.ZAxisMinorTickVisibilityOff() + # self.DrawXGridlinesOn() + # self.DrawYGridlinesOn() + # self.DrawZGridlinesOn() + # self.SetGridLineLocation(self.VTK_GRID_LINES_FURTHEST) - units = kwargs.get('units', 'meters') - self.SetXTitle('East') + # self.XAxisMinorTickVisibilityOff() + # self.YAxisMinorTickVisibilityOff() + # self.ZAxisMinorTickVisibilityOff() + + units = kwargs.get('units', None) + self.SetXTitle('N') self.SetXUnits(units) - self.SetYTitle('North') + self.SetYTitle('E') self.SetYUnits(units) self.SetZTitle('TVD') self.SetZUnits(units) + # self.SetTickLocation(self.VTK_GRID_LINES_FURTHEST) + # self.ForceOpaqueOff() self.SetFlyModeToClosestTriad() + # Try and prevent scientific numbering on axes + self.SetLabelScaling(0, 0, 0, 0) + self.SetXLabelFormat("%.0f") + self.SetYLabelFormat("%.0f") + self.SetZLabelFormat("%.0f") - self.ForceOpaqueOff() + plt.axes_instances[r] = self + plt.renderer.AddActor(self) def plot( From a508c70962334e4e8dd7b84cf61d2db3a925ab7d Mon Sep 17 00:00:00 2001 From: jonnymaserati Date: Sat, 6 Jan 2024 16:52:06 +0100 Subject: [PATCH 2/4] Updated dms2decimal and decimal2dms functions in utils to be more flexible and handle conversions without directions. --- tests/test_utils.py | 28 +++++++++++++++++++++++++++- welleng/utils.py | 41 +++++++++++++++++++++++++++++++++-------- 2 files changed, 60 insertions(+), 9 deletions(-) diff --git a/tests/test_utils.py b/tests/test_utils.py index 971d337..61d1e8e 100644 --- a/tests/test_utils.py +++ b/tests/test_utils.py @@ -61,6 +61,23 @@ def test_decimal2dms(self): ) assert (degrees, minutes, round(seconds, 4), direction) == LAT + dms = decimal2dms( + (LAT[0] + LAT[1] / 60 + LAT[2] / 3600), ndigits=4 + ) + assert np.all(np.equal( + dms, + np.array((np.array(LAT[:3]))) + )) + + dms = decimal2dms(np.array([ + (LAT[0] + LAT[1] / 60 + LAT[2] / 3600), + (LON[0] + LON[1] / 60 + LON[2] / 3600) + ]).reshape((-1, 1)), ndigits=4) + assert np.all(np.equal( + dms, + np.array((np.array(LAT[:3]), np.array(LON[:3]))) + )) + dms = decimal2dms(np.array([ (LAT[0] + LAT[1] / 60 + LAT[2] / 3600, LAT[3]), (LON[0] + LON[1] / 60 + LON[2] / 3600, LON[3]) @@ -80,10 +97,19 @@ def test_dms2decimal(self): ], dtype=object) )) + decimal = dms2decimal((LAT[:3])) + assert decimal == LAT[0] + LAT[1] / 60 + LAT[2] / 3600 + + decimals = dms2decimal((LAT[:3], LON[:3])) + assert np.all(np.equal( + decimals, + np.array((dms2decimal(LAT[:3]), dms2decimal(LON[:3]))) + )) + decimals = dms2decimal((LAT, LON)) assert np.all(np.equal( decimals, - np.array((dms2decimal(LAT), dms2decimal(LON))) + np.array((dms2decimal(LAT), dms2decimal(LON))).reshape(decimals.shape) )) decimals = dms2decimal(((LAT, LON), (LON, LAT))) diff --git a/welleng/utils.py b/welleng/utils.py index 20395ea..3800faa 100644 --- a/welleng/utils.py +++ b/welleng/utils.py @@ -486,7 +486,8 @@ def errors_from_cov(cov, data=False): if data: return { i: { - 'nn': _nn, 'ne': _ne, 'nv': _nv, 'ee': _ee, 'ev': _ev, 'vv': _vv + 'nn': _nn, 'ne': _ne, 'nv': _nv, + 'ee': _ee, 'ev': _ev, 'vv': _vv } for i, (_nn, _ne, _nv, _ee, _ev, _vv) in enumerate(zip(nn, ne, nv, ee, ev, vv)) @@ -660,18 +661,25 @@ def annular_volume(od: float, id: float = None, length: float = None): def _decimal2dms(decimal: tuple, ndigits: int = None) -> tuple: - _decimal, direction = decimal + try: + _decimal, direction = decimal + except (TypeError, ValueError): + _decimal = decimal[0] if isinstance(decimal, np.ndarray) else decimal + direction = None _decimal = float(_decimal) minutes, seconds = divmod(abs(_decimal) * 3600, 60) _, minutes = divmod(minutes, 60) return np.array([ + int(_decimal), + int(minutes), + seconds if ndigits is None else round(seconds, ndigits) + ]) if direction is None else np.array([ int(_decimal), int(minutes), seconds if ndigits is None else round(seconds, ndigits), - direction], - dtype=object - ) + direction + ], dtype=object) def decimal2dms(decimal: tuple | NDArray, ndigits: int = None) -> tuple | NDArray: @@ -680,7 +688,7 @@ def decimal2dms(decimal: tuple | NDArray, ndigits: int = None) -> tuple | NDArra Parameters ---------- decimal : tuple | arraylike - A tuple of (lat, direction) or (lon, direction) or arraylike of + A tuple of (lat, direction) or (lon, direction) or arraylike of ((lat, direction), (lon, direction)) coordinates. ndigits: int (default is None) If specified, rounds the seconds decimal to the desired number of @@ -702,7 +710,10 @@ def decimal2dms(decimal: tuple | NDArray, ndigits: int = None) -> tuple | NDArra [[52 4 43.1868 'N'] [4 17 19.6368 'E']] """ - dms = np.apply_along_axis(_decimal2dms, -1, decimal, ndigits) + try: + dms = np.apply_along_axis(_decimal2dms, -1, decimal, ndigits) + except np.exceptions.AxisError: + dms = _decimal2dms(decimal, ndigits) if dms.shape == (4,): return tuple(dms) @@ -711,10 +722,20 @@ def decimal2dms(decimal: tuple | NDArray, ndigits: int = None) -> tuple | NDArra def _dms2decimal(dms: NDArray, ndigits: int = None) -> NDArray: - degrees, minutes, seconds, direction = dms + try: + degrees, minutes, seconds, direction = dms + except ValueError: + degrees, minutes, seconds = dms + direction = None + decimal = abs(degrees) + minutes / 60 + seconds / 3600 return np.array([ + np.copysign( + decimal if ndigits is None else round(decimal, ndigits), + degrees + ) + ]) if direction is None else np.array([ np.copysign( decimal if ndigits is None else round(decimal, ndigits), degrees @@ -756,5 +777,9 @@ def dms2decimal(dms: tuple | NDArray, ndigits: int = None) -> NDArray: if result.shape == (): return float(result) + elif result.shape == (1,): + return float(result[0]) + elif result.shape[-1] == 1: + return result.reshape(-1) else: return result From ad856db55fdd0b2df9be04e518beb92cba74da8a Mon Sep 17 00:00:00 2001 From: jonnymaserati Date: Sat, 6 Jan 2024 17:26:56 +0100 Subject: [PATCH 3/4] Finessed the utils.decimal2dms a bit more. --- ...ror_mwdrev5_iscwsa_validation_results.xlsx | Bin 21244 -> 21244 bytes tests/test_utils.py | 11 +++++++++-- welleng/utils.py | 11 ++++++++--- 3 files changed, 17 insertions(+), 5 deletions(-) diff --git a/tests/test_data/error_mwdrev5_iscwsa_validation_results.xlsx b/tests/test_data/error_mwdrev5_iscwsa_validation_results.xlsx index 4e134980b19c0ad526d4de01ba6445c9799023af..678c7adcb1d4f2afd96ee9c2ff258ca2c34c14f4 100644 GIT binary patch delta 410 zcmeyfl=06}MxFp~W)=|!1_lm>rEL>=)^aRuQ;Rq;;g83}gYxyh{)Y{C+Ww!KQvG^M z>g0}HaSqnJoSmKud4Z2|ug>CEb$k2IHuI*Hz#WTY?;juE7No5|h5LuhoZ}9<$Brgg zNwro$q_9d^O7gSI4Uu=T29qus%N+17oy6WU>9eSW%1TvFiPo9(d*3X&k`?x= z|AH^mS_W&A3{9tr;swc^_0662WYn2mu)^c2{RErRdQwW-PP`>`4-(druw)oUndrL1* zm0hGKE}%YP!VH#2VwU&BE^Y}-y7KO==jm&W589@9|Ed4?MuWGivv$Gq3)NgXn(N%I zmjy1mpXQ#X&!{AK@@0g|BJYR$`2C7^wJ}%pzW?PAF!d>$%to_Wrap)6=pTEt#$<=a zU-yjkM#lxyWJ+ZmMVQv<)CU)QJl|#+dBX5kt5n-UiEGBw&4fPu`tiX=fBtQ~P46Zc zU6q(3;ZZDq<%?3i)NQtXqO(80SYxuGe=o~jbZaWMl&|1Yn7YDID#EP zbhhUOCLlf8*5^5htL7Wc2;w&TJqB?z{C&XOZ2>m$P@McLz#OE+IM5zMl?1wgs9k}U WAnI!%SWYj<8pO*EasyGDgFFG*Rke%& diff --git a/tests/test_utils.py b/tests/test_utils.py index 61d1e8e..24f981e 100644 --- a/tests/test_utils.py +++ b/tests/test_utils.py @@ -1,5 +1,3 @@ -import inspect -import sys import unittest import numpy as np @@ -78,6 +76,15 @@ def test_decimal2dms(self): np.array((np.array(LAT[:3]), np.array(LON[:3]))) )) + dms = decimal2dms(np.array([ + (LAT[0] + LAT[1] / 60 + LAT[2] / 3600), + (LON[0] + LON[1] / 60 + LON[2] / 3600) + ]), ndigits=4) + assert np.all(np.equal( + dms, + np.array((np.array(LAT[:3]), np.array(LON[:3]))) + )) + dms = decimal2dms(np.array([ (LAT[0] + LAT[1] / 60 + LAT[2] / 3600, LAT[3]), (LON[0] + LON[1] / 60 + LON[2] / 3600, LON[3]) diff --git a/welleng/utils.py b/welleng/utils.py index 3800faa..9f966ba 100644 --- a/welleng/utils.py +++ b/welleng/utils.py @@ -710,15 +710,20 @@ def decimal2dms(decimal: tuple | NDArray, ndigits: int = None) -> tuple | NDArra [[52 4 43.1868 'N'] [4 17 19.6368 'E']] """ + flag = False + _decimal = np.array(decimal) + if _decimal.dtype == np.float64: + _decimal = _decimal.reshape((-1, 1)) + flag = True try: - dms = np.apply_along_axis(_decimal2dms, -1, decimal, ndigits) + dms = np.apply_along_axis(_decimal2dms, -1, _decimal, ndigits) except np.exceptions.AxisError: - dms = _decimal2dms(decimal, ndigits) + dms = _decimal2dms(_decimal, ndigits) if dms.shape == (4,): return tuple(dms) else: - return dms + return dms.reshape((-1, 3)) if flag else dms def _dms2decimal(dms: NDArray, ndigits: int = None) -> NDArray: From 98aed5eb04f095f03efcafd7dbef587bec070cd3 Mon Sep 17 00:00:00 2001 From: jonnymaserati Date: Sat, 6 Jan 2024 17:31:21 +0100 Subject: [PATCH 4/4] Add vedo to default requirements. --- setup.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/setup.py b/setup.py index c38b865..f9a71b4 100644 --- a/setup.py +++ b/setup.py @@ -34,6 +34,7 @@ 'pyproj', # required for getting survey parameters 'PyYAML', 'setuptools', + 'vedo', 'vtk' ]) @@ -43,8 +44,7 @@ 'networkx', 'tabulate', 'trimesh', - 'utm', - 'vedo' + 'utm' ]) # this is the troublesome requirement that needs C dependencies