loopstructural

Constrained implicit geological models

Use LoopStructural when observations constrain continuous geological scalar fields. A scalar value labels an interface; it is not automatically an age, measured depth or elevation. Record the coordinate CRS, shared XYZ units, positive-Z convention, feature order and meaning of each scalar level.

A constrained planar model

This synthetic example uses metre coordinates with a nonzero origin and an asymmetric box. Its scalar field is Z - 3030 + 0.1*(X - 1050) - 0.05*(Y - 2100) in metres. Scalar observations and gradient constraints anchor both level and orientation; a single point alone cannot define a geological surface.

from itertools import product
from LoopStructural import GeologicalModel
import numpy as np
import pandas as pd

origin = np.array([1000., 2000., 3000.])
maximum = np.array([1100., 2200., 3060.])
xyz = np.array(list(product([1010., 1050., 1090.], [2020., 2180.], [3010., 3050.])))
data = pd.DataFrame(xyz, columns=['X', 'Y', 'Z'])
data['feature_name'] = 'strat'
data['val'] = xyz[:, 2] - 3030 + .1 * (xyz[:, 0] - 1050) - .05 * (xyz[:, 1] - 2100)
data[['gx', 'gy', 'gz']] = [.1, -.05, 1.]
model = GeologicalModel(origin, maximum)
model.data = data
model.create_and_add_foliation('strat', interpolatortype='FDI', nelements=1000)
model.update()
query = np.array([[1020., 2040., 3020.], [1080., 2160., 3040.], [1050., 2100., 3030.]])
values = model.evaluate_feature_value('strat', query)

The expected values are approximately −10, +10 and 0. Pass world coordinates to model.evaluate_feature_value(..., scale=True) (the default). Direct feature evaluation has a different local-coordinate contract. GeologicalModel 1.8 accepts the two bounds positionally; do not assume origin= and maximum= constructor keywords from older examples still work.

Evaluate a VTK grid

Continue with the model above. regular_grid returns an N×3 array; disable shuffling and use Fortran ordering for VTK's x-fastest points. Write actual scalar values into point data before saving.