22 Creating Temperature Maps from GemPy Models

Temperature maps are an important tool for Geologists to visualize the temperature distribution of layers at depth. A simple model (Example 1) is created of which temperature maps for the two existing layers are created. The temperature will be calculated as function of the thickness between a given topography and the distance to the layer at depth.

Set File Paths and download Tutorial Data

If you downloaded the latest GemGIS version from the Github repository, append the path so that the package can be imported successfully. Otherwise, it is recommended to install GemGIS via pip install gemgis and import GemGIS using import gemgis as gg. In addition, the file path to the folder where the data is being stored is set. The tutorial data is downloaded using Pooch (https://www.fatiando.org/pooch/latest/index.html) and stored in the specified folder. Use pip install pooch if Pooch is not installed on your system yet.

import gemgis as gg

file_path ='data/22_creating_temperature_maps_from_gempy_models/'

WARNING (theano.configdefaults): g++ not available, if using conda: `conda install m2w64-toolchain`
C:\Users\ale93371\Anaconda3\envs\test_gempy\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.

gg.download_gemgis_data.download_tutorial_data(filename="22_creating_temperature_maps_from_gempy_models.zip", dirpath=file_path)

Loading the data

import geopandas as gpd
import rasterio

interfaces = gpd.read_file(file_path + 'interfaces.shp')
orientations = gpd.read_file(file_path + 'orientations.shp')
extent = [0,972,0,1069, 300, 800]
resolution = [50, 50, 50]

WARNING (theano.configdefaults): g++ not available, if using conda: `conda install m2w64-toolchain`
C:\Users\ale93371\Anaconda3\envs\test_gempy\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.

interfaces.head()

	level_0	level_1	formation	X	Y	Z	geometry
0	0	0	Sand1	0.26	264.86	353.97	POINT (0.256 264.862)
1	0	0	Sand1	10.59	276.73	359.04	POINT (10.593 276.734)
2	0	0	Sand1	17.13	289.09	364.28	POINT (17.135 289.090)
3	0	0	Sand1	19.15	293.31	364.99	POINT (19.150 293.313)
4	0	0	Sand1	27.80	310.57	372.81	POINT (27.795 310.572)

orientations['polarity'] = 1
orientations.head()

	formation	dip	azimuth	X	Y	Z	geometry	polarity
0	Ton	30.50	180.00	96.47	451.56	477.73	POINT (96.471 451.564)	1
1	Ton	30.50	180.00	172.76	661.88	481.73	POINT (172.761 661.877)	1
2	Ton	30.50	180.00	383.07	957.76	444.45	POINT (383.074 957.758)	1
3	Ton	30.50	180.00	592.36	722.70	480.57	POINT (592.356 722.702)	1
4	Ton	30.50	180.00	766.59	348.47	498.96	POINT (766.586 348.469)	1

Creating the GemPy Model

import sys
sys.path.append('../../../../gempy-master')
import gempy as gp

geo_model = gp.create_model('Model1')
geo_model

Model1  2021-01-02 13:32

Initiating the Model

import pandas as pd

gp.init_data(geo_model, extent, resolution,
             surface_points_df = interfaces,
             orientations_df = orientations,
             default_values=True)
geo_model.surfaces

Active grids: ['regular']

	surface	series	order_surfaces	color	id
0	Sand1	Default series	1	#015482	1
1	Ton	Default series	2	#9f0052	2

The vertices and edges are currently NaN values, so no model has been computed so far.

geo_model.surfaces.df

	surface	series	order_surfaces	isBasement	isFault	isActive	hasData	color	vertices	edges	sfai	id
0	Sand1	Default series	1	False	False	True	True	#015482	NaN	NaN	NaN	1
1	Ton	Default series	2	True	False	True	True	#9f0052	NaN	NaN	NaN	2

Mapping Stack to Surfaces

gp.map_stack_to_surfaces(geo_model,
                         {"Strat_Series": ('Sand1', 'Ton')},
                         remove_unused_series=True)
geo_model.add_surfaces('basement')

	surface	series	order_surfaces	color	id
0	Sand1	Strat_Series	1	#015482	1
1	Ton	Strat_Series	2	#9f0052	2
2	basement	Strat_Series	3	#ffbe00	3

Loading Topography

geo_model.set_topography(
    source='gdal', filepath='../../../../gemgis_data/data/22_creating_temperature_maps_from_gempy_models/raster.tif')

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']

Grid Object. Values:
array([[   9.72      ,   10.69      ,  305.        ],
       [   9.72      ,   10.69      ,  315.        ],
       [   9.72      ,   10.69      ,  325.        ],
       ...,
       [ 970.056     , 1059.28181818,  622.0892334 ],
       [ 970.056     , 1063.16909091,  622.06713867],
       [ 970.056     , 1067.05636364,  622.05786133]])

Setting Interpolator

gp.set_interpolator(geo_model,
                    compile_theano=True,
                    theano_optimizer='fast_compile',
                    verbose=[],
                    update_kriging = False
                    )

Compiling theano function...
Level of Optimization:  fast_compile
Device:  cpu
Precision:  float64
Number of faults:  0
Compilation Done!
Kriging values:
                   values
range            1528.90
$C_o$           55655.83
drift equations      [3]

<gempy.core.interpolator.InterpolatorModel at 0x201ed3abfa0>

Computing Model

sol = gp.compute_model(geo_model, compute_mesh=True)

The surfaces DataFrame now contains values for vertices and edges.

geo_model.surfaces.df

	surface	series	order_surfaces	isBasement	isFault	isActive	hasData	color	vertices	edges	sfai	id
0	Sand1	Strat_Series	1	False	False	True	True	#015482	[[29.160000000000004, 194.27877317428587, 305....	[[2, 1, 0], [2, 0, 3], [3, 4, 2], [2, 4, 5], [...	0.26	1
1	Ton	Strat_Series	2	False	False	True	True	#9f0052	[[29.160000000000004, 365.78652999877926, 305....	[[2, 1, 0], [2, 0, 3], [3, 4, 2], [2, 4, 5], [...	0.21	2
2	basement	Strat_Series	3	True	False	True	True	#ffbe00	NaN	NaN	NaN	3

Plotting the 3D Model

gpv = gp.plot_3d(geo_model, image=False, show_topography=True,
                 plotter_type='basic', notebook=True, show_lith=True)

Creating Temperature Maps

Temperature Maps in GemGIS can be created by passing the Digital Elevation Model as a Rasterio Object and the surface in the subsurface for which the temperature is supposed to be calculated.

Creating Depth Maps

By providing a list of surfaces, a dict containing the data for different surfaces is created.

dict_all = gg.visualization.create_depth_maps_from_gempy(geo_model=geo_model,
                                                         surfaces=['Sand1', 'Ton'])

dict_all

{'Sand1': [PolyData (0x201ee61f640)
    N Cells:    4174
    N Points:   2303
    X Bounds:   9.720e+00, 9.623e+02
    Y Bounds:   1.881e+02, 9.491e+02
    Z Bounds:   3.050e+02, 7.250e+02
    N Arrays:   1,
  '#015482'],
 'Ton': [PolyData (0x201fb4bf9a0)
    N Cells:    5111
    N Points:   2739
    X Bounds:   9.720e+00, 9.623e+02
    Y Bounds:   3.578e+02, 1.058e+03
    Z Bounds:   3.050e+02, 7.265e+02
    N Arrays:   1,
  '#9f0052']}

The temperature for the Sand Surface is now being calculated.

mesh = dict_all['Sand1'][0]
mesh.clear_arrays()
mesh

PolyData	Information
N Cells	4174
N Points	2303
X Bounds	9.720e+00, 9.623e+02
Y Bounds	1.881e+02, 9.491e+02
Z Bounds	3.050e+02, 7.250e+02
N Arrays	0

Coordinates of the vertices.

mesh.points

pyvista_ndarray([[ 29.16      , 194.27877317, 305.        ],
                 [  9.72      , 188.12027063, 305.        ],
                 [  9.72      , 203.11      , 314.39133763],
                 ...,
                 [958.49513855, 545.19      , 495.        ],
                 [962.28      , 546.6037711 , 495.        ],
                 [962.28      , 545.19      , 494.12916183]])

Opening the Digital Elevation Model

dem = rasterio.open(file_path + 'raster.tif')
dem.read(1)

array([[482.82904, 485.51953, 488.159  , ..., 618.8612 , 620.4424 ,
        622.05786],
       [481.6521 , 484.32193, 486.93958, ..., 618.8579 , 620.44556,
        622.06714],
       [480.52563, 483.18893, 485.80444, ..., 618.8688 , 620.4622 ,
        622.08923],
       ...,
       [325.49225, 327.21985, 328.94498, ..., 353.6889 , 360.03125,
        366.3984 ],
       [325.0538 , 326.78473, 328.51276, ..., 351.80603, 357.84106,
        363.96167],
       [324.61444, 326.34845, 328.0794 , ..., 350.09247, 355.87598,
        361.78635]], dtype=float32)

Creating the mesh containing the temperature data

mesh = gg.visualization.create_temperature_map(dem=dem,
                                               mesh=mesh,
                                               name= 'Thickness [m]',
                                               apply_threshold=True,
                                               tsurface=10,
                                               gradient=0.03)
mesh

Header	Data Arrays
UnstructuredGrid Information N Cells 3946 N Points 2130 X Bounds 9.720e+00, 9.623e+02 Y Bounds 1.881e+02, 9.491e+02 Z Bounds 3.050e+02, 7.250e+02 N Arrays 2	Name Field Type N Comp Min Max Thickness [m] Points float64 1 9.321e-02 2.020e+02 Temperature [°C] Points float64 1 1.000e+01 1.606e+01

mesh['Thickness [m]']

array([34.28413844, 41.96035767, 50.86764526, ...,  0.58549309,
        0.42822266,  0.42822266])

mesh['Temperature [°C]']

array([11.02852415, 11.25881073, 11.52602936, ..., 10.01756479,
       10.01284668, 10.01284668])

import pyvista as pv
surf = pv.PolyData(geo_model._grid.topography.values)

Plotting the temperature map and the topography. It can be seen that the temperatures at the edges are equal to 10 degrees C. These are the areas where the layer is outcropping at the surface. The temperatures are higher at the center of the surface where the layer has a higher thickness.

sargs = dict(fmt="%.0f", color='black')


p = pv.Plotter(notebook=True)
p.add_mesh(mesh=mesh, scalars='Temperature [°C]', cmap='coolwarm', scalar_bar_args=sargs)
p.add_mesh(surf, opacity=0.025)

p.set_background('white')
p.show_grid(color='black')
p.show()