This chapter shows how to make geological models look more realistic by incorporating topography data: The data can either be stored in a raster file that can be conveniently loaded into gempy. For demonstration purposes GemPy also provides a tool to create a random hilly landscape to limit the models on the surface.
import sys
sys.path.append("../../..")
import gempy as gp
import numpy as np
import matplotlib.pyplot as plt
import os
C:\Users\elisa\Anaconda3\envs\gempy_n\lib\site-packages\dask\config.py:168: YAMLLoadWarning: calling yaml.load() without Loader=... is deprecated, as the default Loader is unsafe. Please read https://msg.pyyaml.org/load for full details. data = yaml.load(f.read()) or {} WARNING (theano.configdefaults): g++ not available, if using conda: `conda install m2w64-toolchain` C:\Users\elisa\Anaconda3\envs\gempy_n\lib\site-packages\theano\configdefaults.py:560: UserWarning: DeprecationWarning: there is no c++ compiler.This is deprecated and with Theano 0.11 a c++ compiler will be mandatory warnings.warn("DeprecationWarning: there is no c++ compiler." WARNING (theano.configdefaults): g++ not detected ! Theano will be unable to execute optimized C-implementations (for both CPU and GPU) and will default to Python implementations. Performance will be severely degraded. To remove this warning, set Theano flags cxx to an empty string. WARNING (theano.tensor.blas): Using NumPy C-API based implementation for BLAS functions.
geo_model = gp.create_model('Single_layer_topo')
data_path= '../..'
gp.init_data(geo_model, extent=[450000, 460000, 70000,80000,-1000,500],resolution = (50,50,50),
path_i = data_path+"/data/input_data/tut-ch1-7/onelayer_interfaces.csv",
path_o = data_path+"/data/input_data/tut-ch1-7/onelayer_orient.csv")
Active grids: ['regular']
Single_layer_topo 2019-12-03 20:01
# use happy spring colors!
geo_model.surfaces.colors.change_colors({'layer1':'#ff8000','basement':'#88cc60'})
surface | series | order_surfaces | isBasement | color | id | |
---|---|---|---|---|---|---|
0 | layer1 | Default series | 1 | False | #ff8000 | 1 |
1 | basement | Basement | 1 | True | #88cc60 | 2 |
%matplotlib inline
gp.map_series_to_surfaces(geo_model, {'series':('layer1','basement')})
surface | series | order_surfaces | isBasement | color | id | |
---|---|---|---|---|---|---|
0 | layer1 | series | 1 | False | #ff8000 | 1 |
1 | basement | series | 2 | True | #88cc60 | 2 |
s = {'s1': ([450000,75000],[460000,75500],[100,100])}
geo_model.grid.create_section_grid(s)
start | stop | resolution | dist | |
---|---|---|---|---|
s1 | [450000, 75000] | [460000, 75500] | [100, 100] | 10012.492197 |
fp = data_path+"/data/input_data/tut-ch1-7/bogota.tif"
geo_model.set_topography(source='gdal',filepath=fp)
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([[450100. , 70100. , -985. ], [450100. , 70100. , -955. ], [450100. , 70100. , -925. ], ..., [460000. , 75500. , 469.6969697 ], [460000. , 75500. , 484.84848485], [460000. , 75500. , 500. ]])
If there is no topography file, but you think that your model with topography would look significantly cooler, you can use gempys function to generate a random topography based on a fractal grid:
geo_model.set_topography(source='random')
[200. 500.] Active grids: ['regular' 'topography' 'sections']
Grid Object. Values: array([[450100. , 70100. , -985. ], [450100. , 70100. , -955. ], [450100. , 70100. , -925. ], ..., [460000. , 75500. , 469.6969697 ], [460000. , 75500. , 484.84848485], [460000. , 75500. , 500. ]])
It has additional keywords to play around with:
geo_model.set_topography(source='random',fd=1.9, d_z=np.array([0,250]), resolution=np.array([200,200]))
Active grids: ['regular' 'topography' 'sections']
Grid Object. Values: array([[450100. , 70100. , -985. ], [450100. , 70100. , -955. ], [450100. , 70100. , -925. ], ..., [460000. , 75500. , 469.6969697 ], [460000. , 75500. , 484.84848485], [460000. , 75500. , 500. ]])
Note that each time this function is called, a new random topography is created. If you particularly like the generated topography or if you have loaded a large file with gdal, you can save the topography object and load it again later:
#save
geo_model.grid.topography.save('test_topo')
saved
#load
geo_model.set_topography(source='saved',filepath='test_topo.npy')
Active grids: ['regular' 'topography' 'sections']
Grid Object. Values: array([[450100. , 70100. , -985. ], [450100. , 70100. , -955. ], [450100. , 70100. , -925. ], ..., [460000. , 75500. , 469.6969697 ], [460000. , 75500. , 484.84848485], [460000. , 75500. , 500. ]])
gp.set_interpolation_data(geo_model,
compile_theano=True,
theano_optimizer='fast_compile')
Compiling theano function... Level of Optimization: fast_compile Device: cpu Precision: float64 Number of faults: 0 Compilation Done!
<gempy.core.interpolator.InterpolatorModel at 0x26b4a7d9e80>
gp.compute_model(geo_model, compute_mesh=False, set_solutions=True)
Lithology ids [2. 2. 2. ... 1. 1. 1.]
Now, the solutions object does also contain the computed geological map. It can be visualized using the plot_map function:
gp.plot.plot_map(geo_model, contour_lines=False, show_data=False)
gp.plot.plot_section(geo_model, cell_number=4, block=geo_model.solutions.lith_block,
direction='y', show_data=True,show_faults=False)
<gempy.plot.visualization_2d.PlotSolution at 0x26b4b7459b0>
gp.plot.plot_section_by_name(geo_model,'s1')
../../..\gempy\plot\visualization_2d.py:208: UserWarning: the orientations are not converted to apparent dip. warnings.warn('the orientations are not converted to apparent dip.')