推荐：用 NSDT设计器 快速搭建可编程3D场景。
在设计任何类型的工程结构时，确定地面以下的东西并有效地将其映射出来是首要也是最重要的部分之一。 地下建模带有很大的误差范围，因为即使是我们今天使用的最先进的地下调查方法也无法完全绘制出更大区域的土壤、岩石、水位等。 无论你使用的是钻孔（SPT 或 CPT）、地震折射（seismic refraction）还是其他一些方法，最终都必须进行大量插值来构建你的地下模型。
Bentley 的 GINT 等软件非常适合创建 2D 和 3D 土壤剖面，其精度实际上仅受所提供数据和工程师经验的限制。 这些地下模型也可以使用 Python 创建，虽然它不像 GINT 那样简单，但它给了我们更多的自由来按照我们认为合适的方式操作我们的模型。 一旦掌握了它，你就可以比某些软件更快地生成良好的土壤剖面，而且绝对比在 CAD 中绘制它更快。
下面是有关如何使用 Python 创建地下3D模型的教程，主要关注开源地质建模库 GemPy。 它能够构建复杂的褶皱结构（folded structures）、断层网络（fault networks）和不整合面（unconformities）的 3D 地质模型。
1、安装配置GemPy
可以使用如下命令安装GemPy：
pip install gempy
我们需要来自项目现场和地下调查的数据，需要CSV 格式，因为我们将使用 pandas 将它们转换为数据框：
- Surface_points →该区域 (x,y,z) 的现场测量表面高程
- Orientations → 为表面点（Gx、Gy、Gz 或方位角）分配地层值。 本质上，它是与表示穿过层的断层的表面点相关联的角度
- Grid → 模型的大小（列数 x 行数 x 层数）
- Surfaces →土壤/岩石层
- Series → 将表面和断层分组为单独的系列
- Faults →要建模的断层的位置和信息（必须单独插入）
对于本教程，我们将包含这 6 个参数的 CSV 文件称为“MyProject.csv”。
使用 GemPy 建模需要先导入必要的包：
from pyvista import set_plot_theme
%matplotlib inline
import gempy as gp
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
2、使用GemPy载入地质数据
首先，我们需要将包含站点信息的 CSV 文件中的数据提取到数据框中，并使用 GemPy 创建模型。
geo_model = gp.create_model('MyProject')
# Importing the data from CSV-files and setting our boundaries
gp.init_data(geo_model, [0, 2000., 0, 2000., 0, 2000.], [100, 100, 100],
path_o= "MyProject_model_orientations.csv",
path_i= " MyProject_model_points.csv",)
这将生成地址模型，我们的数据为 100 行 x 100 列 x 100 层，跨越 2 公里 X 2 公里 X 2 公里的体积，分别是我们的“分辨率”和“网格大小”。 你做得越大，运行模型所需的时间就越长。 该模型同时接收 Surface_point 数据 (path_i) 和将用于插值的方向 (path_o) 并生成我们的 2D 和 3D 图。
使用 .surface() 我们看到用“MyProject”数据制作的 geo_model 有 5 个这样的表面，其中Basement是自动生成的，并充当下界：
geo_model.surfaces
表面层需要分成两个单独的分类系列，地层表面层 (Strat_Series) 和断层 (Fault_Series)。 当涉及到这些表面的顺序时，GemPy 有点挑剔。断层总是需要移动到顶部，相互独立，它们的顺序决定了哪个先穿过表面（构造关系） . Strat_Series 表面的顺序应代表我们正在建模的真实世界位置，地下室始终位于底部。 我们可以通过以下方式实现：
gp.map_stack_to_surfaces(geo_model,
{"Fault_Series": 'Main_Fault',
"Strat_Series": ('Sandstone_2',
'Siltstone', 'Shale', 'Sandstone_1',
'basement')},
remove_unused_series=True)
geo_model.surfaces
我们所做的只是将一个 Fault_Series 放在最上面，但现在模型表面是有序的，Gempy 应该可以毫无问题地处理它们。
3、使用GemPy可视化初始数据
我们可以为初始数据创建 2D 和 3D 图以进行可视化，这也向我们展示了界面和方向。
# 2-D
gp.plotta(geo_data,
direction='y')
# 3-D
gp.plot_3d(geo_model,
plotter_type='basic')
4、使用GemPy插值地质数据
上面的图看起来没什么特别的，尤其是在 3 维空间中。 我们想填满这个网格，基本上使它成为一个 2km x 2km x 2 km 的实心块，就好像你已经撕掉了一个巨大的地球块来检查它一样。 为此，我们将通过对当前数据进行插值来填充绘图。
GemPy 使用称为 Kriging 的方法进行插值，也称为高斯过程回归，这是一种由协方差控制的方法，可在未采样位置提供最佳线性无偏预测。 我们这样插入数据：
interpol_data = gp.InterpolatorData(geo_data, output='geology',
compile_theano=True,
theano_optimizer='fast_compile')
可以根据需要更改克里金参数，但我个人不会乱用它们。
gp.get_data(geo_model, 'kriging')
5、使用GemPy计算地质模型
插值后，我们可以使用 GemPy 创建我们的模型：
gp.compute_model(interpol_data)
生成的模型将生成两个数组，一个用于地质表面，一个用于穿过地表的断层（Faults）。
interpol_data.geo_data_res.get_formations()
结果如下：
Surfaces Array: [Main_Fault, Sandstone_2, Siltstone, Shale, Sandstone_1]
Fault Array:(5, object) [Main_Fault, Sandstone_2, Siltstone, Shale, Sandstone_1]
这两个数组也有 2 个子数组作为块模型的解。
6、使用GemPy可视化地质模型
可以2D或3D的方式可视化GemPy生成的地质模型，并且可以添加地形图。
6.1 2D可视化
使用 Matplotlib，我们可以绘制刚刚生成的模型并创建我们一直追求的土壤剖面。 下面的代码基本上将取出 2km³ 块（由 100 行 x 100 列 x 100 层网格表示）并在 y=50 标记处对其进行切片，从而为我们提供模型中间的横截面。
gp.plotting.plot_section(geo_data,
gp.compute_model(interpol_data)[0],
cell_number=50, direction='y',
plot_data=False)
6.2 3D可视化
我们可以使用 VTK 以 3d 形式可视化我们的模型，VTK 是一个用于 3D 计算机图形、图像处理和科学可视化的开源软件系统。 它已经整合到 GemPy 中，因此无需导入任何新内容。
我们首先需要从我们的表面和函数 get_surface 中获取顶点和三角剖分。 然后我们可以使用这些来绘制我们的 3-D 模型：
ver, sim = gp.get_surfaces(interpol_data, original_scale=True)
gp.plotting.plot_surfaces_3D(geo_data, ver, sim,
plotter_type='basic')
6.3 添加地形图
添加地形总能使模型看起来比只是一些块更好。 地形可以高度自定义，因为你可以使用真实世界的光栅图像，创建自己的表面叠加层，基于海拔的颜色变化，还有很多选择。 我们将根据深度和图层颜色添加地形图。
geo_model.set_topography(d_z=(1000, 2000))
gp.plot_3d(interpol_data, plotter_type='basic',
show_topography=True, show_surfaces=True,
show_lith=True, image=False)
7、结束语
这个教程是为了创建一个基本的地质 3D 模型，你可以从中提取 2D 横截面土壤剖面并将它们用于设计。 它肯定需要结合地下水高程、土壤/岩石参数等内容，使每一层真正变得有用。 最后，我们得到一个交互性很强的模型，可以轻松地从中提取信息。
原文链接：GemPy地质建模实战 — BimAnt
