#import warnings
import sys
sys.path.append("../..")
import gempy as gp
import numpy as np
import matplotlib.pyplot as plt
import os
path_interf = os.pardir+"/data/input_data/AlesModel/2018_interf.csv"
path_orient = os.pardir+"/data/input_data/AlesModel/2018_orient_clust_n_init5_0.csv"
path_dem = os.pardir+"/data/input_data/AlesModel/_cropped_DEM_coarse.tif"
resolution = [100, 100, 100]
extent = np.array([729550.0, 751500.0, 1913500.0, 1923650.0, -1800.0, 800.0])
geo_model = gp.create_model('Alesmodel')
gp.init_data(geo_model, extent = extent, resolution = resolution,
path_i = path_interf,
path_o = path_orient)
Active grids: ['regular']
../../gempy/core/data.py:1742: UserWarning: If pole_vector and orientation are passed pole_vector is used/ warnings.warn('If pole_vector and orientation are passed pole_vector is used/')
Alesmodel 2020-03-05 15:30
sdict = {'section1':([732000, 1916000],[745000,1916000],[200,150])}
geo_model.grid.create_section_grid(sdict)
start | stop | resolution | dist | |
---|---|---|---|---|
section1 | [732000, 1916000] | [745000, 1916000] | [200, 150] | 13000.0 |
print(len(geo_model.orientations.df))
print(len(geo_model.surface_points.df))
print(len(geo_model.surfaces.df))
198 126 7
#sorting of lithologies
gp.map_series_to_surfaces(geo_model,{'fault_left':('fault_left'),
'fault_right':('fault_right'),
'fault_lr':('fault_lr'),
'Trias_Series':('TRIAS','LIAS'),
'Carbon_Series':('CARBO'),
'Basement_Series':('basement')},remove_unused_series=True)
surface | series | order_surfaces | isActive | color | id | |
---|---|---|---|---|---|---|
0 | fault_left | fault_left | 1 | True | #015482 | 1 |
2 | fault_right | fault_right | 1 | True | #ffbe00 | 2 |
1 | fault_lr | fault_lr | 1 | True | #9f0052 | 3 |
3 | TRIAS | Trias_Series | 1 | True | #728f02 | 4 |
4 | LIAS | Trias_Series | 2 | True | #443988 | 5 |
5 | CARBO | Carbon_Series | 1 | True | #ff3f20 | 6 |
6 | basement | Basement_Series | 1 | True | #5DA629 | 7 |
colordict = {'LIAS':'#015482', 'TRIAS': '#9f0052', 'CARBO':'#ffbe00','basement':'#728f02',
'fault_left':'#2a2a2a','fault_right':'#545454', 'fault_lr': '#a5a391'}
geo_model.surfaces.colors.change_colors(colordict)
surface | series | order_surfaces | isActive | color | id | |
---|---|---|---|---|---|---|
0 | fault_left | fault_left | 1 | True | #2a2a2a | 1 |
2 | fault_right | fault_right | 1 | True | #545454 | 2 |
1 | fault_lr | fault_lr | 1 | True | #a5a391 | 3 |
3 | TRIAS | Trias_Series | 1 | True | #9f0052 | 4 |
4 | LIAS | Trias_Series | 2 | True | #015482 | 5 |
5 | CARBO | Carbon_Series | 1 | True | #ffbe00 | 6 |
6 | basement | Basement_Series | 1 | True | #728f02 | 7 |
a = gp.plot.plot_data(geo_model,direction='y')
geo_model.rescaling
values | |
---|---|
rescaling factor | 26804.2 |
centers | [739263.6925, 1917834.9612500002, 402.6439239999999] |
gp.plot.plot_section_traces(geo_model, contour_lines=True, show_all_data=False)
geo_model.set_is_fault(['fault_right', 'fault_left', 'fault_lr'], change_color=True)
Fault colors changed. If you do not like this behavior, set change_color to False.
isFault | isFinite | |
---|---|---|
fault_left | True | False |
fault_right | True | False |
fault_lr | True | False |
Trias_Series | False | False |
Carbon_Series | False | False |
Basement_Series | False | False |
gp.set_interpolation_data(geo_model,
output=['geology'], compile_theano=True,
theano_optimizer='fast_run', dtype='float64',
verbose=[])
Compiling theano function... Level of Optimization: fast_run Device: cpu Precision: float64 Number of faults: 3 Compilation Done!
<gempy.core.interpolator.InterpolatorModel at 0x7fdd68706c50>
geo_model.set_topography(source='gdal', filepath=path_dem)
Cropped raster to geo_model.grid.extent. depending on the size of the raster, this can take a while... storing converted file... Active grids: ['regular' 'topography' 'sections']
Grid Object. Values: array([[ 7.29659750e+05, 1.91355075e+06, -1.78700000e+03], [ 7.29659750e+05, 1.91355075e+06, -1.76100000e+03], [ 7.29659750e+05, 1.91355075e+06, -1.73500000e+03], ..., [ 7.45000000e+05, 1.91600000e+06, 7.65100671e+02], [ 7.45000000e+05, 1.91600000e+06, 7.82550336e+02], [ 7.45000000e+05, 1.91600000e+06, 8.00000000e+02]])
from scipy.spatial import Delaunay
tri=Delaunay(geo_model.grid.topography.values[:,:2])
f = tri.simplices
f
array([[1900, 1899, 1826], [1900, 1826, 1827], [1201, 1129, 1202], ..., [1622, 1549, 1623], [1549, 1477, 1550], [1549, 1476, 1477]], dtype=int32)
geo_model.grid.topography.values_3D.shape, geo_model.grid.topography.values_3D_res.shape, geo_model.grid.topography.values.shape
((34, 73, 3), (100, 100, 3), (2482, 3))
geo_model.surfaces
surface | series | order_surfaces | isActive | color | id | |
---|---|---|---|---|---|---|
0 | fault_left | fault_left | 1 | True | #527682 | 1 |
2 | fault_right | fault_right | 1 | True | #527682 | 2 |
1 | fault_lr | fault_lr | 1 | True | #527682 | 3 |
3 | TRIAS | Trias_Series | 1 | True | #9f0052 | 4 |
4 | LIAS | Trias_Series | 2 | True | #015482 | 5 |
5 | CARBO | Carbon_Series | 1 | True | #ffbe00 | 6 |
6 | basement | Basement_Series | 1 | True | #728f02 | 7 |
_=gp.compute_model(geo_model, compute_mesh=True, compute_mesh_options={'rescale': True})
../../gempy/core/solution.py:280: UserWarning: Attribute error. Using non masked marching cubesmarching_cubes_lewiner() got an unexpected keyword argument 'mask'. warnings.warn('Attribute error. Using non masked marching cubes' + str(e)+'.')
geo_model.interpolator.theano_graph.n_surfaces_per_series.get_value()
array([0, 1, 2, 3, 5, 6], dtype=int32)
geo_model.interpolator.print_theano_shared()
len sereies i [ 0 10 17 21 88 120] len sereies o [ 0 5 6 9 110 198] len sereies w [ 0 28 41 57 433 735] n surfaces per series [0 1 2 3 5 6] n universal eq [3 3 3 3 3] is finite [0 0 0 0 0 0] is erosion [0 0 0 1 0] is onlap [0 0 0 0 0]
gp.plot.plot_section(geo_model, 4, direction='y',show_topo=True, show_data=True, show_faults=True, show_all_data=True)
p1 [729550.0, 1913906.0] p2 [751500.0, 1913906.0]
<gempy.plot.visualization_2d.PlotSolution at 0x7fdd60c2c780>
gp.plot.plot_map(geo_model, show_data=False, contour_lines=False)
#np.save('Ales_vert3', geo_model.solutions.vertices)
#np.save('Ales_edges3', geo_model.solutions.edges)
gp.plot.plot_ar(geo_model)