Smooth & curved surfaces in gempy; nugget?

[shakasaki]

Hello,

I am trying to create smooth surfaces in gempy but i fail to do so with v3. I think it is because the nugget for surface orientations does not seem to change anything. Previously , there was a "smooth" value that actually did redice the weight of the normal vector for each point, so i coud use that to avoid "overly complex" surfaces being created. I show an example below, where I use three points to define a surface, all with normals facing vertically up. If the points deviate from the horizontal plane, even slightly, the resulting surface is nothing but "geological". Even if the deviation is rather small, like 5 m compared to 25 m inter-point spacing.

I recall that in previous versions of gempy i could reduce the weight of orientation points and the surfaces would simply follow the points, but not strictly match the orientation. I think changing the "nugget" in v3 does not really do anything. I tried the code below with nuggets ranging form 10^-4 to 10^4 with no considerable change in the surface.

Any suggestions on how to proceed?

Also, is there an effective way of moving around points and recomputing the surface quickly without having to recreate the model (re-instantiate) ?

surface with all points lying on a horizontal plane

surface with one point slightly higher than the rest

Minimal code example below:

`import gempy as gp
import gempy_viewer as gpv

#define resolution of the model

create a test geometry

points = [[25, 25, 25],
[50, 50, 25], # try changing the z value here slightly to see how it affects the plane
[70, 50, 25]]

normals = [[0, 0, 1],
[0, 0, 1],
[0, 0, 1]]

resolution = (50, 50, 50)

%%

gp_model = gp.create_geomodel(project_name='test nugget',
# These are the boundaries of the model
structural_frame=gp.data.StructuralFrame.initialize_default_structure(),
extent=[0,100,0,100,0,100],)

Add structural element to the geo model

gp_model.structural_frame.structural_elements[0].name = 'plane'
gp.add_surface_points(
geo_model=gp_model,
x=[point[0] for point in points],
y=[point[1] for point in points],
z=[point[2] for point in points],
elements_names=['plane']*len(points),
)

is there a faster way to add points and orientations at the same time?

for index, normal in enumerate(normals):
gp.add_orientations(
geo_model=gp_model,
x=[points[index][0]],
y=[points[index][1]],
z=[points[index][2]],
elements_names=['plane'],
pole_vector=[normal],
nugget=[1000],
)

%%

gp.compute_model(gp_model,)
gpv.plot_3d(
model=gp_model,
show_surfaces=True,
show_data=True,
show_lith=False,
image=False
)
`

[NilsChudalla]

Hi @shakasaki
a few very valid points.

• The shifting issue:
  By default, Gempy is scaling the input data, which can be beneficial in large scale models with structural complexity. In your case, it is actually backfiring. Therefore, the solution to your "non geological models" after shifting is disabling the transform befor computation:

gp_model.update_transform(gp.data.GlobalAnisotropy.NONE)

• The orientations issue:
  I do not have an example that allows me to investigate the nugget effect that much in detail right now at hand. To tweak the impact of orientations, you can also look at the length of your gradient. Maybe this will give you a desired result.

• Effective way to modify points:
  You could try gp.modify_surface_points() and gp.modify_orientations(), as alternative means. This comes with a word of caution: depending on the order in which you entered points, Gempy can reorder the surface points after computation. This means, point x might be point y after computing a model. This makes modifying points with these functions a bit tricky. This might be addressed in a future version of Gempy in some way. Until then, setting up the model in each iteration again might be the safest and most efficient method.

I hope these replies help :)

[shakasaki]

Hi @NilsChudalla you are right to the point! Disabling the transform works, but only if placed in the right place. If I do this before adding points (or orientations) this is again overwritten. It seems to work only if I do it just before the model computation.

You can see below the effect of the nugget. I also attach a short code that shows the nugget comparison for my scale of model (100 units in all 3 dimensions).

Which brings another question.. When trying to access the points to modify them, or even to plot them outside of gempy's plotter (as I do in the attached example) I notice that these seem to be somehow scaled. I found a scaling vector in there, but just wanted to check what the actual procedure is. Where could I access these points for plotting/visualisation/computation before altering them?

Also regarding modifying points.. I want to put a forward modeling scheme which uses gempy into PyMC for probabilistic inference, so any time I can save is precious. That's why I want to modify the points and recompute the model, instead of creating a new model from scratch..

Attached a picture of the results and the code to reproduce it. Orientation arrows are still badly plotted..

gempy_nugget_effect.py