Example: Mapping buried valleys in Kasted, Denmark¶
[1]:
%load_ext autoreload
%autoreload 2
import mpslib as mps
import pyvista as pv
import numpy as np
import matplotlib.pyplot as plt
plt.rc('figure', figsize=(20.0, 10.0))
pv.set_plot_theme("document")
pv.global_theme.jupyter_backend = 'panel' # use this in jupyter lab
#pv.global_theme.jupyter_backend = 'pythreejs' # use this in notebook
Get the training image and conditional data¶
[2]:
# Get the training image
dx=100
TI, TI_fname = mps.trainingimages.kasted(dx=dx)
url_base = 'https://raw.githubusercontent.com/ergosimulation/mpslib/master/data/kasted'
remote_files=['kasted_soft_well.dat', 'kasted_soft_ele.dat', 'kasted_soft_res.dat','kasted_hard_well_consistent.dat' ]
for local_file in remote_files:
mps.trainingimages.get_remote('%s/%s' %(url_base,local_file), local_file)
# Get conditional point data
EAS_well=mps.eas.read(remote_files[0])
EAS_ele=mps.eas.read(remote_files[1])
EAS_res=mps.eas.read(remote_files[2])
EAS_well_hard=mps.eas.read(remote_files[3])
Beginning download of https://raw.githubusercontent.com/ergosimulation/mpslib/master/data/kasted/kasted_ti_dx100.dat to ti_kasted_dx100.dat
Beginning download of https://raw.githubusercontent.com/ergosimulation/mpslib/master/data/kasted/kasted_soft_well.dat to kasted_soft_well.dat
Beginning download of https://raw.githubusercontent.com/ergosimulation/mpslib/master/data/kasted/kasted_soft_ele.dat to kasted_soft_ele.dat
Beginning download of https://raw.githubusercontent.com/ergosimulation/mpslib/master/data/kasted/kasted_soft_res.dat to kasted_soft_res.dat
Beginning download of https://raw.githubusercontent.com/ergosimulation/mpslib/master/data/kasted/kasted_hard_well_consistent.dat to kasted_hard_well_consistent.dat
[3]:
# Get geometry
x_pad = 4*dx
x_min = np.min(EAS_well['D'][:,0])-x_pad
x_max = np.max(EAS_well['D'][:,0])+x_pad
y_min = np.min(EAS_well['D'][:,1])-x_pad
y_max = np.max(EAS_well['D'][:,1])+x_pad
z_min = np.min(EAS_well['D'][:,2])
z_max = np.max(EAS_well['D'][:,2])
x_min = dx*np.floor(x_min/dx)
y_min = dx*np.floor(y_min/dx)
z_min = dx*np.floor(z_min/dx)
nx=np.int16(np.ceil((x_max-x_min)/dx))
ny=np.int16(np.ceil((y_max-y_min)/dx))
nz=np.max([np.int16(np.ceil((z_max-z_min)/dx)),1])
grid_cell_size = np.array([1, 1, 1])*dx
origin = np.array([x_min, y_min, z_min])
simulation_grid_size=np.array([nx, ny, nz])
print(origin)
print(grid_cell_size)
print(simulation_grid_size)
print('[%f,%f] [%f,%f] [%f,%f]' % (x_min, x_max, y_min, y_max, z_min, z_max))
[ 561300. 6225000. 0.]
[100 100 100]
[155 108 1]
[561300.000000,576787.000000] [6225000.000000,6235796.000000] [0.000000,0.000000]
[4]:
# REMOVE ALL NON-INFORMATIVE SOFT DATA
print(EAS_well['D'].shape)
informed_well_soft_data=np.argwhere(np.abs(EAS_well['D'][:,4]-0.5)>0.05)
i_use=informed_well_soft_data.flatten()
EAS_well['D']=EAS_well['D'][i_use,:]
print(EAS_well['D'].shape)
(1254, 5)
(147, 5)
PyVIsta plotting of conditional data¶
[5]:
# Create PyVista meshes
mesh_well=pv.PolyData(EAS_well['D'][:,0:3])
mesh_well["Pv"]=EAS_well['D'][:,3]
mesh_ele=pv.PolyData(EAS_ele['D'][:,0:3])
mesh_ele["Pv"]=EAS_ele['D'][:,3]
mesh_res=pv.PolyData(EAS_res['D'][:,0:3])
mesh_res["Pv"]=EAS_res['D'][:,3]
cmap = plt.cm.get_cmap("viridis", 255)
cmap = plt.cm.get_cmap("PuOr", 255)
#from pyvistaqt import BackgroundPlotter
#l = BackgroundPlotter()
pl = pv.Plotter(shape=(1,3))
pl.subplot(0,0)
pl.add_text('P(channel|resistivity)', font_size=2)
pl.add_mesh(mesh_res, point_size=1, render_points_as_spheres=True, style='points', show_scalar_bar=True, cmap=cmap, clim=[0, 1])
pl.subplot(0,1)
pl.add_text('P(channel|ele)', font_size=30)
pl.add_mesh(mesh_ele, point_size=5, render_points_as_spheres=False, style='points', show_scalar_bar=True, cmap=cmap, clim=[0, 1])
pl.subplot(0,2)
pl.add_text('P(channel|well)', font_size=30)
pl.add_mesh(mesh_well, point_size=4, render_points_as_spheres=False, style='points', show_scalar_bar=True, cmap=cmap, clim=[0, 1])
pl.show_axes()
pl.show_grid()
pl.show()
[6]:
# Plot all soft data
cmap = plt.cm.get_cmap("viridis", 255)
cmap = plt.cm.get_cmap("PuOr", 255)
#from pyvistaqt import BackgroundPlotter
#l = BackgroundPlotter()
pl = pv.Plotter()
pl.add_mesh(mesh_res, point_size=10001, render_points_as_spheres=False, style='points', cmap=cmap, clim=[0, 1])
pl.add_mesh(mesh_ele, point_size=5000, render_points_as_spheres=True, style='points', cmap=cmap, clim=[0, 1])
pl.add_mesh(mesh_well, point_size=80000, render_points_as_spheres=True, style='points',cmap=cmap, clim=[0, 1])
pl.show_axes()
pl.show_grid()
pl.show()
MPSlib in Kasted¶
Setup and run MPSlib¶
[7]:
O=mps.mpslib(method='mps_snesim_tree', verbose_level=-11,
simulation_grid_size=simulation_grid_size,
origin=origin,
grid_cell_size=grid_cell_size)
O.ti=TI
O.par['n_real']=24;
# Set soft data
O.d_soft=EAS_well['D']
# Set Hard
# O.d_hard=EAS_well_hard['D']
[8]:
print(O.sim)
None
[10]:
# Make some plots
# Plot the simulation grid
O.plot_simulation_grid()
# Plot the simulation training image
O.plot_ti()
# plot the hard data
O.plot_hard()
# plot the soft data
O.plot_soft()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Input In [10], in <cell line: 10>()
7 O.plot_ti()
9 # plot the hard data
---> 10 O.plot_hard()
11 # plot the soft data
12 O.plot_soft()
File /mnt/f/PROGRAMMING/mpslib/scikit-mps/mpslib/mpslib.py:1037, in mpslib.plot_hard(self, **kwargs)
1034 import numpy as np
1035 import pyvista as pv
-> 1037 update_xyz()
1039 cmap = kwargs.get('cmap',"viridis")
1041 vmin=0
NameError: name 'update_xyz' is not defined
[ ]:
# Run MpsLib
O.delete_local_files()
#O.run()
O.run_parallel();
Plot MPSlib results¶
[ ]:
O.plot_reals()
O.plot_etype()
Estimation¶
[ ]:
O.par['do_estimation']=1;
O.par['do_entropy']=1;
O.par['n_real']=1
#O.run()
[ ]:
O.plot()
[ ]:
O.plot_soft()
[ ]: