The topography obscures the soil layer

[wdk11]

Hi gempy team，

1. I created a model and added terrain data, and found that the terrain covered the strata. Is this because the terrain forces the overlapping strata to be covered? And I found that when I draw 3D graphics, show_topography=True is one situation, and show_topography=False is another situation, so that the soil layer data I model with gempy will not exceed the ground at least.

2.I want to know how to export mesh grids for terrain like strata

3.I built a model using both ground scatter data and soil scatter data, and set the ground to erode. However, I found that some soil layers crossed over the ground, forming an interpenetration. Is this correct.The red color in the picture represents the ground

4.I can see soil layers emerging from the ground in the 3D model, but I still haven't found any soil layers emerging in the 2D profile I drew. Is this due to the internal algorithm of Gempy. I don't quite understand this situation, I'm sorry, can you explain it

Thank u.I am very puzzled about this.

[wdk11]

orientations-20 - 去除地面线.csv
orientations-20.csv
points-20 - 去除地面线.csv
points-20.csv
地形点.txt

`import numpy as np
import pandas as pd

%%

import gempy as gp
import gempy_viewer as gpv
import torch`
data_path = 'D:/Program Files/jupyter/data
topo_input1 = pd.read_table(data_path+"/地形点.txt")
XX = np.array(topo_input1.iloc[:,0])
YY = np.array(topo_input1.iloc[:,1])
ZZ = np.array(topo_input1.iloc[:,2])
topo_input = np.hstack([XX.reshape(-1,1),YY.reshape(-1,1),ZZ.reshape(-1,1)])
data_path = 'D:/Program Files/jupyter/data'
surface_points = pd.read_csv(data_path+"/points-20 - 去除地面线.csv")
orientations = pd.read_csv(data_path+"/orientations-20 - 去除地面线.csv")
surface_points_table: gp.data.SurfacePointsTable = gp.data.SurfacePointsTable.from_arrays(
x=surface_points['X'].values,
y=surface_points['Y'].values,
z=surface_points['Z'].values,
names=surface_points['formation'].values.astype(str)
)

orientations_table: gp.data.OrientationsTable = gp.data.OrientationsTable.from_arrays(
x=orientations['X'].values,
y=orientations['Y'].values,
z=orientations['Z'].values,
G_x=orientations['G_x'].values,
G_y=orientations['G_y'].values,
G_z=orientations['G_z'].values,
names=orientations['formation'].values.astype(str),
name_id_map=surface_points_table.name_id_map # ! Make sure that ids and names are shared
)

structural_frame: gp.data.StructuralFrame = gp.data.StructuralFrame.from_data_tables(
surface_points=surface_points_table,
orientations=orientations_table
)
geo_model: gp.data.GeoModel = gp.create_geomodel(
project_name='Example',
extent=[604071.3221, 604779.6705, 3005914.7854, 3007065.9975, 1850, 2320],
resolution = [30,30,30],
# resolution=resolution_low,
refinement=6,
structural_frame=structural_frame
)
gp.set_topography_from_arrays(grid=geo_model.grid,xyz_vertices=topo_input)
gp.map_stack_to_surfaces(
gempy_model=geo_model,
mapping_object={
# "First_Series": ('dmx'),
"Second_Series": ('Qcol+dl'),
"Third_Series": ('Qedl'),
"Fouth_Series": ('Qgl'),
# "Basement_Series":("Basement")
}
)geo_model.structural_frame.structural_groups[0].structural_relation = gp.data.StackRelationType.ONLAP
geo_model.structural_frame.structural_groups[1].structural_relation = gp.data.StackRelationType.ONLAP
geo_model.structural_frame.structural_groups[2].structural_relation = gp.data.StackRelationType.ONLAP
gp.set_section_grid(
grid=geo_model.grid,
section_dict={
'section10':([604414.9051,3005914.7854],[604084.8336,3007065.9975],[500,500]),
'section11':([604742.3340,3005914.7854],[604513.4478,3007065.9975],[500,500]),
'section5':([604071.3221,3006536.4907],[604779.6705,3006688.5895],[500,500]),
'section6':([604071.3221,3006438.3689],[604779.6705,3006309.6896],[500,500]),
'section7':([604071.3221,3006199.1698],[604779.6705,3006013.5121],[500,500]),
'section8':([604071.3221,3006587.8083],[604779.6705,3006553.1792],[500,500]),
'section9':([604071.3221,3006027.4361],[604244.7041,3005914.7854],[500,500]),
'section4':([604071.3221,3006861.3935],[604740.7076,3007065.9975],[500,500])
}
)
sol = gp.compute_model(geo_model,
engine_config=gp.data.GemPyEngineConfig(
backend=gp.data.AvailableBackends.PYTORCH,
use_gpu=True,
dtype="float64"
))'
This is my code and data

[NilsChudalla]

Hi @wdk11,
thanks for your questions.

1. Terrain is masking the solution and surfaces should not be affected shape-wise. I must say, that I do not understand or see the issue in the screenshots: surfaces computed by gemy are not exceeding the topography in your example. Also have a look at the scalar field, if the surfaces seem to make no sense to you. Could you clarify this point further?
2. Do you want to have the surfaces from the 3D plotter or simply the vertices?

• 3D object:

plotter = gpv.plot_3d(geo_model, show_lith=False, show=False)
plotter.p.export_obj("filename.obj") 

• vertices:

# vertices have to be backtransformed
back_transformed_vertices = geo_model.input_transform.apply_inverse(geo_model.solutions.dc_meshes[0].vertices)

3. Do I understand that right, that you have the topography as something you want to compute with GemPy? This is critical in your example, as I spot some regions in your model, where points are heavily overlapping. The artefacts you are computing could be a result from that.
4. I do not understand your problem here. It could be similar to the cause I am describing in 3.

When I have some spare time, I will run your model and have a look with the data.

[wdk11]

Hi @NilsChudalla . Thank you very much for your answer. Let me explain further here.
1.The strata I have framed in the picture have a downward trend when there is a terrain surface, but they are upward when there is no terrain surface. So I think if the geological mesh mesh I created is exported to another software, it will cross the terrain surface. If that's the case, the model I created is absolutely wrong. The geological mesh mesh cannot exceed the terrain surface. Therefore, I would like to ask Gempy if the geological and terrain surfaces are interpolated separately when modeling? So when drawing a 3D model, if show_topography=True, the terrain surface will be forced to cover the mesh mesh of the strata, so the 3D model I see is normal.
Simply put, I found that after setting the actual terrain surface, I wanted to ask if it was possible for the strata to cross the terrain surface, and if the terrain surface would forcibly prevent the strata from interpolating on its upper part. If it crosses, how to cut off the area where the strata cross the terrain surface; If it hasn't crossed, is there a problem with the display of my model.

2.I extracted the scatter data of the lower boundaries of various strata from the profile diagram, which may be very close in the vertical direction, but they do not overlap. The following figure shows the model and profile diagram I built. Here, I modeled the terrain surface as a stratum, that is, dmx. Then I found that the model I built, green represents the terrain surface, and other colors represent soil layers. In the figure, various colors are interspersed on the green, which means that the soil layer crosses the terrain surface. In fact, in the data I provided for modeling, the scatter data of dmx is always higher than that of any soil layer.

There is no problem with this situation in the two-dimensional profile, as the top layer is always green. What I want to say is that if I make a cut along the black line of the 3D model in the following figure, there is a yellow area in the middle at the top, so the displayed 2D section should be yellow at the top of the middle area. But in reality, this situation has never occurred before, so I am very puzzled about it. The situation shown in the 3D model does not match the 2D profile.

The above is my confusion.Look forward to your reply, Thank you @NilsChudalla