68. Creating finite fault.

[vincent1990-adongo]

I have a problem creating finite fault using example 68. Anytime I try to run this set of code, I get an error. I do not know exactly where this error is coming from.

gp.init_data(geo_model, [0, 1000, -500, 1500, -600, 0], [50,50,50],
surface_points_df=interfaces,
orientations_df=orientations,
default_values=True)

gp.map_stack_to_surfaces(geo_model,
{
'Fault1': ('Fault1'),
'Fault2': ('Fault2'),
'Strata1': ('Layer1', 'Layer2', 'Layer3'),
},
remove_unused_series=True)
geo_model.add_surfaces('Basement')
geo_model.set_is_fault(['Fault1', 'Fault2'])
geo_model.set_topography(source='random')
gp.set_interpolator(geo_model,
compile_theano=True,
theano_optimizer='fast_compile',
verbose=[],
update_kriging=False
)
geo_model.surfaces

###I encounted this kind of error [output]

Active grids: ['regular']

---

AssertionError Traceback (most recent call last)
Cell In[13], line 1
----> 1 gp.init_data(geo_model, [0, 1000, -500, 1500, -600, 0], [50,50,50],
2 surface_points_df=interfaces,
3 orientations_df=orientations,
4 default_values=True)
6 gp.map_stack_to_surfaces(geo_model,
7 {
8 'Fault1': ('Fault1'),
(...)
11 },
12 remove_unused_series=True)
13 geo_model.add_surfaces('Basement')

File ~\anaconda3\envs\gempy_tutorial\Lib\site-packages\gempy\gempy_api.py:192, in init_data(geo_model, extent, resolution, **kwargs)
190 kwargs['update_surfaces'] = False
191 if 'orientations_df' in kwargs:
--> 192 geo_model.set_orientations(kwargs['orientations_df'], **kwargs)
194 return geo_model

File ~\anaconda3\envs\gempy_tutorial\Lib\site-packages\gempy\core\model.py:1084, in ImplicitCoKriging.set_orientations(self, table, **kwargs)
1081 c = np.array(self._orientations._columns_o_1)
1082 orientations_read = table.assign(
1083 **dict.fromkeys(c[~np.in1d(c, table.columns)], np.nan))
-> 1084 self._orientations.set_orientations(
1085 coord=orientations_read[[coord_x_name, coord_y_name, coord_z_name]],
1086 pole_vector=orientations_read[[g_x_name, g_y_name, g_z_name]].values,
1087 orientation=orientations_read[
1088 [azimuth_name, dip_name, polarity_name]].values,
1089 surface=orientations_read[surface_name])
1091 self.map_geometric_data_df(self._orientations.df)
1092 self._rescaling.rescale_data()

File ~\anaconda3\envs\gempy_tutorial\Lib\site-packages\gempy\core\data_modules\orientations.py:103, in Orientations.set_orientations(self, coord, pole_vector, orientation, surface)
101 # Add nugget effect
102 self.df['smooth'] = 0.01
--> 103 assert ~self.df['surface'].isna().any(), 'Some of the surface passed does not exist in the Formation'
104 'object. %s' % self.df['surface'][self.df['surface'].isna()]

AssertionError: Some of the surface passed does not exist in the Formationobject. 3 NaN
4 NaN
Name: surface, dtype: category
Categories (2, object): ['Layer1', 'Layer2']

How do I resolve this error?
Also, the data 'intersection.csv' and 'orientations.csv' have not been included in your github repository, can have access to that data?

[github-actions]

Hello and welcome! Thanks for posting your first issue in the GemGIS project! Someone from our developers will get back to you. If your question is support related, we may transfer it to the Discussions.

[AlexanderJuestel]

orientations.csv
interfaces.csv

@vincent1990-adongo,
thanks for opening this discussion.

Please find the input data attached. Maybe this will solve the problem already.

[vincent1990-adongo]

Thanks so much Alexander. It was the input data just as you rightly indicated. I have been able to resolve that with your new data. I would like to plot only the finite fault without the original extent of the fault previously plotted in gray. Could you enlighten me on that?

I also tried to repeat the 'clip_fault_of _gempy_model' function to plot finite Fault 2,

mesh = gg.postprocessing.clip_fault_of_gempy_model(geo_model, fault='Fault2',
which='both',
buffer_first=0,
buffer_last=0)
mesh
[output]

{'Fault2': [PolyData (0x22a44b0cc40)
N Cells: 0
N Points: 0
N Strips: 0
X Bounds: 1.000e+299, -1.000e+299
Y Bounds: 1.000e+299, -1.000e+299
Z Bounds: 1.000e+299, -1.000e+299
N Arrays: 1,
'#527682']}
This was the ouput of the input data needed to create the finite fault2. It looks like it has 0 cells and points, and that the mesh cannot be plotted with input 0.

##I had the following error.
ValueError: Empty meshes cannot be plotted. Input mesh has zero points.

However, I do not see the difference between Fault1 and Fault2 apart from their orientation?

[vincent1990-adongo]

Any comments or suggestion @AlexanderJuestel

[AlexanderJuestel]

@vincent1990-adongo It may take a while to answer as we are also following regular jobs. So, a little patience would be appreciated.

[vincent1990-adongo]

Thank you very much. I hope you have understood what I am trying to point out @AlexanderJuestel ?

[AlexanderJuestel]

@vincent1990-adongo I can reproduce your error. I will have a look at what is going on!

[vincent1990-adongo]

@AlexanderJuestel Alright! It is important because to produce a complex geometry, different fault segments must be finite so they can be linked with each other. The finite fault concept must be applicable to all fault segments present.

I will be happy if you can help

[AlexanderJuestel]

@vincent1990-adongo , there is nothing wrong with the code. Depending on the direction of the fault, you may have to use the invert option. The following returns the correct mesh

mesh = clip_fault_of_gempy_model(geo_model,fault='Fault2', 
                                                   which='both', 
                                                   buffer_first=0, 
                                                   buffer_last=0,
                                                    invert_first=False,
                                                     invert_last=True)
mesh

[vincent1990-adongo]

@AlexanderJuestel. Thanks a lot for your response yesterday. The above code worked for the Fault2. I managed to plot both on the same plot figure.

I was trying to plot the finite fault only without the initial fault surfaces during the preprocessing phase. It looks like the finite fault is plotted on the infinite fault.
I am seeking to plot only the finite faults with the horizontal surfaces.

#Also, the fault is finite along the strike only. Is there any possibility of making the fault down the dip finite as well?
That is, there should be a user-controled fault length down the direction of dip.

You may kindly review this when you are free.

[AlexanderJuestel]

Hello @vincent1990-adongo,

just do not plot the original fault, see the code snipped below ;)

p = pv.Plotter(notebook=True)

p.add_mesh(mesh['Fault1'][0], color='red', opacity=0.5)
p.add_mesh(mesh1['Fault1'][0], color=mesh1['Fault1'][1])  #<-- plotting command for the original fault, not needed actually 

p.add_mesh(fault_polydata, color='red', point_size=40, render_points_as_spheres=True)

p.set_background('white')
p.set_scale(1,1,1)
p.show_bounds(color='black', font_size=12)

p.show()

[vincent1990-adongo]

@AlexanderJuestel Thanks a lot!

Any suggestions on plotting the finite fault and the horizon surfaces in one plot view. A truncation along the dip where the two segments intersect.

[AlexanderJuestel]

You can extract the meshes from the GemPy model, see tutorial 18 I believe. Then, plot everything in the same pyvista plot

[vincent1990-adongo]

@AlexanderJuestel. I tried that and it worked!

I have a little problem with the extension down the dip. I would like to control the fault length along the dip direction.

A common intersection line between the two segments. As shown in the images below.

Please, is there a way around that?

[AlexanderJuestel]

You can use fault relations as described on the GemPy documentation page to terminate one fault against the other. Is it that what you need?

[vincent1990-adongo]

@AlexanderJuestel

I think that is what I need but looking into previous tutorials I realized that the fault was infinite hence the interpolation goes beyond the points defined on the fault surfaces. Now the along strike is resolved, the other issue has to do with down-dip.

I will look into the fault relations again and see if I can figure out somthing!

[AlexanderJuestel]

@vincent1990-adongo

just let me know if you have further questions :)

[vincent1990-adongo]

@AlexanderJuestel
My gemPy is really messing up and I dont know why I am facing this error when I want to create a model

geo_model = gp.create_model('Graben_Model')
geo_model

IndexError Traceback (most recent call last)
Cell In[4], line 1
----> 1 geo_model = gp.create_model('Graben_Model')
2 geo_model

File ~\anaconda3\envs\gempy_tutorial\Lib\site-packages\gempy\gempy_api.py:114, in create_model(project_name)
99 def create_model(project_name='default_project') -> Project:
100 """Create a Project object.
101
102 Args:
(...)
112 TODO: Adding saving address
113 """
--> 114 return Project(project_name)

File ~\anaconda3\envs\gempy_tutorial\Lib\site-packages\gempy\core\model.py:1628, in Project.init(self, project_name)
1625 def init(self, project_name='default_project'):
1627 self.meta = MetaData(project_name=project_name)
-> 1628 super().init()

File ~\anaconda3\envs\gempy_tutorial\Lib\site-packages\gempy\core\model.py:85, in ImplicitCoKriging.init(self)
78 self._rescaling = ScalingSystem(self._surface_points, self._orientations,
79 self._grid)
80 self._additional_data = AdditionalData(self._surface_points,
81 self._orientations, self._grid,
82 self._faults,
83 self._surfaces, self._rescaling)
---> 85 self._interpolator = InterpolatorModel(self._surface_points,
86 self._orientations, self._grid,
87 self._surfaces,
88 self._stack, self._faults,
89 self._additional_data)
91 self.solutions = Solution(self._grid, self._surfaces, self._stack)
93 # Previous values of sfai.

File ~\anaconda3\envs\gempy_tutorial\Lib\site-packages\gempy\core\interpolator.py:650, in InterpolatorModel.init(self, surface_points, orientations, grid, surfaces, series, faults, additional_data, **kwargs)
646 def init(self, surface_points: "SurfacePoints", orientations: "Orientations", grid: "Grid",
647 surfaces: "Surfaces", series, faults: "Faults", additional_data: "AdditionalData",
648 **kwargs):
--> 650 super().init(surface_points, orientations, grid, surfaces, series, faults,
651 additional_data, **kwargs)
652 self.len_series_i = np.zeros(1)
653 self.len_series_o = np.zeros(1)

File ~\anaconda3\envs\gempy_tutorial\Lib\site-packages\gempy\core\interpolator.py:66, in Interpolator.init(self, surface_points, orientations, grid, surfaces, series, faults, additional_data, **kwargs)
63 self.aesara_graph = self.create_aesara_graph(additional_data, inplace=False)
64 self.aesara_function = None
---> 66 self._compute_len_series()

File ~\anaconda3\envs\gempy_tutorial\Lib\site-packages\gempy\core\interpolator.py:822, in InterpolatorModel.compute_len_series(self)
817 self.len_series_f = np.atleast_1d(len_series_f.astype(
818 'int32')) # [:self.additional_data.get_additional_data()['values']['Structure', 'number series']]
820 self._old_len_series = self.len_series_i
--> 822 self.len_series_i = self.len_series_i[non_zero]
823 self.len_series_o = self.len_series_o[non_zero]
824 # self.len_series_f = self.len_series_f[non_zero]

IndexError: invalid index to scalar variable.

I have spent the whole day trying to figure this out. It is really disturbing as i have reinstalled GemPy and install severally..

[vincent1990-adongo]

@AlexanderJuestel

When I used the current version of numpy numpy==1.26.1 I run into an error
when i use the lower version mumpy==1.23.5, I get into another error.

[AlexanderJuestel]

@vincent1990-adongo what is your pandas version? There has been a bug since pandas 2.0.2,so it is recommended to use pandas 2.0.1

[vincent1990-adongo]

@AlexanderJuestel I have to check and get back to you tomorrow. I left the lab a while ago. I will appreciate if we could arrange for a meeting via zoom to talk extensively anytime you are less busy. ..

[AlexanderJuestel]

Next week may work for a call

[vincent1990-adongo]

@AlexanderJuestel
Thanks a lot. I appreciate your time and effort. I will be happy to have a chat next week.

I resolved the issue with the numPy error by installing numpy==1.24.1

I am still struggling with the fault intersection along the dip incase you have any suggestions?

I will let you know of anything else!!

[vincent1990-adongo]

@AlexanderJuestel

Checking up on you and trying to know if you have anytime within the week for the short discussion.

I hope to hear from you soon

[AlexanderJuestel]

This week will not work. Maybe next week

[AlexanderJuestel]

@vincent1990-adongo,

how does the week of the 13th of November look for you in terms of a meeting?

[vincent1990-adongo]

@AlexanderJuestel. That will be great! I will appreciate a specific date. I will have my RSP meeting on 16th November and 17th is not bad at all.

I hope to hear from you soon.

[AlexanderJuestel]

Please send a message to [email protected] to talk further.