Question: How to add a fault in Resqml model with grid.

[AbdrahamaneBerete]

Hello, I'm trying to add a fault epc model with grid and one single. Here is my code:
def save_grid_to_resqml(x, y, z, properties, fault_data, resqml_file_path):
"""
Saves a grid with fault and additional properties to RESQML format using resqpy.
"""
model = rq.Model(epc_file=resqml_file_path, new_epc=True, create_basics=True, create_hdf5_ext=True)

crs = rqc.Crs(model)
crs.create_xml()

x_unique = np.unique(x)
y_unique = np.unique(y)
z_unique = -np.unique(z)

extent_kji = (len(z_unique), len(y_unique), len(x_unique))
dxyz = (x_unique[1] - x_unique[0], y_unique[1] - y_unique[0], z_unique[1] - z_unique[0])
origin = (x_unique[0], y_unique[0], z_unique[0])
#print(f"{origin =}")

grid = grr.RegularGrid(model, extent_kji=extent_kji, dxyz=dxyz, origin=origin, crs_uuid=crs.uuid, as_irregular_grid=True)

grid.make_regular_points_cached()

grid.write_hdf5()
grid.create_xml(write_geometry=True)

grid_x, grid_y, grid_z = np.meshgrid(x_unique, y_unique, -z_unique, indexing='ij')
points = np.vstack((x, y, z)).T
grid_points = np.vstack((grid_x.flatten(), grid_y.flatten(), grid_z.flatten())).T
values = np.array(properties, dtype=np.float64)
prop_array = griddata(points, values, grid_points, method='linear').reshape(grid_x.shape)

nan_mask = np.isnan(prop_array)
if np.any(nan_mask):
    nearest_values = griddata(points, values, grid_points, method='nearest').reshape(grid_x.shape)
    prop_array[nan_mask] = nearest_values[nan_mask]

prop_array = prop_array.transpose(2, 1, 0)

single_property = rqp.Property.from_array(
    parent_model=model,
    cached_array=prop_array.astype(float),
    source_info='Generated by script',
    keyword='Struct',
    support_uuid=grid.uuid,
    property_kind='Structural',
    indexable_element='cells',
    uom='m'
)

# Handling fault data
fault_indices = np.where(fault_data['Fault'].values == 1)[0]
fault_points = fault_data[['X', 'Y', 'Z']].values[fault_indices]
point_set = rqs.PointSet(model, points_array=fault_points, crs_uuid=crs.uuid)
point_set.write_hdf5()
point_set.create_xml()

fault_surface = rqs.Surface(model, crs_uuid=crs.uuid)
fault_surface.set_from_point_set(point_set)
fault_surface.write_hdf5()
fault_surface.create_xml(add_as_part=True, add_relationships=True, title="Fault Surface")

fault_feature = rqo.TectonicBoundaryFeature(model, kind='fault', feature_name='Fault1')
fault_feature.create_xml()
fault_interp = rqo.FaultInterpretation(model, tectonic_boundary_feature=fault_feature, is_normal=True)
fault_interp.create_xml(add_as_part=True, add_relationships=True)

if fault_interp.uuid is None:
    fault_interp.uuid = bu.new_uuid()
if fault_surface.uuid is None:
    fault_surface.uuid = bu.new_uuid()

# Create reciprocal relationship
model.create_reciprocal_relationship(fault_interp, 'destinationObject', fault_surface, 'sourceObject')

model.store_epc()

return model, grid

I have a problem with create_reciprocal_relationship: PS C:\Users\ABerete\Documents\code\3DGeoContext.Computation.Utils> & c:/Users/ABerete/Documents/code/3DGeoContext.Computation.Utils/.venv/Scripts/python.exe c:/Users/ABerete/Documents/code/3DGeoContext.Computation.Utils/mss3_utils/scripting/adding_fault_resqml_.py
Traceback (most recent call last):
File "c:\Users\ABerete\Documents\code\3DGeoContext.Computation.Utils\mss3_utils\scripting\adding_fault_resqml_.py", line 142, in
main(input_file_path, fault_file_path, resqml_file_path)
File "c:\Users\ABerete\Documents\code\3DGeoContext.Computation.Utils\mss3_utils\scripting\adding_fault_resqml_.py", line 135, in main
model, grid = save_grid_to_resqml(x, y, z, properties, fault_data, resqml_file_path)
File "c:\Users\ABerete\Documents\code\3DGeoContext.Computation.Utils\mss3_utils\scripting\adding_fault_resqml_.py", line 121, in save_grid_to_resqml
model.create_reciprocal_relationship(fault_interp, 'destinationObject', fault_surface, 'sourceObject')
File "C:\Users\ABerete\Documents\code\3DGeoContext.Computation.Utils.venv\lib\site-packages\resqpy\model_model.py", line 1867, in create_reciprocal_relationship
return m_x._create_reciprocal_relationship(self,
File "C:\Users\ABerete\Documents\code\3DGeoContext.Computation.Utils.venv\lib\site-packages\resqpy\model_xml.py", line 555, in _create_reciprocal_relationship
uuid_a = node_a.attrib['uuid']
AttributeError: 'FaultInterpretation' object has no attribute 'attrib'

[AbdrahamaneBerete]

@andy-beer I would to add a fault in epc model with grid and one single property. Could help me to find a solution for the error?

[andy-beer]

"I have a problem with create_reciprocal_relationship"

Hello,

The Model. create_reciprocal_relationship() method is rather old and is expecting xml tree root nodes rather than the resqpy objects you are passing to it. I recommend using Model.create_reciprocal_relationship_uuids() instead, which has a very similar signature except that it expects uuids instead of the xml nodes. So, try:

# Create reciprocal relationship
model.create_reciprocal_relationship_uuids(fault_interp.uuid, 'destinationObject', fault_surface.uuid, 'sourceObject')

[AbdrahamaneBerete]

Hi @andy-beer, thanks for your help. I did the change, there is no error again. But I can't see the fault when I import the EPC File in petrel. I see the grid with the property. Do have some examples (code with resqpy) on how to add a fault in epc model with grid and a property ??

Kind regards,

[AbdrahamaneBerete]

Hi @andy-beer, thanks for your help. I did the change, there is no error again. But I can't see the fault when I import the EPC File in petrel. I see the grid with the property. Do have some examples (code with resqpy) on how to add a fault in epc model with grid and a property ??

Kind regards,

[AbdrahamaneBerete]

Hi @emmanesbit and @cflynn3 , I hope you're well !!! Please, do you have some examples (code with resqpy) on adding a fault in epc model with grid and a property ??

Kind regards,

[emmanesbit]

Hi @AbdrahamaneBerete. Thanks for reaching out. Looking at your post above, I wonder if the issue might be around Petrel RESQML import capability, rather than anything missing in your code. In the Petrel version I am using, Petrel is able to import RESQML grids and associated properties, but does not import many other object types (though I am able to import these objects in other packages). I can't say for sure whether this is the same for your data, but if you have the option to load into other RESQML-enabled software, it might be worth a check.

[AbdrahamaneBerete]

Hi @emmanesbit , thanks for your answer. I tried to import some original EPC file with fault in my Petrel version and I can see the fault. The fault folder in the model is activated. have you tried to generate an EPC file with fault ??

[andy-beer]

Hello @AbdrahamaneBerete,

So, there are various classes of RESQML objects that have to do with geological faults, including:

• TectonicBoundaryFeature – your code should create that fine
• FaultInterpretation – your code should create one of those as well, and relate it to the feature
• PointSetRepresentation (resqpy PointSet) or TriangulatedSetRepresentation (resqpy Surface) – your code should generate one each of these, and create a 'soft' relationship between the surface and the interpretation
• GridConnectionSetRepresentation (resqpy GridConnectionSet, in the fault subpackage) – identifying faces in a grid which are considered as representing the fault surface
• actual grid geometry

Your code as it stands does not work on those last two points. As your grid has a regular geometry, I will assume that you do not want to represent the fault with a throw by explicit (irregular) geometry, so I won't comment more on the last point here.

To get a grid connection set, you can use functions in the resqpy grid_surface subpackage. The usual function to call is: find_faces_to_represent_surface(), with documentation here:
https://resqpy.readthedocs.io/en/latest/_autosummary/resqpy.grid_surface.find_faces_to_represent_surface.html#resqpy.grid_surface.find_faces_to_represent_surface

That function (which is computationally intensive) will be much more efficient if it 'knows' that your grid is a RegularGrid. As you are using the as_irregular_grid=True option when generating the grid object, if you reload that grid it will appear as a general resqpy Grid, rather than a RegularGrid. I would recommend setting as_irregular_grid=False unless you intend to modify the geometry.

One other small point: your 'soft' relationship between the fault surface and the interpretation might not be recognised by all software. That particular relationship can be embedded in the main xml for the surface (which is better) by:

• creating the feature and interpretation objects before the surface
• after instantiating the Surface object call: fault_surface.set_represented_interpretation_root(fault_interp.root) before creating the xml for the surface

Those enhancements should improve your RESQML dataset. However, Petrel only has very limited RESQML import capabilitiies. As a personal comment, I recommend you don't use Petrel.

[AbdrahamaneBerete]

Thank you so much @andy-beer.

[AbdrahamaneBerete]

I will try to follow your instructions. @andy-beer Please , I have two more questions:

1. which visualization tool can you propose me ?

2. If I want to create a Pillar Grid (irregular grid ) where the spacing dx, dy are constant in x and y directions and dz is not constant in z direction (dz varies between points in z directions), which kind of Class can I use? (GRID Object, Regular Grid, Unstructured GRID )

Kind regards,

[AbdrahamaneBerete]

Hi @andy-beer , I've followed your instructions with this code below:

def save_grid_to_resqml1(x, y, z, properties, fault_data, resqml_file_path):
"""
Saves a grid with fault and additional properties to RESQML format using resqpy.
"""
model = rq.Model(epc_file=resqml_file_path, new_epc=True, create_basics=True, create_hdf5_ext=True)

crs = rqc.Crs(model)
crs.create_xml()

x_unique = np.unique(x)
y_unique = np.unique(y)
z_unique = -np.unique(z)

extent_kji = (len(z_unique), len(y_unique), len(x_unique))
dxyz = (x_unique[1] - x_unique[0], y_unique[1] - y_unique[0], z_unique[1] - z_unique[0])
origin = (x_unique[0], y_unique[0], z_unique[0])


grid = grr.RegularGrid(model, extent_kji=extent_kji, dxyz=dxyz, origin=origin, crs_uuid=crs.uuid, set_points_cached = True, as_irregular_grid=False)
#grid = grr.RegularGrid(model, extent_kji=extent_kji, dxyz=dxyz, origin=origin, crs_uuid=crs.uuid, set_points_cached = True)

grid.make_regular_points_cached()

grid.write_hdf5()
grid.create_xml(write_geometry=True)

grid_x, grid_y, grid_z = np.meshgrid(x_unique, y_unique, -z_unique, indexing='ij')
points = np.vstack((x, y, z)).T
grid_points = np.vstack((grid_x.flatten(), grid_y.flatten(), grid_z.flatten())).T
values = np.array(properties, dtype=np.float64)
prop_array = griddata(points, values, grid_points, method='linear').reshape(grid_x.shape)

nan_mask = np.isnan(prop_array)
if np.any(nan_mask):
    nearest_values = griddata(points, values, grid_points, method='nearest').reshape(grid_x.shape)
    prop_array[nan_mask] = nearest_values[nan_mask]

prop_array = prop_array.transpose(2, 1, 0)

rqp.Property.from_array(
    parent_model=model,
    cached_array=prop_array.astype(float),
    source_info='Generated by script',
    keyword='Struct',
    support_uuid=grid.uuid,
    property_kind='Structural',
    indexable_element='cells',
    uom='m'
)

#Handling fault data

fault_feature = rqo.TectonicBoundaryFeature(parent_model= model, kind='fault', feature_name='Fault 1')
fault_feature.create_xml()
fault_interp = rqo.FaultInterpretation(parent_model= model,  title='Fault 1 Interpretation', tectonic_boundary_feature=fault_feature, is_normal=True)
fault_interp.create_xml(add_as_part=True, add_relationships=True)


fault_indices = np.where(fault_data['Fault'].values == 1)[0]
fault_points = fault_data[['X', 'Y', 'Z']].values[fault_indices]

point_set = rqs.PointSet(model, points_array=fault_points, crs_uuid=crs.uuid)
point_set.write_hdf5()
point_set.create_xml()

fault_surface = rqs.Surface(model, crs_uuid=point_set.uuid)
fault_surface.set_from_point_set(point_set)
fault_surface.set_represented_interpretation_root(fault_interp.root)
#fault_surface.set_model(model)
fault_surface.write_hdf5()
fault_surface.create_xml(add_as_part=True, add_relationships=True, title='fault 1 surface')

grid_connection = rqgs.find_faces_to_represent_surface(grid, fault_surface, 'fault 1', mode = 'regular', feature_type='fault')
grid_connection.write_hdf5()
grid_connection.create_xml()


model.store_epc()

return model, grid

I've visualized the generated EPC file in ResQML editor tool :

I have this error with GridConnection set object : !ENTRY EPC.editor 2 0 2024-08-08 14:32:20.718
!MESSAGE The feature 'count' of 'resqmlv2.impl.ObjGridConnectionSetRepresentationImpl@37d81587{file:/C:/Users/ABerete/Documents/code/3Dmodels/2023_06_12_Larzac_2/grid_resqml_with_fault_3/obj_GridConnectionSetRepresentation_9af40ae9-5581-11ef-bee5-30894ab42fb1.xml#//@gridConnectionSetRepresentation}' contains a bad value