Obtain geological data

[wdk11]

I think the existing tutorials are all about drawing three-dimensional or two-dimensional images. I want to know how to obtain point data after building a geological model. For example, if I input a point, I can obtain the geological information of that point. For example, if I define a cross-section, I can input a point in the cross-section to obtain the geological information of any point in the cross-section.
Overall, I would like to know an example of an interface for obtaining geological information based on points.

[javoha]

Hello @wdk11,

thanks for the question. I do not know how deep you went into GemPy already but of course we are not drawing images. GemPy is based on an implicit interpolator (Universal Co-Kriging) the returns a scalar value corresponding to structural elements (e.g. rock units) everywhere in sapce. Tp create a model this space is dicretized. By default GemPy uses a regular grid for this discretization. The resolution of this grid is set when we create a model.

GemPy also allows to define other grids and multiple grids for a single model. The most common examples include custom cross sections and topography data. How to create these types of custom grids is covered in the first part of our new video tutorial series that I would highly recommend if you did not already watch it.

Apart from these specific grids you can also define custom grids by basically passing any XYZ coordinates to Gempy and calcualte the models solution at these coordinates. Here is a tutorial on how to this.

A shortcut version of this is to use:

gp.compute_model_at(geo_model, points)

where points is a numpy array of shape (n, 3) containing the xyz coordinates of your custom points. It will return the corresponding ID of the element at this location. This result is also stored under geo_model.solutions.raw_arrays.custom.

I hope this answers your question. Feel free to ask any fllow-ups. As a small note I would like to add that GemPy is an Open Source project and we answer issues and help out of personal motivation. We do appreciate people starting these questiosn with a "Hello" or anything similar.

Cheers,
Jan

[wdk11]

Thanks for your help!@javoha

[wdk11]

Hi,javoha @javoha .
Thank you for your help. I apologize, but I still have a few questions. My model computation takes too long, often several hours, and using gp.compute_model_at(geo_model, points) is extremely time-consuming. My specific idea is:

I want to build a high-resolution, large-scale model with a resolution of (500, 500, 500). Later, using interpolation algorithms (e.g., GemPy can export mesh grid data for the bottom surfaces of each stratigraphic layer), I could determine the geological information of a point by checking which two grids it lies between. I want to know if this approach is feasible.

Another question:

Since model computation is too slow, I want to save the entire model for future use. I used the following code:
import pickle
with open('geo_model.pkl', 'wb') as file:
pickle.dump(geo_model, file)

The resulting .pkl file is over 30 GB. Is this reasonable? Are there more efficient ways to store the model? Additionally, are there methods to accelerate computation while maintaining high precision (I tried PyTorch acceleration, but my GPU memory usage remained largely unchanged)?

My code:
data_path = 'C:/Users/Administrator/Desktop'
geo_model: gp.data.GeoModel = gp.create_geomodel(
project_name='Tutorial_ch1_1_Basics',
extent=[562298, 563540, 3032848, 3033348, 1800, 2700],
resolution=[500,500,500], # * Here we define the resolution of the voxels
refinement=4, # * Here we define the number of octree levels. If octree levels are defined, the resolution is ignored.
importer_helper=gp.data.ImporterHelper(
path_to_orientations=data_path + "/orientations.csv",
path_to_surface_points=data_path + "/points-100.csv",
)
)
gp.map_stack_to_surfaces(
gempy_model=geo_model,
mapping_object={
"Upper_Series": ('dmx','Qdl'),
"Lower_Series": ('Qgl')
}
)
sol = gp.compute_model(geo_model,
engine_config=gp.data.GemPyEngineConfig(
backend=gp.data.AvailableBackends.PYTORCH,
use_gpu=True,
dtype="float64"
))

Thank you!