Claudius

import sys, os

# Importing gempy
import gempy as gp
import gempy_viewer as gpv

# Aux imports
import numpy as np
import pandas as pn

Loading data from repository:

With pandas we can do it directly from the web and with the right args we can directly tidy the data in gempy style:

data_path = os.path.abspath('../../data/input_data/Claudius')

reduce_data_by = 30

dfs = []
for letter in 'ABCD':
    dfs.append(
        pn.read_csv(
            filepath_or_buffer=f"{data_path}/{letter}Points.csv",
            sep=';',
            names=['X', 'Y', 'Z', 'surface', 'cutoff'],
            header=0
        )[::reduce_data_by]
    )

# Add fault:
dfs.append(
    pn.read_csv(
        filepath_or_buffer=f"{data_path}/Fault.csv",
        names=['X', 'Y', 'Z', 'surface'],
        header=0,
        sep=','
    )
)

surface_points = pn.concat(dfs, sort=True)
surface_points['surface'] = surface_points['surface'].astype('str')
surface_points.reset_index(inplace=True, drop=False)
surface_points.tail()
index X Y Z cutoff surface
755 88 551099.250977 7.817652e+06 -10466.863281 NaN Claudius_fault
756 89 551160.807495 7.817503e+06 -10356.463867 NaN Claudius_fault
757 90 551131.898438 7.817659e+06 -10383.323242 NaN Claudius_fault
758 91 551164.412476 7.817654e+06 -10299.957031 NaN Claudius_fault
759 92 551197.192139 7.817647e+06 -10216.820312 NaN Claudius_fault


surface_points.dtypes
index        int64
X          float64
Y          float64
Z          float64
cutoff     float64
surface     object
dtype: object

How many points are per surface

surface_points.groupby('surface').count()
index X Y Z cutoff
surface
0 167 167 167 167 167
250 167 167 167 167 167
330 166 166 166 166 166
60 167 167 167 167 167
Claudius_fault 93 93 93 93 0


Now we do the same with the orientations:

dfs = []

for surf in ['0', '330']:
    o = pn.read_csv(
        filepath_or_buffer=f"{data_path}/Dips.csv",
        sep=';',
        names=['X', 'Y', 'Z', 'G_x', 'G_y', 'G_z', '-'],
        header=1
    )

    # Orientation needs to belong to a surface. This is mainly to categorize to which series belong and to
    # use the same color
    o['surface'] = surf
    dfs.append(o)
orientations = pn.concat(dfs, sort=True)
orientations.reset_index(inplace=True, drop=False)

orientations.tail()
index - G_x G_y G_z X Y Z surface
43 19 0.976692 0.187195 0.138165 -0.972558 550989.3105 7817210.527 -9782.967773 330
44 20 0.248187 -0.080561 -0.043063 -0.995819 550939.3105 7821227.309 -9958.425781 330
45 21 0.649480 -0.161328 0.075208 -0.984031 549276.8105 7820682.980 -9985.125977 330
46 22 0.050025 -0.012103 -0.153309 -0.988104 548976.8105 7820345.121 -9974.265625 330
47 23 0.759421 0.369490 -0.187053 -0.910213 549764.3105 7820457.738 -9901.208984 330


orientations.dtypes
index        int64
-          float64
G_x        float64
G_y        float64
G_z        float64
X          float64
Y          float64
Z          float64
surface     object
dtype: object

Data initialization:

Suggested size of the axis-aligned modeling box: Origin: 548800 7816600 -8400 Maximum: 552500 7822000 -11010

Suggested resolution: 100m x 100m x -90m (grid size 38 x 55 x 30)

Number of voxels:

np.array([38, 55, 30]).prod()

surface_points_table: gp.data.SurfacePointsTable = gp.data.SurfacePointsTable.from_arrays(
    x=surface_points['X'].values,
    y=surface_points['Y'].values,
    z=surface_points['Z'].values,
    names=surface_points['surface'].values
)

orientations_table: gp.data.OrientationsTable = gp.data.OrientationsTable.from_arrays(
    x=orientations['X'].values,
    y=orientations['Y'].values,
    z=orientations['Z'].values,
    G_x=orientations['G_x'].values,
    G_y=orientations['G_y'].values,
    G_z=orientations['G_z'].values,
    names=orientations['surface'].values,
    name_id_map=surface_points_table.name_id_map  # ! Make sure that ids and names are shared
)

structural_frame: gp.data.StructuralFrame = gp.data.StructuralFrame.from_data_tables(
    surface_points=surface_points_table,
    orientations=orientations_table
)

geo_model: gp.data.GeoModel = gp.create_geomodel(
    project_name='Claudius',
    extent=[548800, 552500, 7816600, 7822000, -11010, -8400],
    resolution=[38, 55, 30],
    refinement=5,
    structural_frame=structural_frame
)

group_fault = gp.data.StructuralGroup(
    name='Fault1',
    elements=[geo_model.structural_frame.structural_elements.pop(-2)],
    structural_relation=gp.data.StackRelationType.FAULT,
    fault_relations=gp.data.FaultsRelationSpecialCase.OFFSET_ALL
)

geo_model.structural_frame.get_group_by_name("default_formation").elements.pop(-1)

# Insert the fault group into the structural frame:
geo_model.structural_frame.insert_group(0, group_fault)

gp.set_is_fault(
    frame=geo_model.structural_frame,
    fault_groups=[geo_model.structural_frame.get_group_by_name('Fault1')]
)

print(geo_model)
{'grid': <gempy.core.data.grid.Grid object at 0x7f08ac9d3bb0>,
 'interpolation_options': InterpolationOptions(kernel_options={'range': 5, 'c_o': 10, 'uni_degree': 1, 'i_res': 4, 'gi_res': 2, 'number_dimensions': 3, 'kernel_function': <AvailableKernelFunctions.cubic: KernelFunction(base_function=<function cubic_function at 0x7f090b7b1cf0>, derivative_div_r=<function cubic_function_p_div_r at 0x7f090b7b1d80>, second_derivative=<function cubic_function_a at 0x7f090b7b1e10>, consume_sq_distance=False)>, 'compute_condition_number': False, 'kernel_solver': <Solvers.DEFAULT: 1>}, number_octree_levels=5, current_octree_level=0, compute_scalar_gradient=False, mesh_extraction=True, mesh_extraction_masking_options=MeshExtractionMaskingOptions.INTERSECT, mesh_extraction_fancy=True, debug=True, debug_water_tight=False, sigmoid_slope=50000, _number_octree_levels_surface=4),
 'meta': GeoModelMeta(name='Claudius',
                      creation_date=None,
                      last_modification_date=None,
                      owner=None),
 'structural_frame': StructuralFrame(
        structural_groups=[
StructuralGroup(
        name=Fault1,
        structural_relation=StackRelationType.FAULT,
        elements=[
Element(
        name=Claudius_fault,
        color=#527682,
        is_active=True
)
]
),
StructuralGroup(
        name=default_formation,
        structural_relation=StackRelationType.ERODE,
        elements=[
Element(
        name=0,
        color=#015482,
        is_active=True
),
Element(
        name=250,
        color=#9f0052,
        is_active=True
),
Element(
        name=330,
        color=#ffbe00,
        is_active=True
),
Element(
        name=60,
        color=#728f02,
        is_active=True
)
]
)
],
        fault_relations=
[[False,  True],
 [False, False]],
,
 'transform': {'_is_default_transform': False,
 'position': array([ -550658.0605    , -7819213.97425586,     9824.6748045 ]),
 'rotation': array([0., 0., 0.]),
 'scale': array([9.56397917e-05, 9.56397917e-05, 9.56397917e-05])}}

We are going to increase the smoothness (nugget) of the data to increase the conditional number of the matrix:

gp.modify_surface_points(geo_model, nugget=0.01)
Structural Groups: StructuralGroup:
Name:Fault1
Structural Relation:StackRelationType.FAULT
Elements:
StructuralElement:
Name:Claudius_fault

StructuralGroup:
Name:default_formation
Structural Relation:StackRelationType.ERODE
Elements:
StructuralElement:
Name:0

StructuralElement:
Name:250

StructuralElement:
Name:330

StructuralElement:
Name:60
Fault Relations:
Fault1default_fo...
Fault1
default_formation
True
False


Also the original poles are pointing downwards. We can change the direction by calling the following:

gp.modify_orientations(geo_model, polarity=-1)
Structural Groups: StructuralGroup:
Name:Fault1
Structural Relation:StackRelationType.FAULT
Elements:
StructuralElement:
Name:Claudius_fault

StructuralGroup:
Name:default_formation
Structural Relation:StackRelationType.ERODE
Elements:
StructuralElement:
Name:0

StructuralElement:
Name:250

StructuralElement:
Name:330

StructuralElement:
Name:60
Fault Relations:
Fault1default_fo...
Fault1
default_formation
True
False


We need an orientation per series/fault. The faults does not have orientation so the easiest is to create an orientation from the surface points availablle:

element = geo_model.structural_frame.get_element_by_name("Claudius_fault")
new_orientations: gp.data.OrientationsTable = gp.create_orientations_from_surface_points_coords(
    xyz_coords=element.surface_points.xyz
)
gp.add_orientations(
    geo_model=geo_model,
    x=new_orientations.data['X'],
    y=new_orientations.data['Y'],
    z=new_orientations.data['Z'],
    pole_vector=new_orientations.grads,
    elements_names="Claudius_fault"
)
Structural Groups: StructuralGroup:
Name:Fault1
Structural Relation:StackRelationType.FAULT
Elements:
StructuralElement:
Name:Claudius_fault

StructuralGroup:
Name:default_formation
Structural Relation:StackRelationType.ERODE
Elements:
StructuralElement:
Name:0

StructuralElement:
Name:250

StructuralElement:
Name:330

StructuralElement:
Name:60
Fault Relations:
Fault1default_fo...
Fault1
default_formation
True
False


gpv.plot_2d(geo_model, direction='y')
Cell Number: mid Direction: y
<gempy_viewer.modules.plot_2d.visualization_2d.Plot2D object at 0x7f08acd4f0d0>

We will need to separate with surface belong to each series:

gp.map_stack_to_surfaces(
    gempy_model=geo_model,
    mapping_object={
        'Default series': ('0', '60', '250'),
        'Fault': 'Claudius_fault',
        'Uncomformity': '330',
    }
)
Structural Groups: StructuralGroup:
Name:Default series
Structural Relation:StackRelationType.ERODE
Elements:
StructuralElement:
Name:0

StructuralElement:
Name:60

StructuralElement:
Name:250

StructuralGroup:
Name:Fault
Structural Relation:StackRelationType.ERODE
Elements:
StructuralElement:
Name:Claudius_fault

StructuralGroup:
Name:Uncomformity
Structural Relation:StackRelationType.ERODE
Elements:
StructuralElement:
Name:330
Fault Relations:
Default se...FaultUncomformi...
Default series
Fault
Uncomformity
True
False


So far we did not specify which series/faults are actula faults:

gp.set_is_fault(
    frame=geo_model.structural_frame,
    fault_groups=[geo_model.structural_frame.get_group_by_name('Fault')]
)

geo_model.structural_frame
Structural Groups: StructuralGroup:
Name:Default series
Structural Relation:StackRelationType.ERODE
Elements:
StructuralElement:
Name:0

StructuralElement:
Name:60

StructuralElement:
Name:250

StructuralGroup:
Name:Fault
Structural Relation:StackRelationType.FAULT
Elements:
StructuralElement:
Name:Claudius_fault

StructuralGroup:
Name:Uncomformity
Structural Relation:StackRelationType.ERODE
Elements:
StructuralElement:
Name:330
Fault Relations:
Default se...FaultUncomformi...
Default series
Fault
Uncomformity
True
False


geo_model.interpolation_options.kernel_options.range = 1
gp.compute_model(
    geo_model,
    gp.data.GemPyEngineConfig(
        backend=gp.data.AvailableBackends.numpy,
        use_gpu=False,
        dtype='float64'
    )
)
Setting Backend To: AvailableBackends.numpy
/home/leguark/gempy/gempy/core/data/geo_model.py:164: UserWarning: You are using refinement and passing a regular grid. The resolution of the regular grid will be overwritten
  warnings.warn(
Solutions: 5 Octree Levels, 5 DualContouringMeshes


sect = ['mid']

gpv.plot_2d(geo_model, cell_number=sect, series_n=1, show_scalar=True, direction='x')
Cell Number: mid Direction: x
<gempy_viewer.modules.plot_2d.visualization_2d.Plot2D object at 0x7f08ac70d5a0>
gpv.plot_2d(geo_model, cell_number=sect, show_data=True, direction='x')
Cell Number: mid Direction: x
<gempy_viewer.modules.plot_2d.visualization_2d.Plot2D object at 0x7f08f514e260>
gpv.plot_2d(geo_model, cell_number=[28], series_n=0, direction='y', show_scalar=True)
gpv.plot_2d(geo_model, cell_number=[28], series_n=1, direction='y', show_scalar=True)
gpv.plot_2d(geo_model, cell_number=[28], series_n=2, direction='y', show_scalar=True)
  • Cell Number: 28 Direction: y
  • Cell Number: 28 Direction: y
  • Cell Number: 28 Direction: y
<gempy_viewer.modules.plot_2d.visualization_2d.Plot2D object at 0x7f08f2031930>
gpv.plot_2d(geo_model, cell_number=[28], show_data=True, direction='y')
Cell Number: 28 Direction: y
<gempy_viewer.modules.plot_2d.visualization_2d.Plot2D object at 0x7f08f546ca90>
# sphinx_gallery_thumbnail_number = 8
gpv.plot_3d(geo_model, show_lith=True, show_data=True, show_boundaries=True)
Claudius
<gempy_viewer.modules.plot_3d.vista.GemPyToVista object at 0x7f08f2031930>

Total running time of the script: (2 minutes 5.939 seconds)

Gallery generated by Sphinx-Gallery